© 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 10⁹ to 10¹¹ 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 $$\Delta =\frac{a_2-a_1}{a_1}>2\sqrt{3}\frac{R_{\text{H}}}{a_1}\simeq 3.46\frac{R_{\text{H}}}{a_1},$$(1)

where $$R_{\text{H}}=\frac{a_1+a_2}{2}{\left(\frac{m_1+m_2}{3m_0}\right)}^{1/3}$$(2)

is the mutual Hill radius with m₁, m₂ being the planet masses and m₀, 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 ${\left((m_1+m_2)/m_0\right)}^{2/7}$ (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 10⁵ 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 Δ = (a_k+1 − a_k)∕a_k. 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 R_H, $${\log }_{10}\frac{T_{\text{surv}}}{P}=b\frac{\Delta }{R_{\text{H}}}+c,$$(3)

where P is a typical orbital period and b and c are numerical factors¹. 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 $${\log }_{10}\frac{T_{\text{surv}}}{P}=b^{\prime }\Delta {\left(\frac{m_p}{m_0}\right)}^{-1/4}-c^{\prime }-\log \frac{m_p}{m_0},$$(4)

where b′ and c′ are positive numerical coefficients independent of the masses, m_p is the planet mass and m₀ 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 T_surv∕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 R_H. 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 T_surv 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 T_surv seems to have an exponential dependency on the orbital spacing, measured in units of ${(m_p/m_0)}^{1/4}$. The fit is valid over 6 to 8 orders of magnitude for survival times between 100 and almost 10¹⁰ 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 m_p = 10⁻⁵ M_⊙ and orbit a solar-mass star. The inner orbit is at 1 au and the period ratios between adjacent planets is initially close² to P_k+1 ∕P_k = 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⁻¹⁰.

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 snapshots³. 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 P_k∕P_k+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, $\sqrt{C/G}$, 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 P₂ ∕P₃ as a function of the inner pair period ratio P₁∕P₂. 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 m₁, m₂, and m₃ orbiting a star of mass m₀. The canonical positions r_j and momenta ${\tilde{r}}_j$ 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 a_j, the eccentricities e_j, 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): $$\begin{array}{lll}\hfill & {\Lambda }_j=m_j\sqrt{\mu a_j},\hfill & {\lambda }_j,\hfill \\ \hfill & C_j={\Lambda }_j\left(1-\sqrt{1-e_j^2}\right),\hfill & -{\varpi }_j,\hfill \end{array}$$(5)

where $\mu =\mathcal{G}{m}_0$ and $\mathcal{G}$ 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 (C_j, −ϖ_j) can also be replaced by their associated complex variables $$\begin{array}{lll}\hfill & x_j=\sqrt{C_j}{e}^{\iota {\varpi }_j},\hfill & -\iota {\overline{x}}_j,\hfill \end{array}$$(6)

with $\iota =\sqrt{-1}$ (x_j are the canonical momenta and $-\iota {\overline{x}}_j$ the conjugated positions). For small eccentricities, we have $x_j\simeq \sqrt{{\Lambda }_j/2}{e}_j{e}^{\iota {\varpi }_j}$. The system total angular momentum G and AMD C are given by $$G={\sum }_{j=1}^3({\Lambda }_j-C_j)\text{and}C={\sum }_{j=1}^3{C}_j.$$(7)

The Hamiltonian $\mathcal{H}$ describing the dynamics can be split into an integrable part $${\mathcal{H}}_0={\sum }_{j=1}^3\frac{\Vert {\tilde{r}}_j{\Vert }^2}{2m_j}-\frac{\mu m_j}{r_j}=-{\sum }_{j=1}^3\frac{{\mu }^2{m}_j^3}{2{\Lambda }_j},$$(8)

describing the motion on unperturbed Keplerian orbits, and a perturbation, $$\epsilon {\mathcal{H}}_1=-{\sum }_{j<j^{\prime }}\frac{\mathcal{G}{m}_j{m}_{j^{\prime }}}{|r_j-r_{j^{\prime }}|}+\frac{1}{2m_0}{\Vert {\sum }_{j=1}^3{\tilde{r}}_j\Vert }^2,$$(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 $$\epsilon {\mathcal{H}}_1={\sum }_{k,l,\overline{l},}{\mathcal{C}}_{k,l,\overline{l}}(m_j,{\Lambda }_j)\left({\displaystyle \prod _{j=1}^3{x}_j^{l_j}}{\overline{x}}_j^{{\overline{l}}_j}\right)e^{\iota k\cdot \lambda },$$(10)

where $k\in {\mathbb{Z}}^3$, $l,\overline{l}\in {\mathbb{N}}^3$.

Due to the conservation of angular momentum, the coefficient ${\mathcal{C}}_{k,l,\overline{l}}$ must vanish unless the indices $k,l,\overline{l}$ verify the d’Alembert rules (e.g. Morbidelli 2002) and in particular $${\sum }_{j=1}^3{k}_j+l_j-{\overline{l}}_j=0.$$(11)

In the unperturbed case, the system is said to be in a MMR if the mean motions $$n_j={\dot{\lambda }}_j=\frac{\partial {\mathcal{H}}_0}{\partial {\Lambda }_j}=\frac{{\mu }^2{m}_j^3}{{\Lambda }_j^3}$$(12)

verify an equation of the form $${\sum }_{j=1}^3{k}_j{n}_j=0.$$(13)

The sum k = k₁ + k₂ + k₃ is the “order” of the resonance. The sum K = |k₁| + |k₂| + |k₃| is the “index” of the resonance. In the general case, the perturbation $\epsilon {\mathcal{H}}_1$ also influences the resonantdynamics. The terms in Eq. (10) contributing to the resonance are the ones that depend on the combination of mean longitudes k₁λ₁ + k₂λ₂ + k₃λ₃. 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 $\mathcal{H}={\mathcal{H}}_0+\epsilon {\mathcal{H}}_1$ that depend on angles of theform k ⋅λ with k_j≠0 for all j. This means that there are no three-planet resonances at the first order in ε. There are instead of course $\mathcal{O}(\epsilon )$ two-planet MMR terms.

Three-planet MMRs actually emerge in the perturbative Hamiltonian as $\mathcal{O}({\epsilon }^2)$ 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 MMR⁴, that is, we assume that $k{n}_j-(k+1)n_{j^{\prime }}$ 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 $\mathcal{O}({\epsilon }^2)$ 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 $${\mathcal{H}}^{\prime }=H_0({\Lambda }^{\prime })+\epsilon {\overline{\mathcal{H}}}_1({\Lambda }^{\prime },x^{\prime })+{\epsilon }^2{\mathcal{H}}_{2^{\prime }}({\Lambda }^{\prime },{\lambda }^{\prime },x^{\prime })+\text{O}({\epsilon }^3),$$(14)

where $\epsilon {\overline{\mathcal{H}}}_1$ is the average of $\epsilon {\mathcal{H}}_1$ over the mean longitudes λ, and ${\epsilon }^2{\mathcal{H}}_{2^{\prime }}$ 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 ({εχ₁, ⋅}) where {⋅, ⋅} is the Poisson bracket⁵ and εχ₁ is the solution of the homological equation $$\left\{\epsilon {\chi }_1,{\mathcal{H}}_0\right\}+\epsilon {\mathcal{H}}_1=\epsilon {\overline{\mathcal{H}}}_1.$$(15)

More precisely, if we note $h_1^{(k)}=h_1^{(k)}(x,\overline{x}),$ the complex Fourier coefficients of ${\mathcal{H}}_1$ with respect to the mean longitudes, we can write $$\epsilon {\chi }_1=\epsilon {\sum }_{k\ne 0}\frac{h_1^{(k)}}{\iota k\cdot n}{e}^{\iota k\cdot \lambda }.$$(16)

Due to the expression of $\epsilon {\mathcal{H}}_1$ given in Eq. (10), the denominators k ⋅n are of the form $k_j{n}_j+k_{j^{\prime }}{n}_{j^{\prime }}$ 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 ${\mathcal{H}}_1$ have sizes smaller than ε², 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 ${\epsilon }^2{\mathcal{H}}_{2^{\prime }}$ as $${\epsilon }^2{\mathcal{H}}_{2^{\prime }}=\frac{1}{2}\left\{\epsilon {\chi }_1,\epsilon {\mathcal{H}}_1+\epsilon {\overline{\mathcal{H}}}_1\right\}.$$(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 ε² close to the first-order averaged coordinates, and the new Hamiltonian is the average of ${\mathcal{H}}^{\prime }$ 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 ε³ 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 $$p{n}_1-(p+q)n_2+q{n}_3=0.$$(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 three⁶ degrees of freedom (Λ_j, λ_j). The resonance equation defines a plane in the frequency space (n₁, n₂, n₃). Because the gravitational interactions are scale invariant, we can restrict ourselves to a two-dimensional plane corresponding to the period ratios ν₁₂ and ν₂₃ where $${\nu }_{ij}=\frac{P_i}{P_j}=\frac{n_j}{n_i}.$$(19)

Dividing Eq.(18) by n₂, and reorganising terms, one gets $${\nu }_{23}=1-\frac{p}{q}({\nu }_{12}^{-1}-1),$$(20)

that is the equation of a straight line passing through the point (1,1) with slope − p∕q for the period ratios ν₂₁ = P₂∕P₁ > 1 and ν₂₃ = P₂∕P₃ < 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 ν₂₁ for tightly packed planets (i.e. close to 1) leads to a poorer approximation at first order than expanding ν₁₂. We therefore only consider the variables ν₁₂ and ν₂₃ as done in Eq. (20). In the plane (ν₁₂, ν₂₃), 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 2n₁ − 4n₂ + 2n₃ = 0 is the same as the resonance n₁ − 2n₂ + n₃ = 0. The full network is dense in the (ν₁₂, ν₂₃) 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 defined⁷ to take advantage of this property. We define the resonance locator $$\eta =\frac{1-{\nu }_{23}}{{\nu }_{12}^{-1}-{\nu }_{23}}=\frac{{\nu }_{12}(1-{\nu }_{23})}{1-{\nu }_{12}{\nu }_{23}}=\frac{p}{p+q}.$$(21)

The second equality is only valid on resonances. The position along the resonance can be defined by a generalised period ratio separation $$\nu =\frac{(1-{\nu }_{12})(1-{\nu }_{23})}{1-{\nu }_{12}{\nu }_{23}}=\frac{p}{p+q}({\nu }_{12}^{-1}-1)=\frac{q}{p+q}(1-{\nu }_{23}),$$(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 $${\nu }_{12}=\frac{\eta }{\eta +\nu },\text{and}{\nu }_{23}=\frac{1-\eta -\nu }{1-\eta }.$$(23)

The levels of constant ν are hyperbola with horizontal and vertical asymptotes ${(1+\nu )}^{-1}$. 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 $$\begin{array}{l}{\theta }_{\text{res}}=p{\lambda }_1-(p+q){\lambda }_2+q{\lambda }_3,\\ {\theta }_{\Gamma }={\lambda }_2-{\lambda }_3,\\ {\theta }_{{\rm Y}}={\lambda }_3.\end{array}$$(24)

The conjugated momenta are $$\begin{array}{ll}\Theta \hfill & =\frac{1}{p}{\Lambda }_1,\hfill \\ \Gamma \hfill & =\frac{p+q}{p}{\Lambda }_1+{\Lambda }_2,\hfill \\ {\rm Y}\hfill & ={\Lambda }_1+{\Lambda }_2+{\Lambda }_3.\hfill \end{array}$$(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 $$\mathcal{H}={\mathcal{H}}_0(\Theta )+\epsilon {\overline{\mathcal{H}}}_1(\Theta )+{\epsilon }^2{\overline{\mathcal{H}}}_{2,\text{res}}(\Theta ,{\theta }_{\text{res}})+\text{O}({\epsilon }^3),$$(26)

where ${\epsilon }^2{\overline{\mathcal{H}}}_{2,\text{res}}(\Theta ,{\theta }_{\text{res}})$ 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 $\epsilon {\overline{\mathcal{H}}}_1(\Theta )$ as it is small with respect to ${\mathcal{H}}_0$ 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(ε³). 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 Θ = Θ₀ + dΘ where Θ₀ 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 $$\begin{array}{ll}\text{d}{\Lambda }_1\hfill & =p\text{d}\Theta ,\hfill \\ \text{d}{\Lambda }_2\hfill & =-(p+q)\text{d}\Theta ,\hfill \\ \text{d}{\Lambda }_3\hfill & =q\text{d}\Theta ,\hfill \end{array}$$(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 ν₁₂ and d ν₂₃ as a function of dΘ, $$\begin{array}{ll}\text{d}{\nu }_{12}\hfill & =3{\nu }_{12,0}\left(\frac{p}{{\Lambda }_{1,0}}+\frac{p+q}{{\Lambda }_{2,0}}\right)\text{d}\Theta ,\hfill \\ \text{d}{\nu }_{23}\hfill & =-3{\nu }_{23,0}\left(\frac{p+q}{{\Lambda }_{2,0}}+\frac{q}{{\Lambda }_{3,0}}\right)\text{d}\Theta .\hfill \end{array}$$(28)

We can take the ratio of the small variations dν₂₃ and dν₁₂, and using Eq. (20) to replace p and q we obtain the differential equation $$\frac{\text{d}{\nu }_{23}}{\text{d}{\nu }_{12}}=-\frac{{\nu }_{23}}{{\nu }_{12}}\frac{m_2(1-{\nu }_{12}{\nu }_{23})+m_1{\nu }_{12}^{-1/3}(1-{\nu }_{12})}{m_2(1-{\nu }_{12}{\nu }_{23})+m_3{\nu }_{12}{\nu }_{23}^{1/3}(1-{\nu }_{12})},$$(29)

which gives the direction of the change of period ratios anywhere in the plane (ν₁₂, ν₂₃) 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 $\epsilon {\mathcal{H}}_1$ 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 $\epsilon {\mathcal{H}}_1$ that we consider are therefore (Laskar & Robutel 1995; Murray & Dermott 1999) $$\begin{array}{ll}\epsilon {\mathcal{H}}_1=\hfill & {\sum }_{1\le i<j\le 3}{\sum }_{l\in \mathbb{Z}}{W}_{ij}^l{e}^{\iota l({\lambda }_i-{\lambda }_j)}\hfill \\ \hfill & +{\sum }_{1\le i<j\le 3}{\sum }_{l\ge 0}\left(V_{ij,<}^l{x}_i+V_{ij,>}^l{x}_j\right)e^{\iota (l{\lambda }_i-(l+1){\lambda }_j)}+c.c.,\hfill \end{array}$$(30)

where c.c. designates the complex conjugate of the second sum, and $$\begin{array}{l}{W}_{ij}^l=-\frac{m_i{n}_j{\Lambda }_j}{2m_0}{b}_{1/2}^{(l)}({\alpha }_{ij}),\\ V_{ij,<}^l=\frac{m_i{n}_j{\Lambda }_j}{2m_0}\sqrt{\frac{2}{{\Lambda }_i}}\left(l+1+\frac{{\alpha }_{ij}}{2}\frac{\partial }{\partial \alpha }\right)b_{1/2}^{(l)}({\alpha }_{ij}),\\ V_{ij,>}^l=-\frac{m_i{n}_j{\Lambda }_j}{2m_0}\sqrt{\frac{2}{{\Lambda }_j}}\left(l+\frac{1}{2}+\frac{{\alpha }_{ij}}{2}\frac{\partial }{\partial \alpha }\right)b_{1/2}^{(l)}({\alpha }_{ij}),\end{array}$$

where α_ij = a_i∕a_j and $b_s^{(l)}(\alpha )$ 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 $$b_{1/2}^{(l)}(\alpha )\simeq \frac{1}{2}\left|\ln (1-\alpha )\right|e^{-\left|l\right|(1-\alpha )}.$$(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 $$b_{1/2}^{(l)}(\alpha )\simeq \frac{2{\alpha }^{\left|l\right|}}{\sqrt{\pi \left|l\right|(1-{\alpha }^2)}}.$$(35)

In the close planet approximation, α → 1⁻, we get $$b_{1/2}^{(l)}(\alpha )\simeq \sqrt{\frac{2}{\pi \left|l\right|(1-\alpha )}}{e}^{-\left|l\right|(1-\alpha )}.$$(36)

We refer to Appendix B for a detailed discussion.

The exponential dependency on −|l|(1 − α) of the Laplace coefficients has two consequences for $W_{ij}^l$. The interactions between planets 1 and 3 are always negligible with respect to the interaction in adjacent pairs, that is, for a given l, $\left|W_{13}^l\right|\ll \left|W_{j,j+1}^l\right|$. 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 (ε³).

With these clarifications, let us return to the calculation of the perturbative steps, starting from the perturbative Hamiltonian (Eq. (30)). The solution εχ₁ to the corresponding homological equation has for expression⁸ $$\begin{array}{ll}\epsilon {\chi }_1=\hfill & {\sum }_{1\le i<j\le 3}{\sum }_{l\in {\mathbb{Z}}^{*}}\frac{W_{ij}^l}{\iota l(n_i-n_j)}{e}^{\iota l({\lambda }_i-{\lambda }_j)}\hfill \\ \hfill & +{\sum }_{1\le i<j\le 3}{\sum }_{l\ge 0}\frac{\left(V_{ij,<}^l{x}_i+V_{ij,>}^l{x}_j\right)}{\iota (l{n}_i-(l+1)n_j)}{e}^{\iota (l{\lambda }_i-(l+1){\lambda }_j)}+c.c.,\hfill \end{array}$$(37)

where ${\mathbb{Z}}^{*}=\mathbb{Z}\setminus \left\{0\right\}$. 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 $\mathcal{O}({\epsilon }^2)$ terms, which we now calculate explicitly. In essence, we must only calculate the term ${\epsilon }^2{\overline{\mathcal{H}}}_2$ in Eq. (26). Since we would subsequently average over θ_ϒ and θ_Γ, the only terms that remain in ${\epsilon }^2{\overline{\mathcal{H}}}_2$ must depend on the angle (pλ₁ − (p + q)λ₂ + qλ₃) or its opposite. Because of the form of εχ₁ and $\epsilon {\mathcal{H}}_1$, the only terms contributing at zeroth order in eccentricity to the averaged Hamiltonian ${\epsilon }^2{\overline{\mathcal{H}}}_2$ are of the form $$\begin{array}{l}\left\{\frac{W_{12}^p}{\iota p(n_1-n_2)}{e}^{\iota p({\lambda }_1-{\lambda }_2)},W_{23}^q{e}^{-\iota q({\lambda }_2-{\lambda }_3)}\right\},\\ \left\{\frac{W_{23}^q}{-\iota q(n_2-n_3)}{e}^{-\iota q({\lambda }_2-{\lambda }_3)},W_{12}^p{e}^{\iota p({\lambda }_1-{\lambda }_2)}\right\},\\ \left\{\frac{V_{12,>}^p}{\iota (p{n}_1-(p+1)n_2)}{x}_2{e}^{\iota (p{\lambda }_1-(p+1){\lambda }_2)},V_{23,<}^{q-1}{\overline{x}}_2{e}^{-\iota ((q-1){\lambda }_2-q{\lambda }_3)}\right\}\\ \left\{\frac{V_{23,<}^{q-1}}{-\iota ((q-1)n_2-q{n}_3)}{\overline{x}}_2{e}^{-\iota ((q-1){\lambda }_2-q{\lambda }_3)},V_{12,>}^p{x}_2{e}^{\iota (p{\lambda }_1-(p+1){\lambda }_2)}\right\}\end{array}$$(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 Λ₂, λ₂, and x₂ 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 $\iota \frac{\partial \cdot }{\partial x_2}\frac{\partial \cdot }{\partial {\overline{x}}_2}-\iota \frac{\partial \cdot }{\partial {\overline{x}}_2}\frac{\partial \cdot }{\partial x_2}$. This gives $${\epsilon }^2{\overline{\mathcal{H}}}_2={\epsilon }^2{R}_{pq}\cos ({\theta }_{\text{res}}),$$(39)

where $$\begin{array}{ll}{\epsilon }^2{R}_{pq}=\hfill & -\left(\frac{q}{p(n_1-n_2)}+\frac{1}{n_2-n_3}\right)\frac{\partial W_{12}^p}{\partial {\Lambda }_2}{W}_{23}^q\hfill \\ \hfill & +\left(\frac{1}{n_1-n_2}+\frac{p}{q(n_2-n_3)}\right)W_{12}^p\frac{\partial W_{23}^q}{\partial {\Lambda }_2}\hfill \\ \hfill & -\left(\frac{q}{p{(n_1-n_2)}^2}+\frac{p}{q{(n_2-n_3)}^2}\right)W_{12}^p{W}_{23}^q\frac{\partial n_2}{\partial {\Lambda }_2}\hfill \\ \hfill & +\left(\frac{1}{p{n}_1-(p+1)n_2}+\frac{1}{(q-1)n_2-q{n}_3}\right)V_{12,>}^p{V}_{23,<}^{q-1}.\hfill \end{array}$$(40)

Because ${\epsilon }^2{\overline{\mathcal{H}}}_2$ is small with respect to the Keplerian part, we evaluate the action at the nominal resonance value, that is, pn₁ − (p + q)n₂ + qn₃ = 0. Using Eqs. (18), (19), (31)–(33) as well as $$\frac{\partial {\alpha }_{12}}{\partial {\Lambda }_2}=-2\frac{{\alpha }_{12}}{{\Lambda }_2},\frac{\partial {\alpha }_{23}}{\partial {\Lambda }_2}=2\frac{{\alpha }_{23}}{{\Lambda }_2}\text{ and }\frac{\partial n_2}{\partial {\Lambda }_2}=-3\frac{n_2}{{\Lambda }_2},$$(41)

we can express ε²R_pq as $$\begin{array}{ll}{\epsilon }^2{R}_{pq}=\hfill & \frac{m_1{m}_3{n}_2{\Lambda }_2{\alpha }_{23}}{m_0^2\nu }[(1-\eta )b_{1/2}^{(q)}\left(1+{\alpha }_{12}\frac{\partial }{\partial \alpha }\right)b_{1/2}^{(p)}\hfill \\ \hfill & +\eta {\alpha }_{23}{b}_{1/2}^{(p)}\frac{\partial b_{1/2}^{(q)}}{\partial \alpha }\hfill \\ \hfill & +\frac{3\eta (1-\eta )}{2\nu }{b}_{1/2}^{(p)}{b}_{1/2}^{(q)}]\hfill \\ \hfill & +\frac{m_1{m}_3{n}_2{\Lambda }_2{\alpha }_{23}}{m_0^2(\nu -{(p+q)}^{-1})}[\frac{1-\eta }{2}{b}_{1/2}^{(q)}\left(1+{\alpha }_{12}\frac{\partial }{\partial \alpha }\right)b_{1/2}^{(p)}\hfill \\ \hfill & +\left(\frac{\eta }{2}+\frac{1}{4(p+q)}\right){\alpha }_{23}{b}_{1/2}^{(p)}\frac{\partial b_{1/2}^{(q)}}{\partial \alpha }\hfill \\ \hfill & +\frac{{\alpha }_{12}{\alpha }_{23}}{4(p+q)}\frac{\partial b_{1/2}^{(q)}}{\partial \alpha }\frac{\partial b_{1/2}^{(p)}}{\partial \alpha }\hfill \\ \hfill & +\eta (1-\eta )(p+q)b_{1/2}^{(p)}{b}_{1/2}^{(q)}],\hfill \end{array}$$(42)

where the Laplace coefficients depending on p (resp. q) are evaluated at α₁₂ (resp. α₂₃).

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 ν₁₂ = p∕(p + 1) and ν₂₃ = (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 $${\delta }_{ij}=1-{\alpha }_{ij}\simeq \frac{2}{3}(1-{\nu }_{ij}),$$(43)

which is an excellent estimate for period ratios from 0.5 to 1.

The product of the Laplace coefficients and their derivatives in R_pq introduces an exponential factor of the form e^{−pδ₁₂−qδ₂₃} (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: $$p{\delta }_{12}\simeq q{\delta }_{23}.$$(44)

By analogy with the generalised period ratio separation ν, we define a generalised orbital spacing that we note $$\delta =\frac{{\delta }_{12}{\delta }_{23}}{{\delta }_{12}+{\delta }_{23}}\simeq \frac{p}{p+q}{\delta }_{12}\simeq \frac{q}{p+q}{\delta }_{23},$$(45)

where the two last equalities are approximations using Eq. (44). We have $$\nu \simeq \frac{3}{2}\delta .$$(46)

Using the newly defined variable δ and η, the expression of the coefficient ε²R_pq can be simplified to $${\epsilon }^2{R}_{pq}=\frac{m_1{m}_3{n}_2{\Lambda }_2}{3\pi m_0^2}\frac{\eta (1-\eta )}{{\delta }^2}\left(17+\frac{21}{(p+q)\delta }\right)e^{-2(p+q)\delta },$$(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 $\epsilon {\overline{\mathcal{H}}}_1$); 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 $${\mathcal{H}}_0\simeq {\sum }_{j=1}^3\left(-\frac{{\mu }^2{m}_j^3}{2{\Lambda }_{j,0}^2}+n_{j,0}({\Lambda }_j-{\Lambda }_{0,j})-\frac{3n_{j,0}}{2{\Lambda }_{j,0}}{({\Lambda }_j-{\Lambda }_{0,j})}^2\right),$$(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 $$-\frac{{\mathcal{K}}_2}{2}=-\frac{3n_2}{2{\Lambda }_2}{(p+q)}^2\left(1+\frac{m_2}{m_1}\frac{{\eta }^2}{{\alpha }_{12}^2}+\frac{m_2}{m_3}{\alpha }_{23}^2{(1-\eta )}^2\right).$$(49)

We note that ${\mathcal{K}}_2$ 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 $${\mathcal{H}}_{\text{res}}=-\frac{{\mathcal{K}}_2}{2}{(\Theta -{\Theta }_0)}^2+{\epsilon }^2{R}_{pq}\cos {\theta }_{\text{res}}.$$(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) $$\Delta \Theta =2\epsilon \sqrt{\frac{2R_{pq}}{{\mathcal{K}}_2}}.$$(51)

The small oscillation frequency is given by $${\omega }_{pq}=\epsilon \sqrt{{\mathcal{K}}_2{R}_{pq}}=n_2\epsilon MA\frac{\sqrt{\eta (1-\eta ))}}{\delta }(p+q)e^{-(p+q)\delta },$$(52)

where $A=\sqrt{\frac{38}{\pi }}=3.47$ and $$\epsilon M=\frac{\sqrt{m_1{m}_3+m_2{m}_3{\eta }^2{\alpha }_{12}^{-2}+m_1{m}_2{\alpha }_{23}^2{(1-\eta )}^2}}{m_0}$$(53)

is the relevant mass ratio for the studied problem. For equal mass and equal tight spacing we have $\epsilon M=\sqrt{\frac{3}{2}}\frac{m_p}{m_0}$. 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 $$\frac{\text{d}\eta }{\text{d}\Theta }=\frac{{\mathcal{K}}_2}{n_2(p+q)}\frac{\eta (1-\eta )}{\nu }.$$(54)

This means that the width in terms of η can be estimated as $$\Delta {\eta }_{pq}=\frac{4\sqrt{2}\eta (1-\eta ){\omega }_{pq}}{3(p+q)\delta n_2}=6.55\epsilon M\frac{{(\eta (1-\eta ))}^{3/2}}{{\delta }^2}{e}^{-(p+q)\delta }.$$(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 (ν₁₂, ν₂₃) where resonances overlap. Because of the expression of resonance width Δη_pq, we see that the width of the resonances close to a given point (ν₁₂, ν₂₃) 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 locally⁹ 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 index¹⁰. 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 index¹¹ 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 $${\eta }_{p+1,q-1}-{\eta }_{p,q}=\frac{1}{p+q}=\frac{1}{k}.$$(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, $${\rho }_k=k\Delta {\eta }_{pq}=\epsilon M\frac{{(\eta (1-\eta ))}^{3/2}}{{\delta }^2}k{e}^{-k\delta }.$$(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 $$\begin{array}{ll}{\rho }_{\text{tot}}\hfill & =6.55\epsilon M\frac{{(\eta (1-\eta ))}^{3/2}}{{\delta }^2}{\displaystyle {\int }_0^{+\infty }k}{e}^{-k\delta }\text{d}k\hfill \\ \hfill & =6.55\epsilon M\frac{{(\eta (1-\eta ))}^{3/2}}{{\delta }^4},\hfill \end{array}$$(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 δ⁻⁴. 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 $${\delta }_{\text{ov}}={\epsilon }^{1/4}{M}^{1/4}{(6.55)}^{1/4}{(\eta (1-\eta ))}^{3/8}.$$(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 δ $${\rho }_{\text{tot}}={\left(\frac{\delta }{{\delta }_{\text{ov}}}\right)}^{-4}.$$(60)

In the case of equal mass and spacing systems, Eq. (59) becomes $${\delta }_{\text{ov},\text{eq}}=1.16{\left(\frac{m_p}{m_0}\right)}^{1/4}.$$(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 (ν₁₂, ν₂₃) for three equal-mass planets. The mass of each planet is 10⁻⁵ 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 R_pq (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 m_p 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 T_surv.

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) $$D_{p+q}=\frac{{\omega }_{pq}}{2\pi }{\left(\Delta {\eta }_{pq}\right)}^2.$$(62)

This expression is an estimate that is valid¹² 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 k_ov such that the subnetwork composed of the three-planet MMR with index k smaller than k_ov locally covers the space. Using Eqs. (57) and (60) we have $$\frac{{\delta }_{\text{ov}}^4}{{\delta }^2}{\displaystyle {\int }_0^{k_{\text{ov}}}k}{e}^{-k\delta }\text{d}k={\left(\frac{{\delta }_{\text{ov}}}{\delta }\right)}^4\left[1-(k_{\text{ov}}\delta +1)e^{-k_{\text{ov}}\delta }\right]=1.$$(63)

We thus have an implicit definition of k_ov. We also note that the equation depends on k_ovδ rather than k_ov alone. As a result, we define the variable $${\xi }_{\text{ov}}=k_{\text{ov}}\delta ,$$(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) $${\xi }_{\text{ov}}=-1-W_{-1}\left(\frac{{\delta }^4-{\delta }_{\text{ov}}^4}{e{\delta }_{\text{ov}}^4}\right),$$(65)

where W₋₁ is the real valued branch with values smaller than −e defined between − 1 and 0. The function W is the inverse function to z →ze^z. Using bounds on W₋₁ by Chatzigeorgiou (2013), an excellent approximation of ξ_ov is $${\xi }_{\text{ov}}=\sqrt{-2\ln \left(1-\frac{{\delta }^4}{{\delta }_{\text{ov}}^4}\right)}-\frac{2}{3}\ln \left(1-\frac{{\delta }^4}{{\delta }_{\text{ov}}^4}\right).$$(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 D_pq 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 η₀, and becomes unstable once η reaches η₀ + η_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 $$\frac{\text{d}\eta }{\text{d}t}=s(\eta )b(t),$$(67)

where s is related to the local diffusion coefficient and b(t) is a Gaussian white noise with zero average, verifying $$\langle b(t)b(t^{\prime })\rangle =2\delta (t-t^{\prime }).$$(68)

For s = s₀ constant, Eq. (67) is the classical Langevin equation. The associated diffusion equation for the probability density p is $$\frac{\partial p}{\partial t}=s_0^2\frac{{\partial }^{2}p}{\partial {\eta }^2},$$(69)

where the diffusion coefficient is $s_0^2$. By analogy with the case where s is constant, we define $s(\eta )=\sqrt{D_{p+q}}$, where p and q define the largest MMR that contains the point η. Since in the region considered, the resonances overlap, the interval (η₀, η₀ + η_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 $$P\left(\eta \text{isin3-pl}.\text{MMRofindexk}\right)=\{\begin{array}{ll}{\rho }_k\text{d}k,\hfill & k\le k_{\text{ov}},\hfill \\ 0,\hfill & k>k_{\text{ov}}\hfill \end{array}.$$(70)

We note that the considered value of η could also be containedin a higher index resonance than k_ov. 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 $$y(\eta )={\displaystyle {\int }_{{\eta }_0}^{\eta }\frac{\text{d}{\eta }^{\prime }}{s({\eta }^{\prime })}}.$$(71)

Indeed, computing the time derivative of y using Eqs. (67) and (71), we have $$\frac{\text{d}y}{\text{d}t}=\frac{1}{s(\eta )}\frac{\text{d}\eta }{\text{d}t}=b(t),$$(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 $\overline{\eta }$ the average values, we have $$\begin{array}{ll}\overline{y}(\eta )\hfill & ={\displaystyle {\int }_0^{k_{\text{ov}}}{\displaystyle {\int }_{{\eta }_0}^{\overline{\eta }}y}}({\eta }^{\prime })P({\eta }^{\prime })\text{d}{\eta }^{\prime }\text{d}k={\displaystyle {\int }_{{\eta }_0}^{\overline{\eta }}\text{d}}{\eta }^{\prime }{\displaystyle {\int }_0^{k_{\text{ov}}}\frac{{\rho }_k}{\sqrt{D_k}}}\text{d}k\hfill \\ \hfill & =\left(\overline{\eta }-{\eta }_0\right){\displaystyle {\int }_0^{k_{\text{ov}}}\frac{{\rho }_k}{\sqrt{D_k}}}\text{d}k.\hfill \end{array}$$(73)

We see that the variation of ȳ is proportionalto the variation of $\overline{\eta }$, the integral being almost constant¹³ 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 $$D_{\text{eff}}={\left({\displaystyle {\int }_0^{k_{\text{ov}}}\frac{{\rho }_k}{\sqrt{D_k}}}\text{d}k\right)}^{-2}={\left({\displaystyle {\int }_0^{k_{\text{ov}}}\frac{k}{\sqrt{{\omega }_k}}}\text{d}k\right)}^{-2}.$$(74)

Indeed, deriving Eq. (73) gives $$\frac{\text{d}\overline{\eta }}{\text{d}t}=\sqrt{D_{\text{eff}}}b(t).$$(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 article¹⁴.

For the remainder of the discussion, a very good estimate¹⁵ is obtained with the expression $$D_{\text{eff}}\simeq \epsilon MA{n}_2\frac{9\sqrt{\eta (1-\eta )}}{4\pi \sqrt{2}}\frac{{\delta }_{\text{ov}}^6}{{\delta }^4}\left(1-\frac{{\delta }^4}{{\delta }_{\text{ov}}^4}\right){10}^{-\sqrt{-\ln (1-{(\delta /{\delta }_{\text{ov}})}^4)}}.$$(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 u₀ = (η₀ − η₋)∕Δη 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 $$T_0=\frac{\Delta {\eta }^2}{D_{\text{eff}}}.$$(77)

The distribution of log₁₀T_surv∕T₀ is plotted in Fig. 7 as well as the mean and standard deviation of log₁₀T_surv∕T₀. We give the expression of the distribution in Appendix C. As shown in Fig. 7, the mean of log₁₀T_surv∕P₁ can be well approximated as $$\langle {\log }_{10}\frac{T_{\text{surv}}}{P_1}\rangle ={\log }_{10}\frac{T_0}{P_1}+{\log }_{10}\left(\frac{3}{2}{u}_0^2{(1-u_0)}^2\right),$$(78)

where we choose to normalise the survival time T_surv by the innermost planet period P₁ to include T₀ in the right-hand side expression. Similarly, we approximate the standard deviation as a second-order polynomial as $$\sigma \left({\log }_{10}\frac{T_{\text{surv}}}{P_1}\right)=0.9-2.25u_0(1-u_0).$$(79)

We note that the standard deviation does not depend on the value of log₁₀T₀∕P₁ 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 $$\Delta {\eta }_{\max }=\frac{P}{P+Q}-\frac{P-1}{P+Q}=\frac{1}{P+Q}=\nu ,$$(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ν² ≪ ν 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 (ν₁₂, ν₂₃). For our simplified expression, we take Δη = η_max ≃ 3∕2δ. Furthermore, we assume u₀ = 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 $$\begin{array}{lll}\langle {\log }_{10}\frac{T_{\text{surv}}}{P_1}\rangle \hfill & \simeq \hfill & -{\log }_{10}\left(\frac{16\sqrt{2}A\epsilon M\sqrt{\eta (1-\eta )}}{3}\right)\hfill \\ \hfill & \hfill & +{\log }_{10}\left(\frac{{\delta }^6}{{\delta }_{\text{ov}}^6}\frac{1}{1-{(\delta /{\delta }_{\text{ov}})}^4}\right)\hfill \\ \hfill & \hfill & +\sqrt{-\ln \left(1-{\left(\frac{\delta }{{\delta }_{\text{ov}}}\right)}^4\right)}.\hfill \end{array}$$(81)