Nekhoroshev Estimates for the Survival Time of Tightly Packed Planetary Systems

TL;DR

This paper derives a Nekhoroshev-based estimate for the survival time of tightly packed planetary systems, explaining the exponential growth of stability time with interplanetary spacing and planetary mass through high-order resonances.

Contribution

It introduces a new heuristic formula for planetary system survival time based on Nekhoroshev's theorem, improving understanding of long-term stability beyond previous power-law models.

Findings

The proposed formula matches N-body simulation results well.

Survival time increases exponentially with interplanetary spacing and planetary mass.

High-order resonances drive large deviations from initial conditions.

Abstract

$N$-body simulations of non-resonant tightly-packed planetary systems have found that their survival time (i.e. time to first close encounter) grows exponentially with their interplanetary spacing and planetary masses. Although this result has important consequences for the assembly of planetary systems by giants collisions and their long-term evolution, this underlying exponential dependence is not understood from first principles, and previous attempts based on orbital diffusion have only yielded power-law scalings. We propose a different picture, where large deviations of the system from its initial conditions is due to few slowly developing high-order resonances. Thus, we show that the survival time of the system $T$ can be estimated using a heuristic motivated by Nekhoroshev's theorem, and obtain a formula for systems away from overlapping two-body mean-motion resonances as: $T/P=c_1 \frac{a}{\Delta a} \exp \left(c_2 \frac{\Delta a}{a} /\mu^{1/4}\right)$, where $P$ is the average Keplerian period, $a$ is the average semi major axis, $\Delta a\ll a$ is the difference between the semi major axes of neighbouring planets, $\mu$ is the planet to star mass ratio, and $c_1$ and $c_2$ are dimensionless constants. We show that this formula is in good agreement with numerical N-body experiments for $c_1=5 \cdot 10^{-4}$ and $c_2=8$.

Equations

$$H\left(\mathbf{I},\mathbf{\phi }\right)=H_0\left(\mathbf{I}\right)+\epsilon H_1\left(\mathbf{I},\mathbf{\phi }\right)$$

$$|\mathbf{I}\left(t\right)-\mathbf{I}\left(0\right)|<C_1{\epsilon }^{\frac{1}{2N}}$$

$$t<C_3\exp \left[C_2/{\epsilon }^{\frac{1}{2N}}\right]$$

$$H_1\left(\mathbf{I},\mathbf{\phi }\right)=\sum _k{h}_{\mathbf{k}}\left(\mathbf{I}\right)\exp \left(i\mathbf{k}\cdot \mathbf{\phi }\right).$$

$${\dot{\mathbf{I}}}_1=-{\nabla }_{\mathbf{\phi }}H\approx i\epsilon \sum _{\mathbf{k}}{h}_{\mathbf{k}}\left({\mathbf{I}}_0\right)\mathbf{k}\exp \left(i\mathbf{k}\cdot {\mathbf{\phi }}_0\right)$$

$${\mathbf{I}}_1=\epsilon \sum _{\mathbf{k}}{h}_{\mathbf{k}}\left(I_0\right)\mathbf{k}\frac{\exp \left(i\mathbf{k}\cdot {\mathbf{\phi }}_0\right)}{\mathbf{k}\cdot \mathbf{\omega }}.$$

$$H_1\propto \frac{1}{\sqrt{r_1^2+r_2^2-2r_1{r}_2\cos \theta }}\propto \sum _j{b}_{\frac{1}{2}}^j\left(\frac{r_2}{r_1}\right)\cos \left(j\theta \right)$$

$$b_{\frac{1}{2}}^j\left(\alpha \right)=\frac{1}{\pi }{\int }_0^{2\pi }\frac{\cos \left(j\theta \right)d\theta }{\sqrt{1-2\alpha \cos \theta +{\alpha }^2}}$$

$$b_{\frac{1}{2}}^j\left(\alpha \right)\propto \frac{1}{2}\left|\ln \frac{\Delta a}{a}\right|\exp \left[-\frac{\Delta a}{a}|j|\right]$$

$$|{\mathbf{I}}_1|\propto t\cdot h_{\mathbf{j}}\left({\mathbf{I}}_0\right)\propto t\cdot \exp \left[-\frac{\Delta a}{a}|j|\right]$$

$$\ln T\propto \frac{\Delta a}{a}.$$

$$H\left(\mathbf{\Lambda },\mathbf{\Gamma },\mathbf{\lambda },\mathbf{\varpi }\right)=H_0\left(\mathbf{\Lambda },\mathbf{\Gamma }\right)+\epsilon H_1\left(\mathbf{\Lambda },\mathbf{\Gamma },\mathbf{\lambda },\mathbf{\varpi }\right)$$

$$H\left(\mathbf{\Theta },\mathbf{\theta }\right)=H_0\left(\mathbf{\Theta }\right)+\epsilon H_1\left(\mathbf{\Theta },\mathbf{\theta }\right)$$

$${\theta }_1$$

$${\theta }_2$$

$${\theta }_3$$

$${\theta }_4$$

$$\ln T\propto {\epsilon }^{-1/(2N)}\propto {\left(\frac{M}{m}\right)}^{1/4}.$$

$$T=c_1\frac{a}{\Delta a}\sqrt{\frac{a^3}{GM}}\exp \left(c_2\frac{\Delta a}{a}{\left(\frac{M}{m}\right)}^{1/4}\right)$$

$$T/P\approx 5\times 10^{-4}\cdot \left(\frac{a}{\Delta a}\right)\exp \left(8\left(\frac{\Delta a}{a}\right){\left(\frac{M}{m}\right)}^{1/4}\right),$$

$$\log (T/P)\approx -4.4+\log \left(\frac{a}{\Delta a}\right)+3.5\left(\frac{\Delta a}{a}\right){\left(\frac{M}{m}\right)}^{1/4}.$$

$$\delta v_{\perp }\approx \frac{Gm}{\Delta a^2}\sqrt{\frac{a^3}{GM}}$$

$$\delta v_t\approx e\delta v_{\perp }\approx e\frac{Gm}{\Delta a^2}\sqrt{\frac{a^3}{GM}}$$

$$\delta e^2\approx \frac{\delta L}{L}\approx e\delta v_t/v_k\Rightarrow \delta e\approx \frac{m}{M}\frac{a^2}{\Delta a^2}$$

$$T\approx t_s{\left(\frac{e_f}{\delta e}\right)}^2\approx \sqrt{\frac{a^3}{GM}}{\left(\frac{\Delta a}{a}\right)}^5{\left(\frac{M}{m}\right)}^2.$$

$$\delta \vec{e}=\pm \frac{m}{M}\frac{a^2}{\Delta a^2}\hat{\theta }$$

Full text

editorials, notices — miscellaneous — catalogs — surveys

††journal: ApJ

1 Introduction

According to planet formation theory, small rocky cores are formed in the protoplanetary disk. These cores, as shown in computer simulation, are expected to be of similar mass and at a distance of a few Hill radii from one another (Kokubo and Ida,, 1998). After the disk evaporates, these rocky cores begin to collide and merge, eventually forming the terrestrial planets (Agnor et al.,, 1999; Chambers,, 2001). The rate at which these collisions occur is of tremendous importance to planet formation and the survivability of planetary systems like our own. Collisions play a role, for example, in the Nice model (Gomes et al.,, 2005; Tsiganis et al.,, 2005; Morbidelli et al.,, 2005), although there eccentricities are excited by the gas giants, whereas in this work we consider planetary systems with just terrestrial planets.

The stability of planetary systems in general, and of our solar system in particular, is a long standing problem in physics and mathematics (Laskar,, 1989; Sussman and Wisdom,, 1992; Murray and Holman,, 1999), which has been studied extensively, both analytically (Celletti and Chierchia,, 2007; Efthymiopoulos,, 2008; Féjoz,, 2004; Giorgilli et al.,, 2009, 2017; Locatelli and Giorgilli,, 2007; Robutel,, 1995; Sansottera et al.,, 2010) and numerically (Batygin and Laughlin,, 2008; Hayes,, 2007; Laskar et al.,, 1994; Laskar and Gastineau,, 2009). The chaotic evolution of such systems makes their dynamics difficult to analyze. To illustrate this issue, we have plotted the time evolution of the eccentricity and semi-major axis of one planet in a tightly packed four planet system in Figure 1. The planets end up colliding, but rather than increasing gradually, the eccentricity and semi-major axis seem to be bounded for about $10^{8}$ orbits, and then shoot up at a seemingly arbitrary time.

The major analytic progress in this field is the celebrated Kolmogorov-Arnol’d-Moser (KAM) theory (Kolmogorov,, 1954; Arnol’d,, 1963; Möser,, 1962; Pöschel,, 2009). The main result of this theory is that for small enough perturbations to non degenerate, non resonant integrable Hamiltonian system, the motion lies on an invariant torus with fixed frequencies. One of the supplements to the KAM theory is the Nekhoroshev theory (Nekhoroshev,, 1977, 1979; Pöschel,, 1993; Niederman,, 2012), which states that even if a system does not satisfy the condition for KAM stability, the deviations from the unperturbed solution remain bounded for exponentially long times.

In its original formulation, Nekhoroshev’s estimates apply to non degenerate systems, whereas the planetary problem is degenerate. Various authors have shown that a canonical change of coordinates can circumvent the degeneracy (Niederman,, 1996; Guzzo,, 2007), and applied the Nekhoroshev theory to study the long-term orbital stability of bodies in our Solar system (Cellett and Ferrara,, 1996; Morbidelli and Guzzo,, 1997). Pavlović and Guzzo, (2008) also showed that the Nekhoroshev theorem is applicable to the Koronis and Veritas asteroid families. We use a similar method to study long-term evolution of planetary systems, but focused on compact extra-solar systems for which the relative separation is comparable to their Hill radius. Such compact configurations are expected for rocky cores soon after the protoplanetary disk evaporates, while some systems might retain this compactness for billions of years as it has been revealed by the Kepler sample (Pu and Wu,, 2015).

The stability of tight multi-planet systems with more than two planets has been explored primarily using numerical N-body simulations (e.g. Chambers et al.,, 1996; Yoshinaga et al.,, 1999; Zhou et al.,, 2007; Chatterjee and Ford,, 2014; Smith and Lissauer,, 2009; Funk et al.,, 2010). These simulations reveal that the survival time of the system—the time until the first collision happens—increases very steeply (exponentially) with the initial separation between the orbits (see Pu and Wu,, 2015, for a summary of previous numerical results). This behavior, which has been reproduced in multiple independent studies, is not understood theoretically (see attempts by Chambers et al.,, 1996; Zhou et al.,, 2007; Quillen,, 2011).

In this paper we propose that this survival time can be estimated using Nekhoroshev’s theory. For simplicity, we make a number of assumptions. We assume that all the planets are of equal mass move in co-planar, circular, prograde (i.e. all planets orbit in the same direction), regularly spaced orbits. We also assume that the system does not start out in any first order two body resonance (Deck et al.,, 2012).

The plan of the paper is as follows: in section 2 we present Nekhoroshev’s theory and use it to obtain a survival time for tight planetary systems. In section 3 we compare our results to numerical N-body simulations. Finally, in section 4 we discuss our results.

2 Mathematical Model

2.1 Nekhoroshev Estimates

We consider a nearly integrable Hamiltonian system of the form

[TABLE]

with $N$ degrees of freedom, $\varepsilon\ll 1$ measures the strength of the perturbation, and $\bm{I}$ and $\bm{\varphi}$ are action-angle variables. We note that when $\varepsilon=0$ then all $I$ are constant and all angles increase linearly with time. If $H_{0}\left(\bm{I}\right)$ satisfies certain reasonable conditions on its convexity, and if $\varepsilon$ is small enough, then it can be shown that the system remains stable for exponentially long times. Mathematically, it is stated in the following way: the deviation from the unperturbed solution is bounded from above by some power law of the perturbation strength

[TABLE]

for a duration of time that increases exponentially as the strength of the perturbation decreases

[TABLE]

where $C_{1}$ and $C_{2}$ are constants that do not depend on $\varepsilon$ or time, and $C_{3}$ can, at most, be polynomial in $\varepsilon$. This constraint is known as the Nekhoroshev estimate. We emphasise that this condition guarantees that the system remains bounded within a certain time period, but afterwards the system can either remain stable or deviate further from the unperturbed solution. We note that since these are upper bounds, different authors arrive at different power law indices of $\varepsilon$, depending on their assumptions about the behavior of the Hamiltonian (Guzzo,, 2007).

A complete proof of this constraint is given in Lochak and Neishtadt, (1992), and we will not repeat it here. Instead, we only review some of the key ideas behind the proof. First, we express the perturbing function as a Fourier transform in the angle variables

[TABLE]

where $\bm{k}$ represents a list of integers. Expanding Hamilton’s equation for the momenta to first order in the small parameter $\varepsilon$ we obtain

[TABLE]

where $\bm{I}_{0}$ are constants and $\bm{\varphi}_{0}=\bm{\omega}t+\bm{\varphi}_{i}$, $\bm{\varphi}_{i}$ are constants, and $\bm{\omega}=\nabla_{\bm{I}}H_{0}|_{\bm{I}=\bm{I}_{0}}$. Time integration yields

[TABLE]

We immediately notice a problem with this approach. For arbitrarily large values of the integers $\bm{k}$, the product $\bm{k}\cdot\bm{\omega}$ can be brought arbitrarily close to zero. In fact, a number theory result called Dirichlet’s approximation theorem states that the convergence is better than $|\bm{k}|^{1-n}$, where $n$ is the number of integers and $|\bm{k}|=\sum_{l}|k_{l}|$ is the sum of their magnitudes. Terms for which $\bm{k}\cdot\bm{\omega}=0$ are called resonant terms. Non-resonant terms oscillate in time, while resonant terms grow linearly with time.

On the other hand, the Paley-Wiener theorem tells us that at large indices $\bm{k}$, the Fourier coefficients must decline at least exponentially. Therefore, resonant terms eventually grow and invalidate the perturbative solution, but they take at least an exponentially long time to grow. The remaining challenge is to estimate the index limit $|\bm{k}|$.

2.2 Adaptation to Tightly Packed Planetary Systems

The physical system we are interested in is a small number of planets of the same mass $m$ around a star of mass $M\gg m$. The average distance between the star and the planet is $a$ and the initial separation between neighbouring planets’ semi major axes is $\Delta a\ll a$. All planets are assumed to initially move on circular, coplanar orbits (i.e. both positions and velocities are coplanar, so the motion is restricted to a single plane). We want to estimate the average time it takes for the first planet-planet collision to occur, which it translates into a condition on the growth of the action $I$ (either semimajor axes or eccentricities).

In order to illustrate the evolution of such systems and we show an example in Figure 1 for a four planet system. This example shows that the evolution the actions remain bounded for a long time and after $\sim 10^{8}$ orbits, they break loose and diffuse to the absorbing boundaries (i.e., orbit crossing).

Hereafter, we will use Poincaré variables, the action-angle pairs for a given planet are $\{\Lambda\simeq m\sqrt{GMa},\lambda\}$ and $\{\Gamma\simeq m\sqrt{GMa}e^{2}/2,-\varpi\}$ with $\lambda$ and $\varpi$ the mean longitude and the longitude of pericenter, respectively (Murray and Dermott,, 1999).

2.2.1 Scaling with Spacing $\Delta a/a$

From the discussion in the previous section, we know that this timescale will depend on the exponential decline of the Fourier coefficients. Indeed, the gravitational binding energy between two planets is given by

[TABLE]

where $r_{1,2}$ are distances of the planets from the host star, $\theta$ is the angle between the lines connecting the planets to the star, and

[TABLE]

are the Laplace coefficients, which decline exponentially with the index $b_{\frac{1}{2}}^{j}\left(\alpha\right)\propto\exp\left(-\left(\alpha-1\right)|j|\right)$ (Quillen,, 2011). Thus, for nearly circular orbits:

[TABLE]

(though we will neglect the logarithmic term, as it only introduces a small correction) and from Equation (6) we know that the deviation from the unperturbed solution in the resonant case to grow linearly with time roughly as

[TABLE]

The system survives as long as the action variable does not exceed some critical threshold $|\bm{I}_{1}|<I_{\rm crit.}$, and from Equation (10) the survival time $T$ scales as

[TABLE]

2.2.2 Scaling with the Mass Ratio $\mu=m/M$

We now turn to calculate the dependence of the survival time on the mass ratio. Unlike the dependence with the spacing, this requires specifying the nature of the perturbations that lead to growth of the actions.

We shall start by assuming that the evolution of a single planetary orbit is mainly determined by its two closest neighbors (thus neglecting resonances with non adjacent bodies) so the relevant perturbative Hamiltonian involves the gravitational interactions of only three planets. This choice is justified by the fact that the strength of interactions, which is proportional to the Laplace coefficients and decays exponentially with interplanetary separation (Quillen,, 2011).

Away from overlapping two-body resonances, Quillen, (2011) showed that three-body resonances dominate the chaotic dynamics of low-mass multi-planet system. To first-order in the eccentricities, this Hamiltonian can be written as

[TABLE]

where $\bm{\Lambda}=(\Lambda_{1},\Lambda_{2},\Lambda_{3})$, $\bm{\Gamma}=(\Gamma_{1},\Gamma_{2},\Gamma_{3})$, $\bm{\lambda}=(\lambda_{1},\lambda_{2},\lambda_{3})$, $\bm{\varpi}=(\varpi_{1},\varpi_{2},\varpi_{3})$, and the small parameter is $\varepsilon=\left(m/M\right)^{2}$, since three body resonances do not occur in the linear expansion in the mass ratio. After canonical transformations, this Hamiltonian has four degrees of freedom (since two were eliminated by the introduction of resonant angles) and can be expressed as

[TABLE]

where $\bm{\Theta}$ is a function of the Poincaré actions and $\bm{\theta}$ a linear combination of $\bm{\lambda}$ and $\bm{\varpi}$. In particular, including the modes that are zeroth- and first-order in the eccentricities, one possible combination is (Quillen,, 2011)

[TABLE]

Thus, since this Hamiltonian involves $N=4$ degrees of freedom, Nekhoroshev’s theorem predicts the survival time scales as

[TABLE]

Our argument to get scaling with the mass assumes three-body interactions. Alternatively, there might be other set of perturbations involving two-body interactions that can destabilize the system. Since two-body interactions are linear in the masses, $\epsilon=m/M$, a two-degree-of-freedom model would reproduce the same scaling. We remark, however, that our analysis excludes overlapping two-body resonances and these, as clearly shown in simulations by Obertas et al., (2017) with uniform spacing111Uniform spacing leads to nearly uniform period ratio distribution enhancing the appearance of overlapping two-body resonances., produce wild oscillations in the survival times. In turn, we suspect that the weaker, but much denser (abundant set of angles above that can resonate), three-body resonances would map into a very smooth distribution of $T$ with mutual spacing (high density in $\Delta a/a$) as observed in simulations with unequal spacing (e.g. Chambers et al.,, 1996; Faber and Quillen,, 2007).

In summary, our estimates using three-body interactions and the Nekhoroshev’s theorem results in a scaling with the mass as $\ln T\propto\left(M/m\right)^{1/4}$. This scaling differs from the commonly used Hill scaling (i.e., $\ln T\propto\left(M/m\right)^{1/3}$). We note, however, that previous simulations reported better fits with $\left(M/m\right)^{1/4}$ than that of the Hill’s scaling (e.g., Figure 4 of Chambers et al.,, 1996). A more in-depth study on the mass scaling can better test our results.

2.3 The Prefactor

The remaining component for the survival time is the prefactor which contains a time scale. In this context, the relevant time scale is the time between conjunctions of neighboring planets, or the synodic time. Both planets orbit the star at roughly the Keplerian time $t_{k}\approx\sqrt{\frac{a^{3}}{GM}}$, and the relative difference in periods is proportional to $\Delta a/a$, so the synodic time is given by $t_{s}\approx\frac{a}{\Delta a}t_{k}\approx\frac{a}{\Delta a}\sqrt{\frac{a^{3}}{GM}}$.

Finally we are able to write a complete expression for the survival time of tight planetary systems

[TABLE]

where, again, $c_{1}$ and $c_{2}$ are numerical constants that depend neither on $\Delta a/a$ nor $m/M$, and cannot be determined from scaling arguments. In the next section we compare this formula to numerical results, and using those results we are able to calibrate the numerical coefficients $c_{1}$ and $c_{2}$.

3 Comparison to Simulations

In this section we compare our theoretical prediction in Equation (16) to numerical N-body simulations of tight planetary systems. The results of most simulations done in the past are summarized in Pu and Wu, (2015), but here we simply highlight specific works that explore relevant parts for our discussion.

The first set of simulations we consider was performed in (Zhou et al.,, 2007). They simulated nine planets, spaced with a constant Hill parameter $K\approx 2\frac{a_{i+1}-a_{i}}{a_{i+1}+a_{i}}\left(\frac{M}{m}\right)^{1/3}$. They found a very steep dependence of the survival time on the Hill parameter. They fit this dependence to a power law in the Hill parameter, and infer $d\ln T/d\ln K\approx 20$. They also varied the mass ratio $\mu=m/M$, and find $d\ln T/d\ln\mu=-0.27$. They developed an analytic theory for diffusion in eccentricity that reproduces the dependence on the mass ratio $\mu$, but whose dependence on the Hill parameter is too shallow ($K^{5}$ instead of $K^{20}$). A discussion of this theory and its shortcomings is given in the appendix.

A similar suite of simulations was performed by Rice et al., (2018). They explored a wider range of Hill parameters than Zhou et al., (2007), and they fit the numerical result to an exponential relation between the survival time and Hill parameter. A similar fit to an exponential was obtained by Petrovich, (2015) for the case of two eccentric planets. From the results of Rice et al., (2018) we are able to constrain the remaining free parameters $c_{1}$ and $c_{2}$ to obtain a final form for the survival time of tight planetary systems

[TABLE]

or more commonly expressed as

[TABLE]

This form is similar to the expression found by Faber and Quillen, (2007), which has a coefficient of 3.7 accompanying the term $\frac{\Delta a}{a}\left(\frac{M}{m}\right)^{1/4}$ instead of 3.5, and logarithmic dependence on the mass ratio, not the spacing.

A comparison between the numerical results of Rice et al., (2018) and Equation (17) can be seen in figure 2. We note that Rice et al., (2018) argue that the argument in the exponent should be proportional to the Hill parameter $\frac{\Delta a}{a}\left(\frac{M}{m}\right)^{1/3}$, unlike our results. However, we note that some numerical experiments agree well with the $\frac{\Delta a}{a}\left(\frac{M}{m}\right)^{1/4}$ scaling (Chambers,, 2001; Faber and Quillen,, 2007). Similarly, this scaling is born out from theoretical works estimating the regions for the onset of chaos either from three-body resonance overlap (Quillen,, 2011) or two-body resonance overlap for eccentric planets (Hadden and Lithwick,, 2018). A more exhaustive study on the mass scaling of instability times will shed light on this issues.

Finally, Obertas et al., (2017) ran another suite of simulations, with a much higher resolution in the interplanetary spacing compared to all previous works. They found a general trend consistent with the previous studies, and also large and rapid variations in the survival time as a function of $\Delta a/a$. They show that the locations of these oscillations coincide with two-body mean-motion resonances. From the discussion in section 2.1 we can understand why our model fails for systems close to a resonance. In such cases, the Fourier series for the perturbing potential must be truncated at a indices smaller than $|\bm{k}|\propto\varepsilon^{-1/(2N)}$. We note that a resonance is a double edged sword, in the sense that it can rapidly change the parameters of the system, but only within a lower dimensional subspace of parameter space. If a collision is possible within this subspace, then it can happen on timescales much shorter than Nekhorosev’s estimates. On the other hand, if a collision is not possible within this parameter space, then a resonance can prolong the lifetime of the system considerably. Both of these effects can be seen in the results of Obertas et al., (2017).

4 Discussion

In this work we study the survival time (time to first close encounter) of tight planetary systems (mutual spacing $\lesssim 10$ mutual Hill radii). We use Nekhoroshev’s stability limit to derive an analytic formula for the survival time as a function of interplanetary spacing and planet-to-star mass ratio. This formula is defined up to two dimensionless numbers and it reproduce the exponential scaling with interplanetary spacing and mass ratio observed several numerical simulations. By calibrating these constants we arrive at a complete expression for the survival time of tight planetary systems (Equation 17).

We note that most previous attempts to explain the steep dependence of the survival time on the separation focused on estimating the chaotic diffusion coefficients, which only led to polynomial growth of the survival time. To our knowledge, this is the first work attempting to explain the numerical results with an exponential, although the idea has been suggested in previous works (Zhou et al.,, 2007; Quillen,, 2011). In our picture, a high order resonant mode grows linearly with time, but has a tiny prefactor, so it takes a long time for it to exceed low order non resonant modes.

We note that while this work provides a useful tool for estimating the survival time for tight planetary systems, our arguments are heuristic and do not reveal the nature of the chaotic interactions leading to instability. We suggest that three-body resonances are suspect, but have not explored the role of two-body resonances. Based on recent numerical simulations showing large deviations from our results in the vicinity of overlapping mean-motion resonances, and it is also possible that these resonant systems follow a different scaling law, possibly linked to rapid diffusion. Finally, in this work we have assume a highly idealized system in which all planets have equal masses and uniformly spaced in coplanar orbits (i.e. both positions and velocities are coplanar, so the motion is restricted to a single plane). All these effects merit further exploration.

Acknowledgments

We would like to thank Norm Murray and Alysa Obertas for the useful feedback. AY is supported by the Vincent and Beatrice Tremaine Fellowship. C.P. acknowledges support from the Gruber Foundation Fellowship and Jeffrey L. Bishop Fellowship. This work made use of the matplotlib python package (Hunter,, 2007).

Appendix A Diffusion timescale of the eccentricity vector: understanding the power-law scaling using the impulse approximation

In this section we discuss the analytic theory for diffusion in eccentricity developed in (Zhou et al.,, 2007). The original derivation relies on a Hamiltonian formalism, but we can reproduce the results using the impulse approximation. In this approximation, we neglect all interaction between the planets until conjunction, and then assume all exchange of energy and momentum is instantaneous. Again, we assume initially co-planar planets on circular orbits. The mass of the planets is $m$, the mass of the star is $M$, the average semi major axis is $a$, and the separation between the planets is $\Delta a\ll a$.

Both planets move at roughly the Keplerian velocity $v_{k}\approx\sqrt{GM/a}$, and the relative velocity between them is $\Delta v_{k}\approx\frac{\Delta a}{a}\sqrt{\frac{GM}{a}}$. At conjunction, the acceleration between the two planets due to mutual attraction is $Gm/\Delta a^{2}$. The duration of the interaction is roughly given by the time during which the distance between the two planets is $\Delta a$, namely $\Delta t\approx\Delta a/\Delta v_{k}\approx t_{k}\approx\sqrt{\frac{a^{3}}{GM}}$. Hence, due to the interaction, the planets receive a kick velocity normal to the direction of velocity difference

[TABLE]

In the case of a perfectly circular orbit, this kick velocity does not change the angular momentum, since the velocity is parallel to the radial direction. If, however, the orbit has an eccentricity $e$, and conjunction happens at some random angle, then the angle between the kick velocity and the radial direction would be proportional to the eccentricity $e$, and the kick velocity would have a small component in the tangential direction

[TABLE]

Depending on the true anomalies at conjunction, this kick could either be prograge or retrograde. The change in angular momentum is given by $\Delta L\approx a\delta v_{t}$, and the change in eccentricity is

[TABLE]

where $L\approx\sqrt{GMa}$ is the orbital angular momentum per unit mass. For the orbits to cross, the eccentricity has to grow to $e_{f}\approx\Delta a/a$. Since we assume a random walk, the number of conjunction needed is $\left(e_{f}/\delta e\right)^{2}$. Since the time between conjunctions is $t_{s}\approx t_{k}a/\Delta a$, the survival time according to this model is

[TABLE]

Thus we reproduce the expression for survival time obtained by Zhou et al., (2007). This expression does not reproduce the observed exponential dependence of $T$ with the spacing and the mass ratio.

The problem with this model is that it assume a random angle between planets at conjunction, whereas the planets quickly align their elliptical orbits. To see why this is the case, let us use the impulse approximation to calculate how the direction of periapse changes. We already calculate the kick velocity in the radial direction in equation A1. This component is always in the same direction, and does not depend on the eccentricity. We recall that the eccentricity vector is given by $\vec{e}=\frac{\vec{v}\times\vec{L}}{GM}-\hat{r}$, where $\vec{v}$ is the orbital velocity and $\vec{L}$ is the angular momentum. The change in this vector is given by

[TABLE]

where $\hat{\theta}$ represents the tangential direction. In this case the length of the eccentricity vector does not change, only its direction. We note that since the inferior and superior planets get opposite radial kick velocities but have roughly the same angular momentum vector, then the changes in the eccentricity vectors have opposite signs. If both planets had the same eccentricity vector, this would cause them to rotate in opposite ways. On the other hand, if both planets had equal but opposite eccentricity vectors, then this would cause both to rotate in the same way, thus preserving their anti-alignment. This anti-alignment can be maintained even for a large mass ratio between the planets, since in this case the more massive planet would have a much smaller eccentricity, which can be easily rotated by the smaller planet.