FIRST-ORDER RESONANCE OVERLAP AND THE STABILITY OF CLOSE TWO-PLANET SYSTEMS

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

$$\begin{eqnarray} \frac{d X}{d {t}^{\prime }} = 0 && \longrightarrow Y\left (\frac{1}{2}(X^2+Y^2)-\Gamma ^{\prime }\right) = 0 \nonumber \\ &&\longrightarrow Y=0 \nonumber \\ \frac{d Y}{d {t}^{\prime }} = 0 &&\longrightarrow X\left (\frac{1}{2}(X^2+Y^2)-\Gamma ^{\prime }\right) +1 = 0 \nonumber \\ &&\longrightarrow X^3-2\Gamma ^{\prime }X+2 =0. \end{eqnarray} \tag{ 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 (38) is Δ = 32(Γ'³ − 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray} \lambda _2 + m \Delta \lambda + \psi &=& (2 k+1) \pi \nonumber \\ \Delta \lambda &=& \frac{(2k+1)\pi -\psi -\lambda _2}{m}, \end{eqnarray} \tag{ 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 in Figure 2 for initially zero eccentricity orbits.

Zoom In Zoom Out Reset image size

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_sep. 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_sep = 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

$$\begin{eqnarray} X_{\star 1 } & =& -X_3+2/\sqrt{X_3} \nonumber \\ X_{\star 2 } & =& -X_3-2/\sqrt{X_3}, \end{eqnarray} \tag{ 40 }$$

and that

$$\begin{eqnarray} &&\Delta X = X_{\star 1}-X_{\star 2} = 4/\sqrt{X_3}. \end{eqnarray} \tag{ 41 }$$

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

$$\begin{eqnarray} \Delta \Psi ^{\prime } &=& \Psi ^{\prime }_{\star 1} - \Psi ^{\prime }_{\star 2} \nonumber \\ & =&{\frac{1}{2}}(X_{\star 1}-X_{\star 2})\vert X_{\star 1} + X_{\star 2} \vert = 4 \sqrt{X_3} \nonumber \\ \Delta (\delta \Theta ^{\prime })& = &\Delta \Psi ^{\prime }, \end{eqnarray} \tag{ 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 (24). 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 _p, 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 $\Psi ^{\prime }_{\star 1} = X_{ \star 1}^2/2$ passes through zero, X_⋆1 corresponds to the left boundary of the separatrix (at smaller period ratios) in Figure 3.

Zoom In Zoom Out Reset image size

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 separatrix). 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 separatrix can only pass through zero eccentricity when X_⋆1 = 0. Using Equation (40), this implies that X₃ = 4^1/3. Therefore, the width of the resonance at zero eccentricity is $\Delta (\delta \Theta ^{\prime }) = 4 \sqrt{X_3} = 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 (38), we find that the value of Γ' corresponding to zero eccentricity on the separatrix is Γ' = 3/4^1/3 = 1.88988. Since Ψ' = 0 at this point, δΘ' = −1.88988 at zero eccentricity on the separatrix. 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),

$$\begin{eqnarray} s =\sqrt{\frac{\alpha }{\alpha _{{\rm res}}}}& = \frac{1-\delta \Theta ^{\prime } Q m/(\alpha _{{\rm res}}/\zeta +1)}{1+\delta \Theta ^{\prime } Q m/(1+\zeta /\alpha _{{\rm res}})}. \end{eqnarray} \tag{ 43 }$$

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

$$\begin{eqnarray} s^2 & \sim & 1-\delta \Theta ^{\prime } c_0 \epsilon _p^{2/3} + O \left(\epsilon _p^{4/3}\right) \nonumber \\ c_0 & =& 2 \left [\frac{ f_{31}^2 (R^2 +\sqrt{\alpha _{\rm res}} \zeta) (\alpha _{\rm res}+\zeta)}{ 9 \sqrt{\alpha _{\rm res}}(m+1)(1+\zeta)^2} \right]^{1/3}. \end{eqnarray} \tag{ 44 }$$

Therefore, the width of the resonance δα_m/α_res is

$$\begin{eqnarray} \frac{\delta \alpha _m}{\alpha _{{\rm res}}} &=& s^2(\delta \Theta ^{\prime }_1) -s^2(\delta \Theta ^{\prime }_2) \nonumber \\ &=& \Delta (\delta \Theta ^{\prime }) \epsilon _p^{2/3}c_0 + O\left(\epsilon _p^{4/3}\right). \end{eqnarray} \tag{ 45 }$$

Note that if the resonance bounds were symmetric, there would be no contribution at $O(\epsilon _p^{4/3})$, since that term depends on $(\delta \Theta ^{^{\prime }2}_2 - \delta \Theta ^{^{\prime }2}_1)$.

When overlap occurs, α is near unity, so we will replace R² with 1 and α_res with 1. Moreover, we will replace f_{m + 1, 31} with 0.802m + 0.87 ∼ 0.802m for large m. This approximation assumes that 1/m is small. We will also neglect the $\epsilon _p^{4/3}$ terms in this approximation and will justify it afterward. Then the width of the resonance for close orbits (high m) is

$$\begin{eqnarray} \frac{\delta \alpha }{\alpha _{{\rm res}}} &=&\frac{2}{9^{1/3}}\Delta (\delta \Theta ^{\prime }) \epsilon _p^{2/3}\left [\frac{ (0.802m)^2 (1 +\zeta) (1+\zeta)}{ (m+1)(1+\zeta)^2} \right ]^{1/3} \nonumber \\ & \approx & 2\left (\frac{0.802}{3}\right)^{2/3}\Delta (\delta \Theta ^{\prime }) \epsilon _p^{2/3}\left [\frac{ m^2}{ (m+1)} \right ]^{1/3} \nonumber \\ & \approx & 4.18 \epsilon _p^{2/3} m^{1/3}, \end{eqnarray} \tag{ 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

$$\begin{eqnarray} \Delta \alpha & =& \frac{d}{dm}\alpha _{{\rm res}}(m) \Delta m \nonumber \\ & =& \frac{d}{d m} \left (\frac{m}{m+1}\right)^{2/3} (1) \nonumber \\ & =& \alpha _{{\rm res}}(m) \frac{2}{3 m^2} \frac{m}{m+1} \nonumber \\ & \sim & \alpha _{{\rm res}}(m) \frac{2}{3 m^2} \nonumber \\ \frac{\Delta \alpha }{\alpha _{{\rm res}}} &\sim & \frac{2}{3 m^2}. \end{eqnarray} \tag{ 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

$$\begin{eqnarray} \frac{\Delta \alpha }{\alpha _{{\rm res},m}} &\approx & 0.5\times \left (\frac{\delta \alpha }{\alpha _{{\rm res},m }}+\frac{\delta \alpha }{\alpha _{{\rm res},m+1}} \right)\nonumber \\ \frac{2}{3 m^2} & \sim & 4.18 \epsilon _p^{2/3} m^{1/3} \nonumber \\ m & \sim & 0.455 \epsilon _p^{-2/7}. \end{eqnarray} \tag{ 48 }$$

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

$$\begin{eqnarray} \frac{P_2}{P_1} &\lesssim & 1+\frac{1}{m} \nonumber \\ &\lesssim & 1+ 2.2\epsilon _p^{2/7}, \end{eqnarray} \tag{ 49 }$$

or equivalently

$$\begin{eqnarray} &&\frac{a_2-a_1}{a_1}\lesssim 1.46 \epsilon _p^{2/7}, \end{eqnarray} \tag{ 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 (46), we would have found $m\sim 0.51 \epsilon _p^{-2/7}$, or $(a_2-a_1)/a_1 \lesssim 1.3 \epsilon _p^{2/7}$ 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 separatrix, though this effect is probably less important than the first.

In the above calculation, we kept terms of order 1/m² and $\epsilon _p^{2/3} m^{1/3}$, and we neglected terms of order $\epsilon _p^{2/3}/m^{2/3}, 1/m^3$, and $\epsilon _p^{4/3} m^{1/3}$. For Jupiter mass planets, _p ∼ 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 $\epsilon _p^{4/3}$ terms will be smaller than that incurred by neglecting $\epsilon _p^{2/3}{/}m^{2/3}$ and 1/m³ terms. This justifies ignoring the $\epsilon _p^{4/3}$ terms in the above calculation. However, it is clear that the approximation of α_res = 1, R² = 1, and 1/m ≪ 1 becomes less accurate for high-mass planets in the regime where complete resonance overlap occurs. We will study these effects numerically.

We now briefly address developing a resonance overlap criterion as a function of eccentricity (while still assuming that 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 _p = 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 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.

Zoom In Zoom Out Reset image size

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 (38) in this case reduces to X³ − 2Γ'X + 2 ≈ X³ − 2Γ'X = 0. This admits X = 0 and $X=\pm \sqrt{2 \Gamma ^{\prime }}$ as solutions. The resonance center occurs at $X_1\approx -\sqrt{2 \Gamma ^{\prime }}$, the other center at X₂ ≈ 0, and the unstable fixed point at $X_3 \approx \sqrt{2 \Gamma ^{\prime }}$. In this limit, then, the width of the resonance in terms of δΘ' is $\Delta (\delta \Theta ^{\prime }) = 4 \sqrt{X_3} \approx 4(2\Gamma ^{\prime })^{1/4}$.

Since the location of the center of resonance X₁ grows in magnitude as Γ' increases, the forced resonant eccentricity $\sigma _{{\rm res}} \propto \sqrt{\Psi ^{\prime }_1} \propto \vert X_1 \vert \sim \sqrt{2 \Gamma ^{\prime }}$ 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 separatrix with the X-axis (see Figure 1). These two crossings in X, denoted X_⋆1 and X_⋆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 an initial value of σ to be the width of the separatrix assuming σ = σ_⋆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

$$\begin{eqnarray} \Psi ^{\prime }_{\star 1} & =& \frac{X_{\star 1}^2}{2} = \frac{1}{2} (X_3 -2/\sqrt{X_3})^2 \nonumber \\ & \approx & \Gamma ^{\prime } -2(2 \Gamma ^{\prime })^{1/4} +\frac{\sqrt{2}}{\sqrt{\Gamma ^{\prime }}} \end{eqnarray} \tag{ 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

$$\begin{eqnarray} \Delta (\delta \Theta ^{\prime }) & \approx & 4(2g_m R^2 \sigma ^2/Q)^{1/4}. \end{eqnarray} \tag{ 52 }$$

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

$$\begin{eqnarray} m & \approx & 0.482 (\epsilon _p \sigma)^{-1/5} \end{eqnarray} \tag{ 53 }$$

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

$$\begin{eqnarray} \frac{P_2}{P_1} &\lesssim & 1+2.08 (\epsilon _p \sigma)^{1/5} \end{eqnarray} \tag{ 54 }$$

and

$$\begin{eqnarray} &&\frac{a_2-a_1}{a_1} \lesssim 1.38 (\epsilon _p \sigma)^{1/5}. \end{eqnarray} \tag{ 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

$$\begin{eqnarray} &&m \sim 0.482 (\epsilon _p \sigma)^{-1/5} \gtrsim 0.455 \epsilon _p^{-2/7} \end{eqnarray} \tag{ 56 }$$

or

$$\begin{eqnarray} &&\sigma \gtrsim 1.33 \epsilon _p^{3/7}, \end{eqnarray} \tag{ 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 the total planetary mass relative to the host star of _p = 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 those 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 in Figure 4. The curves plotted come from Equation (54) and take the form $\log _{10}\lbrace \sigma \rbrace = 5\log _{10}\lbrace P_2/P_1-1\rbrace +\log _{10}\lbrace C^{-5} \epsilon _p^{-1}\rbrace$, where C = 2.08 (or 2.7; using Mustill & Wyatt 2012). However, the actual slope is closer to 4.6 in the _p = 2 × 10⁻⁵ case, and we found it to be closer to 4.2 in the _p = 2 × 10⁻⁴ case. Therefore, the scaling of P₂/P₁ − 1∝(_pσ)^1/5 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).

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.

Zoom In Zoom Out Reset image size

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_2/P_1 = 1.\bar{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 separatrices 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).

Zoom In Zoom Out Reset image size

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.

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 (39), with e₂(0) = 0 and all angles zero except λ₁, yields the location of the resonant islands in terms of λ₁:

$$\begin{eqnarray} &&\lambda _1 \approx M_1 = -\frac{2k\pi }{m}, \end{eqnarray} \tag{ 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, appears 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 separatrices 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 does 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 Δϖ were initially zero for all of these integrations, the individual values of e₁, e₂, and Δϖ 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 Δϖ = π/2; see the Appendix for details.

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 _p, 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.

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

The resonance widths grow with the total mass of the planets, relative to the host star, as $\epsilon _p^{2/3}$, and so we expect that the overall extent of the chaotic web should grow with _p as well. Figure 9 shows the chaotic structure of the phase space in the period ratio—M₁ plane for three different values of _p (the bottom panel, with _p = 2 × 10⁻⁵, is the same as upper panel of Figure 6). The eccentricity of the inner orbit e₁ was fixed at 0.01, and 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.

Zoom In Zoom Out Reset image size

At the highest masses shown (_p = 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 separatrices of the 3:2 and 5:3 resonances.

The effect of the transformation between real and averaged coordinates grows ∝_p (see Equation (A14)). As a result, the shift at M₁ = 0 is much more extreme in the _p = 2 × 10⁻⁴ case compared to the _p = 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.

How well do the derived resonance overlap criteria given in Equations (50) and (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₁₀{_p}, 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 _p and the semimajor axis of the inner planet were varied. Recall that the criterion predicts, given _p, 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 _p, 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.

Zoom In Zoom Out Reset image size

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.

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 _p. We also show the estimate of Mustill & Wyatt (2012), which lies parallel to the estimate derived in this work—both have slopes of five, but they 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.

Zoom In Zoom Out Reset image size

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).

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 tracks 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.

Zoom In Zoom Out Reset image size

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.

We find that the estimated Lyapunov times 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 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.

Zoom In Zoom Out Reset image size

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 versus σ (for example, Figure 5), other slices of the phase space (particularly in period ratio versus 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 hint 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) or 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.

Zoom In Zoom Out Reset image size

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.

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

$$\begin{eqnarray} &&\frac{a_2-a_1}{a_1} \gtrsim 2.4 \epsilon _p^{1/3}. \end{eqnarray} \tag{ 59 }$$

We have found that for initially circular orbits, first-order resonances will overlap and result in a chaotic phase space if $(a_2-a_1)/a_1\lesssim 1.46\epsilon _p^{2/7}$. 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 _p = (1.46/2.4)²¹. Both boundaries are plotted in Figure 10. For planets with _p 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 _p 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 (_p ≲ 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_2-a_1)/a_1\lesssim 1.46\epsilon _p^{2/7}$) 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 does 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

$$\begin{eqnarray} &&\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 > \frac{p}{a} \vert _{{\rm crit}}, \end{eqnarray} \tag{ 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, respectively (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 _p = 2 × 10⁻⁵. The dotted contour is the Hill boundary; above this curve all orbits fail the Hill criterion, while below it all orbits are Hill stable.

Zoom In Zoom Out Reset image size

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 (1) moderately eccentric orbits of planets of any _p and (2) higher mass systems (_p ≳  few  × 10⁻⁵) 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 ζ.

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?

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 _p = 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.

Zoom In Zoom Out Reset image size

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 (_p ∼ 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 in 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 Equation (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 provides 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.

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 λ_i, 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

$$\begin{eqnarray} H &=& H_0(\mathbf {\Lambda })+\epsilon _1 H_1(\mathbf {\Lambda },\mathbf {x},\mathbf {\lambda },\mathbf {y}) \nonumber \\ H_0 & =& -\frac{\mu _1}{2\Lambda _1^2} -\frac{\mu _2}{2\Lambda _2^2} \nonumber \\ \epsilon _1 & =& \frac{m_1}{m_\star } \nonumber \\ H_1 & =& -\frac{\mu _2}{\Lambda _2^2}\left [ \frac{f_{m+1,27}(\alpha _{{\rm res}})}{\sqrt{\Lambda _1}} (\cos {\theta _m} x_1 -\sin {\theta _m} y_1)\right. \nonumber\\ &&+\frac{f_{m+1,31}(\alpha _{{\rm res}})}{\sqrt{\Lambda _2}} (\cos {\theta _m} x_2 -\sin {\theta _m} y_2) \nonumber \\ &&+ \sum _{k>0,k \ne m} \frac{f_{k+1,27}(\alpha _{{\rm res}})}{\sqrt{\Lambda _1}} (\cos {\theta _k} x_1 -\sin {\theta _k} y_1)\nonumber\\ &&+\frac{f_{k+1,31}(\alpha _{{\rm res}})}{\sqrt{\Lambda _2}} (\cos {\theta _k} x_2 -\sin {\theta _k} y_2)\nonumber\\ &&\left. + \sum _{j>0}f_{j,1}(\alpha _{{\rm res}})\cos {\phi _j} \right], \end{eqnarray} \tag{ A1 }$$

where all Laplace coefficients are evaluated at the resonance location α_res = (m/(m + 1))^2/3. The variables x and y are the same canonical eccentricity polar coordinates introduced in Equation (13). For notational simplicity, we have also defined θ_i = (i + 1)λ₂ − iλ₁ and ϕ_j = j(λ₂ − λ₁).

Again note that f_{2, 31} must include an indirect contribution of −2α_res, and that f_{1, 1} includes an indirect contribution of −α_res. We have ignored two terms in the disturbing function which come from the indirect part. These terms are proportional to e₁α_res 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 will return to this point later.

This expression can be written more simply as

$$\begin{eqnarray} H &=& H_{{\rm Kepler}}(\Lambda _1,\Lambda _2)+ \epsilon _1H_{R,{\rm FO},m }(\mathbf {z})\nonumber\\ &&+ \sum_{k>0,k \ne m} \epsilon _1 H_{{\rm SP},{\rm FO},k}(\mathbf {z})+ \sum _{k>0} \epsilon _1H_{{\rm SP},{\rm circ},k}(\mathbf {z}), \end{eqnarray} \tag{ A2 }$$

where R implies resonant, SP implies short-period terms, FO implies a first-order term, circ implies that 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

$$\begin{eqnarray} &&F_2(\lambda _i,y_i,\tilde{\Lambda }_i,\tilde{x}_i) = \lambda _i \tilde{\Lambda }_i + y_i\tilde{x}_i + \epsilon _1 {f}(\lambda _i,y_i,\tilde{\Lambda }_i,\tilde{x}_i), \end{eqnarray} \tag{ A3 }$$

where ${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

$$\begin{eqnarray} \Lambda _i & =& \frac{\partial F_2}{\partial \lambda _i} = \tilde{\Lambda }_i+\epsilon _1\frac{\partial {f}}{\partial \lambda _i} \nonumber \\ x_i & =& \frac{\partial F_2}{\partial y_i} = \tilde{x}_i+\epsilon _1\frac{\partial {f}}{\partial y_i} \nonumber \\ \tilde{\lambda }_i & =& \frac{\partial F_2}{\partial \tilde{\Lambda }_i} = \lambda _i+\epsilon _1\frac{\partial {f}}{\partial \tilde{\Lambda }_i} \nonumber \\ \tilde{y}_i & = &\frac{\partial F_2}{\partial \tilde{x}_i} = y_i+\epsilon _1\frac{\partial {f}}{\partial \tilde{x}_i}. \end{eqnarray} \tag{ A4 }$$

The new (averaged) Hamiltonian is $\tilde{H} (\tilde{\mathbf {z}}) = H(\mathbf {z}(\tilde{\mathbf {z}}))$ since the generating function is time independent. To first order in ₁,

$$\begin{eqnarray} \tilde{H}(\tilde{\mathbf {z}}) &=& H(\mathbf {z}(\tilde{\mathbf {z}})) \nonumber \\ &=& H_0\left(\tilde{\Lambda }_1+\epsilon _1\frac{\partial {f}}{\partial \lambda _1},\tilde{\Lambda }_2+\epsilon _1\frac{\partial {f}}{\partial \lambda _2}\right) + \epsilon _1 H_{R,{\rm FO},m }(\tilde{\mathbf {z}})\nonumber\\ &&+\sum _{k>0,k \ne m} \epsilon _1 H_{{\rm SP},{\rm FO},k}(\tilde{\mathbf {z}})+ \sum _{k>0} \epsilon _1 H_{{\rm SP},{\rm circ},k}(\tilde{\mathbf {z}}) \nonumber \\ &=& H_0(\tilde{\Lambda }_1, \tilde{\Lambda }_2) + \epsilon _1 \left (\frac{\partial H_0(A,B)}{\partial A} \right \vert _{(A,B) = (\tilde{\Lambda }_1, \tilde{\Lambda }_2)} \frac{\partial {f}}{\partial \lambda _1} (\tilde{\mathbf {z}})\nonumber\\ &&+ \frac{\partial H_0(A,B)}{\partial B} \left \vert _{(A,B) = (\tilde{\Lambda }_1, \tilde{\Lambda }_2)} \frac{\partial {f}}{\partial \lambda _2} (\tilde{\mathbf {z}}) \right)\nonumber \\ && + \epsilon _1H_{R,{\rm FO},m }(\tilde{\mathbf {z}})+\sum _{k>0,k \ne m} \epsilon _1H_{{\rm SP},{\rm FO},k}(\tilde{\mathbf {z}})\nonumber\\ &&+ \sum _{k>0} \epsilon _1 H_{{\rm SP},{\rm circ},k}(\tilde{\mathbf {z}}) \nonumber \\ & =& H_0(\tilde{\Lambda }_1, \tilde{\Lambda }_2) + \epsilon _1 \left(n_{0,1} \frac{\partial {f}}{\partial \lambda _1} +n_{0,2} \frac{\partial {f}}{\partial \lambda _2} \right) +\epsilon _1 H_{R,{\rm FO},m }(\tilde{\mathbf {z}})\nonumber \\ &&+\sum _{k>0,k \ne m} \epsilon _1 H_{{\rm SP},{\rm FO},k}(\tilde{\mathbf {z}})+\sum _{k>0} \epsilon _1H_{{\rm SP},{\rm circ},k}(\tilde{\mathbf {z}}), \end{eqnarray} \tag{ A5 }$$

where n_{0, i} is the unperturbed mean motion of planet i. Note that z and $\tilde{\mathbf {z}}$ are interchangeable in all terms of order ₁.

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

$$\begin{eqnarray} &&\tilde{H} = H_0(\tilde{\Lambda }_1, \tilde{\Lambda }_2) + \epsilon _1 H_{R,{\rm FO},m }(\tilde{\mathbf {z}}), \end{eqnarray} \tag{ A6 }$$

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

$$\begin{eqnarray} &&\left(n_{0,1} \frac{\partial {f}}{\partial \lambda _1} +n_{0,2} \frac{\partial {f}}{\partial \lambda _2}\right) = - \sum _{k>0,k \ne m} H_{{\rm SP},{\rm FO},k}- \sum _{k>0} H_{{\rm SP},{\rm circ},k}.\nonumber\\&& \end{eqnarray} \tag{ A7 }$$

The solution is given as

$$\begin{eqnarray} f &=& -\frac{\mu _2}{\tilde{\Lambda }_2^2}\left [ \sum _{k>0,k \ne m} A_{0,k} \left (\frac{f_{k+1,27}(\alpha _{{\rm res}})}{\sqrt{\tilde{\Lambda }_1}} (\sin {\theta _k} \tilde{x}_1 +\cos {\theta _k} y_1)\right. \right. \nonumber\\ &&\left. +\frac{f_{k+1,31}(\alpha _{{\rm res}})}{\sqrt{\tilde{\Lambda }_2}} (\sin {\theta _k} \tilde{x}_2 +\cos {\theta _k} y_2)\right) \nonumber \\ &&\left. +\sum _{j>0}B_{0,j} f_{j,1}(\alpha _{{\rm res}})\sin {\phi _j} \right]. \end{eqnarray} \tag{ A8 }$$

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 = ∑_iC_icos (β_i). This is true of the combinations x_icos θ − y_isin θ, which can be written as $\sqrt{2 P_i} \cos {(\theta +p_i)}$. The function f is simply the same function, scaled by constant coefficients, but with the cosines switched to sines: $f = \sum _i C^{\prime }_i \sin {(\beta _i)}$. Then,

$$\begin{eqnarray} \frac{\partial {f}}{\partial \lambda _1} & =& \sum _{k>0,k \ne m} A_{0,k} \frac{\partial \theta _k }{\partial \lambda _1} H_{{\rm SP},{\rm FO},k} + \sum _{k>0} B_{0,k}H_{{\rm SP},{\rm circ},k} \frac{\partial \phi _k }{\partial \lambda _1} \nonumber \\ & =& \sum _{k>0,k \ne m} A_{0,k} (-k) H_{{\rm SP},{\rm FO},k} + \sum _{k>0} B_{0,k} H_{{\rm SP},{\rm circ},k} (-k) \nonumber \\ \frac{\partial {f}}{\partial \lambda _2} & =& \sum _{k>0,k \ne m} A_{0,k} \frac{\partial \theta _k }{\partial \lambda _2} H_{{\rm SP},{\rm FO},k} + \sum _{k>0} B_{0,k}H_{{\rm SP},{\rm circ},k} \frac{\partial \phi _k }{\partial \lambda _2} \nonumber \\ & =& \sum _{k>0,k \ne m} A_{0,k} (k+1) H_{{\rm SP},{\rm FO},k} + \sum _{k>0} B_{0,k} H_{{\rm SP},{\rm circ},k} (k).\nonumber\\&& \end{eqnarray} \tag{ A9 }$$

Using the derivatives of f given in Equation (A9) in Equation (A7) then results in an equation for the A_{0, k} and B_{0, k}:

$$\begin{eqnarray} &&\sum_{k>0,k \ne m} (A_{0,k} (-k n_{0,1}+ (k+1) n_{0,2}) +1) H_{{\rm SP},{\rm FO},k}\nonumber\\ &&\qquad+ \sum _{k>0} (B_{0,k}(-k n_{0,1} +k n_{0,2}) +1) H_{{\rm SP},{\rm circ},k} =0 \nonumber \\ &&A_{0,k} = \frac{-1}{(k+1) n_{0,2}-k n_{0,1}} = \frac{1}{(k+1)n_{0,2}(X_k-1)} \nonumber \\ &&B_{0,k} = \frac{-1}{k(n_{0,2}-n_{0,1})} = \frac{1}{k n_{0,2}(P_2/P_1-1)} \nonumber \\ &&X_k = \frac{k P_2}{(k+1) P_1}, \end{eqnarray} \tag{ A10 }$$

where the first equation must hold for each k individually.

To find the transformation between the variables, we use Equation (A4). First, let us determine how to map initial period ratios used in the numerical integrations to "averaged" period ratios. The ${real}$ period ratio of the planets is given by P₂/P₁|_R = (Λ₂/Λ₁)³*(m₁/m₂)³, while the averaged period ratio is given by $P_2/P_1 \left \vert _{A}= (\tilde{\Lambda }_2/\tilde{\Lambda }_1)^{3}*(m_1/m_2)^{3}\right.$. Therefore,

$$\begin{eqnarray} && P_2/P_1 \left\vert_{R} \approx P_2/P_1 \right\vert_{A}\left[1+3 \epsilon_1\frac{(\partial {f}/{\partial \lambda _2})} {\Lambda _2} \right ] \left [1+3 \epsilon _1\frac{(\partial {f}/{\partial \lambda _1})}{\Lambda _1} \right ]^{-1} \nonumber \\ &&\approx P_2/P_1 \left\vert _{A} \left[1-3 \epsilon _1\frac{(\partial {f}/{\partial \lambda _1})}{\Lambda _1} + 3 \epsilon _1 \frac{(\partial {f}/{\partial \lambda _2})}{\Lambda _2}\right]\right. . \end{eqnarray} \tag{ A11 }$$

The required derivatives of f are

$$\begin{eqnarray} \frac{\partial {f}}{\partial \lambda _1} /\Lambda _1 & =& \sum _{k>0,k \ne m} \frac{k }{k+1} \frac{\Lambda _2}{\Lambda _1} \frac{1}{X_k-1}\nonumber\\ &&\times [f_{k,27} e_1 \cos {(\theta _k-\varpi _1)}+f_{k,31} e_2 \cos {(\theta _k-\varpi _2)}]\nonumber\\ &&+ \sum _{k>0} \frac{f_{k,0}}{P_2/P_1-1} \frac{\Lambda _2}{\Lambda _1}\cos {\phi _k} \nonumber \\ \frac{\partial {f}}{\partial \lambda _2} /\Lambda _2 & =& \sum _{k>0,k \ne m} - \frac{1}{X_k-1}\nonumber\\ &&\times [ f_{k,27} e_1 \cos {(\theta _k-\varpi _1)}+ f_{k,31} e_2 \cos {(\theta _k-\varpi _2)}]\nonumber\\ &&+ \sum _{k>0} -\frac{f_{k,0}}{P_2/P_1-1} \cos {\phi _k}, \end{eqnarray} \tag{ A12 }$$

so that

$$\begin{eqnarray} && P_2/P_1 \left \vert _{R} = P_2/P_1 \right \vert _{A} \nonumber \\ && \times\left \lbrace 1-3\epsilon _1 \sum _{k>0,k \ne m} \left [\frac{k }{k+1} \frac{\Lambda _2}{\Lambda _1}+1\right ]\frac{1}{X_k-1}\right. \nonumber\\ && \left.\times [ f_{k,27} e_1 \cos {(\theta _k-\varpi _1)}+ f_{k,31} e_2 \cos {(\theta _k-\varpi _2)}] \right.\nonumber \\ &&\left. - 3\epsilon _1 \sum _{k>0} \frac{f_{k,0}}{P_2/P_1-1} \cos {\phi _k}(1+\Lambda _2/\Lambda _1) \right \rbrace, \end{eqnarray} \tag{ A13 }$$

where Λ₂/Λ₁ = m₂/m₁(P₂/P₁)^1/3. 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_{0, k} and B_{0, k} are small. Though the functions proportional to A_{0, k} also are linear in the eccentricities, the functions f_{k + 1, 27} and f_{k + 1, 31} 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_{0, 2} − qn_{0, 1}. This is why the two indirect period terms in the combinations λ₂ − 2λ₁ and λ₂ are negligible. They appear at first order in eccentricities, but with divisors of n_{0, 2} − 2n_{0, 1} and n_{0, 2}, 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

$$\begin{eqnarray} && P_2/P_1 \left \vert _{R} = P_2/P_1 \right \vert _{A}\times \left \lbrace 1-3\epsilon _1 \sum _{k>0,k \ne m} \left [\frac{k }{k+1} \frac{\Lambda _2}{\Lambda _1}+1\right]\right. \nonumber \\ &&\times\frac{1}{X_k-1} [ f_{k,27} e_1(0) \cos {k M_1}]\nonumber\\ &&\left. - 3 \epsilon_1 \sum _{k>0} \frac{f_{k,0}}{P_2/P_1-1} \cos {k M_1}(1+\Lambda _2/\Lambda _1) \right \rbrace. \end{eqnarray} \tag{ A14 }$$

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 cos kM₁ does since f_{k, 0} > 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 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,

$$\begin{eqnarray} \frac{\partial f}{\partial \tilde{x}_i} &\propto \sin {((k+1)\lambda _2-k\lambda _1)} = 0 \nonumber \\ \frac{\partial f}{\partial y_i} &\propto \cos {((k+1)\lambda _2-k\lambda _1)} = 1. \end{eqnarray} \tag{ A15 }$$

So to first order in ₁, $y_i = \tilde{y}_i$, and

$$\begin{eqnarray} x_1 &= \tilde{x}_1 +\epsilon _1 \frac{\Lambda _2}{\sqrt{\Lambda _1}} \sum _{k>0,k \ne m} \frac{\vert f_{k,27} \vert }{(X_k-1)(k+1)} \nonumber \\ x_2 &= \tilde{x}_2 -\epsilon _1 \sqrt{\Lambda _2} \sum _{k>0,k \ne m} \frac{f_{k,31} }{(X_k-1)(k+1)}. \end{eqnarray} \tag{ A16 }$$

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 (A16) 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_i = \tilde{y}_i$, this implies that both eccentricities and longitudes of periastron are not affected greatly by the transformation. Moreover, this should be true regardless of whether e₂ is nonzero or whether the angles are nonzero because the weighting factor 1/(k + 1) will always appear.

Finally, we check the correction to the mean longitudes λ_i. We do this only for our nominal case, with e₂ = 0 and all other angles equal to zero. In this case, y_i = 0 and all the ϕ_j and θ_k are equal to zero, so all derivatives with respect to $\tilde{\Lambda }_i$ are zero. In other words, $\lambda _i = \tilde{\lambda }_i$ in this case.

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 17 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.

Zoom In Zoom Out Reset image size

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.

Second, 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 decrease. 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.

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 ϖ₁ = π/2. In this case, Δϖ = π/2 and so $\sigma = \sqrt{e_2^2+e_2^2/R^2}$.

For these choices of angles, one can show that the generalized longitude of pericenter $\psi = \arctan {\lbrace s_1/r_1\rbrace } = \arctan {\lbrace R e_1/e_2\rbrace }$. Given the values of e₁ and e₂, ψ ranges from 0° to ∼79°. The value of the resonant angle ϕ, then, is ϕ = −mϖ₁ + ψ = −mπ/2 + ψ. 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 _p). 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 18, but the others do; they are the only regular first-order resonant regions that appear. The predicted separatrices match up well with the boundaries of these regular regions.

Zoom In Zoom Out Reset image size