Issue |
A&A
Volume 641, September 2020
|
|
---|---|---|
Article Number | A176 | |
Number of page(s) | 23 | |
Section | Planets and planetary systems | |
DOI | https://doi.org/10.1051/0004-6361/202038764 | |
Published online | 28 September 2020 |
The path to instability in compact multi-planetary systems
1
Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University,
Box 43,
22100
Lund,
Sweden
e-mail: antoine.petit@astro.lu.se
2
IMCCE, CNRS-UMR8028, Observatoire de Paris, PSL University, Sorbonne Université,
77 Avenue Denfert-Rochereau,
75014
Paris, France
3
Max-Planck-Institut für Astronomie,
Königstuhl 17,
69117
Heidelberg, Germany
Received:
26
June
2020
Accepted:
10
August
2020
The dynamical stability of tightly packed exoplanetary systems remains poorly understood. While a sharp stability boundary exists for a two-planet system, numerical simulations of three-planet systems and higher show that they can experience instability on timescales up to billions of years. Moreover, an exponential trend between the planet orbital separation measured in units of Hill radii and the survival time has been reported. While these findings have been observed in numerous numerical simulations, little is known of the actual mechanism leading to instability. Contrary to a constant diffusion process, planetary systems seem to remain dynamically quiescent for most of their lifetime before a very short unstable phase. In this work, we show how the slow chaotic diffusion due to the overlap of three-body resonances dominates the timescale leading to the instability for initially coplanar and circular orbits. While the last instability phase is related to scattering due to two-planet mean motion resonances (MMRs), for circular orbits the two-planets MMRs are too far separated to destabilise systems initially away from them. The studied mechanism reproduces the qualitative behaviour found in numerical simulations very well. We develop an analytical model to generalise the empirical trend obtained for equal-mass and equally spaced planets to general systems on initially circular orbits. We obtain an analytical estimate of the survival time consistent with numerical simulations over four orders of magnitude for the planet-to-star-mass ratio ε, and 6 to 8 orders of magnitude for the instability time. We also confirm that measuring the orbital spacing in terms of Hill radii is not adapted and that the right spacing unit scales as ε1∕4. We predict that beyond a certain spacing, the three-planet resonances are not overlapped, which results in an increase of the survival time. We confirm these findings with the aid of numerical simulations of three-planet systems with different masses. We finally discuss the extension of our result to more general systems, containing more planets on initially non-circular orbits.
Key words: celestial mechanics / planets and satellites: dynamical evolution and stability
© ESO 2020
1 Introduction
One of the most astonishing results of the Kepler mission has been the discovery of very compact super-Earth multiplanetary systems (Borucki et al. 2011; Fabrycky et al. 2014). These systems, such as Kepler-11 (Lissauer et al. 2011a), can host more than six planets with masses between that of the Earth and Neptune, and with periods of less than 100 days. They have very low mutual inclinations and eccentricities (Johansen et al. 2012; Fang & Margot 2012; Xie et al. 2016) and for the majority, they are not in resonant chains (Lissauer et al. 2011b; Fabrycky et al. 2014). Understanding the orbital properties of these so-called super-Earths or mini-Neptunes is crucial, as it seems that at least 50% of solar-type stars host a close-in planet with a radius comprised between that of Earth and Neptune (Mayor et al. 2011; Petigura et al. 2013; Fressin et al. 2013).
Studies of the Kepler multiplanetary systems have shown that the architecture is most likely sculpted by dynamical stability (Johansen et al. 2012; Pu & Wu 2015). Indeed, it has been shown that the minimum spacing is mass dependent (Weiss et al. 2018), with a lower limit in observed Kepler systems of around 10 Hill radii. As a result, understanding the mechanism leading to the instability of more tightly packed systems is critical to our understanding of planet formation and architecture.
The question of the stability of exoplanetary systems is particularly challenging due to several factors. The observed close-in planets have most likely performed at least 109 to 1011 orbits since their formation, which makes the numerical integration extremely costly if one wants to integrate the system over its whole lifetime. Because of the age of exoplanetary systems, it is often assumed that the systems are stable to constrain the orbital configuration. As a result, in order to understand the architecture of planetary systems, the stability analysis is complementary to observations (Laskar & Petit 2017). Nevertheless, the process is made even more costly because we do not know the exact orbital configuration, let alone the planet masses for systems detected by transits. But even if the present orbital configuration were known perfectly, planetary systems are chaotic, as has been shown for our own Solar System (Laskar 1994; Laskar & Gastineau 2009). As a result, the only approach to a numerical stability analysis is to run several integrations with slight variations of the initial conditions to probe the outcome in a statistical manner. Therefore, for each exoplanetary system, thousands of very costly numerical integrations would need to be run in order to obtain a satisfying understanding of its stability properties. The process could eventually be sped up thanks to the help of machine learning classification (Tamayo et al. 2016, 2020).
Another approach is to rely on analytical stability criteria. Under specific assumptions, it is possible to simplify the dynamics to obtain models accurately describing the behaviour of the system. In particular, one can derive stability criteria that can delineate stable regions from unstable ones where systems will eventually experience close encounters and collisions. Among such analytical criteria, one can cite the Hill stability (Marchal & Bozis 1982; Gladman 1993; Petit et al. 2018) and the overlap of mean motion resonances (MMRs; Wisdom 1980; Deck et al. 2013; Petit et al. 2017; Hadden & Lithwick 2018). For less compact, non-resonant systems, the dynamics are very well approximated by the secular model. In the secular approximation, one averages over the fast motion of the planets on their Keplerian orbits to only consider their long-term deformations. A well-known consequence of this averaging is the conservation of the planet semi-major axes, and thus of the angular momentum deficit (AMD; Laskar 1997, 2000). The AMD gives a dynamically motivated measure of the total eccentricities and mutual inclinations in a planetary system, and thus acts as a dynamical temperature. In particular, if the AMD is low enough, there is no possible orbital rearrangement allowing for close planetary encounters. This concept has been defined as the AMD-stability (Laskar & Petit 2017); it allows for a fast characterisation of the stability of planetary systems away from MMRs, where the secular approximation is valid. Besides the AMD-stability, the AMD has proven to be a versatile tool to understand planet dynamics (e.g. Volk & Malhotra 2020).
However, the transition from the secular regime to regions where the fast interactions between planets shall not be neglected is unclear. This is due to the influence of MMRs which forbid independent averaging over the fast angles of the planets (although it should be noted that an extension of the Lagrange-Laplace secular theory in the vicinity of MMRs is possible Libert & Sansottera 2013; Sansottera & Libert 2019). While theoretical studies in the two-planet case have allowed a sharp limit to be found between the secular and non-secular regions (Hadden & Lithwick 2018; Petit et al. 2018, and references therein), there are no complete studies for three-planet systems and higher. Numerical simulations (Chambers et al. 1996, and references in Sect. 2) have shown a qualitative change in behaviour between two-planet systems and three-planet systems and beyond: multi-planetary systems experience a long quiescent phase wherethe systems are almost secular before a very rapid transition to collisional dynamics. Preliminary analytical studies were proposed by Zhou et al. (2007) and Quillen (2011), but their models did not entirely reproduce the characteristics of the transition zone between long-lived systems and systems where scattering occurs immediately.
The present work attempts to study the mechanism leading to the instability of tightly packed systems. Since the differentstability regime between two-planet and multi-planet systems starts at three planets, we focus on systems composed of three planets. Contrary to previous studies, we do not make any assumptions regarding the masses of the planets (providing that they remain small) and consider unevenly spaced planets. However, we restrict ourselves to initially circular and coplanar systems. Indeed, due to interactions with the protoplanetary disk, compact, close-in systems most likely form in this state due to eccentricity and inclination damping (Lin & Papaloizou 1986). We note that we do not consider planets trapped into resonant chains here and refer to Pichierri & Morbidelli (2020) for an analytical study of stability of resonant chains. In addition, we are interested in systems that should be considered AMD-stable in the sense that no secular interactions can lead to their instability (Petit et al. 2018). Understanding the initially circular systems gives a lower bound for the eccentric ones. By analysing individual simulations, we postulate, as in Quillen (2011), that the instability is driven by the overlap of MMRs between the three planets of each system. Their prominent role comes from the presence of a dense subset of three-planet resonances that covers a large part of the phase space, even for circular orbits. Moreover, the system dynamics in the presence of this subset are not secular, yet they preserve the total AMD, which is a characteristic observed in numerical simulations. This feature is explained in Sect. 4.2. Using estimates of the diffusion rate proposed by Chirikov (1979), we are able to compute an analytical expression for the survival time.
Our analytical approach allows us to determine features in numerical simulations that trace the particular mechanism we study, which leads us to conclude that we isolated the right mechanism for planetary instability. In particular, we confirm that the scaling in terms of Hill radius, widely used in numerical studies (Chambers et al. 1996; Smith & Lissauer 2009; Pu & Wu 2015; Obertas et al. 2017), is not appropriate. By comparing with numerical simulations, we show that our time estimate is valid over four orders of magnitude in mass and almost seven orders of magnitude in survival time.
In the context of exoplanet observations, three-planet resonances are particularly significant as it is possible to assess their dynamical influence from transit data alone (Delisle 2017). They can also be a signpost of the disruption of MMR chains thanks to tidal dissipation (Charalambous et al. 2018; Pichierri et al. 2019). Nevertheless, the interactions between such resonances has not been exhaustively studied.
In Sect. 2, we begin by a review of the works on the problem of tightly packed planetary systems and we perform an in-depth qualitative analysis of the instability. In Sect. 3, we introduce our framework to treat the problem of three-planet MMR. Section 4 contains most of the technical details. We first describe the network of zeroth-order three-planet resonances. We then solve the dynamics for an isolated MMR to finally obtain a criterion delimiting the region where the MMRs overlap. Using the framework developed in Sect. 4, we estimate in Sect. 5 the survival time for a system of three planets, with arbitrary mass distribution and spacing (assuming that the planets are not too massive and tightly packed). We compare our analytical results to numerical simulations in Sect. 6. Finally, we discuss possible extensions to more general systems than three planets on circular and coplanar orbits in Sect. 7. While the analytical derivations make it necessary to define auxiliary variables, we tried where possible to use only variables with a clear physical meaning in the figures to help those readers willing to skip the technical sections.
2 Qualitative description of the instability
The dynamics of tightly packed systems are chaotic, and research on the subject has mainly focused on a qualitative description of their behaviour due to the difficulty of the analytic approach. We review the qualitative description proposed by previous studies and highlight how the instability is triggered.
2.1 Stability in the two-planet case
While the three-body problem is not integrable in general, the problem of the stability of a two planet system is well understood. Most of the stability results come from the existence of a topological boundary in the three-body configuration space leading to the so-called Hill-stability (Marchal & Bozis 1982). In a Hill-stable system, the two planets can never approach one another, which leads to a sharp difference in behaviour. The Hill-stability was popularised by Gladman (1993) for circular orbits, as a minimal distance between orbits normalised by their Hill radius guaranteeing the system’s stability. This stability criterion can be written as
(1)
is the mutual Hill radius with m1, m2 being the planet masses and m0, the star mass. For inclined and eccentric orbits, there exists a critical AMD value depending only on semi-major axis and masses such that a system with a smaller AMD is Hill stable (Petit et al. 2018).
Another stability criterion for a two-planet system can be derived from the overlap of MMR (Wisdom 1980; Deck et al. 2013; Petit et al. 2017; Hadden & Lithwick 2018). While the unperturbed resonant problem is integrable, the interaction between neighbouring MMRs leads to the formation of a chaotic web such that the planets’ orbital elements wander in a random walk fashion. This behaviour is known as the Chirikov (1979) diffusion. For initially circular orbits, the overlap occurs at a distance scaling as (Wisdom 1980). The exponent 2/7 is close to 1/3 but it has been highlighted that there exists a regime where MMRs overlap while the planets are Hill-stable, that is, the system is long-lived while experiencing short-term chaos (Deck et al. 2013; Petit et al. 2018).
This means that a two-planet system is either stable over timescales comparable with the lifetime of the host star or unstable in a very short amount of time (less than 105 revolutions). No such dichotomy is observed for multiplanetary systems. Indeed, a multiplanet system can appear stable if it is numerically integrated over a few million orbits while becoming unstable in less than a billion years.
2.2 Survival time of tightly packed systems
The pioneering work on the stability of tightly packed multi-planetary systems was carried out by Chambers et al. (1996). These latterauthors performed numerical simulations of systems with equal-mass planets on initially equally spaced circular and coplanar orbits (hereafter referred to as EMS systems). The constant orbit spacing is given by Δ = (ak+1 − ak)∕ak. For various planetary masses and numbers of planets, these latter authors recorded the survival time of a system, defined as the integration time before the distance between the two planets becomes smaller than a Hill radius. As shown by Rice et al. (2018), the time between such a close encounter and the proper collision is usually negligible. Chambers et al. (1996) observed that the survival time grows exponentially with the spacing Δ rescaled by the Hill radius RH,
(3)
where P is a typical orbital period and b and c are numerical factors1. Here, b seems to have a small dependency on the mass ratio and the number of planets, and c also seems to depend on the mass ratio. Analysing Fig. 4 from Chambers et al. (1996), a more appropriate scaling seems to be
(4)
where b′ and c′ are positive numerical coefficients independent of the masses, mp is the planet mass and m0 the star mass. We note that such scaling was also chosen by Faber & Quillen (2007).
Subsequent numerical works on the stability of EMS have been carried out. As the computational capacities increased, Smith & Lissauer (2009) and then Obertas et al. (2017) obtained datasets with a much finer distribution of spacing and longer integration times showing systems becoming unstable after almost 10 Gyr. Beyond the trend already observed by Chambers et al. (1996), Obertas et al. (2017) showed that the survival time is reduced in the vicinity of low-ordertwo-planet MMR. Hussain & Tamayo (2020) show that the spreading around the linear trend for log Tsurv∕P is roughly constant and follows a normal distribution with a standard deviation of 0.43 ± 0.16 dex, indicating that the instability emerges from a chaotic diffusion process. Beyond the EMS initial conditions, Pu & Wu (2015) explored theimpact of small variations of the initial conditions by drawing the spacings, eccentricities, and inclinations from distributions and showed that the exact spacing can be replaced by the minimal separation between the orbits. From these studies, the minimal spacing ensuring the stability over a few billion orbits can be estimated to be around 10 Hill radii.
Following Chambers et al. (1996), most of the previously cited studies fit the survival time with curves similar to Eq. (3) because of the natural parallel with the two-planet case. However, there is no generalisation of the Hill stability in the multi-planet case and the mechanism leading to instability has a priori no reason to be related to the Hill scaling RH. The discrepancy between the two proposed mass renormalisations in Eqs. (3) and (4) is easily explained by the fact that most studies only considered a limited mass range and very small difference between the exponents. Zhou et al. (2007) estimated Tsurv as a power-law in the spacing and using Nekhoroshev estimates, Yalinewich & Petrovich (2020) proposed a scaling similar to Eq. (4).
2.3 Phenomenology of the instability
These qualitative and quantitative studies on EMS systems highlight the key features that the tightly packed system instability presents and that an analytical model should explain.
- a.
The survival time Tsurv seems to have an exponential dependency on the orbital spacing, measured in units of
. The fit is valid over 6 to 8 orders of magnitude for survival times between 100 and almost 1010 orbits. The higher end is limited by computational time. However, the physical interest to go beyond is limited as it approaches the lifetime of the central star in most cases.
- b.
Instabilities occur for spacings larger than that leading to two-planet instabilities. As a result, it is an intrinsically multi-planet phenomenon. In addition, Chambers et al. (1996) have shown that the results were unchanged in systems of four or more planets if the planet interactions are limited to their neighbours. Thus, three planets are necessary but also sufficient to reproduce the effect.
- c.
Systems initially on circular orbits, and therefore AMD-stable, can become unstable. The mechanism at play is thus by nature non-secular and involves some kind of MMR overlap despite the fact that two-planet MMRs do not overlap in the range where the instability can occur. However, the AMD does not evolve regularly during the lifetime of the system. Indeed, as shown in Fig. 2a system can experience almost no AMD evolution during most of its lifetime before a rapid increase shortly before instability.
- d.
The survival time distribution suggests that the evolution is driven by a diffusion process (Hussain & Tamayo 2020). The dips close to first-order two-planet MMRs indicate that these latter play a fundamental role in enabling the orbit crossing.
While the stability of EMS systems has been described extensively from numerical simulations, very few works have developed an analytical framework attempting to describe the observed behaviour. In the most elaborate model, Quillen (2011) proposed that the instability is driven by the overlap of zeroth-order three-planet MMRs. Resonances involving more than two planets emerge as the result of the first-order averaging (e.g. Chap. 2, Morbidelli 2002, see also Sect. 3) and are weaker than the two-planet MMR. Quillen (2011) shows that, despite their smaller width, the three-planet MMRs are more numerous and overlap at larger spacing and smaller eccentricities. The ansatz is that the semi-major axes of the planets evolve randomly through the rich network of these three-planet MMRs until a first-order two-planet MMR is encountered, leading to a rapid AMD increase, and close encounters and collision shortly afterwards. Moreover, the main resonances close to circular orbit preserve the total AMD (see Sect. 4), which is consistent with simulations. We highlight the fact that, in the secular dynamics of the Solar System, slow diffusion leads to a region where the system becomes rapidly unstable (Laskar 1994; Batygin et al. 2015).
To illustrate the mechanism leading to instability, we perform the numerical integration of a typical EMS system. The planets have a mass mp = 10−5 M⊙ and orbit a solar-mass star. The inner orbit is at 1 au and the period ratios between adjacent planets is initially close2 to Pk+1 ∕Pk = 1.175. This particular value was chosen in order to observe the instability after roughly a few million orbits of the inner planet whilebeing outside of a two-planet MMR island. The orbits are initially circular and coplanar and the angles drawn randomly. As in previous studies, we run the simulations up to the first close encounter. In the considered case, the integration lasts 3.33 Myr. The system is integrated with the hybrid integrator MERCURIUS (Rein et al. 2019)from the REBOUND code (Rein & Liu 2012) with a time-step of 0.01 yr. The relative energy error is 5 × 10−10.
Figure 1 shows the evolution of the semi-major axis of the three planets. The envelope around the curve corresponds to the extent of the orbits, that is, the position of the periapses and apoapses, and is thus a measure of the eccentricities of the orbits. The curves are smoothed by performing a rolling averaging over the next ten snapshots3. The verticalaxis is discontinuous to highlight the small variations during the large majority of the integration. As already described by previous authors, the system appears quiescent for the majority of its lifetime. Subsequently, after the pair 2-3 crosses the 7:6 resonance, the system becomes unstable in 127 kyr. This figure emphasises the timescale difference between the lifetime of the system and the proper unstable phase that is almost two orders of magnitude shorter. Explaining the lifetime of tightly packed systems should therefore focus on the quiescent phase as the timescale to reach instability is dominated by this phase.
To show the rapid change of behaviour before the close encounter, Fig. 2a shows the evolution of the two adjacent period ratios Pk∕Pk+1 as a function of the time to the close encounter (we note the logarithmic scale). In Fig. 2b, we plot the evolution of the AMD C of the system (see Eq. (7)) rescaled by the total angular momentum G. The plotted quantity, , scales linearly with eccentricity for close to circular orbits. At the moment when the pair 2-3 enters the 7:6 MMR region, the system enters the scattering phase. This is also the moment where the AMD starts to increase. Nevertheless, the initial phase is not secular despite the near conservation of the AMD; indeed we see that the period ratios are not constant but evolve over a long timescale.
The interaction and the position of the system with respect to the network of two-planet MMRs seems critical to the duration of the quiescent phase. However, Fig. 2 merely shows how the instability is triggered and not the mechanism leading to it. The slow evolution of the system is seen much more explicitly in the period ratio plane plotted in Fig. 3. We plot the period ratio of the outer pair P2 ∕P3 as a function of the inner pair period ratio P1∕P2. In this plane, the two planet MMRs are vertical and horizontal dashed lines. We plot the neighbouring first-order MMR, that is, the 7:6 and the 6:5 in black, indicating their approximate extent for the circular orbit in grey. The second-order resonance 13:11 is plotted in red, but it has a null width for circular orbits. We see the system starts outside of the two-planet MMR. However, ontop of the two-planet MMR network, there also exists the network of three-planet MMRs. For circular orbits, the main three-planet MMRs are of zeroth order (see Sect. 4.1). We plot the loci of the largest three-planet MMR in the vicinity of the initial condition: this network is composed of a set of almost parallel lines which run transversally to the two-planet MMR lines. As predicted by Chirikov (1979) theory, the diffusion takes place perpendicularly to the network of the three-planet MMRs, up until the system reaches the two-planet resonances where the trajectory wanders around rapidly.
This qualitative analysis seems to confirm Quillen’s hypothesis. The survival time is dominated by the diffusion along the three-planet MMR network and the system becomes unstable once it reaches the two-planet resonance where chaotic diffusionis faster and rapidly increases the total AMD. The survival time can be estimated by computing the diffusion rate according to Chirikov’s resonance overlap theory (see Sect. 5.1). The scaling law for the survival time obtained by Quillen (2011) is a very steep power-law instead of having an exponential behaviour. In particular, the timescale is overestimated at short separations and underestimated for large ones. Quillen’s result and its difference with numerical simulations can be explained by some simplifications made in the computations leading to an inexact determination of the effective diffusion rate as well as a limit of the three-planet MMR overlap.
In this study, we consider the general, circular, coplanar three-planet problem. We remain in the framework of tightly packed systems but we relax the assumption on the initial equal spacing and equal masses. We show that it is possible to use Chirikov’s theory to explain the observed survival time scaling.
![]() |
Fig. 1 Semi-major axis as a function of time for an example of a three-planet EMS system (see the text for the full initial conditions). The envelope of the curve represents the extent of the orbits. We note the discontinuous vertical axis. We show the time where the first main two-planet MMR is crossed. The system becomes unstable soon afterwards. |
![]() |
Fig. 2 Panel a: period ratio of the adjacent pairs as a function of the time to the close encounter. We note the logarithmic scale. The vertical dashed line is the same as in Fig. 1. The black horizontal dashed lines corresponds to first-order MMRs, the yellow dashed lines to the second-order MMRs. The width of the first-order MMRs is displayed in grey. Panel b: AMD normalised by the total angular momentum as a function of the time to the close encounter. Using the square root gives a typical value of the planet eccentricities. |
![]() |
Fig. 3 Evolution in the period ratios plane. The points are colour-coded according to their time before the close encounter. We note that the system spends almost all of its time very close to its starting location as shown in Figs. 1 and 2. Green oblique lines correspond the loci of the zeroth-order three-planet MMRs (see Sect. 4.1). Chirikov’s diffusion is expected to occur perpendicularly to the network direction, along the dashed orange line. The extent of the adjacent first-order two-planet MMRs, 7:6 and 6:5, is plotted in grey. The width is computed for circular orbits. The second-order resonance 13:11 is plotted in red. Once the system enters the two-planet resonance network, the diffusion is much more rapid. |
3 Problem considered and mean motion resonances
We summarise most of the notations in Table A.1. We consider a system of three planets of masses m1, m2, and m3 orbiting a star of mass m0. The canonical positions rj and momenta are expressed in canonical heliocentric coordinates (Poincaré 1905; Laskar 1991). The initial orbits are assumed to be circular and coplanar. Let the semi-major axes aj, the eccentricities ej, the mean longitudes λj, and the periapses longitude ϖj be the orbital elements defining the orbits. A set of canonical coordinates for the system is given by the modified Delaunay coordinates (Laskar 1991):
(5)
where and
is the gravitational constant. We note that the gravitational parameter μ is the same for all three planets as in Laskar & Petit (2017). This is possible if we consider the so-called democratic-heliocentric formulation of the planetary Hamiltonian (e.g. Morbidelli 2002). The couples of variables (Cj, −ϖj) can also be replaced by their associated complex variables
(6)
with (xj are the canonical momenta and
the conjugated positions). For small eccentricities, we have
. The system total angular momentum G and AMD C are given by
(7)
The Hamiltonian describing the dynamics can be split into an integrable part
(8)
describing the motion on unperturbed Keplerian orbits, and a perturbation,
(9)
describing the planet interactions. Here, ε is a dimensionless parameter of the order of the planet-to star-mass-ratio to reflect the scale difference between the two parts of the Hamiltonian. In terms of Poincaré coordinates, the perturbation part can be written as
(10)
where ,
.
Due to the conservation of angular momentum, the coefficient must vanish unless the indices
verify the d’Alembert rules (e.g. Morbidelli 2002) and in particular
(11)
In the unperturbed case, the system is said to be in a MMR if the mean motions
(12)
verify an equation of the form
(13)
The sum k = k1 + k2 + k3 is the “order” of the resonance. The sum K = |k1| + |k2| + |k3| is the “index” of the resonance. In the general case, the perturbation also influences the resonantdynamics. The terms in Eq. (10) contributing to the resonance are the ones that depend on the combination of mean longitudes k1λ1 + k2λ2 + k3λ3. Because of d’Alembert rules, the leading order term in the perturbation is of order k in eccentricities (k being the resonance order).
We note that because each term in Eq. (9) only contains contributions from two planets, this is also the case for Eq. (10). In particular, there are no terms in the non-averaged Hamiltonian that depend on angles of theform k ⋅λ with kj≠0 for all j. This means that there are no three-planet resonances at the first order in ε. There are instead of course
two-planet MMR terms.
Three-planet MMRs actually emerge in the perturbative Hamiltonian as terms which appear after applying a perturbation step which eliminates the fast angles λ to first order in ε (this step is sometimes referred to as averaging because to first order in ε it is equivalent to averaging out the fast angles from the Hamiltonian). To do so, we assume that the system is “far enough” from the first-order two-planet MMR4, that is, we assume that
is not too small with respect to ε for some integer k (Morbidelli 2002). This is for example the case for the system considered in the previous section. We sketch the main lines of these perturbative steps below, and we carry out the explicit calculation of the relevant terms in the following section.
As we consider systems far enough from two-planet resonances, we can perform one perturbation step and keep track of all terms up to order in the equations. We use the classical approach from the perturbations theory, the Lie series method (Deprit 1969). We refer to Sect. 2.2 of Morbidelli (2002) and references therein for a complete description of the method. The method has already been applied to provide an analytical model of three-body resonances when one of the bodies is a test particle (Nesvorný & Morbidelli 1998). The idea is to introduce a new set of variables (noted with a prime in the following equations) ε-close to the original ones such that in these new variables the transformed Hamiltonian writes
(14)
where is the average of
over the mean longitudes λ, and
is the leading-order term of a series in ε. The transformation can be explicitly constructed as the flow between 0 and 1 of a generating Hamiltonian vector field exp ({εχ1, ⋅}) where {⋅, ⋅} is the Poisson bracket5 and εχ1 is the solution of the homological equation
(15)
More precisely, if we note the complex Fourier coefficients of
with respect to the mean longitudes, we can write
(16)
Due to the expression of given in Eq. (10), the denominators k ⋅n are of the form
and are not “too small” because we assume the system to be far from two-planet MMRs. Thus the formal series (16) is formally well defined; one can stop the summation at indices k of sufficiently high order so that the remaining Fourier terms in
have sizes smaller than ε2, which is ensured by the exponential decay of the Fourier coefficients. Thus, the solution (Eq. (16)) to the homological equation (Eq. (15)) is well defined. We can express the Hamiltonian
as
(17)
The Poisson bracket in Eq. (17) generates terms involving all three mean longitudes. In other words, three-planet MMRs that were not present in the initial Hamiltonian cannot be neglected at second order in averaging. The study of a particular three-planet MMR can be done by a second averaging over the other fast angles, because all other terms will not contribute small divisors and can therefore be eliminated by an additional perturbative step. In practice, this results in another change of coordinates, which are ε2 close to the first-order averaged coordinates, and the new Hamiltonian is the average of with respect to the fast (e.g. non resonant) angles (see following section). Because we do not need to change back to the initial coordinates,hereafter we drop the primes on the coordinates and Hamiltonian. We also drop the terms of order ε3 and greater.
4 The three-planet zeroth-order resonance network
In Fig. 3, we see that the diffusion mainly occurs perpendicularly to the zeroth-order three-planet MMRs. This is expected forclose to circular orbits because the resonant coefficients do not depend on eccentricity at the leading order in eccentricity. In addition, the structure of the network is easier to describe. We make the hypothesis that the zeroth-order three-planet MMRs are sufficient to explain most of the diffusion leading to the instability. This assumption is well supported in Sect. 6, where we compare the analytical prediction of survival times calculated under this hypothesis with the results of numerical simulations. We analyse these MMRs and compute an overlap criterion in this section. We consider the role of additional MMRs in Sect. 7.
4.1 Network description
A zeroth-order three-planet MMR can be described by two integers p and q. The resonance equation is
(18)
Since such resonance does not depend on the longitude of the periapses, the AMD is unaffected by the resonant terms (see below). We therefore restrict ourselves to the three6 degrees of freedom (Λj, λj). The resonance equation defines a plane in the frequency space (n1, n2, n3). Because the gravitational interactions are scale invariant, we can restrict ourselves to a two-dimensional plane corresponding to the period ratios ν12 and ν23 where
(19)
Dividing Eq.(18) by n2, and reorganising terms, one gets
(20)
that is the equation of a straight line passing through the point (1,1) with slope − p∕q for the period ratios ν21 = P2∕P1 > 1 and ν23 = P2∕P3 < 1. While the resonances can be interpreted easily geometrically with these two period ratios, the fact that one of the variables is larger than 1 and the other smaller can be confusing. Moreover, expanding the period ratio ν21 for tightly packed planets (i.e. close to 1) leads to a poorer approximation at first order than expanding ν12. We therefore only consider the variables ν12 and ν23 as done in Eq. (20). In the plane (ν12, ν23), the resonance loci are hyperbolas passing through (1,1); they however behave to a very good approximation as straight lines with slopes p∕q for tightly packed systems. As shown in the following section, the strength of the resonances depends strongly on their index 2(p + q).
Figure 4 shows the loci of the zeroth-order three-planet MMRs such that p + q < 50 (the dashed black lines correspond to two-planet first-order MMR). Because the slope at the point (1, 1) is p∕q, the resonances are not spread uniformly. Indeed, higher index resonances can lie on top of lower index ones if p and q are not coprime, for example where the resonance 2n1 − 4n2 + 2n3 = 0 is the same as the resonance n1 − 2n2 + n3 = 0. The full network is dense in the (ν12, ν23) plane. However, the resonance width may become so small as p + q increases that the resonances do not interact.
The resonance loci do not depend on the MMR index, but only on the ratio p∕(p + q). We choose this specific ratio as it lies between 0 and 1. One can see that p∕(p + q) can be extended as a continuous function in the period ratio plane. An adapted set of coordinates to describe the period ratio plane can be defined7 to take advantage of this property. We define the resonance locator
(21)
The second equality is only valid on resonances. The position along the resonance can be defined by a generalised period ratio separation
(22)
where η is a constant on a specific resonance whereas ν is a hyperbola along which the resonance strength is roughly comparable.
The variables (ν, η) are well adapted to describe the dynamics governed by the three-planet MMRs. We can express the period ratios as a function of these variables
(23)
The levels of constant ν are hyperbola with horizontal and vertical asymptotes . The curve ν = 0.05 is shown in Fig. 4 in orange.
![]() |
Fig. 4 Zeroth-order three-planet MMR loci in the period ratio plane. The colour indicates the index p + q of the resonance. The dashed lines corresponds to first-order two-planet MMRs (the oblique ones correspond to MMRs between planets 1 and 3). The curve ν = 0.05 is displayed in orange (see Eq. (22)). |
4.2 Single zeroth-order three-planet MMR Hamiltonian
Let us consider a specific resonance described by the integers p and q. One can make a linear change of variables to use explicitly the resonant angle and average over the non-resonant ones. Let us define
(24)
The conjugated momenta are
(25)
We call Γ the scaling parameter by analogy with the two-planet case (Michtchenko et al. 2008). Here, ϒ is the circular and coplanar angular momentum and verifies ϒ = G + C. Using the method described in Sect. 3, we can do a formal second-order averaging over θΓ and θϒ because these angles are not resonant. The Hamiltonian takes the form
(26)
where is the Hamiltonian of Eq. (17) averaged over θΓ and θG, and Θ represents all the actions defined in Eq. (25). In Eq. (26), we can neglect
as it is small with respect to
and only accounts for a correction of the mean motions of order ε. It should be noted that Γ and ϒ are integrals of motion of Eq. (26), up to terms of order O(ε3). As a result, the Hamiltonian only has one degree of freedom and is integrable. Another consequence of the conservation of ϒ is that the zeroth-order three-planet MMRs preserve the system AMD. In particular, if a system is only affected by these resonances, initially circular orbits will remain circular. As such behaviour is observed in numerical simulations before the late instability, this result confirms the decisive role of zeroth-order three-planet MMRs in driving the instability.
We consider small variations of the actions around the resonance. Let us denote Θ = Θ0 + dΘ where Θ0 corresponds to the value of Θ such that the system is on the resonance curve (Eq. (20)). Similarly, we have Λk = Λk,0 + dΛk. In turn, we have
(27)
and so the inner and outer planets are moving in the same direction while the middle planet is moving in the opposite direction. At first order, we can express the change in the period ratio d ν12 and d ν23 as a function of dΘ,
(28)
We can take the ratio of the small variations dν23 and dν12, and using Eq. (20) to replace p and q we obtain the differential equation
(29)
which gives the direction of the change of period ratios anywhere in the plane (ν12, ν23) due to the neighbouring resonance. We note that the equation no longer depends explicitly on the integers p and q. Indeed, while the strength of each resonance depends on the resonance index 2(p + q) (see the following section), the resonant motion direction can be extended as a continuous function of the period ratios using the resonance locus Eq. (20).
The solution of Eq. (29) gives the direction of the Chirikov diffusion if the system dynamics were entirely governed by the zeroth-order three-planet MMRs. The differential equation can be integrated numerically given an initial condition, and the solution for the system studied in Sect. 2 is displayed in orange in Fig. 3. We see that the system follows the diffusion direction for the majority of its lifetime, but leaves it as soon as the dynamics are no longer driven by the three-planet MMR network. We note that the problem has been reduced to study the diffusion along a one-dimensional curve.
4.3 Explicit size of the resonance width
From the previous section, we know that the dynamics of a single zeroth-order three-planet resonance are integrable. Provided that these resonances overlap, we also have seen along which curve the motion should take place. However, there is no guarantee that the neighbouring three-planet MMRs interact. If the resonances are well isolated because their width is smaller than their separation, there is no possibility for large-scale chaos. In this case, a system could be influenced by a single resonance and never jump to the other ones. The system will be almost secular and in principle could be considered as long-term stable. Moreover, the diffusion timescale is linked to the resonance strength. It is also possible that while the resonances are overlapped, the diffusion along the network is so slow that it is meaningless for astrophysical applications.
One therefore needs to study the dynamics in detail to evaluate the strength of the three-planet resonances. We therefore carry out in this section a detailed and quantitative derivation of the perturbative steps described above, keeping track of all the relevant terms which contribute to three-planet MMRs. We limit ourselves to a leading-order computation. As a result, we only keep the terms that do not depend on the eccentricity in the final expression. We also neglect the indirect term of the perturbing Hamiltonian as its value only affects the resonances when either p or q is equal to 1. It is instead necessary to keep terms to first order in eccentricity because they contribute to terms independently of the eccentricity at second order in mass. The terms of the perturbing Hamiltonian
that we consider are therefore (Laskar & Robutel 1995; Murray & Dermott 1999)
(30)
where c.c. designates the complex conjugate of the second sum, and
where αij = ai∕aj and are the Laplace coefficients. We refer to Appendix B for a definition and study of the Laplace coefficients and how to approximate them. Here we note that Quillen (2011) used a simplified approximation that confers the advantage of being easy to manipulate
(34)
In the limit of α very close to 1, the asymptotic equivalent of the Laplace coefficients does not depend on l (Laskar & Robutel 1995). However, we find that for a fixed α, for large |l|, the Laplace coefficient is asymptotic to
(35)
In the close planet approximation, α → 1−, we get
(36)
We refer to Appendix B for a detailed discussion.
The exponential dependency on −|l|(1 − α) of the Laplace coefficients has two consequences for . The interactions between planets 1 and 3 are always negligible with respect to the interaction in adjacent pairs, that is, for a given l,
. Similarly, for a given resonance, higher order terms in the resonant angle such as eιNθres for N > 1 are always negligible. It should also be noted that the formal development (Eq. (30)) is possible because the infinite sum can be replaced by a truncated sum in the actual computation due to the exponential decay of the coefficients. As a result, the neglected terms can be moved into the remainder in O (ε3).
With these clarifications, let us return to the calculation of the perturbative steps, starting from the perturbative Hamiltonian (Eq. (30)). The solution εχ1 to the corresponding homological equation has for expression8
(37)
where . One can note that because of the small denominators in both sums, the coordinate transformation is not valid (i.e. not close to the identity) in the vicinity of the co-orbital resonance (first sum) and the first-order two-planet MMRs (second sum). As explained schematically in the previous section, this perturbation step produces
terms, which we now calculate explicitly. In essence, we must only calculate the term
in Eq. (26). Since we would subsequently average over θϒ and θΓ, the only terms that remain in
must depend on the angle (pλ1 − (p + q)λ2 + qλ3) or its opposite. Because of the form of εχ1 and
, the only terms contributing at zeroth order in eccentricity to the averaged Hamiltonian
are of the form
(38)
or their complex conjugates. We note that in all the considered terms, the Poisson bracket can be reduced to the derivations with respect to the coordinates Λ2, λ2, and x2 of planet 2 only, as they are the only ones to appear on both sides. For the two last terms, we only keep the terms where the eccentricity dependency is dropped due to the derivation operator . This gives
(39)
Because is small with respect to the Keplerian part, we evaluate the action at the nominal resonance value, that is, pn1 − (p + q)n2 + qn3 = 0. Using Eqs. (18), (19), (31)–(33) as well as
(41)
where the Laplace coefficients depending on p (resp. q) are evaluated at α12 (resp. α23).
In this expression, the second prefactor can go to infinity for ν = 1∕(p + q). For the resonance defined by p and q, this happens at the intersection of the two-planet MMRs ν12 = p∕(p + 1) and ν23 = (q − 1)∕q. This result is a consequence of the non-validity of the second-order averaging very close to a two-planet MMR. Since we primarily focus on the regions outside of two-planet MMR, we ignore this feature in the following developments. Moreover, for large p + q, the MMR intersections are within the region of the two-planet MMR overlap.
Under the assumptions made so far, the above expression is exact, and can be used for numerical explorations of the size and typical frequency of each three-planet resonance without impediment (see below). To obtain further analytical insight, it is however quite cumbersome, and it does not clearly show which parameters of the planetary system and of the resonance play a role in determining the properties of the resonant motion. We therefore aim to simplify the above expression, keeping always in mind that we are ultimately interested in the diffusion in period ratio space driven by these three-planet resonances, and specifically in the timescale that is needed for large-scale diffusion. It is expected that the resonances with highest index dominate this timescale (see also below); in the remainder of this section we therefore take the limit 1∕(p + q) → 0 and expand around this value. We note moreover that for 1∕(p + q) → 0, the second term in Eq. (42) blows up when ν → 0, in which case also the first term would go to infinity; however ν → 0 only happens when one of the period ratios νi,i+1 ≃ 1: this limit is beyond the scope of the study, and so we can exclude this case.
With these considerations in mind, the above expression can be considerably simplified. To this end, we make the close planet approximation: 1 − αij ≪ 1. We define
(43)
which is an excellent estimate for period ratios from 0.5 to 1.
The product of the Laplace coefficients and their derivatives in Rpq introduces an exponential factor of the form e−pδ12−qδ23 (cf. Eq. (36)) that sets the order of magnitude of the resonance term. We can therefore simplify expression (42) by taking advantage of the resonance relationship. Indeed, for tightly packed systems, and in the vicinity of a resonance defined by p and q, Eq. (20) can be transformed into a relationship on the planet spacings:
(44)
By analogy with the generalised period ratio separation ν, we define a generalised orbital spacing that we note
(45)
where the two last equalities are approximations using Eq. (44). We have
(46)
Using the newly defined variable δ and η, the expression of the coefficient ε2Rpq can be simplified to
(47)
where we only keep the terms up to second order in (p + q)δ. Since the exponential factor depends on (p + q)δ, the resonance mainly matters in the region where (p + q)δ is of order unity, hence we choose to set 17 + 21∕((p + q)δ)) ≃ 38, which is approximately the value taken for (p + q)δ ≃ 1 as it allows us to carry the computations analytically. This value also gives a more accurate estimate for the width of the resonances (see below).
The expression obtained in Eq. (47) for the three-planet resonance perturbation Hamiltonian is remarkable. The resonance strength only depends on its index and not explicitly on p and q. This means that all resonances with the same index can be compared very easily. In other words, the network of zeroth-order three-planet MMRs can be partitioned into subnetworks consisting of resonances with the same index.
To fully describe the resonant dynamics, we now go back to Eq. (26) (we recall that we can safely drop the term ); we now expand the Keplerian part around the resonance centre (Chirikov 1979; Petit et al. 2017). This is more easily done in the original Delaunay variables Λ, and we have
(48)
where the constant terms can be safely dropped. Using Eqs. (18) and (27), the first-order term vanishes and the coefficient of the second-order term has for expression
(49)
We note that only depends on the index of the resonance and weakly on the planet masses. Passing finally to the resonant canonical variables (Eqs. (24) and (25)), the resonant Hamiltonian has for expression
(50)
This is the standard pendulum Hamiltonian. The width of the resonance in the action space is given by the expression (e.g. Ferraz-Mello 2007)
(51)
The small oscillation frequency is given by
(52)
is the relevant mass ratio for the studied problem. For equal mass and equal tight spacing we have . Contrary to the two-planet case, we note that the relevant mass ratio is not reduced to the sum of the planet masses over the mass of the star. It is also interesting to note that the expression remains meaningful even in the case where one of the planets is reduced to a test particle.
The resonances have a clearer geometrical interpretation in the period ratio space than in the action space, particularly when one needs to compare them. We therefore compute the width of the resonances perpendicularly to the network, that is, the width in terms of the variable η. Using Eqs. (21) and (28) and some algebraic manipulation, we have
(54)
This means that the width in terms of η can be estimated as
(55)
We have thus shown that the width of the resonances in the period ratio plane depends exponentially on the MMR index and the prefactor is a continuous function of the period ratios. In particular, it seems important to compare resonances with the same index because of their similar geometry.
4.4 Resonance overlap
We wish to determine the sections of the period ratio space (ν12, ν23) where resonances overlap. Because of the expression of resonance width Δηpq, we see that the width of the resonances close to a given point (ν12, ν23) mainly depends on the index p + q. It is therefore natural to consider the density of the resonances for a fixed value of p + q.
We denote ρk(δ, η), the local filling factor of the zeroth-order three-planet MMRs of index k = p+ q. Here, ρk corresponds to the proportion of the period ratio space occupied by this subnetwork. Let us also define ρtot (δ, η), the filling factor of all zeroth-order three-planet MMRs. The filling factor measures the space locally9 occupied by all the nearby resonances of arbitrary index with respect to the available space in the period ratio plane. If ρtot is larger than 1, then there are enough resonances to locally cover the period ratio plane.
Such a filling factor is introduced by Quillen (2011) for the same problem. However, these latter authors only consider resonances such that |p − q|≃ 1, which leads to the exponential dependency being neglected. We show here that taking into account all the resonances is critical to obtain an accurate diffusion rate and survival time. The idea to count all the resonances without taking care of their precise position in order to obtain Chirikov’s overlap estimate was also used with success for two-planet MMRs of arbitrary order (Hadden & Lithwick 2018).
We have ρtot ≤∑kρk since some resonances are counted multiple times. Indeed, if p and q are not coprime, the resonance lies on top of a resonance of lower index10. Nevertheless, the contribution of a resonance defined by two integers of the form Np, Nq is negligible with respect to the contribution of the resonance p, q because of the exponential decrease. As a result, we consider that the overall resonance filling factor ρtot, is the sum of the subnetwork filling factors ρp+q.
Let us consider the subnetwork of resonances with index11 k = p + q. The distance between two resonances in terms of η is constant. Indeed, let us consider the resonance defined by integers p and q; its upper neighbour is defined by the integers p + 1 and q − 1, hence
(56)
The filling factor ρk for the subnetwork of resonances with index k can be determined by taking the ratio of the resonance width in terms of η with the distance between two neighbouring resonances in η, that is,
(57)
The filling factor ρk thus depends on the subnetwork index k, the generalised orbit spacing δ, the masses, and the resonant locator η.
We approximate the total resonance filling factor ρtot as the sum of the subnetwork ones. We thus have
(58)
where we have replaced the sum by an integral. The computations are also possible using the infinite sums but they result in a more complicated expression with a very limited gain in accuracy. This approximation is also done by Quillen (2011).
As in Quillen (2011), we find that the filling factor depends linearly on the mass ratio and scales as δ−4. However, our expression is valid for an arbitrary spacing and mass distribution, as long as the system is tightly packed. We confirm that the natural spacing rescaling for the problem is not the Hill radius (Eq. (2)), which scales as ε1∕3, but rather a dependency on ε1∕4. In particular, assuming that M and η are constant, we can define a critical spacing value δov such that the zeroth-order three-planet MMR network fills the entire space. Taking ρtot = 1 and solving for δ, one obtains
(59)
Here, δov is a function of the masses and η. We can rewrite the filling factor ρtot as a function of δov as a power law over δ
(60)
In the case of equal mass and spacing systems, Eq. (59) becomes
(61)
We note that δov,eq corresponds to the generalised spacing defined in Eq. (45). The actual orbit spacing is equal to 2δov,eq in this case.The overlap criterion obtained by Quillen (2011) is similar to ours since the exponent 1∕4 makes the numerical factors almost equal.
We plotin Fig. 5, the number of resonances that overlap at a given point in the plane (ν12, ν23) for three equal-mass planets. The mass of each planet is 10−5 M⊙. The image is computed by creating a square grid of 2000 equally spaced period ratios between 0.7 and 1. For each resonance index between 2 and 200, we compute for each point the closest resonance in terms of η. The closest resonance indicator is defined by ηres. The closest point on the resonance, (ν12,res, ν23,res), is found by gradient descent for the function η − ηres. We then compute the width in terms of η using Eq. (55) at the point (ν12,res, ν23,res) and compare it to the distance to the resonance η − ηres. We use the exact expression for Rpq (Eq. (42)). We count the resonances with multiplicities, that is, even if p and q are not coprime. Wesee at the vicinity of the two-planet MMRs that the width of the three-planet resonances increase due to the second term in Eq. (62).
The number of resonances is to first order a proxy for the filling factor ρtot (Eq. (58)). We see in Fig. 5 that the region where the overlap of the three-planet MMR network takes place extends well beyond the Hill-stability limits (Eq. (1)), particularly for comparable spacings between the two neighbouring planet pairs. However, for very unequal spacings (away from the main lower-left to upper-right diagonal) we see that the overlap of only the three-planet MMRs is not enough to account for the instability and the two-planet interactions should be taken into account for the initial diffusion process.
To quantify how far the overlapping region extends, we consider systems of equally spaced planets with equal masses mp and plot in Fig. 6 the minimal spacing given by the Hill-stability limit (Gladman 1993; Petit et al. 2018), the Wisdom (1980) two-planet MMR overlap criterion, and our three-planet MMR overlap criterion (Eq. (61)) as a function of the mutual Hill radius (Eq. (2)). We see that for small masses, the three-planet MMR overlap region goes to orbital spacing of the order 10 Hill radii, comparable to what was observed in previous numerical simulations. It is also worth noting that we are only considering a restricted part of the resonance network. Higher order three-planet resonances can also contribute to the filling factor even for circular orbits. Moreover, for the larger masses, averaging to the second order in the masses as done in Sect. 3 may not be sufficient. Moreover, chaotic diffusion occurs in general well before the full overlap that is computed here (Chirikov 1979; Lichtenberg & Lieberman 1992). Therefore, the phenomenon studied here could also work beyond the predicted limit. Our criterion should be seen as a lower limit.
Nevertheless, we predict that beyond a certain spacing, we should see an increase of the survival time approximately at the limit where the three-planet MMR network is not fully overlapped. The limit should also appear at smaller separations in terms of Hill radius for larger planets (see Sect. 6).
We have shown that there is a region where three-planet MMRs can contribute to a non-secular evolution, resulting in a diffusion in the semi-major axes of the planets, ultimately leading to the instability of the system when a first-order MMR is crossed. The region can be determined quite accurately with the introduction of adapted variables and the computation of a resonance filling factor ρtot (Eq. (58)). We have seen that the resonance index is critical for determining the resonance width and that each equal-index subnetwork can be considered individually. The minimum spacing between planets such that the resonance network does not cover the period ratio space scales as ε1∕4 and not as the Hill radius. However, we see that in some regions, the filling factor is well above 1, that is, only the smaller index MMRs are necessary to cover the space. As a consequence, the diffusion is faster when only wider resonances are involved, which will lead to the observed differences in Tsurv.
![]() |
Fig. 5 Number of resonances covering the period ratios space for equal-mass planets, with masses 10−5 M⊙. At first order, the number of resonances can be compared to the filling factor ρtot (Eq. (58)). We plot the two-planet circular Hill-stability limits (Petit et al. 2018) for both planet pairs in green and the predicted overlap limit for the three-planet MMRs (Eq. (59)) in blue. |
![]() |
Fig. 6 Comparison of the two-planet stability criteria, the Hill stability (Gladman 1993; Petit et al. 2018), the first-order MMR overlap criterion, and the zeroth-order three-planet MMR overlap criterion expressed in units of the mutual Hill radius as a function of the planet-to-star-mass ratio for equal masses and equally spaced planets. |
5 Diffusion timescale
5.1 Chirikov’s diffusion
In the Chirikov model, a large-scale diffusion of the actions (or the frequencies) occurs when the resonances overlap. Qualitatively, perturbations to the integrable resonance Hamiltonian (50) create a stochastic layer in the vicinity of the separatrix. When overlap occurs, the stochastic layers of adjacent resonances merge which allows diffusion along the network. The diffusion rate depends on the resonance width and the period of the resonance. One can estimate the diffusion rate for the resonance locator η in the vicinity of a resonance as (Chirikov 1979)
(62)
This expression is an estimate that is valid12 when resonances overlap and the stochastic layers are well connected. Studies of the Chirikov diffusion have been carried out for simplified Hamiltonian (Giordano & Cincotta 2004) or on astrophysical problems (Cachucho et al. 2010). Cincotta (2002) discusses the link between Chirikov and Arnold diffusion in astronomy and presents a modern description of Chirikov’s theory. The link with Nekhoroshev (1977) theory is proposed in Cincotta et al. (2014). If the space is fully covered by overlapping resonances, the trajectory can be well approximated by a random walk.
The diffusion direction is perpendicular to the resonance in the action space. This means that the diffusion is not isotropic, and to study the trajectory one would need to compute the contribution of every resonance to the diffusion tensor at every point. If the resonance lines do not intersect in the considered region, the diffusion will be well approximated by an unidimensional random walk perpendicular to the resonance network, with a negligible diffusion parallel to the resonances. We can therefore consider a scalar diffusion coefficient given by the specific resonance that dominates the dynamics around the position in the phase space. In particular, the diffusion coefficient is not constant and depends on the closest resonance width. Such a diffusion corresponds to the behaviour observed in Fig. 3.
In Sect. 4.4, we compute an overlap criterion taking into account every zeroth-order three-planet MMR. However, as the resonance index increases, the associated diffusion rate vanishes, such that in the limit where the diffusion is dominated by smaller and smaller resonances, the timescale effectively tends to infinity. However, for δ < δov, not all the resonances are necessary to cover the phase space. We therefore only need to consider the largest ones to compute the survival time.
5.2 Partial resonance overlap
We consider a small region around a point (δ, η) where the zeroth-order three-planet MMR network is locally overlapped. Since ρtot > 1, not all the resonances are necessary to cover the phase space. The largest resonances lead to the fastest diffusion, and therefore we need to only consider the subset of the widest resonances needed to cover the period ratiospace in this given region. As the distance to reach a two-planet MMR is small, the main difference between the size of the resonances is governed by their index. We therefore define an overlap index kov such that the subnetwork composed of the three-planet MMR with index k smaller than kov locally covers the space. Using Eqs. (57) and (60) we have
(63)
We thus have an implicit definition of kov. We also note that the equation depends on kovδ rather than kov alone. As a result, we define the variable
(64)
which is a function of δ∕δov. There is no solution for Eq. (63) in terms of elementary functions. However, an explicit solution can be obtained using the Lambert W function (Corless et al. 1996, the function is also called inverse product log)
(65)
where W−1 is the real valued branch with values smaller than −e defined between − 1 and 0. The function W is the inverse function to z →zez. Using bounds on W−1 by Chatzigeorgiou (2013), an excellent approximation of ξov is
(66)
5.3 Instability timescale
Starting from a point (δ, η) in the period ratio space, we assume that the system wanders along the diffusion direction computed in Sect. 4.2. While not exactly perpendicular to the resonant network in the period ratio space, the motion along the diffusion direction can be well parameterised by η. We therefore monitor the diffusion in terms of η because the resonance width and densities are easy to compute in terms of this variable. Moreover, as seen in Fig. 3, the systems are not too far from the first-order MMRs, and so the distances to cover are short and we can consider the period ratio as almost constant along the trajectory.
Since the diffusion coefficient Dpq mainly depends on p + q, one can associate a diffusion rate to each of the equal resonance index subnetworks. We note ηd, the distance in terms of η to the closest first-order MMR. In order to describe the diffusion process, we adapt the framework developed by Morbidelli & Vergassola (1997) to compute the escape rate of particles from the vicinity of invariant tori. Let us assume that the system starts initially at a position η0, and becomes unstable once η reaches η0 + ηd. Furthermore, we assume that the dynamics behave as a random walk. The position of the system along the resonance network is described by the equation
(67)
where s is related to the local diffusion coefficient and b(t) is a Gaussian white noise with zero average, verifying
(68)
For s = s0 constant, Eq. (67) is the classical Langevin equation. The associated diffusion equation for the probability density p is
(69)
where the diffusion coefficient is . By analogy with the case where s is constant, we define
, where p and q define the largest MMR that contains the point η. Since in the region considered, the resonances overlap, the interval (η0, η0 + ηd) can be partitioned into the different subnetworks. Each value of η can be associated with an index k that corresponds to the widest resonance that contains it. The probability that a given point η is in a resonance of index k is given by
(70)
We note that the considered value of η could also be containedin a higher index resonance than kov. However, the contribution of this higher order resonance to the diffusion is negligible.
Following Morbidelli & Vergassola (1997), Eq. (67) can be solved by introducing the variable
(71)
Indeed, computing the time derivative of y using Eqs. (67) and (71), we have
(72)
which is the Langevin equation with a unit diffusion coefficient. The evolution of y is therefore known and we obtain the evolution of η by inverting Eq. (71). To do so, we need to determine the value of s(η).
Since we are interested in the overall diffusion speed and not the exact diffusion at a given point, we can attribute a probabilistic value to s using Eq. (70). We thus compute the average value of y as a function of η over all the configurations of the resonance network. Noting ȳ and the average values, we have
(73)
We see that the variation of ȳ is proportionalto the variation of , the integral being almost constant13 around a given point of the period ratio space. In particular, we can define an effective diffusion coefficient, taking into account the contribution of all the resonances necessary to locally cover the phase space
(74)
Indeed, deriving Eq. (73) gives
(75)
We refer to Appendix C.1 for the exact expression, which involves several special functions to evaluate the integrals. It is nevertheless straightforward to evaluate the integrals numerically using scipy (Virtanen et al. 2020) for instance. We also provide a jupyter notebook reproducing the figures of this article14.
For the remainder of the discussion, a very good estimate15 is obtained with the expression
(76)
This expression is surprisingly compact, mainly depends on δ and δov, and does not contain an exponential term. This expression is within an order of magnitude of the exact one over the majority of the considered range of δ. We highlight the fact that the diffusion coefficient goes to zero for δ → δov.
We have reduced the initial problem to a simple continuous unidimensional random walk with constant diffusion coefficient. The system can wander up until it reaches a first-order two-planet MMR on the diffusion direction. The survival time distribution is thus given by the well-studied hitting time probability of a Brownian motion with two absorbing boundary conditions (Redner 2001; Borodin & Salminen 2002; Wax & Cohen 2009). We note η− and η+ (η+ > η−), the value of η at the neighbouring first-order two-planet MMRs along the diffusion direction. The interval where the system can wander has for length Δη = η+ − η−. The initial position on this interval can be measured by the quantity u0 = (η0 − η−)∕Δη that is between 0 and 1.
The time can be rescaled such that the considered segment has length unity and the diffusion coefficient is 1 if we take as time unit
(77)
The distribution of log10Tsurv∕T0 is plotted in Fig. 7 as well as the mean and standard deviation of log10Tsurv∕T0. We give the expression of the distribution in Appendix C. As shown in Fig. 7, the mean of log10Tsurv∕P1 can be well approximated as
(78)
where we choose to normalise the survival time Tsurv by the innermost planet period P1 to include T0 in the right-hand side expression. Similarly, we approximate the standard deviation as a second-order polynomial as
(79)
We note that the standard deviation does not depend on the value of log10T0∕P1 itself. Thus the standard deviation does not depend on the order of magnitude of the instability time. This is a remarkable result, as it has been shown in numerical simulations (Hussain & Tamayo 2020) that the standard deviation of the survival time of extremely close initial conditions has the same properties. Besides, Hussain & Tamayo (2020) measured the standard deviation to be 0.43 ± 0.16 dex which isconsistent with the value we obtain for initial conditions not too close initially to two-planets MMR, as can be seen in Fig. 7.
Equation (78) gives the expression of the mean survival time for any initial condition and while it involves terms that are not easily tractable analytically, they are easy to compute numerically given a specific system. In the remaining part of this section, we seek to obtain a simplified expression to show how the mean survival time depends on the spacing and planet masses.
Let us assume that the system initially starts in between the two-planets MMRs P:P–1 and P+1:P for the pair 1-2 and Q: Q–1 and Q+1:Q for the pair2-3. In the example shown in Sect. 2, P = Q = 6. By evaluating η at the edges of the square created by the resonances, one can find that the maximum variation of η without encountering a MMR is
(80)
where the last equality is true for ν evaluated at the intersections of the MMRs P:P−1 and Q+1:Q or P+1:P and Q:Q−1. In practice, the variation of ν in a given resonance rectangle is 2ν2 ≪ ν and so we neglect the variations of ν in the rectangle. Thus, using ν as the characteristic length for the diffusion interval of η is accurate. Here, ηmax describes the largest possible variation of η, hence for any given point (ν12, ν23). For our simplified expression, we take Δη = ηmax ≃ 3∕2δ. Furthermore, we assume u0 = 0.5. Using Eq. (76) to estimate the effective diffusion coefficient, we can get an order of magnitude for the survival time with the expression
(81)
This expression can be decomposed into a prefactor that mainly depends on the planet-to-star-mass ratio and a function that only depends on δ∕δov. At first glance, this expression isnot linear in δ which seemsto contradict the numerical results (e.g. Chambers et al. 1996; Obertas et al. 2017), which hinted at a linear dependency of the logarithm of the survival time on the orbital spacing. However, the function has an inflection point for δ = 0.629δov and can be well approximated by a linear function of δ∕δov. The linear approximation is correct in the regime of interest, that is, for Hill-stable planet pairs not too close to the overlap limit. We have
(82)
To compare with previous numerical studies (Chambers et al. 1996; Faber & Quillen 2007), we compute the estimated survival time for equal-mass and equally spaced planets
(83)
The slope coefficient 3.56 is very close to values obtained in previous works: Faber & Quillen (2007) estimated it at 3.7 and Yalinewich & Petrovich (2020) at 3.4. The prefactor proportional to m0 ∕mp is consistent with numerical simulations (Chambers et al. 1996; Faber & Quillen 2007) and the numerical constant is very close to the one obtained by Faber & Quillen (2007).
To summarise, we estimate the diffusion rate along the zeroth-order three-planet MMR network by only considering the widest resonances, up to the index k = p + q such that they locally cover the period ratio plane. We show that the complex random walk along the resonance network can be represented by a diffusion process with an effective locally constant diffusion coefficient given by Eq. (74). As observed in numericalsimulations (Hussain & Tamayo 2020), we show that the survival time distribution is approximately log-normal and we recover the same standard deviation. Our estimation of the mean survival timescales as the planet separation in units of ε1∕4 and not in units of Hill radii. In particular, we emphasise the importance of considering systems of various masses in such studies as it allows to discriminate between the physical mechanisms driving the dynamics. Moreover, while our estimate is not exponential in planet spacing as fitted in numerical simulations, we show that for the range of times of interest, it can be considered as such. As Quillen (2011), we predict that beyond the overlap limit, the survival time is likely much larger since the Chirikov diffusion is not an efficient process on its own.
![]() |
Fig. 7 Top panel: probability distribution function of log10Tsurv∕T0 for various values of u0, the normalised distance of the initial condition of the system to the two neighbouring first-order two-planet MMRs. We see that for u0 = 0.5, the time distribution is log-normal as observed in Hussain & Tamayo (2020). Bottom panels: mean and standard deviation of the distribution as a function of u0. We also plotapproximate fits to the curves. We note that the mean remains within an order of magnitude of its maximum for almost all values of u0. For u0 close to 0.5, the standard deviation is close to the value 0.43 ± 0.16 measured byHussain & Tamayo (2020), plotted as a green line with its 1σ error. |
6 Comparison with numerical simulations
A large number of numerical studies have recently been performed on the problem of instability of tightly spaced planets (e.g. Obertas et al. 2017). However, the most recent study limited itself to the minimal setup: an equal-mass and spacing (EMS) system with three planets on initially circular and coplanar orbits was performed by Faber & Quillen (2007). While these latter authors considered different mass ratios, the integration time was limited to 106 inner planetorbits with a limited number of points.
6.1 Equally spaced, equal-mass systems
In order to have a fine enough resolution and a longer integration time, we run our own suite of numerical simulations. We use REBOUND (Rein & Liu 2012) and the symplectic integrators WHFAST (Rein & Tamayo 2015) and MERCURIUS (Rein et al. 2019). We initialise three-planet systems on initially circular and coplanar orbits. As in previous studies (e.g. Chambers et al. 1996), the innermost planet semi-major axis is set to 1 au and the two outer-planet semi-major axes are chosen such that the two-planet pairs have an equal semi-major axis ratio. The initial angles are randomly drawn. We do three suites of simulations with planet masses of 10−7 M⊙, 10−5 M⊙, and 10−3 M⊙. In each suite, the range of period ratios considered starts at 2 mutual Hill radii, that is, in the region where even the planet pairs should be considered unstable, and extends beyond the predicted three-planet MMR overlap criterion derived in Sect. 4.4. Each system is integrated until two planets are closer than 1 Hill radius or until the inner planet has performed 109 orbits. For each of these system, we report the final time as the survival time. A value of 1 Gyr should therefore be understood as a lower limit. We use WHFAST and a time-step of one-twentieth (1/20) of the inner orbit for the two smaller mass suites and MERCURIUS with a time-step of one-thirtieth (1/30) of the inner orbit for the largest mass planets. We stopped running simulations at larger separations when it became clear that the systems were stable for 1 Gyr, outside of the first-order two-planet MMRs. Each of the sets of simulations is composed of between 1600 and 1700 systems.
We respectively plot in Figs. 8–10, the survival times for the sets of simulations with planet masses of 10−7 M⊙, 10−5 M⊙, and 10−3 M⊙ and the associated analytical prediction (Eq. (81)), as a function of the initial period ratio. Additionally, we add various other features that help us to understand the patterns that emerge in the survival time curves. We plot the Hill stability limit (Petit et al. 2018), the two-planet first-order MMR overlap criterion (Wisdom 1980), the nominal position of two-planet first-order MMRs p + 1 : p and of two-planet second-order MMRs p + 2 : p. We also plot an estimate of the width of first-order MMR for initially circular orbits (Petit et al. 2017).
The first thing to be noticed from these figures is that the survival time estimate (Eq. (81)) of the logarithm of the survival time is consistent with the numerical simulations in the range where the former is almost linear. The agreement is particularly striking for the intermediate case (10−5 M⊙). We discuss explanations for the discrepancies for the low-mass and Jupiter-mass planets below. Moreover, the slope being correct in all three figures is another indication that the scaling in ε1∕4 is more appropriate than renormalising the spacing by the Hill radius.
We focus on Fig. 9 to describe more precisely the different features that should be pointed out. First, ifone ignores the variations due to the proximity to the resonances, our estimate lies in the middle of the distribution of survival time from period ratios of 1.05 to the end of the region where we consider that the three-planet MMRs overlap, around 1.195. This means that our criterion slightly underestimates the diffusion time, as in Eq. (81) we assume a maximum value for the distance to the two-planet MMR network. Nevertheless, Chirikov diffusion correctly predicts the slope as well as the right order of magnitude for the survival time.
We then notice, as previous studies have (Smith & Lissauer 2009; Obertas et al. 2017), that the substructure on the curve is very well explained by the two-planet first-order MMR. By considering the regions outside of first-order MMR and the regions inside separately, we can see that the survival time when the dynamics are dominated by the two-planet MMR is roughly two to three orders of magnitude lower than outside of their influence. We can explain why the survival time is shorter for period ratios just below the Keplerian resonance, by noting that due to the shape of thefirst-order MMR the unstable fixed point where the first-order MMR separatrices originate is situated on the left of the resonance in the figures. We note that if we had introduced a fluctuation in the initial period ratio, the pattern would be much less clear. We can also note that the second-order MMRs likely play a role in accelerating the diffusion. A similar effect is also clearly observed regarding the 2:1 MMR in Fig. 10as well as in Fig. 8 where the larger apparent spread is due the very dense two-planet MMR network.
Finally, we focus on the region close to where we predict that the three-planet MMR stops overlapping (period ratios of about 1.19) in Fig. 9. We see that the linear trend followed in the range 1.05 to 1.16 no longer holds due to some systems surviving longer than expected, in particular beyond 1 Gyr. Moreover, the spread of the survival times increases instead of staying constant as shown by Hussain & Tamayo (2020). One can note that a similar feature is also visible in the results of Obertas et al. (2017), although the increase in the spread is less visible therein, most likely because these latter authors consider five-planet systems instead of three. In particular, outside of the 6:5 and the 11:9 resonance, it appears that systems live much longer than one would have expected from extrapolating the linear trend fitted in previous numerical studies. The same behaviour is also observed in Fig. 10. However, the region where short-lived and long-lived systems coexist is much larger because of the larger size of high-order two-planet MMRs that are not taken into account in this analysis.
These two observations, namely the longer survival times and the increased spreading around the overlap limit, are consequences of our analytical derivation. Indeed, beyond the overlap region, the Chirikov diffusion alone can no longer drive the instability over homogeneous regions of the phase space. This does not mean that systems beyond our approximate overlap limit will live indefinitely. However, the instability in these systems is driven by an alternative mechanism to the Chirikov diffusion considered here. In particular, we have neglected the diffusion parallel to the resonances, that is the Arnold diffusion (Cincotta 2002). However, these alternate pathways to instability are most likely much slower.
We now look more in detail at the apparent discrepancies between the estimated survival time and the two extreme cases. In Fig. 8, we see that while the analytical curve shows good agreement beyond period ratios of 1.035, the numerical simulations have a much shorter survival time at very close separations. We postulate that two unaccounted-for effects play a significant role in this regime. First, for very small bodies, the two-planet MMR overlap limit extends beyond the Hill-stability limit (Deck et al. 2013; Petit et al. 2018). In particular, Chirikov diffusion is also possible along the two-planet MMR network. Since the first-order two-planet MMRs are of comparable size, which depends on roughly ε2∕3 for circular orbits, this network has a faster diffusion timescale. Besides, in the derivation of the resonant coefficient expression in Sect. 4.3, we neglect the interactions between the inner and the outermost planets. While these interactions are orders of magnitude smaller for larger masses and separations, they may be taken into account for very close planets. However, the region where these long-distance interactions matter is of little interest and this effect is most likely smaller than the already mentioned two-planet MMR overlap. We also notice in Fig. 8 that we are not able to resolve the difference of behaviour in the non-overlapping region. Indeed, the survival time depends on ε−1 and as a result the transition away from the Chirikov regime occurs for times of about 1010 orbits.
In the case of Jupiter-mass planets (10−3 M⊙), we see that the survival times starts to spread from 104 orbits to more than 109 around a period ratio of 1.7 whereas we estimate the overlap region to extend up to 1.85. Moreover, our analytical estimate is too low in the region where the spread of survival time is consistent with a uniform diffusion process (up until period ratio of 1.7). We should first point out that the analytical results in this mass regime are at most an extrapolation. Indeed, the perturbation theory used in Sect. 3 diverges for large perturbations. Moreover, the close planet approximation is no longer valid for large spacing. In particular, in our approximation of the Laplace coefficients (Eqs. (35) and (36)), we replaced αl with e−l(1−α) to carry the computations analytically. Since 1 − α < −ln(α), we overestimate the resonance width and thus the extent of the overlap region. In order to get a closer estimate for the particular case of equally spaced planets, one can replace 1 − α with − ln(α) in the expression of δ (Eq. (45)) and then in the expression of the resonance density (Eq. (58)) in order to compute the period ratio where overlap occurs. After some algebraic manipulation, we estimate the actual limit to the overlap region to be situated at
(84)
We note that we still rely on formulas obtained with the close planet approximation and as such, this result remains an attempt to understand a discrepancy between the analytical curve and the numerical simulations. Figures 9 and 10 show the estimate (Eq. (84)) as a red dotted line. We see that they lie almost exactly where some systems start surviving beyond 109 orbits.
![]() |
Fig. 8 Survival time for a three-planet EMS system for planet masses of 10−7 M⊙ as a functionof the initial period ratio. The red curve corresponds to the survival time estimate (Eq. (81)), the blue vertical line to the Hill-stability limit (Petit et al. 2018), the orange vertical line to the two-planet MMR overlap criterion (Wisdom 1980). The dashed black (resp. green) lines are the two-planet first(resp. second)-order MMR. The light orange rectangles show an estimate of the width of the two-planet MMR (Petit et al. 2017). |
![]() |
Fig. 9 Survival time for a three EMS planet system for planet masses of 10−5 M⊙ as a function of the initial period ratio. See Fig. 8 for a detailed caption. The red dotted line corresponds to an alternate estimate position for the limit of the overlapped region (see text and Eq. (84)). |
![]() |
Fig. 10 Survival time for a three EMS planet system for planet masses of 10−3 M⊙ as a function of the initial period ratio. See Fig. 8 for a detailed caption. The red dotted line corresponds to an alternate estimate position for the limit of the overlapped region (see text and Eq. (84)). |
6.2 A non-equally spaced case
In order to demonstrate that our result is valid beyond the EMS case, we show an example of a configuration where the spacing between the inner and the outer pairs is different. We initialise systems of three planets of equal mass (10−5 M⊙) on initially circular and coplanar orbits. For a given period ratio ν12 for the inner pair, we set the outer period ratio to
(85)
where is the golden ratio and was chosen to avoid going through the intersections of two-planet MMRs, reducing their importance into the numerical results. This leads to an outer pair spacing in period ratio that is ϕ times larger than for the inner pair. The rest of the setup is similar to the EMS cases.
The survival time as a function of the initial inner period ratio is plotted in Fig. 11. The dashed red curve corresponds to the survival time for the EMS configuration shown in Fig. 9. We see that our model also successfully predicts the survival time in this case. In particular, we show that the EMS estimate underestimates the survival time by roughly an order of magnitude. The overlap limit is also reached at a tighter separation for the inner pair than in the equal spacing case. Because there is no configuration where both planets are initialised close to a two-planet MMR, we observe less features in the numerical results.
![]() |
Fig. 11 Survival time for non-equally spaced systems described in Sect. 6.2 as a function of the initial inner-pair period ratio. The outer-pair period ratio is given at the top of the figure. The dashed survival-time estimate corresponds to the EMS case. See Fig. 8 for a detailed caption. |
6.3 Application to non-equal-mass planets
The relevant mass ratio introduced in our model, εM (Eq. (53)), shows little variation in the case where one body has a significantly smaller mass than the two others. Indeed, its value is maximised for three equal-mass planets for a given maximum planet mass mp. In this case, for equal spacing we get εM ≃ 1.22mp∕m0. As already noted, εM is non-zero as long as two planets are massive. Assuming two equal-mass planets and a test particle, εM = mp∕m0 if the test particle is in the middle and scales as ηmp∕m0 (resp. (1 − η)mp∕m0) if the test particle is the inner (resp. the outer) body. In this case, εM = mp∕m0 varies by less than a factor three between these extreme cases for initially equally spaced bodies. To confirm this numerically, we plot on Figs. 12 and 13 the survival time of a test particle and two equal-mass planets in the case where the test particle is the inner and the middle body, respectively. The massive planets have a mass of 10−5 M⊙ and the remaining setup is similar to the one described in Sect. 6. We plot with a dashed line the predicted time if the three planets were of equal mass.
The agreement is excellent in the case where the test particle is placed between the two planets (Fig. 13). We recover the same features that were discussed in Sect. 6, including the corrected overlap limit (Eq. (84)). We also note that the predicted timescale is slightly longer than in the EMS case and this point seems confirmed by the numerical simulations. When the inner planet is replaced by a test particle (Fig. 12), we note that the estimate remains good before the predicted overlap but it seems that the overlap limit is underestimated in this case. We hypothesise that taking into account the interactions between the inner and the outer body may solve this discrepancy.
These numerical simulations confirm that diffusion along the three-planet MMR network is the main mechanism driving the instability for tightly packed systems. In particular we highlight the very strong change of behaviour occurring at the limit of the fully overlapped region. Moreover, the survival time estimate given by Chirikov diffusion accurately predicts the numerical simulation results over a wide range of survival times and the relevant range in orbital separations.
![]() |
Fig. 12 Survival time for a system composed of an inner test particle and two planets of equal mass (10−5 M⊙) with equal spacing as a function of the initial period ratio. See Fig. 8 for a detailed caption. The full red curve is theprediction using the survival time estimate (Eq. (81)) and the dashed line corresponds to the same estimate for the EMS case shown in Fig. 9. The red dotted line corresponds to an alternate estimate position for the limit of the overlapped region (see text and Eq. (84)). |
7 Beyond three planets on circular orbits
7.1 Systems with four or more planets
As already seen in previous numerical studies, increasing the number of planets beyond three does not fundamentally change the survival timescale. Chambers et al. (1996) show that while there is a slight change in the slope of the survival time between systems of five planets and systems of three, the survival time is mostly unchanged by the addition of other planets into the systems16.
It is therefore natural to try to extend our results to systems composed of more than three planets. Unfortunately, contrary to the three-planet case, it is not possible to reduce the dynamics to a unidimensional Chirikov diffusion. Indeed, the resonance network cannot be projected into a two-dimensional plane as done in Sect. 4. One solution is to consider triplets of adjacent planets, and assume that this triplet is perturbed by the additional planets. The influence of the other planets can be seen as a change in the period ratio ν12 and ν23 due to the resonances with the adjacent planets. Assuming that the planet spacings are comparable, the perturbation of the period ratio is of the same order of magnitude as the one induced by the three-planet MMR from the considered triplet.
As a result, we can modify the resonance density ρk (Eq. (57)) by including a multiplicative factor K representing the influence of the other planets. This is similar to assuming that the network is composed of K times more resonances than previously accounted for. The planets are mainly influenced by their direct neighbours. Taking a conservative approach, we consider that both the inner and outer neighbours of the triplet increase the number of resonances affecting the three-planet subsystem by 50%. We therefore take K = 2 as a reasonable guess. The overlap spacing of a system of five or more planets is given by K1∕4δov. The survival time is also affected because, while ρk has changed, the resonance width has not. As a result, the effective diffusion coefficient (Eq. (74)) becomes K−2Deff, where the change of δov should be accounted for.
We can compare this expression to previous numerical results. We choose the simulations from Obertas et al. (2017) as they run systems composed of five Earth-mass planets with high resolution in terms of period ratio. Moreover, their simulations have been run up to 1010 orbits. We plot in Fig. 14, the survival time from Obertas et al. (2017) for five-Earth-mass-planet EMS systems as a function of the initial period ratio as well as the three-planet survival time estimate (Eq. (81)) and its extrapolation to a five-planet system. As expected, the three-planet theoretical survival time slightly overestimates the computed survival time, and more importantly, the MMR overlap criterion fails to account for the continuation of the trend beyond period ratios of 1.14. On the other hand, the extrapolation to five planets (with K = 2) goes almost to the region where the survival time starts to increase faster (around 1.165). We conclude that our approach can successfully account for the difference between three-planet systems and systems of four or more planets.
![]() |
Fig. 14 Survival time of five-Earth-mass-planet EMS systems from Obertas et al. (2017) as a function of the period ratio. The red curve corresponds to the survival time of a three-Earth-mass-planet system (Eq. (81)) and the dashed line is the extrapolation to a five-planet system. |
7.2 Eccentric planets
It is tempting to generalise our results to systems with eccentric and inclined planets. As noted by Pu & Wu (2015), systems where planets are eccentric have similar survival times to systems with circular orbits if the spacing between the planets is measured by the distance between the apoapsis of the inner planet and the periapsis of the outer planet. In principle, an overlap criterion could be computed by taking into account the full three-planet MMR network with MMRs of arbitrary order in the same way as was done by Hadden & Lithwick (2018) for the two-planet case. Another similarity with the two-planet case is the fact that the dynamics of an isolated first-order three-planet MMR are integrable using the same strategy as in the two-planet case (such a result will be presented in a separate article, Petit, in prep.).
In practice, the structure of the full network is more complicated than the structure of a zeroth-order resonance network and the size of the resonances depends on the individual planet eccentricities. Moreover, in the case of eccentric orbits, the diffusion is not restricted to a one-dimensional direction. Indeed one has to take into account the diffusion along the eccentricity degrees of freedom contrarily to the diffusion for circular orbits, because the zeroth-order MMRs conserve the total AMD. An extension to the eccentric case will be the goal of future works.
8 Conclusion
We analyse the mechanism driving the instability time of closely packed planetary systems. Extending previous work by Quillen (2011), we use an integrable model to compute the dynamics of three zeroth-order three-planet MMRs for systems on initially circular orbits, with arbitrary planet mass distribution and spacing. We then compute the region where these resonances overlap (Eq. (59)), as this is the region where large-scale diffusion can occur (Chirikov 1979). We find that this region extends past the limit of overlap of two-planet MMR and the spacing scales as ε1∕4, where ε measures the planet-to-star mass ratio.
Inside the region of overlapping three-planet MMRs, the dynamics are not secular, despite the near-conservation of the AMD, and the period ratios can diffuse up until they reach one of the larger two-planet MMRs, then leading to rapid instability. We derive an estimated diffusion coefficient by considering only the necessary resonances and as a result, estimate the survival time (Eq. (81)). Although in general, Chirikov diffusion leads to survival times that follow power-laws (Quillen 2011), our expression is well approximated by an exponential curve, as is reported in numerical simulations. Moreover, we predict and observe on numerical simulations a change of behaviour in the region where three-planet MMRs are not overlapping. Beyond the overlap limit, the dynamics cannot be well represented by a relatively uniform diffusion mechanism, and while some other mechanism may destabilise the system, the phenomenon is expected to be much slower. The stability time therefore depends much more on the initial conditions because other mechanisms such as Arnold diffusion may be necessary to allow the planet to reach the instability.
We compare our results with numerical simulations and find excellent agreement with our analytical estimate. Moreover, we discuss how apparent discrepancies can be explained. We also discuss how this result can be extended to systems containing more planets or on eccentric orbits. Moreover, we show that the classical fit where the instability time is an exponential function of the spacing measured in Hill radius fails to capture the physical mechanism at play. In particular, we see that for very small bodies, three-body resonances can drive the instabilities over distances that are much larger than single two-planet interactions. The tools necessary to compute the time estimates and reproduce the figures are made available online17.
In this paper, we focus on systems initially outside of the influence of two-planet MMR. A similar approach could be applied to the vicinity of a two-planet MMR in order to track the system through the rapid final instability phase. Such works are necessary to understand the creation of AMD during scattering events, leading to planet collisions and ejections.
Acknowledgements
We thank the anonymous reviewer for constructive remarks that improved the manuscript. A.P. thanks A. Morbidelli and J. Laskar for useful discussions on the model and M. Pain for providing the references for the random walk exit time. We thank A. Obertas and D. Tamayo for allowing us to reproduce their data, and D. Ragozzine, D. Fabrycky and D. T. for their feedback on the preprint. This work is supported by the Royal Physiographic Society of Lund through the Fund of the Walter Gyllenberg Foundation (number 40730). M.D. and A.P. are supported by the project grant 2014.0017 “IMPACT” from the Knut and Alice Wallenberg Foundation. A.J. and A.P. are supported by the European Research Council under ERC Consolidator Grant agreement 724687-PLANETESYS, the Swedish Research Council (grant 2018-04867), and the Knut and Alice Wallenberg Foundation (grants 2014.0017 and 2017.0287). G.P. thanks the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. This research was made possible by the open-source projects rebound (Rein & Liu 2012), Jupyter (Kluyver et al. 2016), iPython (Perez & Granger 2007), numpy (van der Walt et al. 2011), scipy (Virtanen et al. 2020), pandas (Wes McKinney 2010), and matplotlib (Hunter 2007).
Appendix A Notations summary
We present in Table A.1 a summary of the notations used in this article.
Summary of the main notations used throughout the article.
Appendix B Laplace coefficients
The Laplacecoefficients appear naturally in the study of planetary systems through the development of the perturbation part in the three-body problem. The coefficient corresponds to the lth Fourier coefficient of the function
, i.e.
(B.1)
There are recurrence relations between them and we refer to Laskar & Robutel (1995) for a complete description.
One of the challenges of analytical studies of planet dynamics comes from the estimation of the Laplace coefficients. Indeed, due to the third Keplerlaw, α and l are often tied toeach other. For example, in the study of first-order MMR, it is necessary to compute an approximation for large l of the coefficient for
(Petit et al. 2017). In other words, the order in which the limits in terms of α and l are taken is relevant.
Laskar & Robutel (1995) give an alternative expression for the Laplace coefficients, in terms of hypergeometeric functions
(B.2)
where Γ is the Gamma function and2F1 is the Gaussian hypergeometric function. Laskar & Robutel (1995) use expression (B.2) to show that for α → 1, the Laplace coefficients are independent of l. However, we are interested in an estimate where we fix α and make l take larger and larger values. We cannot therefore use the equivalent they proposed.
In this article, we particularly focus on . For s = 1∕2, Eq. (B.2) becomes
(B.3)
In the limit of large l, the ratio of Γ functions is equivalent to l−1∕2 (it is worth noting that the estimate is already good for l = 1). We can therefore focus on estimating the hypergeometric function. We find that for l → ∞, converges to a value depending on α2 that we notef1∕2(α2). Taking the limit l large into the differential equation verified by the hypergeometric function (Olver et al. 2010)18, we find that f1∕2 is solution of
(B.4)
with f1∕2(0) = 1. As a result, we have
(B.5)
f1∕2 (α2) approximates extremely well the hypergeometric function as shown in Fig. B.1 where we plot for different values of l as a function of α. We note that α is plotted inlogarithmic scale centred on 1 to show where the curve starts to differ. We observe a fast convergence.
We therefore approximate the Laplace coefficients as
(B.6)
The approximation is very good (within 10%) for almost all values of l. We plot in Fig. B.2, the ratio of the exact Laplace coefficient and its estimate as a function of l for different values of α. In order to compare with Quillen (2011) estimate we plot with dashed line the ratio . We use αl instead of e−l(1−α) in Quillen’s expression to avoid an unfair comparison since the difference in the exponential would dominate the difference in the prefactors. It is critical to properly estimate the prefactor. Indeed, because we integrate over the resonance index in Sect. 4.4, the prefactor contributes significantly to the resonance density, and later to the estimate of the survival time. As a result, Quillen (2011) estimate fails to fit the Laplace coefficient and as a result, overestimates the resonance width. In their work, this effect is partially compensated by an underestimation of the Laplace coefficient derivative.
![]() |
Fig. B.1
|
![]() |
Fig. B.2 Ratio of the exact Laplace coefficients |
Appendix C Effective diffusion coefficient estimation
C.1 Exact diffusion coefficient
To compute the effective diffusion coefficient Deff (Eq. (62)), one need to solve the integral
(C.1)
The integral can be evaluated in term of the special Dawson function, which gives for Deff the expression
(C.2)
where ξov = kovδ is given by Eq. (65) and D+ is the Dawson function defined as
(C.3)
Using Eq. (63), we can replace e−ξov to obtain
(C.4)
As mentioned in the main text, we find that ξov(ξov + 1)F(ξov) is extremelywell approximated by
(C.6)
Indeed, the relative difference is below 1% for δ < 0.96 and within a factor of two overall. This expression was found by chance after trying to improve an estimation based on a development around zero in terms of δ∕δov. We plot in Fig. C.1 the exact expression of ξov(ξov + 1)F(ξov) as a function of δ∕δov as well as its estimate. As can be seen, the two curves lie on top of each other.
![]() |
Fig. C.1 ξov(ξov + 1)F(ξov) as a functionof δ∕δov, its development around zero, and the estimate (Eq. (C.6)). |
We also express the numerical factor as a function of δov,
(C.7)
Combining these terms, we obtain an expression of Deff that depends on δ and δov:
(C.8)
C.2 Exit time distribution
We provide here the distribution of the survival time. As in the main text, the interval where the system can wander has for length Δη = η+ − η−. The initial position on this interval can be measured bythe quantity u0 = (η0 − η−)∕Δη that is between 0 and 1. The distribution of the log of the survival time log10Tsurv∕T0 is given by the expression (Borodin & Salminen 2002, Eq. (3.0.2))
(C.9)
where and
.
References
- Batygin, K., Morbidelli, A., & Holman, M. J. 2015, ApJ, 799, 120 [NASA ADS] [CrossRef] [Google Scholar]
- Borodin, A. N., & Salminen, P. 2002, Handbook of Brownian Motion - Facts and Formulae, Probability and Its Applications (Basel: Birkhäuser Basel) [Google Scholar]
- Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19 [NASA ADS] [CrossRef] [Google Scholar]
- Cachucho, F., Cincotta, P. M., & Ferraz-Mello, S. 2010, Celest. Mech. Dyn. Astron., 108, 35 [CrossRef] [Google Scholar]
- Chambers, J., Wetherill, G., & Boss, A. 1996, Icarus, 119, 261 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
- Charalambous, C., Martí, J. G., Beaugé, C., & Ramos, X. S. 2018, MNRAS, 477, 1414 [CrossRef] [Google Scholar]
- Chatzigeorgiou, I. 2013, IEEE Commun. Lett., 17, 1505 [CrossRef] [Google Scholar]
- Chirikov, B. V. 1979, Phys. Rep., 52, 263 [NASA ADS] [CrossRef] [Google Scholar]
- Cincotta, P. M. 2002, New Astron. Rev., 46, 13 [NASA ADS] [CrossRef] [Google Scholar]
- Cincotta, P. M., Efthymiopoulos, C., Giordano, C. M., & Mestre, M. F. 2014, Phys. D Nonlinear Phenomena, 266, 49 [CrossRef] [Google Scholar]
- Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J., & Knuth, D. E. 1996, Adv. Comput. Math., 5, 329 [CrossRef] [MathSciNet] [Google Scholar]
- Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, 129 [NASA ADS] [CrossRef] [Google Scholar]
- Delisle, J.-B. 2017, A&A, 605, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Deprit, A. 1969, Celest. Mech., 1, 12 [NASA ADS] [CrossRef] [Google Scholar]
- Faber, P., & Quillen, A. C. 2007, MNRAS, 382, 1823 [NASA ADS] [CrossRef] [Google Scholar]
- Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146 [Google Scholar]
- Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92 [NASA ADS] [CrossRef] [Google Scholar]
- Ferraz-Mello, S. 2007, Canonical Perturbation Theories, Astrophysics and Space Science Library, (New York, NY: Springer New York), 345 [Google Scholar]
- Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81 [NASA ADS] [CrossRef] [Google Scholar]
- Giordano, C. M., & Cincotta, P. M. 2004, A&A, 423, 745 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Gladman, B. 1993, Icarus, 106, 247 [NASA ADS] [CrossRef] [Google Scholar]
- Hadden, S., & Lithwick, Y. 2018, AJ, 156, 95 [NASA ADS] [CrossRef] [Google Scholar]
- Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
- Hussain, N., & Tamayo, D. 2020, MNRAS, 491, 5258 [CrossRef] [Google Scholar]
- Johansen, A., Davies, M. B., Church, R. P., & Holmelin, V. 2012, ApJ, 758, 39 [NASA ADS] [CrossRef] [Google Scholar]
- Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Schmidt (IOS Press), 87–90 [Google Scholar]
- Laskar, J. 1991, in Predictability, Stability, and Chaos in N-Body Dynamical Systems SE - 7, ed. A. Roy (New York: Springer US), 93 [CrossRef] [Google Scholar]
- Laskar, J. 1994, A&A, 287, L9 [NASA ADS] [Google Scholar]
- Laskar, J. 1997, A&A, 317, L75 [NASA ADS] [Google Scholar]
- Laskar, J. 2000, Phys. Rev. Lett., 84, 3240 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
- Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
- Laskar, J., & Petit, A. C. 2017, A&A, 605, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
- Libert, A.-S., & Sansottera, M. 2013, Celest. Mech. Dyn. Astron., 117, 149 [NASA ADS] [CrossRef] [Google Scholar]
- Lichtenberg, A. J., & Lieberman, M. A. 1992, Regular and Chaotic Dynamics, eds. F. John, J. E. Marsden, & L. Sirovich (New York, NY: Springer New York), 38 [CrossRef] [Google Scholar]
- Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
- Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011a, Nature, 470, 53 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
- Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011b, ApJS, 197, 8 [Google Scholar]
- Marchal, C., & Bozis, G. 1982, Celest. Mech., 26, 311 [Google Scholar]
- Mayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv e-prints [arXiv:1109.2497] [Google Scholar]
- Michtchenko, T. A., Beaugé, C., & Ferraz-Mello, S. 2008, MNRAS, 387, 747 [NASA ADS] [CrossRef] [Google Scholar]
- Morbidelli, A. 2002, Modern Celestial Mechanics: Aspects of Solar System Dynamics (Oxfordshire: Taylor & Francis) [Google Scholar]
- Morbidelli, A., & Vergassola, M. 1997, J. Stat. Phys., 89, 549 [CrossRef] [Google Scholar]
- Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
- Nekhoroshev, N. N. 1977, Russ. Math. Surv., 32, 1 [NASA ADS] [CrossRef] [Google Scholar]
- Nesvorný, D., & Morbidelli, A. 1998, Celest. Mech. Dyn. Astron., 71, 243 [NASA ADS] [CrossRef] [Google Scholar]
- Obertas, A., Van Laerhoven, C., & Tamayo, D. 2017, Icarus, 293, 52 [CrossRef] [Google Scholar]
- Olver, F. W., Lozier, D. W., Boisvert, R. F., & Clark, C. W. 2010, NIST Handbook of Mathematical Functions, 1st edn. (USA: Cambridge University Press) [Google Scholar]
- Perez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21 [Google Scholar]
- Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proc. Natl.l Acad. Sci., 110, 19273 [Google Scholar]
- Petit, A. C., Laskar, J., & Boué, G. 2017, A&A, 607, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Petit, A. C., Laskar, J., & Boué, G. 2018, A&A, 617, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pichierri, G., & Morbidelli, A. 2020, MNRAS, 494, 4950 [CrossRef] [Google Scholar]
- Pichierri, G., Batygin, K., & Morbidelli, A. 2019, A&A, 625, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Poincaré, H. 1905, Leçons De Mécanique Céleste, Tome I (Gauthier-Villars. Paris) [Google Scholar]
- Pu, B., & Wu, Y. 2015, ApJ, 807, 44 [NASA ADS] [CrossRef] [Google Scholar]
- Quillen, A. C. 2011, MNRAS, 418, 1043 [CrossRef] [Google Scholar]
- Redner, S. 2001, A Guide to First-Passage Processes (Cambridge, UK: Cambridge University Press) [CrossRef] [Google Scholar]
- Rein, H., & Liu, S.-F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [NASA ADS] [CrossRef] [Google Scholar]
- Rein, H., Hernandez, D. M., Tamayo, D., et al. 2019, MNRAS, 485, 5490 [NASA ADS] [CrossRef] [Google Scholar]
- Rice, D. R., Rasio, F. A., & Steffen, J. H. 2018, MNRAS, 481, 2205 [CrossRef] [Google Scholar]
- Sansottera, M., & Libert, A.-S. 2019, Celest. Mech. Dyn. Astron., 131, 38 [NASA ADS] [CrossRef] [Google Scholar]
- Smith, A. W., & Lissauer, J. J. 2009, Icarus, 201, 381 [NASA ADS] [CrossRef] [Google Scholar]
- Tamayo, D., Silburt, A., Valencia, D., et al. 2016, ApJ, 832, L22 [NASA ADS] [CrossRef] [Google Scholar]
- Tamayo, D., Cranmer, M., Hadden, S., et al. 2020, Proc. Natl. Acad. Sci., 202001258 [Google Scholar]
- van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
- Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
- Volk, K., & Malhotra, R. 2020, AJ, 160, 98 [CrossRef] [Google Scholar]
- Wax, N., & Cohen, M. H. 2009, Phys. Today, 8, 19 [CrossRef] [Google Scholar]
- Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [NASA ADS] [CrossRef] [Google Scholar]
- Wes McKinney 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & Jarrod Millman, 51 [Google Scholar]
- Wisdom, J. 1980, AJ, 85, 1122 [NASA ADS] [CrossRef] [Google Scholar]
- Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proc. Natl. Acad. Sci., 113, 11431 [NASA ADS] [CrossRef] [Google Scholar]
- Yalinewich, A., & Petrovich, C. 2020, ApJ, 892, L11 [CrossRef] [Google Scholar]
- Zhou, J.-L., Lin, D. N. C., & Sun, Y.-S. 2007, ApJ, 666, 423 [NASA ADS] [CrossRef] [Google Scholar]
This method is a proxy to highlight the evolution of the averaged canonical coordinates defined in Sect. 3, as it removes the rapid oscillations due to the fast dynamics onto the orbits. This is sufficient here as our goal is to illustrate the mechanism at play.
We simplify the general case because the first-order MMRs are the only ones with a non-zero width for close to circular orbits. In general, we should require that is small with respect to the largest coefficient corresponding to this particular order in the sum Eq. (10). But because of their dependence on eccentricity, such a coefficient is negligible for higher order resonances.
cf. Eq. (16); note in the equation that the sum over l does not contain the secular term l = 0.
Here and later, by locally, we mean a region large enough to contain resonances of different indexes such that the exact resonance position is not relevant, but small enough such that δ and η do not vary significantly. A good example is a rectangle delimited by adjacent two-planet first-order MMRs.
A more accurate expression can be derived by analysing the perturbations of Hamiltonian (50), we refer to Chirikov (1979); Cincotta (2002).
We note that Pichierri & Morbidelli (2020) show that adding more planets into a resonant chain changes its stability. Further studies on this topic are required.
All Tables
All Figures
![]() |
Fig. 1 Semi-major axis as a function of time for an example of a three-planet EMS system (see the text for the full initial conditions). The envelope of the curve represents the extent of the orbits. We note the discontinuous vertical axis. We show the time where the first main two-planet MMR is crossed. The system becomes unstable soon afterwards. |
In the text |
![]() |
Fig. 2 Panel a: period ratio of the adjacent pairs as a function of the time to the close encounter. We note the logarithmic scale. The vertical dashed line is the same as in Fig. 1. The black horizontal dashed lines corresponds to first-order MMRs, the yellow dashed lines to the second-order MMRs. The width of the first-order MMRs is displayed in grey. Panel b: AMD normalised by the total angular momentum as a function of the time to the close encounter. Using the square root gives a typical value of the planet eccentricities. |
In the text |
![]() |
Fig. 3 Evolution in the period ratios plane. The points are colour-coded according to their time before the close encounter. We note that the system spends almost all of its time very close to its starting location as shown in Figs. 1 and 2. Green oblique lines correspond the loci of the zeroth-order three-planet MMRs (see Sect. 4.1). Chirikov’s diffusion is expected to occur perpendicularly to the network direction, along the dashed orange line. The extent of the adjacent first-order two-planet MMRs, 7:6 and 6:5, is plotted in grey. The width is computed for circular orbits. The second-order resonance 13:11 is plotted in red. Once the system enters the two-planet resonance network, the diffusion is much more rapid. |
In the text |
![]() |
Fig. 4 Zeroth-order three-planet MMR loci in the period ratio plane. The colour indicates the index p + q of the resonance. The dashed lines corresponds to first-order two-planet MMRs (the oblique ones correspond to MMRs between planets 1 and 3). The curve ν = 0.05 is displayed in orange (see Eq. (22)). |
In the text |
![]() |
Fig. 5 Number of resonances covering the period ratios space for equal-mass planets, with masses 10−5 M⊙. At first order, the number of resonances can be compared to the filling factor ρtot (Eq. (58)). We plot the two-planet circular Hill-stability limits (Petit et al. 2018) for both planet pairs in green and the predicted overlap limit for the three-planet MMRs (Eq. (59)) in blue. |
In the text |
![]() |
Fig. 6 Comparison of the two-planet stability criteria, the Hill stability (Gladman 1993; Petit et al. 2018), the first-order MMR overlap criterion, and the zeroth-order three-planet MMR overlap criterion expressed in units of the mutual Hill radius as a function of the planet-to-star-mass ratio for equal masses and equally spaced planets. |
In the text |
![]() |
Fig. 7 Top panel: probability distribution function of log10Tsurv∕T0 for various values of u0, the normalised distance of the initial condition of the system to the two neighbouring first-order two-planet MMRs. We see that for u0 = 0.5, the time distribution is log-normal as observed in Hussain & Tamayo (2020). Bottom panels: mean and standard deviation of the distribution as a function of u0. We also plotapproximate fits to the curves. We note that the mean remains within an order of magnitude of its maximum for almost all values of u0. For u0 close to 0.5, the standard deviation is close to the value 0.43 ± 0.16 measured byHussain & Tamayo (2020), plotted as a green line with its 1σ error. |
In the text |
![]() |
Fig. 8 Survival time for a three-planet EMS system for planet masses of 10−7 M⊙ as a functionof the initial period ratio. The red curve corresponds to the survival time estimate (Eq. (81)), the blue vertical line to the Hill-stability limit (Petit et al. 2018), the orange vertical line to the two-planet MMR overlap criterion (Wisdom 1980). The dashed black (resp. green) lines are the two-planet first(resp. second)-order MMR. The light orange rectangles show an estimate of the width of the two-planet MMR (Petit et al. 2017). |
In the text |
![]() |
Fig. 9 Survival time for a three EMS planet system for planet masses of 10−5 M⊙ as a function of the initial period ratio. See Fig. 8 for a detailed caption. The red dotted line corresponds to an alternate estimate position for the limit of the overlapped region (see text and Eq. (84)). |
In the text |
![]() |
Fig. 10 Survival time for a three EMS planet system for planet masses of 10−3 M⊙ as a function of the initial period ratio. See Fig. 8 for a detailed caption. The red dotted line corresponds to an alternate estimate position for the limit of the overlapped region (see text and Eq. (84)). |
In the text |
![]() |
Fig. 11 Survival time for non-equally spaced systems described in Sect. 6.2 as a function of the initial inner-pair period ratio. The outer-pair period ratio is given at the top of the figure. The dashed survival-time estimate corresponds to the EMS case. See Fig. 8 for a detailed caption. |
In the text |
![]() |
Fig. 12 Survival time for a system composed of an inner test particle and two planets of equal mass (10−5 M⊙) with equal spacing as a function of the initial period ratio. See Fig. 8 for a detailed caption. The full red curve is theprediction using the survival time estimate (Eq. (81)) and the dashed line corresponds to the same estimate for the EMS case shown in Fig. 9. The red dotted line corresponds to an alternate estimate position for the limit of the overlapped region (see text and Eq. (84)). |
In the text |
![]() |
Fig. 13 Same as Fig. 12 but with the test particle being placed between the two planets.. |
In the text |
![]() |
Fig. 14 Survival time of five-Earth-mass-planet EMS systems from Obertas et al. (2017) as a function of the period ratio. The red curve corresponds to the survival time of a three-Earth-mass-planet system (Eq. (81)) and the dashed line is the extrapolation to a five-planet system. |
In the text |
![]() |
Fig. B.1
|
In the text |
![]() |
Fig. B.2 Ratio of the exact Laplace coefficients |
In the text |
![]() |
Fig. C.1 ξov(ξov + 1)F(ξov) as a functionof δ∕δov, its development around zero, and the estimate (Eq. (C.6)). |
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.