text:
# A generalization of the Lagrangian points: studies of resonance for highly eccentric orbits

## 1 Introduction

The three-body problem, or the dynamics of three masses due to their mutual gravitational influences, has a number of well-known special cases. One of these, the circular planar restricted case, requires that the primary and secondary bodies, $`m_1`$ and $`m_2`$, follow circular orbits about their common center of mass and that the third body be a massless test particle moving in the massive bodies’ orbit plane. These conditions simplify the three-body problem enough to produce an integral of the motion: the Jacobi constant $`C_J=2(E\mathrm{\Omega }h)`$ where $`E`$ is the particle’s energy111We refer to the test particle as ‘the particle’ and to its energy per unit mass and angular momentum per unit mass as its ‘energy’ and ‘angular momentum’., $`h`$ is the particle’s angular momentum, and $`\mathrm{\Omega }`$ is the massive bodies’ constant angular velocity.

Still, the circular planar restricted case has important applications to the dynamics of our solar system. Many of the orbits of major planets about the sun are nearly circular and are roughly confined to a plane; the same goes for many of the orbits of large moons about their planets. Common examples of applications for the circular planar restricted case include the effects of Jupiter on the asteroid belt; of Neptune on the Kuiper belt; of moons on planetary rings; and of giant planets on comets.

This paper describes a study of this problem in the regime where $`m_2m_1`$, the test particle’s eccentricity is large, and the Jacobi constant is greater than but close to $`3`$ in the standard system of units where $`G=1`$, the primary-secondary separation is 1, $`1=m_1+m_2m_1`$, and, therefore, $`\mathrm{\Omega }=1`$. Since values of $`C_J`$ near 3 correspond to test particles on circular orbits close to the secondary, this special regime includes particles originally in circular orbits around a star close enough to a planet for the planet to perturb them into very eccentric orbits. Our interest in this regime arises from an intent to investigate the paths through which small particles are perturbed by a planet until they escape or are captured. This problem was studied by Ford et al. (2001) and Rasio & Ford (1996) via numerical simulations of three massive bodies in three dimensions. Due to this motivation we use ‘star’ and ‘planet’ to refer to the primary and secondary in the remainder of this paper.

In §2 we derive to first order in $`\mu =m_2/(m_1+m_2)=m_2`$ the energy kick received by a particle in a highly eccentric orbit with semimajor axis $`a1`$ at each periapse passage. We show that since the interaction is localized at periapse, this energy kick is essentially independent of $`a`$ and depends only on the periapse distance and the azimuth difference between the planet and particle at periapse. In §3, §4, and §5 we use these energy kicks to find ‘fixed’ particle orbits and describe motion near them. These ‘fixed’ orbits are located at planet-particle mean-motion resonances. When observed stroboscopically at periapse only, they appear as fixed points just like the well-known Lagrangian points. We use both a continuous approximation and a discrete mapping in a derivation of the particle’s motion around resonances, the resonance widths, and the libration periods. When these librations are observed stroboscopically, they likewise become analogies of the well-known tadpole and horseshoe orbits. In §6 we discuss types of chaos for large-eccentricity orbits, and in §7 we summarize and discuss our findings.

## 2  Energy kick to first order in $`𝝁`$

Let $`\mathrm{\Delta }E`$ the changes in the particle’s energy between its consecutive apoapse passages. In our units, where the angular velocity of the star-planet is set to unity, the change in angular momentum is also $`\mathrm{\Delta }E`$ 222Since the angular momentum is always perpendicular to the orbit plane, only one of its components is nonzero. We therefore treat the torque and the change in angular momentum as scalars equal to the components of the vector torque and the vector change in angular momentum which are perpendicular to the orbit plane.. Therefore, it can be calculated by integrating the torque exerted on the particle:

$$\mathrm{\Delta }E=\frac{V}{f}|_rdt$$
(1)
where $`V`$ is the gravitational potential produces by the planet and star and $`f`$ is the particle’s azimuth in inertial space.

To begin with, we estimate the energy kicks to first order in $`\mu `$. We express $`\mathrm{\Delta }E`$ as $`\mathrm{\Delta }E=\mu \mathrm{\Delta }E^1+𝒪(\mu ^2)`$. To evaluate $`\mathrm{\Delta }E^1`$ we calculate the torque assuming that the particle moves on a Keplerian trajectory around the star, with its focus at the center of mass. The effect of the deviation of the trajectory from that description on $`\mathrm{\Delta }E`$ is of order $`\mu ^2`$ or higher.

Since we are considering only the time elapsed between two consecutive apoapses, we choose coordinates such that the time $`t=0`$ when the particle is at periapse and the direction of periapse is along the positive $`x`$-axis. The planet and star are in uniform circular motion, so we can write $`V=V(\theta ,r)`$ where $`\theta `$ is the angle between the planet and the particle and $`r`$ is the particle’s distance from the origin. This gives

$`\mathrm{\Delta }E={\displaystyle \frac{V}{\theta }}|_rdt.`$ (2)
$`V`$ is given explicitly by

$$V=V_{\mathrm{planet}}+V_{\mathrm{star}}=\frac{\mu }{|\stackrel{}{r}\stackrel{}{r}_{\mathrm{planet}}|}+\frac{1}{|\stackrel{}{r}\stackrel{}{r}_{\mathrm{star}}|};$$
(3)
to first order in $`\mu `$, this gives

$$V=\frac{1}{r}+\mu \left(\frac{1}{(r^2+12r\mathrm{cos}\theta )^{1/2}}\frac{\mathrm{cos}\theta }{r^2}\right).$$
(4)
Let $`\varphi `$ be the angle between the planet and the particle at periapse333Thus defined, $`\varphi `$ is the usual resonant argument measured at periapse only., so that $`\theta =\varphi +tf`$. Then the derivative with respect to $`\theta `$ at fixed $`r`$ which appears in Eq. (2) can be replaced by a derivative with respect to $`\varphi `$. To first order in $`\mu `$ the particle trajectory $`r(t)`$ can be assumed fixed and independent of $`\varphi `$, so we can move the $`\varphi `$ derivative outside the integral of Eq. (2). Using the first order expression for $`V`$ we get

$$\mathrm{\Delta }E^1=\frac{dU^1}{d\varphi }$$
(5)
where the effective potential $`U^1`$ is given by

$$U^1=\left[\frac{1}{\sqrt{r^2+12r\mathrm{cos}\theta }}\frac{1}{\sqrt{r^2+1+2r\mathrm{cos}(tf)}}+\frac{\mathrm{cos}\theta \mathrm{cos}(tf)}{r^2}\right]𝑑t.$$
(6)
The integral is performed over one Keplerian orbit of the particle.

In this expression for $`U^1`$, the first term in the brackets is the ‘direct’ term; it represents the planet’s contribution. The second term does not contribute to $`\mathrm{\Delta }E^1`$; it keeps $`U^1`$ from diverging when $`a\mathrm{}`$ and is obtained from the first term by substituting $`\varphi =\pi `$. The third term is the ‘indirect’ term; it represents interactions with the star. $`\mathrm{\Delta }E^1`$ and its effective potential $`U^1`$ are functions of $`\mu `$, $`\varphi `$, and the particle trajectory shape, which determines $`r`$ and $`f`$ as a function of $`t`$. Note that up to a constant, the effective potential $`U^1`$ is simply the time-integrated potential over the trajectory of the particle.

When the apoapse distance $`a(1+e)`$ is much larger than both 1 and the periapse distance $`r_p=a(1e)`$, the perturbing effects of the star and planet on the particle near periapse dominate over perturbing effects on the particle elsewhere in its orbit. In this regime, the entire energy kick $`\mathrm{\Delta }E`$ occurring between consecutive apoapse passages can be thought of as a discrete event associated with a particular periapse passage. In the limit as $`a`$ diverges due to energy kicks but $`C_J`$ remains constant, $`e1`$, $`r_p`$ approaches a constant, and except near apoapse the entire trajectory approaches a parabola independent of $`a`$: $`a\mathrm{}`$. If the particle is outside the planet’s Hill sphere,

$$\underset{a\mathrm{}}{lim}r_p=\underset{a\mathrm{}}{lim}\left[\frac{1}{2}\left(C_J+\frac{1}{a}\right)\frac{1}{\sqrt{1+e}}\right]^2=\frac{C_J^{\mathrm{\hspace{0.33em}2}}}{8}$$
(7)
$$r=2r_p/(1+\mathrm{cos}f)$$
(8)
$$\frac{df}{dt}=\frac{(1+\mathrm{cos}f)^2}{(2r_p)^{3/2}}$$
(9)
$$t=(2r_p)^{3/2}\frac{1}{6}\mathrm{tan}\frac{f}{2}\left(3+\mathrm{tan}^2\frac{f}{2}\right).$$
(10)
For particles that start close to the orbit of the planet, the periapse distance is therefore $`r_p=9/8`$.

Given this asymptotic form for the orbit, we can calculate the asymptotic forms of $`U^1(a,e,\varphi )U^1(r_p,\varphi )`$ and $`\mathrm{\Delta }E^1(a,e,\varphi )\mathrm{\Delta }E^1(r_p,\varphi )`$ in the large-$`a`$ limit. For $`C_J=3`$, the computed values of $`U^1`$ and its derivative $`\mathrm{\Delta }E^1`$ as a function of $`\varphi `$ are shown in Figures 1 and 2. Near $`\varphi =0`$, $`\mathrm{\Delta }E^1`$ is dominated by the direct particle-planet interaction because the minimum planet-particle distance is much less than the star-particle distance. When $`\varphi =0`$, $`\mathrm{\Delta }E^1=0`$ because of symmetry. When $`\varphi <0`$ but $`|\varphi |1`$, the planet lags the particle for most of the time the particle spends near periapse, so $`\mathrm{\Delta }E^1<0`$. Similarly, when $`\varphi >0`$ and $`|\varphi |1`$, $`\mathrm{\Delta }E^1>0`$.

When $`||\varphi |\pi |1`$, the indirect contribution $`\mathrm{\Delta }E_{\mathrm{ind}}^1`$ due to the star’s reflex motion dominates because the star passes closer to the particle. From the $`\varphi `$-dependent part of star’s contribution to the integral in Eq. (6), $`\mathrm{\Delta }E_{\mathrm{ind}}^1`$ is a sinusoidal function of $`\varphi `$:

$$\frac{}{\varphi }_{\mathrm{}}^{\mathrm{}}\frac{\mathrm{cos}\theta }{r^2}𝑑t=\frac{\mathrm{sin}\varphi }{\sqrt{2r_p}}_\pi ^\pi \mathrm{cos}(tf)𝑑f2.0\mathrm{sin}\varphi $$
(11)
where in the last step we use $`r_p=9/8`$ as an example in evaluating the coefficient.

The integral in Eq. (11) seems to suggest that star-particle interactions over all intervals in $`f`$ should contribute significantly to $`\mathrm{\Delta }E_{\mathrm{ind}}^1`$. However, as Eq. (10) shows, $`|t|`$ increases much faster than $`|f|`$ as $`|f|`$ approaches $`\pi `$. As a result, oscillations in $`\mathrm{cos}(tf)`$ kill contributions to the integral at $`|f|`$ near $`\pi `$, and the star-particle interaction is important only near periapse.

We can also get the total contribution of terms that are second order or above in $`\mu `$: this is just the difference between values of $`\mathrm{\Delta }E`$ found by numerical integration of the equations of motion and values of $`\mu \mathrm{\Delta }E^1`$ given by Eq. (6) (see Figure 1).

## 3 First-order resonances

Resonances occur when the particle completes $`p`$ orbits in exactly the time needed for the planet to complete $`p+q`$ orbits for some integers $`p`$, $`q`$. This situation is known as a $`p:p+q`$ resonance. In the standard treatment of these resonances, both orbits in question are usually nearly circular and a significant interaction occurs every time the bodies are at conjunction, i.e. whenever their azimuths coincide. This happens once every resonant cycle if $`q=1`$, so $`q=1`$ resonances are usually termed ‘first order’ resonances. During a conjunction between a test particle and planet in orbits with low eccentricity $`e1`$, the torque exerted on the particle while the particle precedes the planet almost cancels the torque exerted while the particle lags the planet; the residual is of order $`e`$. When $`q>1`$, $`q`$ conjunctions occur during each resonant cycle. Because they occur in different positions in inertial space, their effects tend to cancel each other, leaving a residual torque of order $`e^q`$. Since the interaction strength decreases exponentially with increasing $`q`$ as $`e^q`$, resonances in the standard treatment are usually classified by $`q`$ value. Accordingly, a $`p:p+q`$ resonance is called a ‘$`q`$th-order’ resonance regardless of the value of $`p`$.

However, the high eccentricities of orbits in the large-$`a`$ regime discussed here make the standard definition of resonance order meaningless. Since $`e1`$, resonances of different order under the standard definition have comparable significance because $`e^q1`$. Also, encounters at periapse are physically far more important than conjunctions at other points in the particle’s orbit. We therefore redefine the ‘orders of resonance’ to focus on interactions at periapse. If the planet completes an integer number of orbits in the time it takes the particle to orbit exactly once, then we say the particle is in a ‘first-order’ resonance. In general, we say the particle is in an ‘$`p`$th-order’ resonance if the planet completes an integer number of orbits in the time it takes the particle to orbit $`p`$ times: then there are $`p`$ interactions within one resonant cycle. In terms of the standard resonance treatment, we say a $`p:p+q`$ resonance in the large-$`a`$ regime is ‘$`p`$th-order’ regardless of the value of $`q`$. In both the large- and small-eccentricity cases, the order of the resonance is given by the number of significant interactions within a single resonant cycle.

In the following we show that this revised definition does indeed make sense. We calculate the widths of resonances of various orders in the large-$`a`$ limit and show that with this new definition, their widths decay exponentially with the order of the resonance. We discuss in detail the first-order or $`1:N`$ resonances.

According to this definition of resonance orders, $`\varphi `$ should be constant in time if we consider a particle exactly at a first-order resonance of semimajor axis $`a_{\mathrm{res}}`$ and if we ignore the effects of energy kicks. A particle close to resonance with, say, semimajor axis $`a=a_{\mathrm{res}}+\mathrm{\Delta }a`$ should drift in $`\varphi `$ over time at a constant rate, again ignoring energy kicks. The amount of drift per orbit of the test particle is just the difference between its orbital period $`2\pi a^{3/2}`$ and the resonant one $`2\pi a_{\mathrm{res}}^{3/2}`$. We can express this drift as

$$\frac{d\varphi }{dn}=3\pi a^{1/2}\mathrm{\Delta }a.$$
(12)
Here the $`dn`$ acts as a place-holder indicating that the $`d\varphi `$ is associated with a single particle orbit. The differential is justified assuming $`\mathrm{\Delta }aa^{1/2}`$ so that the change in $`\varphi `$ per particle orbit is much less than $`\pi `$. We refer to this differential form as the continuous approximation.

Energy kicks cause the semimajor axis to evolve in time. To first order in $`\mu `$ we have

$$\frac{d(\mathrm{\Delta }a)}{dn}=2a^2\frac{dE}{dn}=2a^2\mu \frac{dU^1}{d\varphi }.$$
(13)
To justify the differentials here we require that $`\mu `$ be small enough for the change in $`\mathrm{\Delta }a`$ due to a single kick to be much less than the typical $`\mathrm{\Delta }a`$. We differentiate Eq. (12) and substitute Eq. (13) to get

$$\frac{d^2\varphi }{dn^2}=6\pi a^{5/2}\mu \frac{dU^1}{d\varphi }.$$
(14)
This shows that $`\varphi `$ simply evolves as a particle moving in the potential $`U^1(\varphi )`$.

### 3.1 Generalized Lagrangian points

Since $`U^1`$ has four extrema at the four zeroes of $`\mathrm{\Delta }E^1`$, there are four fixed points. These fixed points are simply periodic orbits of the particle. Two of the fixed points are unstable since they correspond to maxima of the potential $`\theta =0,\pi `$. The other two are stable since they correspond to potential minima at $`\theta =\pm 1.21`$. The existence of two extrema at $`\theta =0,\pi `$ is guaranteed by symmetry arguments. The two additional extrema at $`\varphi =\pm 1.21`$ occur where the energy kicks from the planet and star cancel each other exactly. These extrema therefore appear only when the indirect term—or, equivalently, the star’s motion—is taken into account.

This discussion suggests an analogy between the five well known Lagrangian points and the new fixed points. The two stable points correspond to the stable Lagrangian points $`L_4`$ and $`L_5`$, which also appear only when the motion of the star, i.e. the indirect term, is taken into account. The unstable fixed point at $`\varphi =\pi `$ is the analogue of $`L_3`$; the one at $`\varphi =0`$ corresponds to $`L_1`$ and $`L_2`$, which merge in this generalization. For a given resonance $`1:N`$, $`N=a^{3/2}`$, we therefore denote the fixed points by $`L_{12}^N`$, $`L_3^N`$, $`L_4^N`$, and $`L_5^N`$. The positions of these new fixed points in comparison to their standard Lagrangian counterparts is summarized in Table 1.

### 3.2 Generalized tadpoles

The analogy is more obvious when motion around the fixed points is investigated. Small amplitude motion around the stable fixed point $`L_4^N`$ and $`L_5^N`$ can be approximated by expanding $`U^1`$ around its minimum. This results in a harmonic oscillator equation:

$$\frac{d^2\varphi }{dn^2}=6\pi a^{5/2}\mu \left(\frac{d^2U^1}{d\varphi ^2}\right)|_{\varphi =\varphi _{4,5}^N}(\varphi \varphi _{\mathrm{res}}).$$
(15)
The small amplitude libration period around either $`L_4^N`$ or $`L_5^N`$ is therefore given by

$$K=\frac{T_{\mathrm{libration}}}{2\pi a^{3/2}}=\left(\frac{3}{2\pi }\frac{d^2U^1}{d\varphi ^2}|_{\varphi =\varphi _{4,5}^N}\right)^{1/2}a^{5/4}\mu ^{1/2}=0.79a^{5/4}\mu ^{1/2}$$
(16)
where in the last step we use $`r_p=9/8`$ in the large-$`a`$ limit to get $`d^2U^1/d\varphi ^23.3`$ at $`\varphi =\varphi _{4,5}`$. Note that $`K`$ gives the number of periapse crossing per libration period. In our units, where $`2\pi `$ is the period of the massive bodies, the libration period is then $`2\pi a^{3/2}K`$.

Since Eq. (14) describes motion under the influence of a fixed potential, we can write down the conservation of energy equation by multiplying Eq. (14) by $`\frac{d\varphi }{dn}`$ and integrating with respect to $`n`$:

$$\frac{1}{2}\left(\frac{d\varphi }{dn}\right)^2+6\pi a^{5/2}\mu U^1=\mathrm{constant}.$$
(17)
The constant of integration is the ‘energy’ associated with the movement of the orbit in $`\varphi `$ and $`a`$. Since the potential is finite, it can only support a finite particle ‘speed’ in libration around $`L_4^N`$ or $`L_5^N`$. The ‘speed’ is directly related to the deviation of the semimajor axis from the resonance via Eq. (12), so the maximal width of these librations in $`a`$ is given by

$$\mathrm{\Delta }a_{\mathrm{max}}=\left(\frac{4}{3\pi }\right)^{1/2}\mu ^{1/2}a^{3/4}[U^1(\pi )U^1(\varphi _4)]^{1/2}0.78a^{3/4}\mu ^{1/2}.$$
(18)
These librations around the fixed points $`L_4^N`$ or $`L_5^N`$ are analogues of the well-known tadpole orbits. Note that the maximal widths of both the standard and generalized tadpole orbits scale as $`\mu ^{1/2}`$ (see Table 1). The similarity is more apparent if we treat the $`(a,\varphi )`$ parameters, which describe the orbit of the particle, as polar coordinates as shown in Figure 3. Seen in this way, $`(a,\varphi )`$ are analogous but not identical the polar coordinates of the particle in the rotating frame: $`a`$ is the semimajor axis, not the radius, and $`\varphi `$ is the azimuth of the test particle in the rotating frame only at periapse passage. Then the fundamental difference between the $`(a,\varphi )`$ plane and the rotating frame is that while generalized Lagrangian points and the motion around them exist in a surface of section made up of discrete points representing periapse passages, the standard rotating frame with the standard Lagrangian points is made up of continuous trajectories. Therefore, while the standard Lagrangian points are fixed points in the rotating frame, the generalized points represent periodic orbits in that frame.

The generalized tadpoles are equivalent to ‘asymmetric librations’—trajectories whose resonant argument librates about a value other than 0 or $`\pi `$. In this context, $`L_4^N`$ and $`L_5^N`$ correspond to ‘asymmetric periodic orbits’ whose resonant argument is constant but not equal to 0 or $`\pi `$. Our discussion above gives a simple physical argument for the existence of asymmetric librations in all stable $`1:N`$ exterior resonances. Again, note that the existence of these asymmetric librations and asymmetric periodic orbits follows from analysis of $`U^1`$ only when both the direct and indirect terms are accounted for.

### 3.3 Generalized horseshoes

As the energy of the particle moving under the $`U^1`$ potential increases beyond that of the maximal tadpole orbit, it overcomes the lower potential barrier at $`\varphi =\pi `$. As long as its energy is still below the higher barrier at $`\varphi =0`$, the particle will librate around both the $`L_4^N`$ and $`L_5^N`$ points, avoiding only a narrow range in $`\varphi `$ around $`\varphi =0`$. These trajectories are the generalized horseshoe orbits. Using the same method as we used for the tadpoles, we calculate their widths in the continuous approximation to be

$$\mathrm{\Delta }a_{\mathrm{max}}=\left(\frac{4}{3\pi }\right)^{1/2}\mu ^{1/2}a^{3/4}[U^1(0)U^1(\varphi _4)]^{1/2}=1.8a^{3/4}\mu ^{1/2}.$$
(19)
The width of the maximal standard horseshoe does not follow this $`\mu ^{1/2}`$ pattern, since the standard horseshoe case differs qualitatively from its generalized version. For the standard horseshoe, the angular momentum change is concentrated near the horseshoe’s two ends. The close approach of the particle to the planet there increases the strength of the interaction beyond $`\mu `$. As a result, the width of the horseshoe scales as $`\mu ^{1/3}`$ rather than $`\mu ^{1/2}`$. For a generalized horseshoe, the librating particle never gets closer to the planet than $`r_p1`$.

In Figures 4 and 5 we show libration around $`L_{4,5}^4`$ and $`L_{4,5}^{10}`$. In order to focus on the motion close to these points we plot $`a`$ and $`\varphi `$ as Cartesian rather than polar coordinates. In these plots, the librations in the surfaces of section appear to be ‘warped’ when compared to the continuous approximations calculated using the pendulum-like Eq. (14). This ‘warping’ is due to the discrete nature of the motion in the surfaces of section. As a trajectory moves from $`\mathrm{\Delta }a=0`$ toward larger positive $`\mathrm{\Delta }a`$ values, for example, the energy kicks stay positive and $`\mathrm{\Delta }a`$ should keep increasing until the trajectory reaches a $`\varphi `$ value corresponding to a zero in the $`\mathrm{\Delta }E^1`$ vs. $`\varphi `$ curve. Within the continuous approximation, we expect the trajectory to begin moving back toward $`\mathrm{\Delta }a=0`$ at exactly this $`\varphi `$ because $`\mathrm{\Delta }E^1`$ changes sign. A discrete trajectory will ‘overshoot’ the nominal $`\varphi `$ where $`\mathrm{\Delta }E^1=0`$ since a positive energy kick will carry the trajectory past this $`\varphi `$ before the first negative kick is applied. As a result, the libration trajectories in the surfaces of section tend to become warped in the direction in which orbits move when librating. A quantitative discussion of this feature is given in the next section.

## 4 The eccentric mapping

The ‘warping’ noted above suggests that the discrete nature of the surface of section is essential to understanding some feature of the motion in the $`(a,\varphi )`$ plane. To study this, we define a mapping from the $`(a,\varphi )`$ plane to itself. Beginning at an arbitrary point, this mapping produces an infinite series $`(a^i,\varphi ^i)`$ of points visited by the test particle in the $`(a,\varphi )`$ plane. Except perhaps in the few lowest-$`N`$ resonances, we can build an excellent approximation to the correct mapping by applying the first-order kicks in the large-$`a`$ limit:

$$\frac{1}{2a^{i+1}}=\frac{1}{2a^i}+\mu \mathrm{\Delta }E^1(\varphi ^i)$$
(20)
$$\varphi ^{i+1}=\varphi ^i+2\pi (a^{i+1})^{3/2}$$
(21)
where the new value of $`\varphi `$ is calculated modulo $`2\pi `$, i.e. brought back into the interval $`(\pi ,\pi )`$ by adding an integer multiple of $`2\pi `$. Note that the $`a`$-value used to find $`\varphi ^{i+1}`$ itself has index $`i+1`$; physically, this corresponds to the large-$`a`$ limit assumption that each energy kick is a discrete event associated with a given periapse passage. Applying this mapping for several initial values in the $`(a,\varphi )`$ plane results in the right panels of Figures 4 and 5. The close resemblance between trajectories generated with the mapping and with numerical orbit integrations demonstrate this mapping’s accuracy.

It turns out that the warping of the small amplitude tadpoles can be understood completely in terms of the mapping. Close to the fixed points $`L_4^N`$ and $`L_5^N`$, we define $`\mathrm{\Delta }a^i=a^ia_{res}`$ and $`\mathrm{\Delta }\varphi ^i=\varphi ^i\theta _{res}`$ so that the mapping becomes

$$\mathrm{\Delta }a^{i+1}=\mathrm{\Delta }a^i2a_{\mathrm{res}}^2\mu \mathrm{\Delta }\varphi ^i\frac{d^2U^1}{d\varphi ^2}|_{\varphi =\varphi _{4,5}^N}$$
(22)
$$\mathrm{\Delta }\varphi ^{i+1}=\mathrm{\Delta }\varphi ^i+3\pi a_{\mathrm{res}}^{1/2}\mathrm{\Delta }a^{i+1}.$$
(23)
Since these are linear recursive equations, they can be solved analytically by standard techniques. We seek a solution of the form $`(\mathrm{\Delta }a^i,\mathrm{\Delta }\varphi ^i)=(A,\mathrm{\Phi })\alpha ^i`$. (The reader should not confuse indices with exponents—superscripts above Greek letters other than $`\varphi `$ are exponents.) Substituting in the recursive equations, and seeking a non-trivial solution we obtain

$$(\alpha 1)^2+6\pi \mu a_{\mathrm{res}}^{5/2}\frac{d^2U^1}{d\varphi ^2}|_{\varphi =\varphi _{4,5}^N}\alpha =0.$$
(24)
Note that the dimensionless parameter in this equation is simply $`(2\pi /K)^2`$ where $`K`$ is the number of periapse passages per libration in the continuous approximation as given by Eq. (16). If we denote the solutions as $`\alpha _1`$ and $`\alpha _2`$, it is clear from the above equation that their product $`\alpha _1\alpha _2`$ equals $`1`$. Since we are interested in potential minima, $`K^2>0`$.

For $`K\pi `$, the two roots are complex conjugates and each has unity norm. The fixed point is therefore an elliptical point in the discrete mapping as well as in the continuous approximation. The two values of $`\alpha `$ are given by

$$\alpha _{1,2}=12(\pi /K)^2\pm 2\sqrt{1(\pi /K)^2}(\pi /K)i.$$
(25)
The number of periapse passages per libration is given by

$$K^{\mathrm{map}}=\frac{2\pi }{\mathrm{arg}(\alpha )}=2\pi \left[\mathrm{arctan}\left(\frac{2\sqrt{1(\pi /K)^2}(\pi /K)}{12(\pi /K)^2}\right)\right]^1.$$
(26)
Clearly, as $`K\mathrm{}`$, $`K^{\mathrm{map}}/K1`$. This is expected since the continuous approximation is justified in this limit. Using the two values of $`\alpha `$, we can find the eigenvectors:

$$(A,\mathrm{\Phi })=(\mu a_{\mathrm{res}}^2\frac{d^2U^1}{d\varphi ^2},(\pi /K)^2\pm \sqrt{1(\pi /K)^2}(\pi /K)i).$$
(27)
Since the eigenvectors determine the axes of the ellipses representing small librations about the fixed points, the similar shapes and orientations of the smallest librations in the middle and right-hand panels in each of Figures 4 and 5 confirm the eccentric mapping’s accuracy.

For $`a=10^{2/3}`$, $`\mu =10^4`$, the continuous approximation gives $`11.6`$ orbits per tadpole libration (Eq. 16). The eccentric mapping gives $`11.4`$ orbits (Eq. (26)). This is close to the the 10.7 orbits per libration observed for very small librations about the fixed points444The largest tadpole libration shown in Figure 5 breaks into 14 islands. This is an example of the Poincaré-Birkhoff fixed-point theorem. It indicates that this tadpole’s libration period is 14 orbits. This lengthening of the period is expected as the trajectory grows toward the separatrix passing through $`\varphi =\pi `$.. The negative power of $`a`$ in Eq. 16 implies that as $`a`$ increases, the number of periapse passages per tadpole libration period will decrease and the trajectory shapes will become increasingly warped.

In fact, when $`a`$ grows so large that $`K`$ falls below $`\pi `$, the tadpoles are destroyed. For $`K<\pi `$, the roots of Eq. (24) are real and distinct; therefore one of them is larger than unity. Then the fixed point is not stable despite being at a potential minimum. Our quantity $`K`$ is closely related to the residue $`R`$ discussed by Greene (1979): $`R=1(\pi /K)^2`$.

The warping of the tadpoles which, at its extreme, leads to destruction of the resonances is absent in the continuous approximation. However, it can be understood as perturbations from nearby resonances. Interactions between neighboring first-order resonances should become large enough to destroy these resonances when the resonances begin to overlap. Eq. (19) implies that as $`a`$ increases, the resonances widen in $`a`$ while the distance between them decreases. Then we can find a condition on $`\mu `$ and $`a`$ for resonance overlap by dividing half the distance between consecutive first order resonances by the width $`\mathrm{\Delta }a_{\mathrm{max}}`$ of each resonance as given by equation (19). In the large-$`a`$ limit, the distance between resonances is given by $`\frac{2}{3}a^{1/2}`$ so we obtain

$$\frac{\mathrm{resonance}\mathrm{separation}}{2\mathrm{\Delta }a_{\mathrm{max}}}=\left(\frac{\pi }{12}\right)^{1/2}a^{5/4}\mu ^{1/2}[U^1(0)U^1(\varphi _4)]^{1/2}=0.18a^{5/4}\mu ^{1/2}.$$
(28)
This is proportional to $`K^2`$: in the large-$`a`$ limit with $`r_p=9/8`$, the right-hand side is $`0.28K^2`$ and first-order resonances overlap when $`K<1.9`$. In this case, therefore, first-order resonances are destroyed before they formally overlap. With $`\mu =10^3`$, first-order resonances are destroyed above $`a4.02`$. Then if we neglect higher-order effects in $`\mu `$ and differences between the potentials for $`a4`$ and the large-$`a`$ limit, we expect just eight stable first-order resonances.

## 5 Higher-order resonances

As defined in §2, higher-order resonances are the $`p:p+q`$ resonances with $`p>1`$. These resonances are located at $`a_{\mathrm{res}}=(N/p)^{2/3}`$ where $`N=p+q`$ is an integer relatively prime to $`p`$. In analogy to our treatment of first order resonances, we note that if we neglect energy kicks, a particle exactly at resonance should move in $`\varphi `$ by $`2\pi `$ during each resonant cycle and by $`2\pi q/p`$ between consecutive periapse passages. The stationary points of this resonance should therefore occur at regular intervals of $`2\pi /p`$ in $`\varphi `$.

To study motion near but not at resonance, we include energy kicks. For a particle close to resonance, we can follow its trajectory by treating each resonant cycle as $`p`$ applications of the eccentric mapping, one for each periapse passage in the cycle:

$$\left(\begin{array}{c}\mathrm{\Delta }a^{j+1}\\ \mathrm{\Delta }\varphi ^{j+1}\end{array}\right)=\underset{i=0}{\overset{p1}{}}\left(\begin{array}{cc}1& 2a_{\mathrm{res}}^2\mu \frac{d^2U^1}{d\varphi ^2}|_{\varphi =\varphi _p^N+2\pi i/p}\\ 3\pi a_{\mathrm{res}}^{1/2}& 16\pi a_{\mathrm{res}}^{5/2}\mu \frac{d^2U^1}{d\varphi ^2}|_{\varphi =\varphi _p^N+2\pi i/p}\end{array}\right)\left(\begin{array}{c}\mathrm{\Delta }a^j\\ \mathrm{\Delta }\varphi ^j\end{array}\right).$$
(29)
As before, $`\mathrm{\Delta }a^j=a^ja_{\mathrm{res}}`$ and $`\mathrm{\Delta }\varphi ^j=\varphi ^j\varphi _p^N`$ where $`\varphi _p^N`$ corresponds to the nearest fixed point in the resonance. The condition under which the linearization in $`\frac{dU^1}{d\varphi }`$ is valid is now $`\mathrm{\Delta }a^ja^{1/2}/p^2`$ instead of $`\mathrm{\Delta }a^ja^{1/2}`$ because the number of energy kicks per resonant cycle is $`p`$ instead of 1 and because the scale in $`\varphi `$ over which the potential changes is now $`\pi /p`$ instead of $`\pi `$. The condition under which linearization in $`\mu `$ is valid also changes because the largest term linear in $`\mu `$ that appears in the mapping matrix is $`a_{\mathrm{res}}^{5/2}\mu \frac{d^2U^1}{d\varphi ^2}`$. Though $`\mu `$ itself is small, cross-terms of order $`\mu ^2`$ and higher are now important unless $`a_{\mathrm{res}}^{5/2}\mu \frac{d^2U^1}{d\varphi ^2}1`$. This stronger condition is equivalent to $`K1`$, so the higher-order resonance treatment does not offer any simplifying advantages over the eccentric mapping discussed in §4 unless $`K`$ is large.

Since $`\mathrm{\Delta }\varphi `$ changes very little between consecutive periapse passages in this $`K1`$ regime, we can use a variant of the continuous approximation where we neglect the effects of drift in $`\varphi `$ within a single resonant cycle. Then we can treat the particle’s motion in terms of the net energy kick over an entire resonant cycle rather than a single particle orbit. The net energy kick is just the sum of $`p`$ energy kicks spaced $`2\pi /p`$ apart in $`\varphi `$, so the particle appears to move in the potential

$$U_p^1=\underset{k=0}{\overset{p1}{}}U^1(\varphi 2\pi k/p).$$
(30)
Note that effects of the star’s reflex motion do not contribute to $`U_p^1`$ if $`p>1`$: the indirect term in $`U^1`$ is exactly sinusoidal and the sum of $`p`$ identical sine curves spaced $`2\pi /p`$ apart in phase is 0, so $`U_{p,\mathrm{ind}}^1=0`$. Since the part of $`U^1`$ due to the planet’s direct contribution has just one maximum and one minimum at $`\varphi =0,\pi `$ respectively, $`U_p^1`$ has $`p`$ identical maxima and minima (see Figure 6. Then a trajectory librating in one of the minima of $`U_n^1`$ should appear as a series of ‘islands’ spaced evenly in $`\varphi `$ in the $`(a,\varphi )`$ plane. As a result, no asymmetric librations are possible in higher-order resonances. Our result that, among exterior resonances, only $`1:N`$ resonances show asymmetric librations is consistent with work done by Frangakis (1973). He analyzed expressions for the time-averaged direct and indirect terms of the disturbing function to find that asymmetric librations can exist only in $`p:p+q`$ resonances where $`p=\pm 1`$.

Because the $`p`$ energy kicks received by the particle during each resonant cycle are spaced evenly by $`2\pi /p`$ in $`\varphi `$, we expect that the kicks will partially cancel over each resonant cycle and that this cancellation will improve exponentially as $`p`$ increases. We therefore expect the amplitude of $`U_p^1`$ to decrease exponentially with $`p`$. As Figure 7 shows, this exponential decay is observed numerically: a best-fit line in log-log space gives amplitude $`1.20^p`$.

## 6  Chaos in the large-$`𝒂`$ limit

We discuss just a few of the types and regions of chaos which arise when $`a`$ is large. We first discuss ‘global’ chaos, which consists of chaotic regions that span a few resonance widths or more. We then give a few examples of ‘local’ chaos—chaos confined to regions within a single resonance—and compare the structure seen in trajectories in the $`(a,\varphi )`$ plane on local and global scales.

On large scales in $`a`$, chaotic regions arise where there is overlap between neighboring resonances or instability due to a small winding number as discussed in §4. In regions of the $`(a,\varphi )`$ plane where first-order and/or higher-order resonances overlap even partially, we expect to see contiguous ‘globally’ chaotic regions which span large ranges in $`a`$. Any remaining stable regions within resonances will appear as ‘islands’ of stable librations. Particles can undergo large changes in $`a`$ only if they move in these chaotic regions, so such regions provide the only channels through which initially bound particles can escape from the star-planet system.

If $`a`$ is large enough, $`K`$ falls below $`\pi `$ and these ‘islands’ disappear as discussed in §4. For a given value of $`\mu `$, we see from Eq. 16 that this occurs when

$$a>\mu ^{2/5}\left(\frac{3\pi }{2}\frac{d^2U^1}{d\varphi ^2}\right)^{2/5}=0.33\mu ^{2/5}$$
(31)
where, again, the numerical example corresponds to $`r_p=9/8`$. The condition $`K<\pi `$ for resonance destruction in higher-order resonances does not follow simply from Eq. 28: effects of order $`\mu ^2`$ or higher may make the winding number expressions for higher-order resonances differ from the first-order resonance case in Eq. 16. However, numerical experiments suggest that the $`a`$ at which higher-order resonances become unstable is comparable to but less than that given by 31.

If $`\mu `$ is small enough, there should be regions in $`a`$ where the resonances do not overlap. In these regions we expect to see stable trajectories which circulate around the resonances instead of librating in them. We find numerically that stable circulating trajectories exist for $`\mu `$ values up to at least $`\mu =5\times 10^6`$; an example is shown in Figure 8. Greene (1979) suggests that as $`\mu `$ increases, the last stable circulating trajectory should have semimajor axis $`a`$ such that $`a^{3/2}`$ is the golden ratio $`(1+\sqrt{5})/2`$. Our situation differs qualitatively from Greene’s in that our potential depends on its linear coordinate, the semimajor axis, while Greene’s potential, which is given by the standard map, is independent of its linear coordinate $`r`$. Specifically, when $`a`$ is not much larger than 1, $`r_p11`$; this leads to a larger maximum energy kick and potential well depth than is expected for $`r_p=9/8`$, so the resonances are wider and more prone to overlap for a given $`a`$ close to 1 than we would expect in the large-$`a`$ limit. However, as $`a`$ increases the resonance spacing decreases as discussed in §4. These competing effects suggest that the last stable circulating trajectories—those which, in a sense, are ‘farthest’ from any resonances—should lie neither near $`a=1`$ nor at $`a1`$. Also, effects of order $`\mu ^2`$ and higher which are present in our situation have no analogue in Greene’s analysis of the standard map. So it is unsurprising that the last stable circulating trajectories which we found numerically have $`a^{3/2}`$ unrelated to $`(1+\sqrt{5})/2`$.

Continuity and uniqueness imply that in a system with two degrees of freedom, stable trajectories in a two-dimensional surface of section cannot be crossed. In the planar restricted three-body problem, therefore, stable circulating trajectories divide the $`(a,\varphi )`$ plane into separated regions in $`a`$. This implies that for any $`\mu <5\times 10^6`$, chaotic and regular trajectories which start close enough to the planet are confined to a set range in $`a`$. Then the particles associated with these trajectories can never escape from the star-planet system.

This bounding of chaotic regions by stable trajectories also leads to confinement of chaos on very small scales in $`a`$. Regions of small-scale ‘local’ chaos arise from unstable fixed points which must be saddle points due to the area-preserving nature of the eccentric mapping; the separatrices associated with the saddle points are chaotic. If stable continuous trajectories exist near a saddle point, they act as boundaries to the chaotic separatrix. Prominent examples of these separatrices include those dividing tadpole and horseshoe analogue trajectories within individual first-order resonances. These regions are bounded by the largest stable tadpole and smallest stable horseshoe, so their maximum range in $`a`$ is at most the resonance width. One of these is shown in Figure 9.

The existence of similar separatrices on all scales smaller than a single resonance width follows from the Poincaré-Birkhoff theorem, which states that for small enough $`\mu `$, a trajectory with rational winding number $`K`$ is associated with equal numbers of alternating stable and unstable fixed points. According to the KAM theorem, some continuous trajectories—that is, trajectories with irrational $`K`$—should also be stable as long as $`\mu `$ is small enough. If a trajectory with rational $`K`$ is bounded on either side by stable continuous trajectories with irrational $`K`$, then the chaos associated with the unstable fixed points is confined to the region bounded by the continuous trajectories. As for the stable fixed points, they are associated with their own librating trajectories; the tadpole analogue made up of islands shown in Figure 5 gives an example of such librations. We expect some librations like these to have rational winding numbers and, therefore, their own sets of unstable fixed points and confined chaos on an even smaller scale. In principle, this argument can be applied repeatedly within a single resonance to unearth similar chaotic regions on scales as small as desired.

We can treat the entire $`(a,\varphi )`$ plane as an extension of this self-similarity to the largest possible scales. If we plot $`(a,\varphi )`$ as polar coordinates, a $`p:p+q`$ resonance trajectory appears to ‘wind’ around the point $`a=0`$ with rational winding number $`p/(p+q)`$. Also, the corresponding resonance is associated with $`p`$ stable and $`p`$ unstable fixed points when $`p>1`$ and $`2`$ stable and $`2`$ unstable fixed points when $`p=1`$. This provides a striking visual analogy to the librations seen within a single resonance.

## 7 Discussion and conclusions

Using simple physical reasoning instead of explicit analysis of terms in the disturbing function, we have developed a framework for studying particles perturbed into exterior high-eccentricity orbits in the circular planar restricted three-body problem. We have found that, to first order in $`\mu `$, these orbits move in $`(a,\varphi )`$ phase space according to a potential with maxima at $`\varphi =0,\pi `$ separated by symmetrical minima. In the special case of resonance orbits, movement in this potential translates into behavior governed by a modified pendulum equation. Previous pendulum-analogue analyses of this problem have usually been formulated via the disturbing function and the continuous resonant argument (Winter & Murray, 1997a; Dermott & Murray, 1983).

Our analysis, specifically that of mapping, is most similar to that of Malyshkin & Tremaine (1999). They consider the evolution of high-eccentricity comet-like orbits in the low-inclination circular restricted three-body problem by integrating numerically to find the energy kick as a function of the resonant angle at periapse and then using this energy kick to create a mapping which takes one periapse passage to the next. However, Malyshkin & Tremaine (1999) are interested in particle orbits that cross the orbit of the secondary, so the form of their energy kick is qualitatively different from ours. In particular, while a small nonzero orbital inclination would barely affect our energy kick function, it could drastically change the shape of the overall energy kick function in the case of planet-crossing orbits. Partly because of this, they do not discuss their energy kick it in terms of a potential. Also, they focus on the chaotic diffusion of the particle toward escape or capture rather than on motion in resonances.

For $`1:N`$ resonance orbits—that is, those we call first-order resonance orbits—the shape of the potential generates analogues of the Lagrangian points for $`N>1`$ resonances. The potential similarly leads to two kinds of libration analogous to the horseshoe and tadpole orbits seen in a $`1:1`$ resonance. $`p:N`$ resonances—that is, those we call higher-order resonances—show only one kind of libration: when the winding number is large, the sum relating the higher-order resonance potentials to the first-order resonance potential eliminates the indirect term responsible for the tadpole analogues.

Several authors discuss the tadpole analogues’ presence or absence in mean-motion resonances in general under the name ‘asymmetric librations’; Nesvorný & Roig (2001) are the only others we know of to refer to the $`1:N`$ resonance librations as ‘tadpoles’ and ‘horseshoes’, though they do not elaborate on this analogy. Some authors have used analytical studies of the Hamiltonian and the disturbing function to set conditions for the existence of asymmetric resonances (Bruno, 1994; Frangakis, 1973; Message, 1970). In particular, Frangakis (1973) analyzed the time-averaged direct and indirect parts of the disturbing function to deduce that only what we call first-order resonances should show asymmetric librations. Bruno (1994) also found analytically that asymmetric librations only exist in what we call first-order exterior resonances. We confirm this and provide a simple physical explanation.

Others have used numerical methods to confirm the existence of asymmetric librations for particular $`1:N`$ resonances and ranges in eccentricity (see, for example, Winter & Murray, 1997b; Beaugé, 1994; Message & Taylor, 1978; Frangakis, 1973; Message, 1958). Some of these also compare their numerical results to expressions for the Hamiltonian correct to first or second order in eccentricity. Although the agreement is generally good for what we call first-order resonances, the Hamiltonian expressions for what we call higher-order resonances tend to predict spurious asymmetric librations. We believe these are due to extra extrema introduced into the potential by the $`\mathrm{cos}(2\mathrm{resonant}\mathrm{angle})`$ term in the Hamiltonian. In some of the more recent studies involving asymmetric librations (Nesvorný & Roig, 2001; Malhotra, 1996) the discussion is framed in terms of the dynamics of the classical Kuiper Belt and so is confined mostly to what we call low-$`N`$ first-order and low-$`p`$ higher-order resonances in the low- to moderate-eccentricity regime.

We find a limit on $`a`$ for stable first-order resonances. Overlap between the resonances creates chaotic regions of $`(a,\varphi )`$ phase space; for semimajor axes larger than some $`a\mu ^{2/5}`$, the resonance centers are overlapped and no stable librations are possible. This is the high-eccentricity analogue of the well-known chaotic criterion $`|a1|\mu ^{2/7}`$ found by Wisdom (1980) for the circular planar restricted three-body problem in the low eccentricity case. We use the Chirikoff criterion for resonance overlap to estimate the location of the onset of chaos. For sufficiently narrow resonances, or small enough $`\mu `$, there exist regions in $`(a,\varphi )`$ space which lie outside all of the resonances but which are not chaotic. In the planar problem we consider, particles interior to the circulating trajectories in these regions are never able to escape from the star-planet system.

The basic framework for the behavior of high-eccentricity orbits and the properties of chaotic regions in $`(a,\varphi )`$ space can be applied to the orbital evolution of small bodies in the solar system. Objects in the Kuiper Belt, for example, are believed to have arrived there via interactions with Neptune (see, for example, Malhotra et al., 2000); we can apply this framework to study their trajectories. Many of these objects are known to be in resonances (see Chiang et al., 2003, for a recent compilation). The mass ratio between Neptune and the Sun is $`\mu _N=4.4\times 10^5`$. Since this is above the critical $`\mu 5\times 10^6`$, Kuiper Belt objects with $`C_J=3`$ are either librating around a resonance or moving chaotically. The latter could, in principle, be ejected as there is no stable circulation for that value of $`\mu =\mu _N`$ and $`C_J=3`$. However, the known Kuiper Belt objects span a range in $`C_J`$ roughly $`2.6<C_J<3.2`$. In the planar problem, $`\mu =\mu _N`$ and for example $`C_J=3.1`$, stable circulations exist, and protect some of these objects from escape. To study the ultimate fate of such Kuiper Belt objects, the effect of inclination must be understood.

Similarly, we might apply this framework to the scattering of small planetesimals by giant protoplanets and could provide insight to numerical integrations such as those of Rasio & Ford (1996); Ford et al. (2001). Studies like this require an investigation of the way in which the energy kicks move orbits through the ‘global chaos’ region surrounding the resonances in the $`(a,\varphi )`$ plane. Although the antisymmetry of $`\mathrm{\Delta }E^1(\varphi )`$ about $`\varphi =0`$ suggests that the orbits should random walk through phase space, effects of nearby resonances (e.g. Malyshkin & Tremaine, 1999) and terms of higher order in $`\mu `$ become important on timescales long enough for escape to become possible. The importance of second order effects may be understood as follows. Since the amount of extra energy needed to escape is $`1/a_{\mathrm{init}}1`$ and the energy kick per orbit is $`\mu `$, we expect that the average number of kicks needed to escape is $`\mu ^2`$. Note that unlike the first order kicks, the $`𝒪(\mu ^2)`$ energy kicks do not average to 0 over the interval $`(\pi ,\pi ]`$ in $`\varphi `$. This is also apparent from Figure 1. Therefore, with $`\mu ^2`$ kicks, the sum of $`𝒪(\mu ^2)`$ effects produced by the energy kicks will be of order unity—that is, of size comparable to the total first order effect.

MP is supported by an NSF Graduate Research Fellowship and RS holds a Sherma Fairchild senior research fellowship. We thank Peter Goldreich for useful discussions.