First order resonance overlap and the stability of close two planet systems

Katherine M. Deck¹ ³ , Matthew Payne² , Matthew J. Holman²

Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA 02139

Corresponding author: kdeck@mit.edu

Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138

Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138

Abstract

Motivated by the population of observed multi-planet systems with orbital period ratios 1<P₂/P₁ ≲ 2, we study the long-term stability of packed two planet systems. The Hamiltonian for two massive planets on nearly circular and nearly coplanar orbits near a first order mean motion resonance can be reduced to a one degree of freedom problem (Sessin & Ferraz-Mello, 1984; Wisdom, 1986; Henrard et al., 1986). Using this analytically tractable Hamiltonian, we apply the resonance overlap criterion to predict the onset of large scale chaotic motion in close two planet systems. The reduced Hamiltonian has only a weak dependence on the planetary mass ratio m₁/m₂, and hence the overlap criterion is independent of the planetary mass ratio at lowest order. Numerical integrations confirm that the planetary mass ratio has little effect on the structure of the chaotic phase space for close orbits in the low eccentricity (e≲ 0.1) regime. We show numerically that orbits in the chaotic web produced primarily by first order resonance overlap eventually experience large scale erratic variation in semimajor axes and are therefore Lagrange unstable. This is also true of the orbits in this overlap region which satisfy the Hill criterion. As a result, we can use the first order resonance overlap criterion as an effective stability criterion for pairs of observed planets. We show that for low mass (≲ 10M_⊕) planetary systems with initially circular orbits the period ratio at which complete overlap occurs and widespread chaos results lies in a region of parameter space which is Hill stable. Our work indicates that a resonance overlap criterion which would apply for initially eccentric orbits likely needs to take into account second order resonances. Finally, we address the connection found in previous work between the Hill stability criterion and numerically determined Lagrange instability boundaries in the context of resonance overlap.

Subject headings:

celestial mechanics - chaos - planets and satellites: dynamical evolution and stability †

1. INTRODUCTION

Observational results from the Kepler mission have revealed a large population of multi-planet systems with nearly circular and nearly coplanar orbits (Batalha et al., 2013; Fabrycky et al., 2012; Fang & Margot, 2012). The distribution of the orbital separation between adjacent pairs of planets, expressed either in terms of the observed period ratios or in terms of the inferred “dynamical spacing” (mutual Hill radii), encodes a significant amount of information regarding how these systems formed and dynamically evolved (Fabrycky et al., 2012; Fang & Margot, 2013). One interesting question regarding the orbital separations of planets is whether or not the inferred dynamical spacing distribution is dictated by orbital stability requirements or by planetary formation processes. In order to address this, it is important to understand in a theoretical manner the long-term stability of close two-planet systems.

Addressing analytically the question of how small an orbital separation two planets can have while remaining “long-lived” is challenging, as even systems of two planets can be chaotic, and chaotic orbits cannot be described in terms of simple functions. In general, to make analytic progress, researchers must turn to more global methods, rather than the study of individual initial conditions.

It has been proven that the conservation of angular momentum and energy can constrain the motion in two planet systems in a way such that crossing orbits never develop (Marchal & Bozis, 1982; Milani & Nobili, 1983). Orbits which do not cross will never lead to collisions between planets or strong gravitational scattering. As a result, one can use the integrals of the motion to prove that a two planet system is collisionally (Hill) stable. However, planetary orbits which fail this “Hill criterion” and evolve in the region of phase space where crossing orbits are possible are not guaranteed to be unstable; satisfying the criterion is a sufficient but not necessary condition for collisional stability.

Additionally, the Hill criterion gives no information about the long-term behavior of orbits that satisfy the criterion. Even though two bodies may never suffer strong encounters, repeated interactions can lead to a net transfer of angular momentum and energy between the two bodies that result in large, erratic, variations in the orbits. In particular, the semimajor axes can change considerably, and in some cases this leads to ejection of the outer body or a collision between the inner and central bodies. These orbits are referred to as Lagrange (or Laplace) unstable, even though they are stable in the Hill sense. Our definition of Lagrange instability is not restricted to cases in which ejection of the outer planet or collisions between the central body and the inner planet occur.

Systems which are protected from this type of behavior have more constrained orbits that we call Lagrange long-lived. It has been found numerically that orbits which marginally satisfy the Hill criterion exhibit Lagrange instabilities, but that orbits which satisfy the Hill criterion by a larger amount are generally Lagrange long-lived (Barnes & Greenberg, 2006, 2007; Mudryk & Wu, 2006; Deck et al., 2012). These studies found that the transition between Lagrange unstable orbits and Lagrange long-lived orbits, as the “distance” from the Hill boundary grows, occurs very abruptly, over a narrow range in the orbital parameters. A theoretical explanation for this behavior is lacking, though the specific location of the transition point in terms of the planetary mass and the orbital parameters has been explored numerically (Kopparapu & Barnes, 2010; Veras & Mustill, 2013).

Gladman (1993) numerically studied the effectiveness of the Hill criterion as a stability criterion and found that for low mass planets on initially circular orbits, the Hill criterion was effectively necessary for stability, while for more eccentric cases there seemed to be many orbits which failed the criterion but appeared to be long-lived. He also found that there existed a chaotic region of phase space which extended past the Hill boundary. A limited number of numerical experiments with equal mass planets indicated that the width of this region scaled as εₚ²⁄⁷, where εₚ is the total mass of the planets, relative to the star. The specific power of 2/7 is motivated theoretically by analytic work on first-order resonance overlap in the circular restricted three body problem (Wisdom, 1980; Duncan et al., 1989). However, it is not known how this result changes in the case of two massive planets on low eccentricity orbits, with arbitrary mass ratio, or if first order resonance overlap accounts for the chaotic zone discovered by Gladman.

In this work, we study the stability of two-planet systems by extending the resonance overlap criterion to the full planetary problem. We seek to explain Gladman’s numerical results and investigate the intriguing correlation between Hill stable and Lagrange long-lived orbits and the role resonance overlap plays in this relationship.

We also approach the problem “globally” in that we do not evolve individual orbits for very long times. Instead, we focus on predicting where the chaotic regions of phase space lie. A single chaotic orbit will explore densely the entire extent of the chaotic zone available to it. As a result, if a chaotic zone encompasses regions of phase space where Lagrange-type instabilities can set in, any chaotic orbit within that region is Lagrange unstable. An example of this would be a chaotic zone which extended over a large range of period ratios. The problem of long-term stability in planetary systems can therefore be approached by identifying the regions of phase space where widespread chaos will be present, rather than actually following the long-term evolution of single trajectories.

In order to identify chaotic regions of phase space we make use of the resonance overlap criterion of Chirikov (1979); see also Walker & Ford (1969). This heuristic criterion provides an intuitive way to predict where in phase space global chaos appears¹, and has been used successfully to estimate regions of chaos in the restricted three body problem, in close three planet systems, in widely spaced two planet systems, and in explaining ejection of planets orbiting in binary star systems; see for example Lecar et al. (2001); Mardling (2008); Mudryk & Wu (2006); Quillen (2011).

He or she who “desires to reach a cherished islet of stability in the violent (and stochastic!) sea of nonlinear oscillations should not rely upon the beacons only”, where “the beacons” are the rigorous mathematical theorems of nonlinear dynamics (Chirikov, 1979).

We apply the resonance overlap criterion to a system with two massive planets on close orbits (P₂/P₁ ≲ 2) which are nearly circular ( e≲ 0.1) and nearly coplanar. In order to do so, we must first identify the dominant resonances in this regime. These are the first order mean motion resonances, which are important when the planetary period ratio P₂/P₁ ∼(m+1)/m≤2, where m is a positive integer. The distance in period ratio between the first order resonances shrinks as the planetary period ratio grows closer to unity, but the widths of the resonances (in terms of period ratio) do not shrink as quickly. Consequently, the first order resonances must overlap as the period ratio shrinks to unity, leading to chaotic motion.

In Section 2 we show how to reduce the first order resonance problem, at first order in the planetary eccentricities, to a one degree of freedom problem. We derive the resonance overlap criterion in Section 3. In Section 4, we compare the predicted resonance widths and overlap criterion to results from numerical integrations. In Section 5 we discuss the implications of resonance overlap for the long-term stability of close two planet systems.

2. ANALYTIC REDUCTION OF THE HAMILTONIAN

Here we reduce the full Hamiltonian for two planets orbiting near a first order mean motion resonance with nearly circular and nearly coplanar orbits to a one-degree of freedom system with a single free parameter (Sessin & Ferraz-Mello, 1984). We follow a sequence of canonical transformations originally used by Wisdom (1986) and Henrard et al. (1986) (and recently by Batygin & Morbidelli (2013)). Although this reduction of the Hamiltonian is not original work of ours, we include it here for consistency of notation throughout the text and to familiarize the reader.

The Hamiltonian governing the dynamics of a system of two planets of mass mᵢ orbiting a much more massive star of mass m_, written in Jacobi coordinates rᵢ and momenta pᵢ, takes the form

$${H} ={H}_{0}+{H}_{1}\qquad {H}_{0} ={H}_{Kepler1}+{H}_{Kepler2}\qquad {H}_{Kepler1} =\frac{{p}_{1}^{2}}{2\tilde{m}_{1}}-\frac{G\tilde{M}_{1}\tilde{m}_{1}}{\big| \mathbf{{r}_{1}}\big\rvert }\qquad {H}_{Kepler2} =\frac{{p}_{2}^{2}}{2\tilde{m}_{2}}-\frac{G\tilde{M}_{2}\tilde{m}_{2}}{\big| \mathbf{{r}_{2}}\big\rvert }$$ $${H}_{1} =-G{m}_{1}{m}_{2}\left(\frac{1}{\big| \mathbf{{r}_{1}}-\mathbf{{r}_{2}}\big\rvert }-\frac{\mathbf{{r}_{2}}\cdot\mathbf{{r}_{1}}}{\big| \mathbf{{r}_{2}}\big\rvert ^{3}}\right)+O(\epsilon^{2}),$$

(1)

where M_2=m_(m_+m_1+m_2)/(m_+m_1), M_1=(m_+m_1), ε=max(m_i/m_), and m̃ᵢ=mᵢ+O(ε) denote Jacobi masses. A subscript of 1 refers to the inner planet, a subscript of 2 refers to the outer planet. The interaction potential between the planets takes the form of the disturbing function (Murray & Dermott, 1999). We set M_i=m_ and also ignore the difference between Jacobi and physical masses. In the Keplerian piece H₀, this approximation corresponds to a change in the mean motions of the planets by order ε, but this is negligible.

Instead of Cartesian coordinates and momenta, we use the Poincare canonical variables

$$\Lambda_{i} =m_{i}\sqrt{Gm_{\star}a_{i}}\qquad P_{i} =m_{i}\sqrt{Gm_{\star}a_{i}}(1-\sqrt{1-e_{i}^{2}})\qquad G_{i} =m_{i}\sqrt{Gm_{\star}a_{i}(1-e_{i}^{2})}(1-\cos{I_{i}})$$

(2)

and their conjugate angles λ_i=_i+M_i, p_i=-_i, and gᵢ=-Ωᵢ. Here a, e, I, λ, , and M denote the (Jacobi) semimajor axis, eccentricity, inclination, mean longitude, longitude of ascending node, longitude of periastron and mean anomaly of the osculating orbit. When expressed in terms of these variables, the interaction Hamiltonian H₁ takes the form of a power series in eccentricities and inclinations, with each term proportional to a periodic function of a linear combination of the orbital angles (except the term proportional to cos(jλ₂-jλ₁) with j=0).

For systems far from mean motion resonances, these periodic terms are either short period (they average out on timescales comparable to the orbital periods) or secular (varying on long timescales; these terms are independent of the fast angles λᵢ). When the mean motions of the planet form a rational ratio, n₂/n₁ ∼ m/(m+q), with m and q positive integers, the system is near a q-th order mean motion resonance and all periodic terms in the interaction Hamiltonian which are functions of the combination (m+q)λ₂-mλ₁ will be slowly varying (relative to the mean longitudes of the planets λᵢ).

We are interested in the long-term dynamics of two planet systems, and so we will remove all of the short period terms by averaging over the mean longitudes of the planets. There is a canonical transformation between the full Hamiltonian and the averaged Hamiltonian, part of which is derived in the Appendix, which can be important for comparing the analytic and numerical results (see Section 4), though it is not necessary to perform this step explicitly for the following analysis. After averaging, only secular and resonant terms remain in the Hamiltonian.

The secular terms appear only at second order in the eccentricities and the inclinations or higher. The q-th order resonant terms appear at order q in eccentricities and at order q or higher in the inclinations. The first order mean motion resonance terms, with q=1, or n₁/n₂=P₂/P₁ ∼(m+1)/m, therefore appear at first order in eccentricities (though they appear at second order in inclinations due to symmetries of the problem). As a result, they are the most important resonances to include when considering resonance overlap of nearly circular and nearly coplanar orbits - other mean motion resonances have smaller widths and contribute less to overlap. Throughout the rest of the paper, an m with no subscript refers to a particular first order resonance near the period ratio P₂/P₁ ∼(m+1)/m.

Truncating the expansion of H₁ at first order in I removes all (I,Ω) dependence (and the secular terms), so that the pair (Gᵢ,gᵢ) do not appear in the Hamiltonian. Because the Hamiltonian containing only terms at first order in e or I is the same as the coplanar Hamiltonian at first order in e, any results derived below apply to slightly non-coplanar orbits (as long as I is sufficiently small such that I² terms are negligible).

After performing these steps the resulting Hamiltonian is:

$$H =H_{0}(\mathbf{\Lambda})+\epsilon_{1}H_{1}(\mathbf{\Lambda},\mathbf{P},\mathbf{\lambda},\mathbf{p}),$$

(3)

where

$$H_{0} =-\frac{\mu_{1}}{2\Lambda_{1}^{2}}-\frac{\mu_{2}}{2\Lambda_{2}^{2}},$$

(4)

$$\epsilon_{1} =\frac{m_{1}}{m_{\star}},$$

(5)

and

$$H_{1} =-\frac{\mu_{2}}{\Lambda_{2}^{2}}\bigg[ f_{0,1}(\alpha)+\qquad f_{m+1,27}(\alpha)\sqrt{\frac{2P_{1}}{\Lambda_{1}}}\cos{((m+1)\lambda_{2}-m\lambda_{1}+p_{1})}+\qquad f_{m+1,31}(\alpha)\sqrt{\frac{2P_{2}}{\Lambda_{2}}}\cos{((m+1)\lambda_{2}-m\lambda_{1}+p_{2})}\bigg]$$

(6)

and μ_i=G^2m_^2m_i^3, α=a₁/a₂, √(2P/Λ)=e(1+O(e²)), and fₘ₊₁,₂₇, fₘ₊₁,₃₁ and f₀,₁ are functions of Laplace coefficients (Murray & Dermott, 1999, p. 539-556). These can be written as

$$f_{j,1}(\alpha) =\frac{1}{2}A_{j}(\alpha)\qquad f_{j,27}(\alpha) =\frac{1}{2}(-2j-\alpha\frac{d}{d\alpha})A_{j}(\alpha)\qquad f_{j,31}(\alpha) =\frac{1}{2}(-1+2j+\alpha\frac{d}{d\alpha})A_{j-1}(\alpha)\qquad A_{j}(\alpha) =\frac{1}{\pi}\int_{0}^{2\pi}\frac{\cos{(j\phi)}}{\sqrt{1-2\alpha\cos{\phi}+\alpha^{2}}}d\phi.$$

(7)

We could equivalently use ε₂ instead of ε₁ as our small parameter, the overall coefficient of H₁ would just change to -μ₁/Λ₂². Note that when m=1, f₂,₃₁ must be modified to include a contribution from the indirect terms of -2α.

We approximate these functions as constants, evaluated at the nominal resonance location, αᵣₑₛ≡(m/(m+1))²⁄³. These values are listed in Table 1 for the first order resonances with 1≤m≤13.

Table 1The functions fₘ₊₁,₂₇ and fₘ₊₁,₃₁ evaluated at the nominal resonance location, with α=αᵣₑₛ(m). As m→∞, fₘ₊₁,₂₇ →-fₘ₊₁,₃₁.

We find numerically that the functions fₘ₊₁,₂₇ and fₘ₊₁,₃₁, for m≥2, evaluated at αᵣₑₛ(m) and treated as functions of m, are well fit by straight lines. Using a range of m from 2 to 150 yields fits of

fₘ₊₁,₂₇ ≈-0.46-0.802m fₘ₊₁,₃₁ ≈ 0.87+0.802m

(8)

The maximum fractional deviation in the values for fₘ₊₁,₂₇ and fₘ₊₁,₃₁ using the best fit line is 1.9% and 0.4%, respectively. It is clear that the ratio of the two quickly approaches -1 as α→ 1. This is expected for close orbits as in the limit α→ 1, f₂₇ and f₃₁ are dominated by the derivative in their expressions (Quillen, 2011). Since the coefficient of the derivative contribution is the same (up to a sign) in both expressions, f₂₇ ∼-f₃₁ for close orbits (high m).

Only one linear combination of λ’s appears in the Hamiltonian (3), implying that we can reduce the number of degrees of freedom by one. Let

$$\theta =(m+1)\lambda_{2}-m\lambda_{1}\qquad \theta_{1} =\lambda_{1}$$

(9)

be the new angles, and Θ and Θ₁ their corresponding actions. The new actions are found using the generating function

$$F_{2}=\Theta((m+1)\lambda_{2}-m\lambda_{1})+\Theta_{1}\lambda_{1}$$

(10)

as follows:

$$\Lambda_{1} =\frac{\partial F_{2}}{\partial\lambda_{1}}=\Theta_{1}-m\Theta\qquad \Lambda_{2} =\frac{\partial F_{2}}{\partial\lambda_{2}}=(m+1)\Theta,$$

(11)

which we can invert as

$$\Theta =\Lambda_{2}/(m+1)\qquad \Theta_{1} =\frac{m}{m+1}\Lambda_{2}+\Lambda_{1}.$$

(12)

We’ll use the canonical polar variables:

$$x_{i} =\sqrt{2P_{i}}\cos{(p_{i})}\qquad y_{i} =\sqrt{2P_{i}}\sin{(p_{i})},$$

(13)

where xᵢ are momenta and yᵢ are coordinates ( {yᵢ,xᵢ}=1, with {·,·} denoting a Poisson Bracket). This sequence of time independent canonical transformations leads to the Hamiltonian

$$H =H_{0}(\Theta_{1}-m\Theta,(m+1)\Theta)\qquad {}-\epsilon_{1}\frac{\mu_{2}}{((m+1)\Theta)^{2}}\bigg[ f_{0,1}+\qquad \frac{f_{m+1,27}}{\sqrt{\Theta_{1}-m\Theta}}(x_{1}\cos{\theta}-y_{1}\sin{\theta})+\qquad \frac{f_{m+1,31}}{\sqrt{(m+1)\Theta}}(x_{2}\cos{\theta}-y_{2}\sin{\theta})\bigg] .$$

(14)

Note that Θ₁ is an integral of the flow generated by Hamiltonian (2).

We express all actions in units of Θ₁, and the Hamiltonian in units of μ₂/Θ₁². We choose a new independent variable t̂=μ₂/Θ₁³ t, so that when using Ĥ=H/(μ₂/Θ₁²) as the Hamiltonian function Hamilton’s equations will still hold. Hats denote unitless variables: x̂ᵢ=xᵢ/√(Θ₁),ŷᵢ=yᵢ/√(Θ₁), and Θ/Θ₁=Θ̂. We define μ₁/μ₂=(m₁/m₂)³≡ζ³.

In these units, the Hamiltonian takes the form

$$\hat{H} =-\frac{\zeta^{3}}{2(1-m\hat{\Theta})^{2}}-\frac{1}{2((m+1)\hat{\Theta})^{2}}\qquad {}-\epsilon_{1}\bigg[ \frac{f_{0,1}}{((m+1)\hat{\Theta})^{2}}+$$

(15)

where we have defined

$$\delta_{1} =\frac{1}{((m+1)\hat{\Theta})^{2}}\frac{f_{m+1,27}}{\sqrt{1-m\hat{\Theta}}}\qquad \delta_{2} =\frac{1}{((m+1)\hat{\Theta})^{2}}\frac{f_{m+1,31}}{\sqrt{(m+1)\hat{\Theta}}}.$$

(16)

The center of the resonance corresponds approximately to Θ̂=Θ̅, where Θ̅ is implicitly defined by the expression

$$\frac{d\theta}{d\hat{t}} =\frac{\partial\hat{H}}{\partial\hat{\Theta}}\approx\frac{\partial\hat{H}_{0}}{\partial\hat{\Theta}}=0,$$

or

$$\frac{-m\zeta^{3}}{(1-m\bar{\Theta})^{3}}+ \frac{1}{(m+1)^{2}\bar{\Theta}^{3}}=0,$$

(17)

where we have ignored the order ε₁ terms. Including the effects of the f₀,₁ term in solving for Θ̅ is straightforward, but it would only shift the nominal resonance location by order ε₁, and we neglect it. Solving for Θ̅ yields

$$\bar{\Theta} =\frac{\alpha_{res}}{m(\alpha_{res}+\zeta)}\qquad \hat{\Lambda}_{2,res} =(m+1)\bar{\Theta}\qquad \hat{\Lambda}_{1,res} =1-m\bar{\Theta}.$$

(18)

For orbits near the resonance, δΘ=Θ̂-Θ̅ is small, and we can expand the Hamiltonian around Θ̅. In this case, the effective Hamiltonian K=Ĥ-Ĥ₀(Θ̅)+ε₁ f₀,₁/(m+1)²/Θ̅² is

$$K =\frac{\beta}{2}\delta\Theta^{2}-\epsilon_{1}\times\qquad \bigg[ (\bar{\delta}_{1}\hat{x}_{1}+\bar{\delta}_{2}\hat{x}_{2})\cos{\theta}-(\bar{\delta}_{1}\hat{y}_{1}+\bar{\delta}_{2}\hat{y}_{2})\sin{\theta}\bigg] ,$$

(19)

where

$$\bar{\delta}_{1} =\frac{1}{((m+1)\bar{\Theta})^{2}}\frac{f_{m+1,27}}{\sqrt{1-m\bar{\Theta}}}\qquad \bar{\delta}_{2} =\frac{1}{((m+1)\bar{\Theta})^{2}}\frac{f_{m+1,31}}{\sqrt{(m+1)\bar{\Theta}}}\qquad \beta=\frac{\partial^{2}\hat{H}_{0}}{\partial\hat{\Theta}^{2}}\bigg| _{\hat{\Theta}=\bar{\Theta}} =-3\bigg( \frac{m^{2}\zeta^{3}}{(1-m\bar{\Theta})^{4}}+\frac{1}{(m+1)^{2}\bar{\Theta}^{4}}\bigg)\qquad {}=-3\frac{m^{2}(\alpha_{res}+\zeta)^{5}}{\zeta\alpha_{res}}.$$

(20)

Note that δΘ is the canonical momentum conjugate to θ, and can be rewritten as

$$\delta\Theta =\frac{\zeta}{m}(1-s)\frac{1}{\alpha_{res}+\zeta+s\cdot\zeta(1+\zeta/\alpha_{res})}$$

(21)

with s²=α/αᵣₑₛ. We see that when δΘ=0, or Θ=Θ̅, α=αᵣₑₛ.

We have neglected terms of order δΘ³, ε₁ δΘ, and higher in the Hamiltonian (2). If the first two neglected terms are of roughly the same magnitude, this would imply that δΘ∼√(ε₁) and that both terms in K are the same magnitude. Note that the expansion about the resonance location is actually necessary to do - in order for the following sequence of canonical transformations to work out, the functions δ₁ and δ₂ must be constants (we require Θ̂ be replaced by Θ̅).

Then we define the new canonical set

$$r_{1} =\frac{\bar{\delta}_{1}\hat{x}_{1}+\bar{\delta}_{2}\hat{x}_{2}}{\bar{\delta}}\qquad s_{1} =\frac{\bar{\delta}_{1}\hat{y}_{1}+\bar{\delta}_{2}\hat{y}_{2}}{\bar{\delta}}\qquad r_{2} =\frac{\bar{\delta}_{2}\hat{x}_{1}-\bar{\delta}_{1}\hat{x}_{2}}{\bar{\delta}}\qquad s_{2} =\frac{\bar{\delta}_{2}\hat{y}_{1}-\bar{\delta}_{1}\hat{y}_{2}}{\bar{\delta}}.$$

(22)

The rᵢ are momenta, and sᵢ their canonical conjugate coordinates, and δ̅=√(δ̅₁²+δ̅₂²). Substituting Equation (2) into the Hamiltonian (2) yields:

$$K =\frac{\beta}{2}\delta\Theta^{2}-\epsilon_{1}\bar{\delta}\times\bigg[ r_{1}\cos{\theta}-s_{1}\sin{\theta}\bigg]$$

(23)

Note that r₂ and s₂ are conserved. Next, let r₁=√(2Ψ)cosψ and s₁=√(2Ψ) sinψ. Ψ is proportional to e²:

$$\Psi =\frac{1}{2}(r_{1}^{2}+s_{1}^{2})\qquad {}=\frac{\zeta}{\alpha_{res}+\zeta}\frac{1}{2(R^{2}+\zeta\sqrt{\alpha_{res}})}\times\qquad \bigg[ R^{2}e_{1}^{2}+e_{2}^{2}-2Re_{1}e_{2}\cos{(\Delta\varpi)}\bigg]\qquad {}=g_{m}R^{2}\bigg[ e_{1}^{2}+e_{2}^{2}/R^{2}-2/Re_{1}e_{2}\cos{(\Delta\varpi)}\bigg]$$

(24)

with R=| fₘ₊₁,₂₇/fₘ₊₁,₃₁ | and gₘ implicitly defined. If we define: e _1 e_1(_1,_1) and the same for e₂, Ψ∝(e₁-(1)/(R)e₂)², so it is a measure of the total eccentricity, and we will refer to the quantity

$$\sigma \equiv\sqrt{e_{1}^{2}+e_{2}^{2}/R^{2}-2/Re_{1}e_{2}\cos{(\Delta\varpi)}}=\bigg[ \frac{\Psi}{g_{m}R^{2}}\bigg] ^{1/2}$$

(25)

as the weighted eccentricity of the system. The angle = s_1/r_1 is a generalized longitude of pericenter. In the special case when e₂=0, =π-_1, and when e₁=0, =-_2. These new “eccentricity” vector components (rᵢ,sᵢ) are obtained from a linear transformation of the original (xᵢ,yᵢ). As pointed out by Batygin & Morbidelli (2013), the transformation is a rotation - it preserves length, and indeed you can show that Ψ+Ψ₂=(P₁+P₂)/Θ₁, where Ψ=1/2(r₁²+s₁²) and Ψ₂=1/2(r₂²+s₂²). Therefore Ψ+Ψ₂ is proportional to the angular momentum deficit.

In terms of Ψ and ψ, the Hamiltonian K becomes

$$K =\frac{\beta}{2}\delta\Theta^{2}-\epsilon_{1}\bar{\delta}\sqrt{2\Psi}\cos{(\theta+\psi)}.$$

(26)

And now we will choose new angles ϕ=θ+ψ and γ=-θ. Using the generating function

F₂=Φ(θ+ψ)-Γθ,

(27)

the new set of actions is

Ψ=Φ δΘ=-Γ+Φ.

(28)

Then the Hamiltonian becomes

$$K =\frac{1}{2}{\beta}(\Phi-\Gamma)^{2}-\epsilon_{1}\bar{\delta}\sqrt{2\Phi}\cos{\phi},$$

(29)

with Γ an integral of the motion and Θ̅ a constant. At this point, we are essentially finished. The three conserved quantities Γ, Θ₁, and r₂ make the original four degree of freedom problem a one degree of freedom problem². The total angular momentum is a function of these conserved quantities. However, in order to simplify the analysis, it is worthwhile to do a bit more algebra in order to reduce the Hamiltonian to a one parameter function.

The conserved quantity s₂ is not in involution with r₂ as {s₂,r₂}=1, so there are only three independent conserved integrals of the motion.

To that end, we rescale the variables. Let actions be divided by the unitless parameter Q; angles are unchanged. Let primes mark the rescaled variables. We also rescale both the Hamiltonian K and the time t as H′=K/a and t′=t̂/a. After this step, the Hamiltonian is

$$H^{\prime} =K/a\qquad {}=-\frac{1}{2}|\beta|Q^{2}\frac{1}{a}(\Phi^{\prime}-\Gamma^{\prime})^{2}-\epsilon_{1}\bar{\delta}\sqrt{Q}\frac{1}{a}\sqrt{2\Phi^{\prime}}\cos{\phi}.$$

(30)

If we choose a=Q²|β|, and Q=(ε₁δ̅/|β|)²⁄³, we are only left with one free parameter (Γ′). It can be shown that

$$Q =\epsilon_{p}^{2/3}\zeta\bigg( \frac{f_{m+1,31}^{2}}{9(m+1)(1+\zeta)^{2}}\bigg) ^{1/3}\frac{\alpha_{res}^{5/6}}{m}\bigg( \frac{R^{2}+\zeta\sqrt{\alpha_{res}}}{(\alpha_{res}+\zeta)^{5}}\bigg) ^{1/3},$$

(31)

where we have replaced ε₁=ζεₚ/(1+ζ) with the total mass of the planets εₚ=ε₁+ε₂.

This rescaling results in a Hamiltonian function

$$H^{\prime} =-\frac{1}{2}(\Phi^{\prime}-\Gamma^{\prime})^{2}-\sqrt{2\Phi^{\prime}}\cos{\phi}.$$

(32)

We perform a final canonical transformation

$$X =\sqrt{2\Phi^{\prime}}\cos{\phi}\qquad Y =\sqrt{2\Phi^{\prime}}\sin{\phi},$$

(33)

so that - at last!

$$H^{\prime} =-\frac{1}{2}\bigg[ \frac{1}{2}(X^{2}+Y^{2})-\Gamma^{\prime}\bigg] ^{2}-X$$

(34)

This is the one degree of freedom Hamiltonian with a single free parameter that we have been seeking. It is interesting that the conserved quantities s₂ and r₂ do not appear at all, not even as constant parameters. We note that the Hamiltonian (34), which applies for arbitrary mass ratio ζ, has the same functional form as the Hamiltonian for the motion of a test particle near a first order resonance with a planet (ζ=0 or ζ→∞). It must be true that in either limit the Hamiltonian (34) reduces to the restricted case.

However, we have found something stronger to be true: the Hamiltonian (34) is approximately independent of the mass ratio ζ. To show this, we write s=√(α/αᵣₑₛ)≈ 1+Δα/(2αᵣₑₛ), and expand Equation (21) in powers of Δα/αᵣₑₛ,

$$\delta\Theta =-\frac{\alpha_{res}\zeta}{m(\alpha_{res}+\zeta)^{2}}\frac{1}{2}\frac{\Delta\alpha}{\alpha_{res}}+\ldots$$

(35)

so that

$$\Gamma^{\prime} =\frac{\Psi-\delta\Theta}{Q}\qquad {}=\frac{1}{Q}\times\bigg( g_{m}R^{2}\sigma^{2}+\frac{\zeta}{m(\alpha_{res}+\zeta)^{2}}\frac{1}{2}\Delta\alpha\bigg)\qquad {}=\frac{1}{2}\bigg( \frac{9(m+1)}{\epsilon_{p}^{2}f_{m+1,31}^{2}\alpha_{res}^{5/2}}\bigg) ^{1/3}\bigg( \frac{(1+\zeta)^{2}(\alpha_{res}+\zeta)^{2}}{(R^{2}+\zeta\sqrt{\alpha_{res}})}\bigg) ^{1/3}\times\qquad \bigg[ \frac{mR^{2}\sigma^{2}}{R^{2}+\zeta\sqrt{\alpha_{res}}}+\frac{\Delta\alpha}{\alpha_{res}+\zeta}\bigg] .$$

(36)

We are considering close orbits, where αᵣₑₛ ≈ 1,R² ≈ 1, and σe_1^2+e_2^2-2e_1e_2Δ. In this limit

$$\Gamma^{\prime} \sim\frac{1}{2}\bigg( \frac{9(m+1)}{\epsilon_{p}^{2}f_{m+1,31}^{2}}\bigg) ^{1/3}\bigg( \frac{(1+\zeta)^{2}(1+\zeta)^{2}}{(1+\zeta)}\bigg) ^{1/3}\times\qquad \bigg[ \frac{m\sigma^{2}}{1+\zeta}+\Delta\alpha\frac{1}{1+\zeta}\bigg]\qquad {}=\frac{1}{2}\bigg( \frac{9(m+1)}{\epsilon_{p}^{2}f_{m+1,31}^{2}}\bigg) ^{1/3}\times\bigg[ m\sigma^{2}+\Delta\alpha\bigg] .$$

(37)

The mass ratio dependence of the Hamiltonian drops out for close orbits near first order mean motion resonances. Although the function form of the Hamiltonian (34) is identical to that for a test particle and a planet near a first order mean motion resonance, we would not have a priori expected the coefficients to work out in such a way. Note that this does not imply that the actual motion is independent of mass ratio - for example, the relative amplitude of the eccentricity oscillation of the two planets is still proportional to the mass ratio of the planets.

Also, the mass ratio will certainly be a factor for the 2:1 mean motion resonance, where α is further from unity and R=2.78 (and hence the approximations used to get to expression (2) do not hold). The significant deviation of R from 1 is because of the indirect contribution to the disturbing function for this commensurability.

2.1. Fixed Point Analysis and Conditions for Resonance

The analysis of the level curves and fixed points of the Hamiltonian (34) is the same as in the case of the second fundamental model for resonance (also called an Andoyer Hamiltonian with index 1, see for example Henrard & Lamaitre (1983); Ferraz-Mello (2007)). To summarize, the fixed points satisfy the equations:

$$\frac{dX}{d{t}^{\prime}}=0 \longrightarrow Y\bigg( \frac{1}{2}(X^{2}+Y^{2})-\Gamma^{\prime}\bigg) =0\qquad \longrightarrow Y=0\qquad \frac{dY}{d{t}^{\prime}}=0 \longrightarrow X\bigg( \frac{1}{2}(X^{2}+Y^{2})-\Gamma^{\prime}\bigg) +1=0\qquad \longrightarrow X^{3}-2\Gamma^{\prime}X+2=0.$$

(38)

When Y=0 and X>0, the fixed point corresponds to ϕ=0. When Y=0 and X<0, the fixed point corresponds to ϕ=π. The discriminant for the cubic equation (2.1) is Δ=32(Γ^′ 3-27/8), but only the sign of Δ matters for the number of real roots. If Δ>0, there are three real roots; when Δ<0, there is only a single real root. This is the only bifurcation for this system. The transition corresponds to Γ′=3/2. When there are three real roots, the root with the largest value of X is a saddle point, and the other two are centers (elliptic fixed points). Two homoclinic orbits, together comprising the separatrix, are the level curve connecting the unstable fixed point X₃ to itself, as shown in Figure 1 in red. When the unstable fixed point and the stable fixed point at X₂ merge at the value of Γ′=3/2, a single center remains (at X₁), and there is no longer a separatrix. This Hamiltonian is only a function of Y², so there is symmetry in the contours of H(X,Y) about the Y=0 line.

Figure 1.— Contours of the Hamiltonian function (34) for a value of Γ′=1.7. The separatrix is marked in red. The fixed point with the largest value of X, X₃, is the unstable fixed point, whereas the other two are centers. Note that some resonant orbits correspond to circulation of the resonant angle (resonant contours which enclose the origin, which is marked with a small black dot). It is only when Γ′≥1.88988 that the seperatrix encloses only librating orbits. The dashed blue line illustrates the resonance width.

We clarify here the distinction between an orbit with an oscillating resonant angle and an orbit “in resonance”. An orbit in resonance evolves within the resonant region, which is the area located between the two homoclinic orbits, corresponding to oscillations around X₁ (Henrard & Lamaitre, 1983; Ferraz-Mello, 2007). But depending on the value of Γ′, not all of the enclosed resonant orbits correspond to oscillatory motion of ϕ, and there also exists oscillatory motion of ϕ outside of the resonant region (and also when there is no separatrix at all). Whether or not the resonant angle circulates or oscillates depends on the choice of the origin; the intrinsic dynamics do not depend on coordinate choice.

The fixed point X₁ at the center of the resonance region corresponds to ϕ=π, or

$$\lambda_{2}+m\Delta\lambda+\psi =(2k+1)\pi\qquad \Delta\lambda =\frac{(2k+1)\pi-\psi-\lambda_{2}}{m},$$

(39)

where k is an integer. This implies that there should be m values of Δλ corresponding to a single center in the (X,Y). In terms of the original angles, there are m resonant islands for the m:m+1 resonance. This is exhibited Figure 2 for initially zero eccentricity orbits.

Figure 2.— Seperatrices of the first order mean motion resonances at zero eccentricity as a function of period ratio and initial Δλ. All angles but λ₁ are set to zero, so using Equation (2.1), Δλ=(2k+1)π/m corresponds to ϕ=π. There are m separate resonance centers, marked in red, in Δλ at the m:m+1 resonance. Note that the separatrices do not close on one side because one of the boundaries of the separatrix does not extend through all values of ϕ. The masses of the planets are set to ε₁=ε₂=10⁻⁵. At period ratios close to one, the seperatrices associated with first order resonances of different m overlap.

3. RESONANCE OVERLAP

The motion determined by an integrable 2n degree of freedom Hamiltonian reduces to that of n angles ϕᵢ changing at constant rates ωᵢ. When the unperturbed frequencies ωᵢ are commensurate, satisfying Σᵢ₌₀ⁿ⁻¹ aᵢ ωᵢ ≈ 0, where aᵢ is an integer and at least two of the aᵢ are nonzero, the system is near a resonance (with Σᵢ₌₀ⁿ⁻¹ aᵢ=ϕₐ as the resonant angle). When the system is weakly perturbed (in such a way as to be non-integrable), if only one such integer vector a=(a₀,aₙ₋₁) exists, the evolution of orbits near the resonance is approximately integrable. However, if there exists another (independent) integer vector b such that Σᵢ₌₀ⁿ⁻¹ bᵢ ωᵢ ≈ 0, with a corresponding resonant angle ϕ_b, the motion is more complicated. The resonance overlap criterion states that for a weakly perturbed (and non-integrable) Hamiltonian system, if more than one linearly independent combination of the angles, treated individually and without interaction, are simultaneously resonant, the motion will be chaotic. The region in which a particular angle is resonant is within the boundaries of the seperatrix, and so we must first determine the widths of the first order resonances, in terms of their seperatrices, in order to determine if resonance overlap occurs.

The resonance overlap criterion has been applied to the first order mean motion resonances in the case of an eccentric test particle and a circular planet orbiting a much more massive central body (Wisdom, 1980). It was found that the orbit of a test particle with initially zero eccentricity will be chaotic if its semi major axis a satisfies | a-aₚₗ | /aₚₗ ≲ 1.3εₚₗ²⁄⁷ (where the subscript pl stands for planet) unless the test particle is protected by the 1:1 resonance. An alternative analytic derivation confirmed this, though accompanying numerical integrations indicated that the coefficient is slightly larger (∼ 1.5, Duncan et al. (1989)). This result was extended to slightly eccentric particle orbits by Mustill & Wyatt (2012), who found the overlap region satisfied | a-aₚₗ | /aₚₗ ≲ 1.8(εₚₗ e)¹⁄⁵ when e≳ 0.2εₚₗ³⁄⁷. We are interested here in determining how these thresholds translate in the case of two massive planets on eccentric orbits. Our approach closely follows that of Wisdom (1980).

3.1. Determining the Widths of the Resonances

Our goal is to determine the boundary of the resonant region, defined by the separatrix, as a function of orbital elements. The resonant region only appears when Γ′>3/2, so we only consider this case. For a given value of Γ′>3/2, we first determine the position of the unstable fixed point X₃. The value of the Hamiltonian along the separatrix is H′(X₃,0)=Hₛₑₚ. We will quantify the “width” of the resonant region as the maximal distance between the inner and outer curve of the separatrix. This maximal width is the distance in X between the two crossings of the separatrix with the X axis, marked with a blue dotted line in Figure 1. The crossing points are the solutions of the equation H(X,0)-Hₛₑₚ=0 (where one of the roots is X₃).

It can be shown (see for example Ferraz-Mello (2007)) that the crossing points X_ 1 and X_ 2 have values of

$$X_{\star 1} =-X_{3}+2/\sqrt{X_{3}}\qquad X_{\star 2} =-X_{3}-2/\sqrt{X_{3}},$$

(40)

and that

$$\Delta X=X_{\star 1}-X_{\star 2} =4/\sqrt{X_{3}}.$$

(41)

This width in terms of X is related to the width in terms of Ψ′ and δΘ′ as

$$\Delta\Psi^{\prime} =\Psi^{\prime}_{\star 1}-\Psi^{\prime}_{\star 2}\qquad {}=\frac{1}{2}(X_{\star 1}-X_{\star 2})\big| X_{\star 1}+X_{\star 2}\big\rvert =4\sqrt{X_{3}}\qquad \Delta(\delta\Theta^{\prime}) =\Delta\Psi^{\prime},$$

(42)

where we have used the conserved quantity Γ′=Ψ′-δΘ′. The boundaries of the resonance in terms of the scaled quantities (δΘ′,Ψ′) are the same for every resonance (but there is dependence on m contained in the scaling factor Q).

The boundaries in terms of Ψ′ and δΘ′ correspond uniquely to a width in terms of the weighted eccentricity σ and the semimajor axis ratio of the planets, through Equations (21) and (2). The weighted eccentricity is a weak function of the semimajor axis ratio, as R depends on the Laplace coefficients.

The resulting boundaries for a single first order resonance are shown in the period ratio - σ plane in Figure 3. To make this plot, only ε₁, ε₂, and the particular resonance integer m are required. Three instances of the resonance boundaries are plotted corresponding to the cases of ζ=10⁻⁴,ζ=1, and ζ=10⁴, confirming that the mass ratio of the planets ζ is relatively unimportant for determining resonance widths (and therefore that we only require εₚ, not both ε₁ and ε₂). Note that orbits inside the resonance have different libration amplitudes about the fixed point and do not in general explore the full width of the resonance. Since ^′_ 1=X_ 1^2/2 passes through zero, X_ 1 corresponds to the left boundary of the seperatrix (at smaller period ratios) in Figure 3.

Figure 3.— Left: Boundaries of the m=3, 4:3 mean motion resonance as a function of the period ratio of the planets and their weighted eccentricity σ (Equation 25). The two sets of dashed lines correspond to cases of an interior or exterior test particle with a massive planet, while the solid line shows the resonance bounds for the case of equal mass planets. In all cases the total mass of the planets is 2× 10⁻⁵ relative to mass of the central body. Right: The difference in period ratio and σ of the resonance boundaries between the case of ζ=10⁻⁴ and ζ=10⁴. Plotted is log₁₀| P₂/P₁(ζ=10⁴)-P₂/P₁(ζ=10⁻⁴))| /[(m+1)/m], i.e. the difference in period ratio between the two cases, relative to the period ratio at the nominal resonance location, and log₁₀[σ(ζ=10⁴)-σ(ζ=10⁻⁴)] for several of the first order resonances. Although ζ changes by 8 orders of magnitude, the resonance widths only change by a small amount.

3.2. A Resonance Overlap Criterion for Initially Circular Orbits

We now turn to the question of resonance overlap. If the overlap of two resonances is moderate, only the high amplitude libration orbits (oscillating around the resonant fixed point) will be chaotic (since only the high amplitude libration orbits come very near the seperatrix). For a particular value of σ, however, if the overlap is extreme enough that even the resonance centers are in the overlapped region, all resonant motion at higher values of σ will be chaotic. If the resonances overlap so much that even initially zero eccentricity orbits are chaotic, then approximately all orbits with higher eccentricities will be chaotic as well. As a result, the period ratio at which the first order resonances overlap, even for initially zero eccentricity orbits, is a critical period ratio. Essentially all systems of planets with orbits more tightly packed than this critical period ratio will be chaotic, regardless of the eccentricities.

The seperatrix can only pass through zero eccentricity when X_ 1=0. Using Equation (3.1), this implies that X₃=4¹⁄³. Therefore the width of the resonance at zero eccentricity is Δ(δΘ′)=4√(X₃)=5.04. Since we are working in scaled variables, this width is the same for all m. If we substitute this value of X₃ into Equation (2.1), we find that the value of Γ′ corresponding to zero eccentricity on the seperatrix is Γ′=3/4¹⁄³=1.88988. Since Ψ′=0 at this point, δΘ′=-1.88988 at zero eccentricity on the seperatrix. In the original application of the first order resonance overlap criterion to the restricted three body problem, Wisdom defined an effective resonance width as symmetric about the resonance center, i.e. as Δ(δΘ′)=1.88988-(-1.88988)=3.78. This is a very slight underestimate.

A width in terms of δΘ′ is related to a width in terms of α by solving for s in Equation (21),

$$s=\sqrt{\frac{\alpha}{\alpha_{res}}} =\frac{1-\delta\Theta^{\prime}Qm/(\alpha_{res}/\zeta+1)}{1+\delta\Theta^{\prime}Qm/(1+\zeta/\alpha_{res})}.$$

(43)

We Taylor expand the expression for s in powers of εₚ. The result is

$$s^{2} \sim 1-\delta\Theta^{\prime}c_{0}\epsilon_{p}^{2/3}+O(\epsilon_{p}^{4/3}))\qquad c_{0} =2\bigg[ \frac{f_{31}^{2}(R^{2}+\sqrt{\alpha_{res}}\zeta)(\alpha_{res}+\zeta)}{9\sqrt{\alpha_{res}}(m+1)(1+\zeta)^{2}}\bigg] ^{1/3}.$$

(44)

Therefore the width of the resonance δαₘ/αᵣₑₛ is

$$\frac{\delta\alpha_{m}}{\alpha_{res}} =s^{2}(\delta\Theta^{\prime}_{1})-s^{2}(\delta\Theta^{\prime}_{2})\qquad {}=\Delta(\delta\Theta^{\prime})\epsilon_{p}^{2/3}c_{0}+O(\epsilon_{p}^{4/3}).$$

(45)

Note that if the resonance bounds were symmetric there is no contribution at O(εₚ⁴⁄³), since that term depends on (δΘ^′2₂-δΘ^′2₁).

When overlap occurs, α is near unity, so we will replace R² with 1 and αᵣₑₛ with 1. Moreover, we’ll replace fₘ₊₁,₃₁ with 0.802m+0.87∼ 0.802m for large m. This approximation assumes 1/m is small. We will also neglect the εₚ⁴⁄³ terms in this approximation, and will justify it afterwards. Then the width of the resonance for close orbits (high m) is

$$\frac{\delta\alpha}{\alpha_{res}} =\frac{2}{9^{1/3}}\Delta(\delta\Theta^{\prime})\epsilon_{p}^{2/3}\bigg[ \frac{(0.802m)^{2}(1+\zeta)(1+\zeta)}{(m+1)(1+\zeta)^{2}}\bigg] ^{1/3}\qquad {} \approx 2\bigg( \frac{0.802}{3}\bigg) ^{2/3}\Delta(\delta\Theta^{\prime})\epsilon_{p}^{2/3}\bigg[ \frac{m^{2}}{(m+1)}\bigg] ^{1/3}\qquad {} \approx 4.18\epsilon_{p}^{2/3}m^{1/3}.$$

(46)

where in the last line we have used Δ(δΘ′)=5.04. The mass ratio ζ can only appear in the higher order terms, confirming the weak dependence of the widths on the mass ratio.

The distance between two neighboring resonances is

$$\Delta\alpha =\frac{d}{dm}\alpha_{res}(m)\Delta m\qquad {}=\frac{d}{dm}\bigg( \frac{m}{m+1}\bigg) ^{2/3}(1)\qquad {}=\alpha_{res}(m)\frac{2}{3m^{2}}\frac{m}{m+1}\qquad {} \sim\alpha_{res}(m)\frac{2}{3m^{2}}$$ $$\frac{\Delta\alpha}{\alpha_{res}} \sim\frac{2}{3m^{2}}.$$

(47)

When the distance between two neighboring resonances is equal to the sum of their half widths, the resonance overlap criterion is satisfied. This occurs when

$$\frac{\Delta\alpha}{\alpha_{res,m}} \approx 0.5\times\bigg( \frac{\delta\alpha}{\alpha_{res,m}}+\frac{\delta\alpha}{\alpha_{res,m+1}}\bigg)\qquad \frac{2}{3m^{2}} \sim 4.18\epsilon_{p}^{2/3}m^{1/3}\qquad m \sim 0.455\epsilon_{p}^{-2/7}.$$

(48)

Hence all orbits should be chaotic if their averaged period ratio satisfies

$$\frac{P_{2}}{P_{1}} \lesssim 1+\frac{1}{m}\qquad {} \lesssim 1+2.2\epsilon_{p}^{2/7},$$

(49)

or equivalently

$$\frac{a_{2}-a_{1}}{a_{1}}\lesssim 1.46\epsilon_{p}^{2/7},$$

(50)

and the planets are not protected by the 1:1 resonance. We note that this criterion is for averaged coordinates, but in most cases the transformation between the averaged and true Hamiltonian is negligible. This also must hold for the restricted case; indeed, this result has the same function form as that derived by Wisdom (1980), albeit with a ∼ 10% larger coefficient. The numerical coefficient of 1.46 is dependent on our specific definition of the width of the resonance. If we had carried through the calculation using Δ(δΘ′)=3.78, as Wisdom did, in Equation (3.2), we would have found m∼ 0.51εₚ⁻²⁄⁷, or (a₂-a₁)/a₁ ≲ 1.3εₚ²⁄⁷ in agreement with the Wisdom (1980) result.

This criterion should be interpreted as a minimum criterion for widespread chaos. We have neglected, for example, the effect of higher order resonances. The chaotic separatrices of these resonances serve to link two neighboring first order resonances before they would overlap if there were not higher order resonances present. We have also ignored the finite extent of the chaotic zone around the seperatrix, though this effect is probably less important than the first.

In the above calculation, we kept terms of order 1/m² and εₚ²⁄³ m¹⁄³, and we neglected terms of order εₚ²⁄³/m²⁄³,1/m³, and εₚ⁴⁄³ m¹⁄³. For Jupiter mass planets, εₚ ∼ 10⁻³, and the criterion predicts overlap at m∼ 3, and the neglected terms are of order 5× 10⁻³,3× 10⁻², and 10⁻⁴, respectively. For smaller mass planets, the error incurred by neglecting the εₚ⁴⁄³ terms will be smaller than the that incurred by neglecting εₚ²⁄³/m²⁄³ and 1/m³ terms. This justifies ignoring the εₚ⁴⁄³ terms in the above calculation. However, it is clear that the approximation of αᵣₑₛ=1, R²=1, and 1/m≪1 become less accurate for high mass planets in the regime where complete resonance overlap occurs. We will study these effects numerically.

3.3. A Resonance Overlap Criterion for Initially Eccentric Orbits

We now briefly address developing a resonance overlap criterion as a function of eccentricity (while still assuming eccentricities are small, so that our approximation of the Hamiltonian holds). The resonant widths grow with eccentricity, and so we expect resonance overlap to occur at a period farther from unity when eccentricities are nonzero. Figure 4 shows the separatrices of several first order mean motion resonances in the case of εₚ=2× 10⁻⁵; a resonance overlap criterion for initially eccentric orbits would pass through the crossing points (shown in red) of neighboring resonances. We remind the reader that the the theory does not apply for real systems with large eccentricities; however any analytic overlap criterion derived from the theory should be able to predict where overlap occurs even in a regime where the theory does not apply, and so we consider large eccentricities here as an exercise.

Figure 4.— The separatrices of the first order mean motion resonances for arbitrary eccentricities (in black) when εₚ=2× 10⁻⁵. Note that although the seperatrices are plotted for all σ<1, they do not apply in reality to high eccentricity systems! However, they still can be used to test an overlap criterion as a function of σ. The positions where neighboring first order resonances first overlap are marked by red dots. In both plots, we show the resonance overlap criterion for initially circular orbits (in purple), the estimate derived here for arbitrary eccentricites (in blue), and the Mustill & Wyatt (2012) criterion (in green). The curves are an approximation to the red dots.

An approximate resonance overlap criterion as a function of eccentricity has been obtained for the restricted three body problem (Mustill & Wyatt, 2012), and so we expect to be able to recover the results obtained in that case. To reproduce their result, we look at the large Γ′ limit (which we will show amounts to a large eccentricity limit).

As Γ′ increases, the value of the unstable fixed point X₃ grows. The equation for the fixed points (2.1) in this case reduces to X³-2Γ′ X+2≈ X³-2Γ′ X=0. This admits X=0 and X=±√(2Γ′) as solutions. The resonance center occurs at X₁ ≈-√(2Γ′), the other center at X₂ ≈ 0, and the unstable fixed point at X₃ ≈√(2Γ′). In this limit, then, the width of the resonance in terms of δΘ′ is Δ(δΘ′)=4√(X₃)≈ 4(2Γ′)¹⁄⁴.

Since the location of the center of resonance X₁ grows in magnitude as Γ′ increases, the forced resonant eccentricity σᵣₑₛ ∝√(Ψ′₁)∝| X₁ | ∼√(2Γ′) grows as well. Hence a large Γ′ limit amounts to a high σ limit. How large the eccentricities must be to use this approximation we will quantify momentarily.

Recall that the width of the resonance is determined by the two crossings of the seperatrix with the X axis (see Figure 1). These two crossings in X, denoted X_ 1 and $X_{\star 2}$, correspond to a minimum and maximum value of Ψ′, respectively, for an orbit in resonance for a particular value of Γ′. These equivalently denote a minimum and maximum value of σ for the orbit. We will define the width of the resonance at a initial value of σ to be the width of the separatrix assuming $\sigma=\sigma_{\star 1}$ (the minimum value of σ for the resonant orbit). This is the same definition of width as we used above for initially circular orbits.

Now, we need to express the width Δ(δΘ′), which is a function of Γ′, in terms of $\Psi^{\prime}_{\star 1}$. This involves (approximately) inverting the function

$$\Psi^{\prime}_{\star 1} =\frac{X_{\star 1}^{2}}{2}=\frac{1}{2}(X_{3}-2/\sqrt{X_{3}})^{2}\qquad {} \approx\Gamma^{\prime}-2(2\Gamma^{\prime})^{1/4}+\frac{\sqrt{2}}{\sqrt{\Gamma^{\prime}}}$$

(51)

for the relationship $\Gamma^{\prime}=\Gamma^{\prime}(\Psi^{\prime}_{\star 1})$. To obtain the Mustill & Wyatt (2012) result, we take the limit of large Γ′ and use $\Gamma^{\prime}=\Psi^{\prime}_{\star 1}$. Therefore, the width of the resonance is

$$\Delta(\delta\Theta^{\prime}) \approx 4(2g_{m}R^{2}\sigma^{2}/Q)^{1/4}$$

(52)

Following the same analysis as above, we find the critical m for overlap to be

m ≈ 0.482(εₚ σ)⁻¹⁄⁵

(53)

leading to the critical period ratio and separation where the first order resonances overlap

$$\frac{P_{2}}{P_{1}} \lesssim 1+2.08(\epsilon_{p}\sigma)^{1/5}$$

(54)

and

$$\frac{a_{2}-a_{1}}{a_{1}} \lesssim 1.38(\epsilon_{p}\sigma)^{1/5}$$

(55)

Note that these results do not apply in the case of zero σ because of the approximations made (if so, they would incorrectly imply that the widths of the resonances go to zero at zero σ). When

$$m\sim 0.482(\epsilon_{p}\sigma)^{-1/5}\gtrsim 0.455\epsilon_{p}^{-2/7}$$

(56)

or

σ≳ 1.33εₚ³⁄⁷

(57)

we expect that the overlap criterion in Equation (54) should approximately hold. This suggests that systems with σ≲(0.013,0.035,0.093) and total planetary mass relative to the host star of εₚ=2×(10⁻⁵,10⁻⁴,10⁻³), respectively, are adequately described by the resonance overlap criterion for initially circular orbits.

Equations (54), (55), and (57) have the same functional form as that derived for the restricted three body problem by Mustill & Wyatt (2012), as expected, but with different coefficients. Both the formula derived here and that of Mustill & Wyatt (2012) are correct to within a factor of a few, as demonstrated by Figure 4. The curves plotted come from Equation (54), and take the form log₁₀{σ}=5log₁₀{P₂/P₁-1} +log₁₀{C⁻⁵ εₚ⁻¹}, where C=2.08 (or 2.7, using Mustill & Wyatt (2012)). However, the actual slope is closer to 4.6 in the εₚ=2× 10⁻⁵ case, and we found it to be closer to 4.2 in the εₚ=2× 10⁻⁴ case. Therefore, the scaling of (P₂)/(P₁)-1∝(εₚ σ)¹⁄⁵ is only approximate as well.

There may be a tractable way of obtaining a more accurate first order resonance overlap criterion in the case of nonzero eccentricities. However, such a formula may not be very useful in practice. As we will show in Section 4, the effects of higher order resonances become important even for small eccentricities, and hence the approximation of only considering the first order resonances in deriving an overlap criterion is no longer as valid (though it remains valid at ∼ zero eccentricity, where the widths of the higher order resonances are zero).

4. COMPARISON TO NUMERICAL INTEGRATION

Our model of first order mean motion resonances makes clear predictions for the locations of first order mean motion seperatrices and for the location of the chaos resulting from first order resonance overlap. We numerically evolved suites of initial conditions with a range of period ratios, eccentricities, mass ratios, and total mass of the planets, relative to the star using a Wisdom-Holman symplectic integrator (Wisdom & Holman, 1991). We integrate the tangent equations (Lichtenberg & Lieberman, 1992) concurrently with the equations of motion to determine whether or not the resulting orbits were chaotic. The initial tangent vector used was randomized and normalized to unity. At the end of the integration, the final length of the tangent vector, d, is reported. If an orbit is chaotic, log{d} grows exponentially in time, with a characteristic e-folding time equal to the Lyapunov time. If the orbit is regular, we expect polynomial growth in time: log{d}=plog{t_final}, where p is the polynomial exponent. The estimated Lyapunov time T_Ly, then, is T_Ly=t_final/log{d}, where t_final is the length of the integration.

After an ensemble of integrations, we expect an approximately bimodal distribution of Lyapunov times. The estimated Lyapunov times of regular orbits should be sharply peaked at approximately 10-100 times shorter than the integration time (assuming p≲ 5). The estimated Lyapunov time of chaotic orbits will fall to shorter times than this peak. Long enough integrations clearly differentiate between the two, though we cannot prove that an orbit is regular through numerical integration alone. Our integrations were 10⁶ days long unless noted otherwise. This is ∼ 3,000 orbits of the outer planet, the orbital period of which was fixed at one year. We tested how the integration time could affect our results by comparing the output from integrations of 10⁶, 10⁷ and 10⁸ days and found that 10⁶ days provided reliable results for the orbits studied; see the Appendix for more details.

A symplectic corrector (Wisdom et al., 1996) identical to the one used in Mercury (Chambers & Migliorini, 1997) was implemented to improve the accuracy of the integrations. With εₚ ∼ 10⁻⁵, the typical fractional energy error was smaller than 10⁻⁸ for orbits that do not suffer close encounters. We may not be resolving close encounters with high fidelity, but, orbits that experience close encounters and strong gravitational scatterings are certainly chaotic, and hence for our purposes it is less important to resolve the close encounters accurately than to identify that the orbits are chaotic (by measuring a short Lyapunov time).

Symplectic mappings are known to introduce chaotic behavior when commensurabilities occur between the timestep frequency 2π/Δ t and orbital frequencies of the physical problem (Wisdom & Holman, 1992). In practice, this chaotic behavior is confined to an exponentially small region of phase space as long as the shortest timescale is resolved by a factor of ∼ 20 (Rauch & Holman, 1999). Our time step of ∼ 1 day, for orbits with a minimum period of half of a year, ensures that we do not need to worry about integrator chaos.

4.1. Resonance Widths as a Function of Orbital Elements

The analytic model for first order resonances can be used to predict widths as a function of σ and period ratio as in Figure 3, or, for a fixed σ, the widths as a function of period ratio and initial orbital angles as in Figure 2. In Figure 5 we show the results of our numerical integrations as a function of σ and period ratio for three different values of the initial mean anomaly of the inner planet, M₁(0). In particular, the masses of the planets are held fixed at ε₁=ε₂=10⁻⁵, e₂(0)=0, and all other orbital angles are fixed at zero. The coloring indicates the measured Lyapunov time of the orbits, in units of the initial inner planet’s period; darker colors reflect shorter times. In each case, we have overplotted the predicted resonance widths based on the analytic model.

Figure 5.— Chaotic structure of phase space as a function of period ratio of the planets and the weighted eccentricity σ. In these plots, e₂(0)=0, so σ=e₁. Darker colors indicate shorter Lyapunov times. The three panels show the results for different initial values of the mean anomaly of the inner planet, M₁. The color scale is cut off at a value of T_Ly=3P₁, and orbits with shorter Lyapunov times than this are in black. Any orbits with Lyapunov times longer than the upper value of the color scale are in yellow. The seperatrices for each resonance as predicted by the analytic theory are shown in grey.

The upper panel shows the case of M₁(0)=0. A regular (yellow) region appears for every first order resonance until the first order resonances are completely overlapped, at period ratios near unity, where a large chaotic zone appears that extends to zero eccentricity. However, it is clear that higher order resonances are important in explaining much of the chaotic behavior we see at higher eccentricities. For example, even though the 5:4 at P₂/P₁=1.25 (m=4) and the 4:3 at P₂/P₁=1.3 (m=3) first order resonances are far from overlapping, one can see both the chaotic separatrices of and regular regions associated with the two third order resonances between them (though the second order resonance between them appears to be chaotic). As the period ratio nears unity, overlap between second order and first order resonances starts to occur, and hence the region around the first order seperatrices becomes chaotic even though the first order resonances themselves do not overlap with each other. It is clear that the chaos is not confined to a region where the first order resonances overlap; this is why we believe that a first order resonance overlap criterion as a function of eccentricity may not be very applicable in practice.

In the lower two panels, the same picture emerges, but not every first order resonance appears as a yellow region - in some cases, the first order resonances are entirely chaotic when they were entirely regular at M₁(0)=0. The reason for this is demonstrated in Figure 6, which shows the chaotic regions of the phase space as a function of period ratio and M₁, with all other initial orbital elements held fixed. The plots in the figure, which correspond to nonzero eccentricities, agree with the predictions of the model as to the number of and location of resonance islands for a particular value of m (and should be compared with Figure 2).

Recall that the resonant overlap criterion only predicts the extent of the chaotic zone where ∼all orbits will be chaotic. The chaotic zone at larger period ratios contains many regular islands; as a result, we refer to it as the chaotic web. These regular islands correspond to first order mean motion resonances that are only weakly overlapped with neighboring first (or second) order resonances. The resonant islands at a particular period ratio are smaller in the lower panel (where e₁ is larger) compared to the upper panel. This is because the resonance widths grow with eccentricity and consequently the overlap becomes more extreme.

Figure 6.— Map of the chaotic structure as a function of M₁ and period ratio, with all else held fixed. Yellow regions embedded in the chaotic web corresponds to resonant islands. The islands in a particular chain appear disjoint because there is an unstable fixed point in between each one. Chaos appears first near the unstable fixed points. Note also that all other orbital angles are initially zero, so the initial value of ϕ(0) depends only on M₁.

It is only in the case of M₁=0 that islands associated with every m get sampled (think of taking a vertical slice in Figure 6 at M₁=0). Other values of M₁ do not pass through as many resonant islands. Using Equation (2.1), with e₂(0)=0 and all angles zero except λ₁ yields the location of the resonant islands in terms of λ₁:

$$\lambda_{1}\approx M_{1} =-\frac{2k\pi}{m},$$

(58)

with k an integer ranging from 0 to m-1. We see that the case M₁=0 corresponds to the center of a resonant island for every m (k=0). The center of the islands for the 6:7 resonance at P₂/P₁ ∼ 1.1667 (m=6), for example, appear at M₁=-(k)/(3)π=(0,π/3,2π/3,π,4π/3,5π/3). A value of M₁=π/2 falls directly in between two of these islands, right at the unstable fixed point. It makes sense, then, that the 6:7 resonance appears entirely chaotic when M₁=π/2, while it appears regular for M₁=0. The dependence of the location of the resonant islands on m and M₁ explains many of the features in Figure 5. It is clear that the case of M₁=0 (and all other angles equal to zero) is a special case; even the small change to M₁(0)=π/10 leads to significant changes in the structure of the chaotic zone.

Another interesting feature of Figure 6 is that the locations of the resonances at M₁(0)=0 are consistently shifted to smaller period ratios compared to neighbors in the same resonant chain. However, the analytic model developed in Section 2 predicts no such behavior. This effect is caused by the transformation between real coordinates (those used in the numerical integration) and averaged variables (those used in the analytic development). The averaging step was necessary to remove the short period terms in the disturbing function. In the Appendix, we derive the canonical transformation between the two sets of variables and show that in the case of small eccentricities, the shift is in the direction observed and it is maximized when M₁=0. In Figure 5 we have corrected the theoretical curves in the case of M₁=0 to account for this transformation (but only using the zero eccentricity terms). There is still a slight deviation between the analytic curves and the numerically determined seperatrices at period ratios near unity and nearly zero eccentricities. This discrepancy could presumably be removed if we carried the transformation to second order in the masses.

In summary, the initial value of the angles do not greatly affect the overall extent of the chaotic web in terms of period ratio (as shown in Figure 6). The locations of the resonant islands in this web do depend on the orbital angles, but in a way that is predicted by the analytic theory and confirmed numerically. Although e₂ and $\Delta\varpi$ were initially zero for all of these integrations, the individual values of e₁, e₂, and $\Delta\varpi$ do not matter. Only the weighted eccentricity σ (and the generalized longitude of pericenter ψ) matter. We confirmed that the predicted boundaries of the resonance in terms of σ matched well with those observed numerically in the case of e₁ variable, e₂=0.02, and $\Delta\varpi=\pi/2$: see the Appendix for details.

4.2. Planetary Mass Ratio Dependence

One of the main results of the analytic work is that the widths of first order mean motion resonances should be independent of the mass ratio of the planets, for a fixed εₚ, especially for resonances with higher values of m. This implies that the width of the chaotic zone due to overlap of first order mean motion resonances is independent of ζ. We confirm numerically that the structure of the chaotic zone does not change significantly over 12 orders of magnitude in ζ; the results are shown in Figure 7 for both ∼zero eccentricity (e₁=10⁻³,e₂=0) and nonzero eccentricity (e₁=0.05,e₂=0). These plots can be thought of as horizontal slices at a particular σ in Figure 5.

At very small eccentricities, the chaotic zone is due almost entirely to first order mean motion resonance overlap, and so it makes sense that we see no significant dependence on ζ. However, it is surprising that the mass ratio is relatively unimportant even at higher eccentricities and at larger period ratios, where the chaos we are observing is due to second and third order mean motion resonance overlap. These results suggest that the widths of higher order resonances, in the case of two massive planets on close orbits, can be approximated analytically using the restricted three body problem. We defer an investigation of this to future work.

The independence of the structure of the chaotic region on the mass ratio is even more striking when we consider a wide range of eccentricities. Figure 8 shows the chaotic regions of phase space for the same initial conditions used to create the bottom panel in Figure 5. The only parameter that has changed is the mass ratio of the planets ζ, which increased from 1 to 10. The total mass of the planets is unchanged. Figure 8 is almost indistinguishable from the bottom panel of Figure 5. Tests with ζ=100 and ζ=0.1 showed similar behavior. It is clear that the mass ratio of the planets matters very little in determining the structure of the chaotic zone in this regime (P₂/P₁ ≲ 1.5-2,σ≲ 0.1).

Figure 7.— Map of the chaotic regions of phase space as a function of period ratio and mass ratio of the planets. The color scale reflects the numerically estimated Lyapunov time, relative to the initial period of the inner planet. The effect of the mass ratio of the planets on the chaotic zones is negligible, even at higher eccentricities. These integrations used initially circular planets with εₚ=10⁻⁵. Only the initial period ratio and mass ratio of the planets is varied. All orbital angles are set to zero (the initial orbital angles chosen will affect the location and widths of these bands of chaos).

Figure 8.— Chaotic structure of phase space as a function of period ratio of the planets and the weighted eccentricity σ. In these plots, e₂(0)=0, so σ=e₁. We set εₚ=2× 10⁻⁵ and ζ=10, with all angles set to zero initially except M₁=π/2. This plot is completely analogous to and should be compared to the bottom panel of Figure 5; the only difference is the value of ζ used. Changing the planetary mass ratio makes no significant changes to the structure of the chaotic region.

4.3. Effect of the Total Mass of the Planets Relative to the Mass of the Host Star

The resonance widths grow with the total mass of the planets, relative to the host star, as εₚ²⁄³, and so we expect that the overall extent of the chaotic web should grow with εₚ as well. Figure 9 shows the chaotic structure of the phase space in the period ratio - M₁ plane for three different values of εₚ (the bottom panel, with εₚ=2× 10⁻⁵) is the same as upper panel of Figure 6). The eccentricity of the inner orbit e₁ was fixed at 0.01, e₂ was fixed at 0. The discussion accompanying Figure 6 explained why resonant islands appear at all, and how their sizes depend on the eccentricities. Here we see how they change as the total mass of the planets grows.

Figure 9.— The effect of the total mass of the planets relative to the mass of the host star, εₚ, on the chaotic structure of phase space for close orbits. Top: εₚ=2× 10⁻³, Middle: εₚ=2× 10⁻⁴, Bottom: εₚ=2× 10⁻⁵. For more massive planets, relative to the star, the first order resonances are more overlapped, and less of the phase space at low eccentricity and period ratios near unity is regular. In each plot, the thick neon green line corresponds approximately to the Hill boundary. For low mass planets there are more resonant islands.

At the highest masses shown (εₚ=2× 10⁻³), there are no resonant islands at all, since the widths of the resonances are large enough to completely overlap. Presumably the abrupt transition at period ratios of ∼ 1.45 to a mostly regular phase space happens because the 3:2 resonance does not overlap much with the second order resonance above it (the 5:3) resonance, while it does with the second order resonance below it (the 7:5) since the 7:5 is closer. The planets would have to be more massive than Jupiter or more eccentric in order to merge the chaotic seperatrices of the 3:2 and the 5:3 resonance.

The effect of the transformation between real and averaged coordinates grows ∝εₚ (see Equation (A-1)). As a result, the shift at M₁=0 is much more extreme in the εₚ=2× 10⁻⁴ case compared to the εₚ=2× 10⁻⁵ case.

Finally we point out that for low mass planets on low eccentricity orbits, the Hill boundary (shown in green) lies within the chaotic web. The implications of this will be discussed in Section 5.

4.4. Resonance Overlap Scaling

How well do the derived resonance overlap criteria given in Equation (50) and Equation (55) predict the size of chaotic zone in practice? In Figure 10 we show, in terms of the initial value of log₁₀{a₂/a₁-1} and log₁₀{εₚ}, the location of initially circular orbits we determined to be chaotic numerically. Only the total mass of the planets relative to the mass of the host star εₚ and the semimajor axis of the inner planet were varied. Recall that the criterion predicts, given εₚ, the extent of the chaotic zone where there are no regular orbits (aside from those protected by the 1:1). The derived formula for the boundary of the chaotic zone closely tracks the numerically determined edge across a wide range of εₚ, regardless of the initial values of the orbital angles. When the mass of the planets is larger, the simple scaling appears to underestimate slightly the extent of first order resonance overlap for initially zero eccentricity orbits. A slight discrepancy was expected here due to higher order terms, as discussed in Section 3. Interestingly, the Hill boundary appears to track the extent of the chaotic zone for more massive planets, at least when M₁(0) is further from zero.

For intermediate mass planets, with initial orbital configurations M₁≠M₂, a regular region appears corresponding to the 1:1 mean motion resonance. The 1:1 region will be explored in future work (Payne et al, in prep).

Figure 10.— Comparison of the size of the chaotic zone for initially circular orbits as predicted by the resonance overlap criterion (Equation (50), green line) to that determined by numerical integrations. The black points represent the location of chaotic orbits in terms of the initial value of a₁ and the total mass of the planets relative to the mass of the host star, εₚ. In red is the Hill boundary (Equation (59)); the resonance overlap criterion is a stricter criterion for low mass planets. Only a₁ and εₚ were varied within each panel, with M₁ varied between panels. All other angles were fixed at zero. We have not applied the correction between numerical coordinates and averaged coordinates in this figure. The correction is small and will mainly affect the results when M₁=0.

Figure 11.— Comparison of the size of the chaotic zone for initially eccentric orbits as predicted by the eccentric resonance overlap criterion (Equation (55), purple line) to that determined by the analysis of Mustill & Wyatt (2012) (the green line) and the results of numerical integrations. The estimate of the size of the chaotic zone predicted by the initially circular resonance overlap criterion (Equation (50), cyan line) is also shown. It is clear that first order resonance overlap alone cannot account for all of the chaos shown, particularily at lower masses. Within each panel, only the initial orbital separation of the planets, e₁, and M₁ vary. M₁ was chosen randomly. All other angles were fixed at zero, and e₂=0. The different panels correspond to different values for εₚ, ζ was fixed at 1.

Figure 11 shows the boundary of the chaotic zone as predicted by the overlap criterion for initially eccentric orbits (Equation (55)) on the log₁₀{σ} -log₁₀{a₂/a₁-1} plane for three different values of εₚ. We also show the estimate of Mustill & Wyatt (2012) (which lies parallel to the estimate derived in this work, both have slopes of 5 but differ by a constant coefficient). These lines apply in the large σ regime, after they cross the overlap boundary predicted for initially circular orbits (the vertical line). The Mustill & Wyatt (2012) coefficient does fit better than the one derived in Section 3.3. However, it is clear that the first order resonance scaling as a function of σ is approximate. It tracks the boundary of the chaotic zone only for a small range of σ. Once σ is large enough, second and third order resonances become important, and the slope of the edge of the chaotic zone on the log₁₀{σ} -log₁₀{a₂/a₁-1} plane changes significantly from the estimated value of 5.

5. DISCUSSION

5.1. Chaotic Dynamics and the Lifetime of Chaotic Orbits

In order to use the resonance overlap criterion as a stability criterion, we need to ensure that the orbits in the chaotic zone are indeed Lagrange unstable. Paardekooper et al. (2013) found that in a region of phase space similar to that inhabited by the Kepler 36 system (period ratio close to 6:7, e₁=0,e₂<0.04) orbits that were chaotic also eventually experienced at least a 10% variation in semimajor axis (compared to the initial value) and hence were Lagrange unstable. A smaller subset of these initial conditions led to ejections of the outer planet. Hence we might expect that orbits in the chaotic zone are indeed Lagrange unstable on short timescales. Gladman also found that the majority of the chaotic orbits were Lagrange unstable, but that near the edge of the chaotic web orbits appeared to be relatively quiescent for many thousands of Lyapunov times, despite having very short Lyapunov times (on the order of the synodic period of the planets) (Gladman, 1993).

Figure 12.— The location of chaotic orbits (in black triangles, as determined by integrations of 10⁶ days) compared to those that show Lagrange instabilities within 10⁹ days (in colored circles). The initial conditions used here are the same as those used to create the lower panel of Figure 6. The location of unstable orbits closely track the chaotic region of phase space - underneath nearly every colored point is a black triangle. Note that regions with no points plotted were integrated and found to be both regular and long-lived. For reference, period ratios larger than ∼ 1.148 are Hill stable, so the chaotic web extends past the Hill boundary. The large black dot corresponds to the initial condition shown in Figure 14 (M₁(0)∼ 0,P₂/P₁ ∼ 1.17).

To determine the stability of the orbits in the chaotic zone, we studied the evolution of a subset of our initial conditions. The majority of the initial conditions were integrated for 10⁸ days. Chaotic orbits which did not show instability on this timescale were integrated for 10⁹ days. Energy was conserved in these integrations to within a part in 10¹⁰ for orbits not resulting in strong scattering events. Orbits where the semimajor axes changed by more than 5% from their initial values were classified as Lagrange unstable; the first instance that | a(t)-a(0)| /a(0)>0.05 is satisfied is the instability time. We find that many of the chaotic orbits are indeed Lagrange unstable on short timescales, regardless of whether they satisfy the Hill criterion. Figure 12 shows the results of a test using the same initial conditions as the bottom panel of Figure 6. The location of the unstable orbits (in colored points) closely track the structure of the chaotic web (in black triangles). A comparison of Figure 12 and Figure 6 illustrates how similar the two are; underneath almost every colored circle is a black triangle. Note that the resonance overlap criterion does not immediately apply to this case, as the eccentricities are nonzero.

All unstable orbits are chaotic, but it is also clear that a small fraction of the chaotic orbits do not show instabilities during the 10⁹ day integrations. As expected based on Gladman’s work, these orbits are grouped near the edge of the chaotic web; all satisfy the Hill criterion. Since the timescale to exhibit instability is consistently longer nearer the edge of the chaotic web, we predict that chaotic orbits that have not shown instabilities in 10⁹ days will reveal erratic motion in longer integrations (though if the instability timescale is longer than ∼ 10¹² days the system is effectively long-lived). The dependence of an instability timescale on the initial location of an orbit in the chaotic zone has been observed (Murray & Holman, 1997; Deck et al., 2012) and interpreted as a dependence of the chaotic diffusion coefficient on the orbital elements. Note that chaotic orbits near the edge of the chaotic web (in period ratio) are those that satisfy the Hill criterion by the largest amount, a quality that has been found to correlate with the timescale to exhibit Lagrange instabilities (Deck et al., 2012). All of the initial conditions leading to ejections or extremely wide two planet systems within 10⁹ days fail the Hill criterion.

Figure 13.— The distribution of the ratio of the Lagrange instability time to the Lyapunov time of the unstable orbits. Orbits which fail the Hill criterion exhibit instabilities in significantly fewer Lyapunov times compared to those which satisfy the criterion. There is not a clear relationship between Lyapunov time and instability time.

We find that the estimated Lyapunov time of these orbits are 0.3-30 times the initial synodic period of the pair, similar to Gladman’s result. Note that both the estimated Lyapunov time and the synodic period of Lagrange unstable orbits change can change as orbital instabilities set in. The distribution of the ratio of the instability time to the Lyapunov time is shown in Figure 13 for the unstable initial conditions, showing that indeed the Lyapunov time is orders of magnitude shorter than the instability time.

How do trajectories evolve within this chaotic web once instabilities set in? Although the chaotic zones at different period ratios look disjoint in the plots of period ratio vs. σ (for example, Figure 5), other slices of the phase space (particularly in period ratio vs M₁, e.g. Figure 6) suggest that these zones are part of a larger connected chaotic region. The similar Lyapunov times across the chaotic regions also hints that the same mechanism is responsible for the chaos. In either case, however, we are looking at a projection of the chaotic zone onto a two dimensional subspace and so we cannot claim based on these figures that the chaotic web is indeed connected, making up a single chaotic zone.

Figure 14 shows the time evolution of the period ratio of a single Hill stable initial condition. The trajectory indeed wanders throughout the chaotic web, which extends over a large range of period ratios. The distribution of the period ratio is what one would expect based on resonance overlap - the planets spend very little time near the first order mean motion resonances with low m (P₂/P₁=1.5,1.33,1.25,1.2) nor near the second order resonances with low m (P₂/P₁=1.4,1.28). This is because there are large regular islands at these period ratios. The chaotic web has fractionally less volume of phase space near these first and second order resonances, and the orbit spends proportionally less time here. Taken together, these observations strongly suggest that the chaotic web associated with resonance overlap is all connected, as it appears to be in Figure 6. Note that as the period ratio grows, the eccentricities of the planets increase as well. At higher eccentricities, though, the chaotic web extends to higher values of the period ratio. This is how an orbit beginning at P₂/P₁ ∼ 1.17 and low eccentricity can reach period ratios as large as 1.7.

Even if a pair of planets will always remain bound to their host star while the star is on the main sequence, we expect the multi-planet systems we observe to be essentially stationary, in that their period ratios should not be varying erratically on short timescales. As a result, we predict that we should not observe systems of two planets with period ratios smaller than the critical period ratio determined by the overlap criterion, even though these systems may be Hill stable.

Our work suggests that the Lagrange instability time of Hill stable close orbits is much shorter than the timescales for ejections. This is in agreement with recent numerical studies of planetary ejection in Hill stable systems (Veras & Mustill, 2013). In the region of parameter space studied here, the reason for this may be that the chaotic zone associated with first order mean motion resonances can only extend to P₂/P₁ ∼ 2, and diffusion to period ratios larger than this in a sparse chaotic web will be very slow. It is important to extend the definition of instability to include erratic variations in the semimajor axes especially if using long-term stability as a check for whether or not a given planetary system is stable enough that we would be likely to observe it in its current configuration.

It is interesting that the chaotic zone caused by first order resonance overlap is in fact divided by the Hill boundary for low mass planets (as explained in the following section, and demonstrated in the lower panels of Figure 9). The Hill boundary corresponds to an invariant surface in the phase space. Orbits on one side cannot cross to the other. The chaotic web is a single chaotic zone in the sense that the mechanism for chaos is the same across it, but it is probably more correct to think of it as two separate chaotic zones, since chaotic orbits cannot explore both sides of the Hill boundary.

Figure 14.— Evolution of the period ratio from a single integration of two planets in a Hill stable, but Lagrange unstable, orbital configuration. The distribution of the period ratio is what one expects from resonance overlap of the first and second order mean motion resonances. The location of this initial condition in the chaotic web is shown in Figure 12 and Figure 15. The maximum fractional deviation in the energy for this integration is max| (E(t)-E(0))/E(0)| ∼ 3× 10⁻¹⁰, and the estimated Lyapunov time is 21,144 days, or ∼ 67P₁(0).

5.2. A Relationship Between the Hill Criterion, Lagrange Instabilities, and Resonance Overlap

In the previous section, we demonstrated that orbits in the chaotic web caused by first order resonance overlap are unstable, and therefore we can use resonance overlap criteria to determine long-term stability. We now compare the resonance overlap criterion to the Hill criterion and more generally consider the effectiveness of the Hill criterion in the context of the resonance overlap.

In the case of initially circular orbits, two planets are Hill stable if their initial semimajor axes satisfy

$$\frac{a_{2}-a_{1}}{a_{1}}\gtrsim 2.4\epsilon_{p}^{1/3}.$$

(59)

We have found that for initially circular orbits, first order resonances will overlap and result in a chaotic phase space if (a₂-a₁)/a₁ ≲ 1.46εₚ²⁄⁷. This resonance overlap criterion is absolute in that orbits which fail it, regardless of the eccentricities, are effectively guaranteed to be chaotic (and Lagrange unstable) if they are not protected by the 1:1 resonance.

As Gladman pointed out, these two boundaries cross, and our work predicts the crossing to occur at a value of εₚ=(1.46/2.4)²¹. Both boundaries are plotted in Figure 10. For planets with εₚ greater than 3× 10⁻⁵, corresponding to about 10 Earth masses around a solar mass star, the Hill criterion is stricter. However, for smaller planets, the region of complete resonance overlap (where all orbits are expected to be chaotic) extends to larger period ratios P₂/P₁ than the Hill boundary.

Recall, however, that the resonance overlap criterion as calculated here is a minimum criterion for chaos as we have not included the effects of higher order resonances in between the first order resonances. Since the coefficient is taken to the 21st power, even a small increase in the coefficient significantly increases the value of εₚ where the Hill boundary and the overlap boundary are equal.

Gladman found that for initially circular orbits the Hill criterion was effectively a necessary condition for collisional stability of low mass (εₚ ≲ 10⁻⁵, arbitrary ζ>1) planets - though it is only formally sufficient - but that at higher eccentricities many orbits which failed the criterion were seemingly long-lived. We suggest that the Hill criterion is an effectively necessary condition for stability of low mass planets on initially circular orbits because the region of complete resonance overlap (that which is predicted by the criterion (a₂-a₁)/a₁ ≲ 1.46εₚ²⁄⁷) extends past the Hill boundary and there are no regular regions within it (ignoring the 1:1). Hence, initially circular orbits which fail the Hill criterion are almost guaranteed to be chaotic, and, as we have seen in the previous section, unstable. For initially circular orbits, the Hill criterion do not depend on ζ and so this result should be independent of the planetary mass ratio.

For orbits with arbitrary eccentricities, the criterion for Hill stability can be written as

$$\frac{p}{a}=\frac{-2(m_{\star}+m_{1}+m_{2})}{G^{2}(m_{\star}m_{1}+m_{2}m_{1}+m_{\star}m_{2})^{3}}L^{2}H \gt \frac{p}{a}\lvert_{crit},$$

(60)

where L is the total angular momentum, H is the total orbital energy, (p/a)|_crit is a function only of the masses of the bodies, and p and a are the generalized semi-latus rectum and semimajor axis of the problem (Marchal & Bozis, 1982; Gladman, 1993). Figure 15 shows contours of (p/a)/(p/a)|_crit and the underlying chaotic zone as a function of period ratio and σ for two equal mass planets with εₚ=2× 10⁻⁵. The dotted contour is the Hill boundary, above this curve all orbits fail the Hill criterion, below it all orbits are Hill stable.

At larger eccentricities there are still regular regions of phase space which do not satisfy the Hill criterion but will be long-lived. We posit that the Hill criterion is not effectively necessary for stability of (i) moderately eccentric orbits of planets of any εₚ and (ii) higher mass systems ($\epsilon_{p}\gtrsimfew\times 10^{-5}$) on initially circular orbits since in these cases the region of complete resonance overlap does not encompass the Hill boundary.

As we have seen in Section 4.2, the mass ratio of the planets does not greatly affect the structure of the chaotic zone and so we show the Hill boundary ((p/a)/(p/a)|_crit=1) in the cases of ζ=0.2 and ζ=5 in Figure 15 as well. The Hill criterion depends strongly on the planetary mass ratio ζ.

Figure 15.— Various contours of (p/a)/(p/a)|_crit, with ζ=1, overplotted in grey onto the bottom panel of Figure 5. The dotted line corresponds to the Hill boundary, where (p/a)/(p/a)|_crit=1. Initial conditions above this line fail the Hill criterion, those below it satisfy the criterion. The contours correspond to values of (p/a)/(p/a)|_crit=0.9992,0.9995,0.9998,1.0,1.0002,1.0005,1.0008 from top to bottom. The large black dot corresponds to the initial condition shown in Figure 14. Though the Kepler 36 system does not fit directly onto this plot, that system has a similar total planetary mass relative to the mass of the host star εₚ, with σ≲ 0.04 and P₂/P₁ ∼ 1.17. We also show the critical Hill contour in the case of εₚ=2× 10⁻⁵ and ζ=0.2 and ζ=5, to show how the Hill contours depend on the mass ratio ζ.

We now turn to one of the motivating questions for this work: why is the proximity of an orbit to the Hill boundary (how close (p/a)/(p/a)|_crit is to 1) seemingly so important for determining whether or not a given planetary system is Lagrange long-lived? And why is the transition between Lagrange unstable and Lagrange long-lived orbits such a sharp function of this distance from the Hill boundary?

Figure 16.— A study of how the chaotic fraction of the phase space changes as a function of (p/a)/(p/a)|_crit for equal mass planets. Over 400,000 initial conditions were used in making this plot, taken from the initial conditions used in the three panels of Figure 5 (εₚ=2× 10⁻⁵,ζ=1). A wide range of period ratios and eccentricities were used, as well as different initial orbital angles. Despite this, we see a consistent trend: a sharp transition from an almost entirely chaotic phase space to an entirely regular phase space. The exception is at the highest eccentricities, where regular regions associated with resonances are more important.

In the case of comparable mass planets, the contours of (p/a)/(p/a)|_crit lie approximately parallel to the edge of the chaotic zone, and a small change in the parameter (p/a)/(p/a)|_crit moves orbits entirely out of the chaotic region. To explore this further, we determined the fraction of the orbits in each bin of (p/a)/(p/a)|_crit which were chaotic using all of the initial conditions in each panel of Figure 5. These orbits had εₚ=2× 10⁻⁵ and ζ=1. The results are shown in Figure 16 for three different eccentricity groups. It is clear that at higher eccentricities the fraction of the phase space which is chaotic decreases much less steeply as a function of (p/a)/(p/a)|_crit than it does for lower eccentricities.

We have shown that chaotic orbits are typically Lagrange unstable, and therefore the sharp transition between an almost entirely chaotic phase space and an almost entirely regular one can naturally explain the observed sharp transition between Lagrange unstable and Lagrange long-lived orbits with (p/a)/(p/a)|_crit - in the case of low eccentricities (σ≲ 0.05) and comparable mass planets. For reference, a stability analysis of the Kepler 36 system (εₚ ∼ 4× 10⁻⁵,e₁ ∼ 0.02,e₂ ∼ 0.01,ζ=1.8) has shown that the orbits must satisfy (p/a)/(p/a)|_crit ≳ 1.0007 in order to be long-lived.

The curves shown in Figure 16 will depend on the planetary mass ratio ζ (because the contours of (p/a)/(p/a)|_crit depend strongly on the planetary mass ratio, see Figure 15). However, since in the limit of zero eccentricity the Hill boundary is independent of ζ (see Equation (59)), we would predict a similar transition for very low eccentricity orbits regardless of ζ. However, we would not expect a sharp transition to occur in close two planet systems with eccentricities 0.1≳σ≳ 0.05, especially for planetary mass ratios far from unity. This is because at larger eccentricities there are regions which fail the Hill criterion yet remain long-lived, as demonstrated by Figure 16, and because contours of (p/a)/(p/a)|_crit do not necessarily track the edge of the chaotic zone when $\zeta\not\approx 1$.

We note that in the extreme limit of the elliptic restricted three body problem, (p/a)|_crit → 1 and (p/a)→ 1-e², where e is the eccentricity of the massive planet. All orbits fail the Hill criterion (see eqn. (50), Marchal & Bozis (1982), or Milani & Nobili (1983)), though this does not imply that all orbits are unstable. In this case, the Hill contours begin to look horizontal on the (P₂/P₁,σ) plane as there is no P₂/P₁ dependence. As a result, there would be no sharp transition in the Lagrange stability of orbits as a function of (p/a)/(p/a)|_crit.

We cannot address the abrupt transition between Lagrange unstable and Lagrange long-lived orbits observed for two planet systems in areas of parameter space different from those studied here, with higher eccentricities (e≳ 0.1) and P₂/P₁ ≳ 2 (Barnes & Greenberg, 2006; Kopparapu & Barnes, 2010; Veras & Mustill, 2013). A more complete understanding of a connection between resonance overlap (and Lagrange instability) and the Hill criterion will require more work.

We hypothesize that instability of close orbits with nearly circular orbits (σ≲ 0.1) is caused by first order resonance overlap, either between first order resonances directly, or indirectly between first and second order resonances. However, the first order resonances are crucial to the overlap because they effectively connect regions of phase space with P₂/P₁ ∼ 1 to regions of phase space with P₂/P₁ ∼ 1.5-2; the first order resonance provide the backbone structure of the chaotic web. An analytic derivation of the boundary of the chaotic zone resulting from first and second order resonance overlap, as a function of eccentricity and period ratio, is beyond the scope of this paper. However, we predict that a derivation for the extent of this chaotic zone would serve as an effective Lagrange criterion for stability of systems of two planets. This is left to future work.

6. CONCLUSION

The Hamiltonian governing the dynamics of two planets near first order mean motion resonance can be reduced to a one degree of freedom system when the orbital eccentricities and inclinations are small (Sessin & Ferraz-Mello, 1984; Wisdom, 1986; Henrard et al., 1986). We have used the Hamiltonian in this reduced form to study the structure of the first order resonances and to derive their widths as a function of orbital elements of the planets. These widths are predicted to be independent of the planetary mass ratio ζ=m₁/m₂. Using these analytically determined widths, we have derived a resonance overlap criterion for two massive planets in the cases of initially circular orbits and initially eccentric orbits. For the former case, we find that the first order mean motion resonances will overlap when (a₂-a₁)/a₁<1.46εₚ²⁄⁷, where εₚ is the total mass of the planets, relative to the mass of the host star. This implies that approximately all orbits failing this criterion (even those that are eccentric) will be chaotic unless they are protected by the 1:1 resonance.

These analytic results were compared extensively to numerical integrations, which confirmed the general predictions of the analytic theory as to the locations and widths of the first order resonances. These integrations demonstrated that the derived overlap criterion for initially circular orbits closely tracks the edge of the numerically determined chaotic zone at zero eccentricity. However, it is clear that chaotic structure of phase space at higher eccentricities and larger period ratios is partially due to higher order mean motion resonances between the first order resonances. As a result, the analytic estimate of the width of the chaotic zone caused by first order resonance overlap as a function of eccentricities may not be very applicable in practice.

Our numerical investigations show that the chaotic structure caused by first (and higher) order resonance overlap forms a “chaotic web” in phase space, with regular regions appearing only deep in low order resonances. We have found that the structure of this chaotic web is indeed independent of the mass ratio of the planets ζ. Though this is what we expected based on the first order resonance theory, we find it surprising as it suggests that the widths of second order resonances are approximately independent of ζ in the regime P₂/P₁ ≲ 2,σ≲ 0.1. Numerical integrations show that (i) chaotic orbits within the chaotic web are generally Lagrange unstable on short timescales, (ii) the timescale to exhibit instability grows the closer the orbit is to the boundary of the chaotic web, and (iii) the timescale to exhibit instability is apparently not related to the Lyapunov time of the orbits. Finally, we have demonstrated that chaotic orbits explore the entire extent of the chaotic web (and therefore the chaotic web is likely a single connected chaotic zone).

These studies demonstrate that resonance overlap criteria can be used as stability criteria. We have shown that for systems with εₚ ≲ 10M_⊕ and initially circular orbits the resonance overlap criterion is a more restrictive stability criterion than the Hill criterion. This led us to an explanation for Gladman’s numerical result that the Hill criterion is effectively a necessary criterion for stability of low mass initially circular two planet systems. More generally, our work suggests that orbital instabilities of close planets are generated by overlap of first (and higher) order mean motion resonances. Therefore we believe that a more extensive resonance overlap criterion which incorporated first and second order resonances could be used as a Lagrange stability criterion for observed pairs of planets (and possibly be a better indicator of stability than the Hill criterion).

Lastly, we have taken some initial steps towards understanding the observed relationship between the Hill criterion and measured Lagrange instability boundaries. In the special case of approximately equal mass planets, contours of the distance from the Hill boundary appear to lie parallel to the edge of the chaotic zone caused by resonance overlap. This naturally produces a sharp transition from an entirely chaotic (and unstable) phase space to an almost entirely long-lived (and regular) phase space for low eccentricity orbits (e≲ 0.05). However, it is clear there is still much to be done to understand the connection between Lagrange instability and the Hill criterion for more eccentric orbits and for systems of two planets with very unequal masses.

K.M.D. acknowledges support from an NSF Graduate Research Fellowship. M.J.P. acknowledges supports from the NASA Kepler Participating Scientists Program and from the NASA Origins of Solar Systems Program. We would also like to thank the anonymous referee for reading the text so closely and for many helpful suggestions.

References

• Barnes & Greenberg (2006)

  Barnes, R., & Greenberg, R. 2006, ApJ, 647, L163

• Barnes & Greenberg (2007)

  —. 2007, ApJ, 665, L67

• Batalha et al. (2013)

  Batalha, N. M., et al. 2013, ApJS, 204, 24

• Batygin & Morbidelli (2013)

  Batygin, K., & Morbidelli, A. 2013, A&A, in press. arXiv:1305.6513

• Chambers & Migliorini (1997)

  Chambers, J. E., & Migliorini, F. 1997, in Bulletin of the American Astronomical Society, Vol. 29, AAS/Division for Planetary Sciences Meeting Abstracts #29, 1024

• Chirikov (1979)

  Chirikov, B. V. 1979, Phys. Rep., 52, 263

• Deck et al. (2012)

  Deck, K. M., Holman, M. J., Agol, E., Carter, J. A., Lissauer, J. J., Ragozzine, D., & Winn, J. N. 2012, ApJ, 755, L21

• Duncan et al. (1989)

  Duncan, M., Quinn, T., & Tremaine, S. 1989, Icarus, 82, 402

• Fabrycky et al. (2012)

  Fabrycky, D. C., et al. 2012, ArXiv e-prints

• Fang & Margot (2012)

  Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92

• Fang & Margot (2013)

  —. 2013, ApJ, 767, 115

• Ferraz-Mello (2007)

  Ferraz-Mello, S. 2007, Astrophysics and Space Science Library, Vol. 345, Canonical Perturbation Theories - Degenerate Systems and Resonance (Springer)

• Gladman (1993)

  Gladman, B. 1993, Icarus, 106, 247

• Henrard & Lamaitre (1983)

  Henrard, J., & Lamaitre, A. 1983, Celestial Mechanics, 30, 197

• Henrard et al. (1986)

  Henrard, J., Milani, A., Murray, C. D., & Lemaitre, A. 1986, Celestial Mechanics, 38, 335

• Kopparapu & Barnes (2010)

  Kopparapu, R. K., & Barnes, R. 2010, ApJ, 716, 1336

• Lecar et al. (2001)

  Lecar, M., Franklin, F. A., Holman, M. J., & Murray, N. J. 2001, ARA&A, 39, 581

• Lichtenberg & Lieberman (1992)

  Lichtenberg, A. J., & Lieberman, M. A. 1992, Regular and chaotic dynamics (Springer-Verlag)

• Marchal & Bozis (1982)

  Marchal, C., & Bozis, G. 1982, Celestial Mechanics, 26, 311

• Mardling (2008)

  Mardling, R. A. 2008, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 760, The Cambridge N-Body Lectures, ed. S. J. Aarseth, C. A. Tout, & R. A. Mardling, 59

• Milani & Nobili (1983)

  Milani, A., & Nobili, A. M. 1983, Celestial Mechanics, 31, 213

• Mudryk & Wu (2006)

  Mudryk, L. R., & Wu, Y. 2006, ApJ, 639, 423

• Murray & Dermott (1999)

  Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge University Press)

• Murray & Holman (1997)

  Murray, N., & Holman, M. 1997, AJ, 114, 1246

• Mustill & Wyatt (2012)

  Mustill, A. J., & Wyatt, M. C. 2012, MNRAS, 419, 3074

• Paardekooper et al. (2013)

  Paardekooper, S.-J., Rein, H., & Kley, W. 2013, MNRAS, submitted. arXiv:1304.4762

• Quillen (2011)

  Quillen, A. C. 2011, MNRAS, 418, 1043

• Rauch & Holman (1999)

  Rauch, K. P., & Holman, M. 1999, AJ, 117, 1087

• Sessin & Ferraz-Mello (1984)

  Sessin, W., & Ferraz-Mello, S. 1984, Celestial Mechanics, 32, 307

• Veras & Mustill (2013)

  Veras, D., & Mustill, A. J. 2013, MNRAS, in press. arXiv:1305.5540

• Walker & Ford (1969)

  Walker, G. H., & Ford, J. 1969, Physical Review, 188, 416

• Wisdom (1980)

  Wisdom, J. 1980, AJ, 85, 1122

• Wisdom (1986)

  —. 1986, Celestial Mechanics, 38, 175

• Wisdom & Holman (1991)

  Wisdom, J., & Holman, M. 1991, AJ, 102, 1528

• Wisdom & Holman (1992)

  —. 1992, AJ, 104, 2022

• Wisdom et al. (1996)

  Wisdom, J., Holman, M., & Touma, J. 1996, Fields Institute Communications, 10, 217

APPENDIX

A-1. Transformation of Coordinates

The numerical integrations evolve the equations of motion for the full, exact, Hamiltonian, while the theory predicts the extent of the resonances in terms of a different set of canonical variables. In order to compare the two sets of results, one must perform the sequence of canonical transformations between the two set of variables. This requires removing the short period terms via a generating function so that they do not appear in the averaged Hamiltonian. It is equivalent to averaging the full Hamiltonian over the fast angle λᵢ, at least at first order in ε, because the resulting Hamiltonian has the same functional form.

The deviation between the averaged and full trajectories is of order ε, and hence the canonical transformation between the two sets of variables should be very near the identity transformation. In most cases, it is entirely negligible. However, that is not always true when the two orbits are close to each other.

The full Hamiltonian near the m:m+1 resonance can be written to first order in eccentricities as

$$H =H_{0}(\mathbf{\Lambda})+\epsilon_{1}H_{1}(\mathbf{\Lambda},\mathbf{x},\mathbf{\lambda},\mathbf{y})\qquad H_{0} =-\frac{\mu_{1}}{2\Lambda_{1}^{2}}-\frac{\mu_{2}}{2\Lambda_{2}^{2}}\qquad \epsilon_{1} =\frac{m_{1}}{m_{\star}}\qquad H_{1} =-\frac{\mu_{2}}{\Lambda_{2}^{2}}\bigg[ \frac{f_{m+1,27}(\alpha_{res})}{\sqrt{\Lambda_{1}}}(\cos{\theta_{m}}x_{1}-\sin{\theta_{m}}y_{1})+\frac{f_{m+1,31}(\alpha_{res})}{\sqrt{\Lambda_{2}}}(\cos{\theta_{m}}x_{2}-\sin{\theta_{m}}y_{2})+$$ $$\sum_{k\gt 0,k\neq m}\frac{f_{k+1,27}(\alpha_{res})}{\sqrt{\Lambda_{1}}}(\cos{\theta_{k}}x_{1}-\sin{\theta_{k}}y_{1})+\frac{f_{k+1,31}(\alpha_{res})}{\sqrt{\Lambda_{2}}}(\cos{\theta_{k}}x_{2}-\sin{\theta_{k}}y_{2})+\sum_{j\gt 0}f_{j,1}(\alpha_{res})\cos{\phi_{j}}\bigg]$$

(A.1)

where all Laplace coefficients are evaluated at the resonance location αᵣₑₛ=(m/(m+1))²⁄³. The variables x and y are the same canonical eccentricity polar coordinates introduced in Equation (2). For notational simplicity, we have also defined θᵢ=(i+1)λ₂-iλ₁ and ϕⱼ=j(λ₂-λ₁).

Again note that f₂,₃₁ must include an indirect contribution of -2αᵣₑₛ, and that f₁,₁ includes an indirect contribution of -αᵣₑₛ. We have ignored two terms in the disturbing function which come from the indirect part. These terms are proportional to e₁ αᵣₑₛ and are independent of the Laplace coefficients. They depend on the combinations of the fast angles λ₂-2λ₁ and λ₂. We do not prove it now, but these terms will be negligible to the transformation and so we ignore them. We’ll return to this point later.

This expression can be written more simply as

$$H =H_{Kepler}(\Lambda_{1},\Lambda_{2})+\epsilon_{1}H_{R,FO,m}(\mathbf{z})+\sum_{k\gt 0,k\neq m}\epsilon_{1}H_{SP,FO,k}(\mathbf{z})+\sum_{k\gt 0}\epsilon_{1}H_{SP,circ,k}(\mathbf{z})$$

(A.2)

where R implies resonant, SP implies short period terms, FO implies a first order term, circ implies the term appears at zero eccentricity, and z denotes the phase space variables.

We seek a generating function of Type 2 which determines a near-identity transformation that removes the short period terms. We will write this generating function as

$$F_{2}(\lambda_{i},y_{i},\tilde{\Lambda}_{i},\tilde{x}_{i}) =\lambda_{i}\tilde{\Lambda}_{i}+y_{i}\tilde{x}_{i}+\epsilon_{1}\it{f}(\lambda_{i},y_{i},\tilde{\Lambda}_{i},\tilde{x}_{i})$$

(A.3)

where $\it{f}$ is yet to be determined.

Here the variables with tildes are the averaged canonical coordinates and momenta, while those without are the exact physical phase space coordinates. The old and new coordinates are related as

$$\Lambda_{i} =\frac{\partial F_{2}}{\partial\lambda_{i}}=\tilde{\Lambda}_{i}+\epsilon_{1}\frac{\partial\it{f}}{\partial\lambda_{i}}\qquad x_{i} =\frac{\partial F_{2}}{\partial y_{i}}=\tilde{x}_{i}+\epsilon_{1}\frac{\partial\it{f}}{\partial y_{i}}\qquad \tilde{\lambda}_{i} =\frac{\partial F_{2}}{\partial\tilde{\Lambda}_{i}}=\lambda_{i}+\epsilon_{1}\frac{\partial\it{f}}{\partial\tilde{\Lambda}_{i}}\qquad \tilde{y}_{i} =\frac{\partial F_{2}}{\partial\tilde{x}_{i}}=y_{i}+\epsilon_{1}\frac{\partial\it{f}}{\partial\tilde{x}_{i}}$$

(A.4)

The new (averaged) Hamiltonian is H̃(z̃)=H(z(z̃)) since the generating function is time independent. To first order in ε₁,

H( z ) =H( z ( z )) =H_0(_1+ε_1fλ_1,_2+ε_1fλ_2)+ε_1H_R,FO,m( z )+Σ_k> 0,k≠ mε_1H_SP,FO,k( z )+Σ_k> 0ε_1H_SP,circ,k( z ) =H_0(_1,_2)+ε_1( H_0(A,B) A _(A,B)=(_1,_2)fλ_1( z )+ H_0(A,B) B _(A,B)=(_1,_2)fλ_2( z )) + ε_1H_R,FO,m( z )+Σ_k> 0,k≠ mε_1H_SP,FO,k( z )+Σ_k> 0ε_1H_SP,circ,k( z )

$${}=H_{0}(\tilde{\Lambda}_{1},\tilde{\Lambda}_{2})+\epsilon_{1}(n_{0,1}\frac{\partial\it{f}}{\partial\lambda_{1}}+n_{0,2}\frac{\partial\it{f}}{\partial\lambda_{2}})+\epsilon_{1}H_{R,FO,m}(\tilde{\mathbf{z}})+\qquad \sum_{k\gt 0,k\neq m}\epsilon_{1}H_{SP,FO,k}(\tilde{\mathbf{z}})+\sum_{k\gt 0}\epsilon_{1}H_{SP,circ,k}(\tilde{\mathbf{z}})$$

(A.5)

where n₀,ᵢ is the unperturbed mean motion of planet i. Note that z and z̃ are interchangeable in all terms of order ε₁.

Now, we would like this averaged Hamiltonian to be equal to

$$\tilde{H} =H_{0}(\tilde{\Lambda}_{1},\tilde{\Lambda}_{2})+\epsilon_{1}H_{R,FO,m}(\tilde{\mathbf{z}})$$

(A.6)

i.e. with all short period terms removed. We must choose the function f to satisfy the relation:

$$+n_{0,2}\frac{\partial\it{f}}{\partial\lambda_{2}})=-\sum_{k\gt 0,k\neq m}H_{SP,FO,k}-\sum_{k\gt 0}H_{SP,circ,k}$$

(A.7)

The solution is given as

$$f =-\frac{\mu_{2}}{\tilde{\Lambda}_{2}^{2}}\bigg[ \sum_{k\gt 0,k\neq m}A_{0,k}\bigg( \frac{f_{k+1,27}(\alpha_{res})}{\sqrt{\tilde{\Lambda}_{1}}}(\sin{\theta_{k}}\tilde{x}_{1}+\cos{\theta_{k}}y_{1})+\frac{f_{k+1,31}(\alpha_{res})}{\sqrt{\tilde{\Lambda}_{2}}}(\sin{\theta_{k}}\tilde{x}_{2}+\cos{\theta_{k}}y_{2})\bigg) +\qquad \sum_{j\gt 0}B_{0,j}f_{j,1}(\alpha_{res})\sin{\phi_{j}}\bigg]$$

(A.8)

One can write this down almost immediately because the short period Hamiltonian can be written entirely in terms of cosine functions, i.e. in the form of H_SP=Σᵢ Cᵢ cos(βᵢ). This is true of the combinations xᵢ cosθ-yᵢ sinθ which can be written as √(2Pᵢ)cos(θ+pᵢ). The function f is simply the same function, scaled by constant coefficients, but with the cosines switched to sines: f=Σᵢ C′ᵢ sin(βᵢ). Then,

fλ_1 =Σ_k> 0,k≠ mA_0,kθ_kλ_1H_SP,FO,k+Σ_k> 0B_0,kH_SP,circ,kϕ_kλ_1 =Σ_k> 0,k≠ mA_0,k(-k)H_SP,FO,k+Σ_k> 0B_0,kH_SP,circ,k(-k) fλ_2 =Σ_k> 0,k≠ mA_0,kθ_kλ_2H_SP,FO,k+Σ_k> 0B_0,kH_SP,circ,kϕ_kλ_2 =Σ_k> 0,k≠ mA_0,k(k+1)H_SP,FO,k+Σ_k> 0B_0,kH_SP,circ,k(k)

(A.9)

Using the derivatives of f given in Equation (A-1) in Equation (A.7) then results in an equation for the A₀,ₖ and B₀,ₖ:

$$\sum_{k\gt 0,k\neq m}\qquad A_{0,k} =\frac{-1}{(k+1)n_{0,2}-kn_{0,1}}=\frac{1}{(k+1)n_{0,2}(X_{k}-1)}\qquad B_{0,k} =\frac{-1}{k(n_{0,2}-n_{0,1})}=\frac{1}{kn_{0,2}(P_{2}/P_{1}-1)}\qquad X_{k} =\frac{kP_{2}}{(k+1)P_{1}}$$

(A.10)

where the first equation must hold for each k individually.

To find the transformation between the variables, we use Equation (A-1). First, let us determine how to map initial period ratios used in the numerical integrations to “averaged” period ratios. The $\it{real}$ period ratio of the planets is given by $P_{2}/P_{1}\bigg\lvert _{R}=(\Lambda_{2}/\Lambda_{1})^{3}\ast (m_{1}/m_{2})^{3}$, while the averaged period ratio is given by $P_{2}/P_{1}\bigg\lvert _{A}=(\tilde{\Lambda}_{2}/\tilde{\Lambda}_{1})^{3}\ast (m_{1}/m_{2})^{3}$. Therefore,

$$P_{2}/P_{1}\bigg\lvert _{R} \approx P_{2}/P_{1}\bigg\lvert _{A}\bigg[ 1+3\epsilon_{1}\frac{\partial\it{f}}{\partial\lambda_{2}}/\Lambda_{2}\bigg] \bigg[ 1+3\epsilon_{1}\frac{\partial\it{f}}{\partial\lambda_{1}}/\Lambda_{1}\bigg] ^{-1}\qquad {} \approx P_{2}/P_{1}\bigg\lvert _{A}\bigg[ 1-3\epsilon_{1}\frac{\partial\it{f}}{\partial\lambda_{1}}/\Lambda_{1}+3\epsilon_{1}\frac{\partial\it{f}}{\partial\lambda_{2}}/\Lambda_{2}\bigg]$$

(A.11)

The required derivatives of f are:

fλ_1/_1 =Σ_k> 0,k≠ m(k)/(k+1)_2_11X_k-1[ f_k,27e_1(θ_k-_1)+f_k,31e_2(θ_k-_2)] +Σ_k> 0f_k,0P_2/P_1-1_2_1ϕ_k fλ_2/_2 =Σ_k> 0,k≠ m-1X_k-1[ f_k,27e_1(θ_k-_1)+f_k,31e_2(θ_k-_2)] +Σ_k> 0-f_k,0P_2/P_1-1ϕ_k

(A.12)

so that

$$P_{2}/P_{1}\bigg\lvert _{R} =P_{2}/P_{1}\bigg\lvert _{A}\times\qquad \bigg\lbrace 1-3\epsilon_{1}\sum_{k\gt 0,k\neq m}\bigg[ \frac{k}{k+1}\frac{\Lambda_{2}}{\Lambda_{1}}+1\bigg] \frac{1}{X_{k}-1}\bigg[ f_{k,27}e_{1}\cos{(\theta_{k}-\varpi_{1})}+f_{k,31}e_{2}\cos{(\theta_{k}-\varpi_{2})}\bigg]\qquad {}-3\epsilon_{1}\sum_{k\gt 0}\frac{f_{k,0}}{P_{2}/P_{1}-1}\cos{\phi_{k}}(1+\Lambda_{2}/\Lambda_{1})\bigg\rbrace$$

(A.13)

where Λ₂/Λ₁=m₂/m₁(P₂/P₁)¹⁄³. Again, we can neglect the difference between real and averaged coordinates at order ε₁.

There are two reasons why this shift is not negligible. The first is that the period ratio of the planets we are considering is near unity, so the divisors in the coefficients A₀,ₖ and B₀,ₖ are small. Though the functions proportional to A₀,ₖ also are linear in the eccentricities, the functions fₖ₊₁,₂₇ and fₖ₊₁,₃₁ are larger when evaluated at α near unity.

In order to remove a short period term in the Hamiltonian which is periodic in the angle pλ₂-qλ₁, the transformation must include a term proportional to the coefficient of that term and divided by pn₀,₂-qn₀,₁. This is why the two indirect terms period in the combinations λ₂-2λ₁ and λ₂ are negligible. They appear at first order in eccentricities, but with divisors of n₀,₂-2n₀,₁ and n₀,₂, respectively, neither of which are small. Moreover, they do not depend on the Laplace coefficients, which makes them even smaller compared to the terms in the transformation coming from the direct terms of first order in eccentricities (and low k).

We only require the transformation from real coordinates to average coordinates at the initial time of the numerical integrations. Hence, for all parameters we substitute their initial values to determine the mapping between real and averaged variables. For example, a numerical simulation using e₂=0 and all of the angles equal to zero except for λ₁=M₁ would require a transformation of

$$P_{2}/P_{1}\bigg\lvert _{R} =P_{2}/P_{1}\bigg\lvert _{A}\times\qquad \bigg\lbrace 1-3\epsilon_{1}\sum_{k\gt 0,k\neq m}\bigg[ \frac{k}{k+1}\frac{\Lambda_{2}}{\Lambda_{1}}+1\bigg] \frac{1}{X_{k}-1}\bigg[ f_{k,27}e_{1}(0)\cos{kM_{1}}\bigg] -3\epsilon_{1}\sum_{k\gt 0}\frac{f_{k,0}}{P_{2}/P_{1}-1}\cos{kM_{1}}(1+\Lambda_{2}/\Lambda_{1})\bigg\rbrace$$

(A.14)

As the period ratio of the planets grows closer to unity, the convergence rates of the two sums are much slower. When e₁ is small, the most important contribution is from the second sum. Each term will have the same sign as coskM₁ does since fₖ,₀>0 and P₂/P₁>1. It is only when M₁=0 that these terms all have the same sign and add coherently. In this case, the period ratio in “real” variables is smaller than the period ratio in averaged variables by a maximum amount. When M₁ is nonzero, the shift is much smaller since terms will not add coherently. This behavior is reflected in Figures 5, 6, and 9. This incoherent adding of terms occurs even for the sum that appears linearly in e₁(0), and so the shift is should be small regardless of e₁(0) if M₁≠0.

The other transformations are straightforward to derive given the generating function. Of interest to us in particular is the transformation of the eccentricities, since we plot the widths of the resonances as a function of eccentricities. We will just give the result here for the case where e₂=0 and all angles equal to zero. In this special case,

$$\frac{\partial f}{\partial\tilde{x}_{i}} \propto\sin{((k+1)\lambda_{2}-k\lambda_{1})}=0\qquad \frac{\partial f}{\partial y_{i}} \propto\cos{((k+1)\lambda_{2}-k\lambda_{1})}=1$$

(A.15)

So to first order in ε₁, yᵢ=ỹᵢ, and

$$x_{1} =\tilde{x}_{1}+\epsilon_{1}\frac{\Lambda_{2}}{\sqrt{\Lambda_{1}}}\sum_{k\gt 0,k\neq m}\frac{\lvert f_{k,27}\rvert}{(X_{k}-1)(k+1)}\qquad x_{2} =\tilde{x}_{2}-\epsilon_{1}\sqrt{\Lambda_{2}}\sum_{k\gt 0,k\neq m}\frac{f_{k,31}}{(X_{k}-1)(k+1)}$$

(A.16)

Note that these expressions are already of order (e), so no e appears in the sum.

This correction is negligible, however. Although the terms in the sums in Equation (A-1) look very similar to those in the correction to the period ratio, they are weighted by a factor of 1/(k+1), and this makes the correction much smaller. Taken together with the relationship yᵢ=ỹᵢ, this implies that both eccentricities and longitudes of periastron are not affected greatly by the transformation. Moreover, this should be true regardless of if e₂ is nonzero or if the angles are nonzero because the weighting factor 1/(k+1) will always appear.

Finally, we check the correction to the mean longitudes λᵢ. We do this only for our nominal case, with e₂=0 and all other angles equal to zero. In this case, yᵢ=0 and all the ϕⱼ and θₖ are equal to zero, so all derivatives with respect to Λ̃ᵢ are zero. In other words, λᵢ=λ̃ᵢ in this case.

A-2. Effect of Integration Time on the Estimated Lyapunov Times

Although numerical integrations can determine if an orbit is chaotic, they cannot prove that an orbit is quasiperiodic. Longer integrations provide better estimates of the Lyapunov time and give a better idea of what orbits could be quasiperiodic. We performed a set of numerical integrations only varying the total integration time order to check the convergence of our estimate of the Lyapunov time. In Figure A.1 we show the structure of the chaotic zone, in terms of estimated Lyapunov times, for integration lengths of 10⁶, 10⁷, and 10⁸ days. The three panels agree very well, and give us confidence that 10⁶ days is an adequate amount of time to reliably estimate Lyapunov times for these types of orbits.

There are two minor differences between the three sets. First, we see that some of the orbits around the second order resonances are chaotic, though the shorter integration did not flag them as such. This is not uncommon; the tangent vector corresponding to a chaotic orbit exponentially grows at different rates, depending on the region of the chaotic zone. For example, the local estimate of the Lyapunov time can be very long when the orbit is near invariant curves associated with resonances.

Secondly, the estimate of the Lyapunov time slightly shifts to longer times during longer integrations. Our hypothesis is that this is caused by instabilities which set in on timescales of 10⁷ or 10⁸ days (but not 10⁶ days). For example, orbits which undergo strong scattering events can lead to a system of a single bound planet and a second planet on a hyperbolic orbit. Though the orbits would still be formally chaotic, the motion is closer to regular after the scattering (the perturbations are weaker). We tracked the log of the norm of the tangent vector for a system which underwent scattering and indeed found that the tangent vector stopped growing exponentially after the scattering event, and hence longer integrations of this orbit will lead to longer and longer estimates of the Lyapunov time.

Even in cases without planetary ejection, Lagrange instabilities cause the planets to become more and more widely spaced. Again, the perturbations between the planets decreases. We also tracked the log of the norm of the tangent vector for a system undergoing Lagrange instabilities, and found that the estimated Lyapunov time grew as the period ratio of the pair grew further from unity. The longer integrations indicate a longer Lyapunov time because a larger fraction of the integration time is spent at larger period ratios.

Figure A.1.— Effect of the integration time on the estimated Lyapunov times. Top: 10⁶ day integration, with e₂ initially zero, all angles set to zero except M₁=π/2, ζ=1, and εₚ=2× 10⁻⁵. Only the chaotic orbits are shown in color. White regions contain quasiperiodic orbits. Middle: the same initial conditions, but evolved only for 10⁷ days. Bottom: the same initial conditions, but evolved for 10⁸ days. 10⁶ days was our standard integration length.

A-3. Results for a Case with e₂ Nonzero

All of our tests used e₂=0, but the analytic theory predicts that only the weighted eccentricity, and not e₁ and e₂ individually, matters for the structure of the first order resonances at low eccentricity. We tested this numerically for a single case with e₂=0.02. These integrations had all angles initially set to zero except $\varpi_{1}=\pi/2$. In this case, $\Delta\varpi=\pi/2$ and so σ=√(e₂²+e₂²/R²).

For these choices of angles, one can show that the generalized longitude of pericenter $\psi=\arctan{\lbrace s_{1}/r_{1}\rbrace }=\arctan{\lbrace Re_{1}/e_{2}\rbrace }$. Given the values of e₁ and e₂, ψ ranges from 0 to $\sim 79^{\circ}$. The value of the resonant angle ϕ, then, is $\phi=-m\varpi_{1}+\psi=-m\pi/2+\psi$. Recall that ϕ=π corresponds to the center of the resonance, and so if the resonance overlap criterion is not satisfied resonances with ϕ=π may appear as regular regions (whether or not they do also depends on e₁ and εₚ). When ϕ is closer to zero the orbits are more likely to be chaotic (especially for period ratios closer to unity). Given the range of ψ, it is clear that only resonances with m=2,6,…, or m=3,7,… should appear regular. These resonances are the 3:2, the 4:3, the 7:6, and the 8:7. The 3:2 does not lie within the range of period ratios shown in Figure A.2, but the others do, they are the only regular first order resonant regions that appear. The predicted seperatrices match up well with the boundaries of these regular regions.

Figure A.2.— The chaotic structure of phase space when e₂ is nonzero and fixed at 0.02. It does not matter what e₁ and e₂ are individually, only the weighted eccentricity σ matters in determining the widths of the first order resonances. This figure shows the chaotic regions of phase space for a set of integrations where e₁ varied from 0.0 to 0.1. Besides e₁, only the period ratio of the planets varied. The masses of the planets were set as ε₁=ε₂=10⁻⁵. These result shown are what we expect if σ is indeed the relevant parameter.