NEPTUNE'S WILD DAYS: CONSTRAINTS FROM THE ECCENTRICITY DISTRIBUTION OF THE CLASSICAL KUIPER BELT

We wish to use inclinations to separate the cold and hot classicals and then examine the eccentricity distributions of these two populations. Traditionally, the observed cold and hot objects have been separated using one inclination cutoff. However, because of the overlap between the hot and cold components in the bimodal inclination distribution, a single cutoff will necessarily result in the misclassification of hot objects as cold and vice versa. For example, if the classical population follows the model KBO inclination distribution derived by Gulbis et al. (2010) and we were to distinguish between the cold and hot populations using an inclination cutoff icut = 4°, 11% of objects with i < 4° would actually be hot objects and 15% of objects with i > 4° would actually be cold objects. Thus, using icut = 4°, 11% of the objects classified as cold would be "contaminated," and 15% of those classified as hot would be contaminated. In Figure 2, we plot the "contaminated" fraction over a range of values of icut for the cold and hot populations based on three models of the debiased inclination distribution (Brown 2001; Gulbis et al. 2010; Volk & Malhotra 2011). For all three models, less than 10% of the cold classicals are contaminated for icut < 2°, while less than 3% of the hot classicals are contaminated for icut > 6°.

Zoom In Zoom Out Reset image size

Therefore, instead of using a single icut, we divide the classicals into a likely cold population (i < 2°), a likely hot population (i > 6°), and an ambiguous population (2° < i < 6°). We then examine the eccentricity distributions of the likely cold and likely hot populations, which are "uncontaminated" samples. We use the uncontaminated eccentricity distributions to identify major features that models must match. In Appendix A, we confirm that our results are consistent if we probabilistically include the ambiguous population.

We wish to identify features of the eccentricity distribution that are sculpted during Neptune's wild days, not by the long-term stability of the region under the influence of the modern solar system planetary configuration or by observational bias. First, we compare the eccentricities of observed likely cold (i < 2°) and likely hot (i > 6°) objects to the survival map of Lykawka & Mukai (2005), generated from a 4 Gyr simulation. Lykawka & Mukai (2005) generated initial conditions for test particles uniformly filling a cube of (a, e, i) in the classical region: 41.375 AU < a < 48.125 AU, 0 < e < 0.3, and 0° < i < 30°. They then performed a 4 Gyr numerical integration including the test particles and the four giant planets (starting on their modern orbits). Then they computed the survival rate of test particles in bins of (a, e) and (a, i). In this work, we consider only the eccentricity survival map. This map bins over all inclinations.² After 4 Gyr of evolution under the influence of the planets in their current configuration, KBOs with an initially uniform eccentricity distribution would be distributed according to this survival map.

However, rather than following the survival map, the observed cold and hot objects exhibit major distinct features. In Figure 3, we plot the sample of observed classical objects from the Minor Planet Center (MPC) and the Canada–France Ecliptic Plane Survey (CFEPS; Kavelaars et al. 2009) on top of the Lykawka & Mukai (2005) stability map. The cold objects are confined to very low eccentricities. From 42.5 to 44 AU, the cold objects appear to be confined to e < 0.05. From 44 to 45 AU, the cold objects appear confined below e < 0.1. This confinement of the cold classicals to below the survival limit implies that they were not excited above these levels because if they had been, we would still observe objects at higher eccentricities. Similarly, Kavelaars et al. (2009) found that classical objects with i < 45 are restricted to 42.5 AU < a < 45 AU and e < 0.1. In contrast, hot objects occupy the upper portion of the survival region and appear uniformly distributed in a from 42 to 47.5 AU. Suggestively, they also appear to be distributed roughly along a scattering line, as if they were scattered into the classical region but did not have time to evolve to low eccentricities before Neptune's eccentricity damped.

Zoom In Zoom Out Reset image size

We use the following major qualitative features to place constraints on Neptune's dynamical history. We consider these criteria "conservative" because they allow for dynamical histories at the very edge of consistency with the observations.

Cold population: confined to low eccentricities of e < 0.1 in the region from 42.5 to 45 AU. In the region between 42.5 and 45 AU, the cold objects have eccentricities well below the distribution that follows the survival map. Therefore, Neptune cannot excite the cold classical objects in this region above e = 0.1. (We choose this value to be conservative in ruling out regions of parameter space and to match Kavelaars et al. (2009), but it appears that cold objects with semimajor axes less than 44 AU are confined below e < 0.05, a tighter constraint.) We indicate this threshold as a solid yellow line in Figure 3.

Hot population: delivered to the upper survival region with q > 34 AU out to 47.5 AU. The observed hot objects occupy the upper portion of the survival region. Therefore, a consistent dynamical history should allow some objects to reach this region. It is not necessary for the transported objects to reach very low eccentricities, only low enough to survive under the current solar system planetary configuration. We set the criterion that the hot classicals must be delivered to periapse q > 34 AU (dashed line in Figure 3) from 42 to 47.5 AU, the edge of the observed population.

In determining which major features serve as constraints on the dynamical history of the solar system, we address several complications.

• 1.  
  The inclinations of objects vary over time (Volk & Malhotra 2011).
• 2.  
  The inclination cutoff between the hot and cold classicals is model dependent.
• 3.  
  Proper elements are more robust than the observed instantaneous elements.
• 4.  
  The features in the eccentricity distributions might be the result of random chance or small number statistics.
• 5.  
  The eccentricity distributions may be impacted by observational bias.

The first complication is addressed by Volk & Malhotra (2011). They find that, at any given time, only 5% of objects will be inconsistent with their original inclination-based classification of hot versus cold. Therefore, we expect the major qualitative features we identify to hold despite variations in the inclinations of some objects.

The second complication is that the cutoff inclination between the hot and cold classicals depends on the parameters and form of the model for the bimodal inclination distribution. The three models we consider (Brown 2001; Gulbis et al. 2010; Volk & Malhotra 2011) each use the functional form of sin i multiplied by a Gaussian but use different widths and cold/hot fractions. They also use different planes for the inclination: Brown (2001) defines the inclination with respect to the ecliptic plane, Gulbis et al. (2010) with respect to the mean plane of the Kuiper Belt, and Volk & Malhotra (2011) with respect to the invariable plane. However, despite these differences, i < 2° and i > 6° are robust cutoffs for establishing an uncontaminated cold and hot population, respectively, for each of the three models. In regard to the functional form of the model, Volk & Malhotra (2011) find that the high-inclination component is not well described by a Gaussian, and Fabrycky & Winn (2009) argue that the most robust generic functional form for a distribution of inclinations is a Fischer distribution. However, the discrepancies between different functional forms are strongest for classifying objects in the intermediate, overlapping portion of the bimodal inclination distribution (K. Volk 2011, private communication), so we argue that our approach of definitively classifying only the "uncontaminated" low- and high-inclination objects is robust.

Regarding the third complication, the stability map of Lykawka & Mukai (2005) is formulated in terms of instantaneous eccentricity and inclination, but the most robust, non-varying formulation of the orbital elements are the proper, or free, elements. However, none of the model inclination distributions are formulated in terms of the proper inclination, nor is the stability map of Lykawka & Mukai (2005) formulated in the proper elements. To compare "apples to apples," we use the instantaneous orbital elements in the plots in this section. We use the proper elements in Appendix A and find that the observational features we identify (Section 2.2) still hold.

Addressing the fourth complication, in Appendix A, we confirm that the confinement of the cold population to low eccentricities is statistically significant. For the hot population, we only impose the constraint that the models must deliver them to the long-term stable region (Section 2.2); we will demonstrate that this loosely formulated restriction ends up imposing strong constraints on Neptune's dynamical history.

2.3.1. Ruling Out Observational Bias through Statistical Tests

We would not expect observational bias (the final complication) to cause the cold classicals to appear to be confined to low eccentricities; KBO searches are more likely to preferentially observe high-eccentricity (i.e., small periapse) objects of a given semimajor axis. However, to ensure that the features on which we base our constraints (Section 2.2) are not created by observational bias, we perform the following test to see if observational bias could generate them:

• 1.  
  We begin by generating a simulated sample of objects uniformly distributed in (a, e). We set the inclinations to follow the unbiased inclination distribution of the classicals, as modeled by Gulbis et al. (2010).
• 2.  
  Then we use the stability map of Lykawka & Mukai (2005) to transform this simulated sample following a uniform eccentricity distribution into a sample following the eccentricity distribution shaped by the four giant planets under their current configuration. We call this process "filtering." For each simulated object, we obtain a predicted survival rate from the (a, e) stability map of Lykawka & Mukai (2005). Then we select a uniform random number between 0 and 1. If the randomly selected number is less than the predicted survival rate, we include the object in our sample. Note that since the survival rates of Lykawka & Mukai (2005) are given as a range (e.g., 10%–20%, 90%–100%), we repeated this entire test (i.e., steps 2–4) three times, once using the minimum of each range, once using the mean of each range, and once using the maximum of each range. As expected, the resulting eccentricity distribution had higher (lower) eccentricities when we used the maximum (minimum) of each range, but the major features we identified still held. The simulated population in Figure 4 uses the mean. Next, we transformed the "survival-rate-filtered" sample from step 2 into an observed sample.
• 3.  
  We randomly assign each object an H magnitude³ between 6 and 8.
• 4.  
  Then we apply the L7 Survey Simulator for the well-characterized CFEPS (Kavelaars et al. 2009).
• 5.  
  We compare the final simulated distribution to the subset of objects that were detected by CFEPS (Kavelaars et al. 2009; Figure 4). Figure 4 is analogous to Figure 3. It includes the simulated distribution (circles) and only the subset of KBOs observed by CFEPS. The simulated distribution is not confined to e < 0.1 from 42.5 to 45 AU, confirming that this feature of the observed eccentricity distribution does not result from observational bias. Note also that simulated hot objects are found at lower eccentricities than observed.

Zoom In Zoom Out Reset image size

A single origin for the hot and cold populations seems unlikely. If the cold and hot classicals formed together in the classical region, where they are observed today, it is difficult to imagine a process that would excite the hot population while leaving the cold population confined to low eccentricities. Hahn & Malhotra (2005) proposed a scenario in which the classical region has been pre-excited. However, this scenario does not produce a population of cold classicals confined to low eccentricities. Moreover, if both the cold and hot classicals were transported from the inner disk, it seems unlikely that a common deposition process would place the cold classicals solely at low eccentricities. Levison & Morbidelli (2003) and Levison et al. (2008) propose scenarios in which both the hot and cold classicals are transported from the inner disk.

In the scenario of Levison & Morbidelli (2003), the cold classicals were pushed outward by the 2:1 resonance and dropped during stochastic migration, while the hot classicals scattered off of Neptune. The feasibility of this mechanism depends on the size distribution of planetesimals because the migration needs to be stochastic in order to drop objects from resonance. When Neptune scatters a planetesimal inward and the planetesimal is ejected by Jupiter, Neptune experiences a net gain in angular momentum and migrates outward. If the planetesimals are small, this is a smooth process, but if they are large, it is a jumpy, stochastic process, in which KBOs can be dropped from resonance. See Murray-Clay & Chiang (2006) for a detailed exploration of stochastic migration; they conclude that planetesimal-driven migration cannot generate the necessary stochasticity unless a large fraction of planetesimals formed very large. This constraint merits a fresh look in light of new planetesimal formation models (see Chiang & Youdin 2010, and references therein). However, even if extreme planetesimal properties allowed this mechanism to work, objects dropped from the 2:1 resonance would have a range of eccentricities, not confined solely to low eccentricities. Therefore, this mechanism holds more potential for producing the hot population than the cold population.

In the scenario of Levison et al. (2008), the cold classicals are objects that, like the hot classicals, were scattered into the Kuiper Belt by an eccentric Neptune but, unlike the hot classicals, evolve down to low eccentricities in regions near resonances. However, this mechanism would create a range of eccentricities for the cold classicals and thus would have trouble producing the confined eccentricities of the region of 42.5 AU < a < 45 AU (Figure 3). They do find some correlation between a particle's final inclination and its initial semimajor axis in one of their simulations (Levison et al. 2008, Figure 11, panel (b)), which may be able to partially account for a difference in physical properties between low- and high-inclination objects. However, transporting the cold classicals from the inner disk is not consistent with the finding by Parker & Kavelaars (2010) that wide binaries—of which the cold population contains a number—cannot survive transportation from the inner disk to the classical region.

3.1.1. The Eccentricity Distribution Was Not Sculpted Solely by a Different Stability Threshold in the Past

Thus, throughout the rest of this paper, we consider the general scenario—also discussed in Morbidelli et al. (2008)—in which the hot objects are transported to the classical region from the inner disk and the cold objects form in situ in the classical region. The cold objects must not be dynamically excited, as quantified by the criterion we established in Section 2.2. In Figure 5, we show a conceptualization of this model. This general scenario encompasses previous models and allows Neptune to undergo any potential combination of high eccentricity, migration, and/or apsidal precession with a range of initial eccentricities and semimajor axes. If our constraints do not rule out all of parameter space for this model, it may be possible to produce both the hot and cold classical populations. Otherwise, a major physical process is missing from current models.

Zoom In Zoom Out Reset image size

For now, we can think of Neptune as having a high eccentricity at one location after undergoing planet–planet scattering. Our results hold in more complicated scenarios as well, as we will describe in Section 5.4.

We have chosen features (Section 2.2) that are not shaped by the long-term survival of KBOs under the current planetary configuration. Therefore, we can focus on modeling the processes that affect these features during the interval of Neptune's wild days instead of treating the entire 4 Gyr. We model these processes analytically in Section 4 and validate our analytical expressions using numerical integrations.

Unless otherwise specified, the integrations are performed as follows. We perform the integrations using the Mercury 6.2 hybrid symplectic integrator (Chambers 1999), an N-body code that allows massless test particles. We employ an accuracy parameter of 10⁻¹² and a step size of 200 days and impose (if applicable) the migration and damping of Neptune's orbit through user-defined forces and velocities. Migration and damping follow the following forms:

$$\begin{eqnarray} a_{\rm N}&=& (a_{\rm N})_f + ((a_{\rm N})_{0} - (a_{\rm N})_{f}) \exp {(-t/\tau _{a_{\rm N}})} \nonumber \\ e_{\rm N}&=& (e_{\rm N})_0 \exp {(-t/\tau _{e_{\rm N}})} , \end{eqnarray} \tag{ 1 }$$

where a_N is the semimajor axis of Neptune at time t, (a_N)₀ is the initial semimajor axis, and (a_N)_f = 30.1 AU is the final semimajor axis. At time t, the eccentricity of Neptune is e_N; Neptune's initial eccentricity is (e_N)₀. The forced evolution of Neptune's orbit is implemented through modifications to Mercury 6.2, described in detail in the Appendix of Wolff et al. (2012). The migration and damping are parameterized by timescales $\tau _{a_{\rm N}}$ and $\tau _{e_{\rm N}}$, respectively, which we specify in the text in the applicable cases. When specified, Neptune is forced to undergo apsidal precession using an artificial stellar oblateness force, built into Mercury 6.2, that we modified to apply only to Neptune, parameterized by a J₂ coefficient chosen to produce the correct precession rate. Unless otherwise noted, Neptune is the only planet included in the integration. The KBOs are modeled as 600 massless test particles with initial a evenly spaced between 40 and 60 AU and initial e = i = 0. The migration, damping, and apsidal precession are not applied to the KBOs, only to Neptune.

The purpose of this paper is to explore whether the generalized scenario described here can even work, i.e., whether it is ever possible to transport the hot classicals from the inner disk to the classical region without disrupting an in situ cold population. Yet alternative scenarios exist that do not fit within this framework, such as the additional planet beyond Pluto proposed by Lykawka & Mukai (2008) and others. Furthermore, we will describe in the conclusion how additional constraints could rule out the generalized model we consider. In that case, development of alternative scenarios would be necessary. Obviously, the constraints we will place do not necessary hold for a scenario that is not encompassed by our general model.

In the generalized model we explore, Neptune may be scattered outward from the inner solar system onto a highly eccentric orbit. After this occurs, Neptune's new orbit crosses the orbits of some planetesimals in the inner disk (see Figure 5), which scatter off the planet. This mechanism can potentially deliver hot objects from the inner disk into the classical region.

The region into which Neptune can scatter objects is defined by the planet's semimajor axis a_N and eccentricity e_N. Neptune can scatter objects outward to periapses q = a(1 − e) between Neptune's periapse r_{p, N} = (a_N − a_H)(1 − e_N) and apoapse r_{a, N} = (a_N + a_H)(1 + e_N). In Figure 6, we show examples of the region into which Neptune can scatter KBOs for two sets of parameters (a_N, e_N). We have adjusted r_{p, N} and r_{a, N} to include the Hill sphere radius, a_H (∼1 AU), the distance from Neptune at which Neptune's gravity overcomes the Sun's tidal gravity. Particles that enter Neptune's Hill sphere will be scattered and, if they are scattered outward into the classical region, will reach—under the approximation that the scattering location becomes the particle's new periapse—a given semimajor axis a with eccentricities between 1 − r_{p, N}/a and 1 − r_{a, N}/a.

Zoom In Zoom Out Reset image size

4.1.1. Implications of Scattering for the Cold Classicals

Secular forcing has a strong effect on the behavior of KBOs during the potential period when Neptune's eccentricity is high. (See Murray & Dermott 2000 for a pedagogical description of first-order secular theory outside of mean motion resonance.) On a timescale of a million years, a cold object is excited to a higher eccentricity via Neptune's secular forcing. (The direct forcing from the other planets on the KBO is negligible and the other planets effectively only affect the KBOs via Neptune, as we will demonstrate in Section 4.2.3.) A hot object that is scattered into the classical region also experiences secular forcing, which can decrease its eccentricity so that its orbit no longer crosses Neptune's. Hot or cold, an object's eccentricity is a vector combination of its forced eccentricity—set by Neptune's eccentricity, Neptune's semimajor axis relative to the object's, and Neptune's apsidal precession rate—and of the object's free eccentricity, which is set by its initial condition before Neptune is scattered to a high eccentricity. The object's free eccentricity precesses about the forced eccentricity at the secular frequency g_KBO. Therefore, as we will demonstrate, a cold object has a well-defined excitation time and amplitude, and a hot object will have a minimum eccentricity it can reach after being scattered into the classical region. Thus, while Neptune's eccentricity is high, secular forcing potentially is an important mechanism for exciting the cold objects and stabilizing the hot ones. As we will show, when Neptune's eccentricity damps quickly, the orbits of the KBOs are "frozen" near the eccentricities they reached through secular evolution.

4.2.1. Basic Secular Evolution

First, we define expressions for the secular evolution of a test particle under the influence of a planetary system containing only Neptune and the Sun. The components of the particle's eccentricity vector are h = esin ϖ and k = ecos ϖ, where ϖ is the particle's longitude of periapse. Secular forcing by Neptune causes h and k to evolve as (to first order in e and e_N)

$$\begin{eqnarray} h &=& e_{\rm free} \sin (g_{\rm KBO}t + \beta) + e_{\rm forced} \sin (\varpi _{\rm N}) \nonumber \\ k &=& e_{\rm free} \cos (g_{\rm KBO}t + \beta) + e_{\rm forced} \cos (\varpi _{\rm N}), \end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray} e_{\rm forced} &=& \frac{b_{3/2}^{(2)}(\alpha)}{ b_{3/2}^{(1)}(\alpha) } e_{\rm N}\nonumber \\ \alpha &=& \frac{a_{\rm N}}{a} \\ g_{\rm KBO}&=& \alpha b_{3/2}^{(1)}(\alpha) \frac{m_{\rm N}}{m_\odot } \frac{n}{4}. \nonumber \end{eqnarray} \tag{ 3 }$$

The constants e_free and β are determined from the initial conditions, and the particle's forced eccentricity is e_forced. Here, ϖ_N is the longitude of periapse of Neptune, e_N is the eccentricity of Neptune, and α is the ratio of Neptune's semimajor axis to that of the particle, all of which are assumed to be constant. The functions b are standard Laplace coefficients (see Murray & Dermott 2000). The secular frequency of the KBO is g_KBO, m_N is the mass of Neptune, m_☉ is the mass of the Sun, n = (Gm_☉/a³)^1/2 is the particle's mean motion, and G is the universal gravitational constant.

Consider a cold object with e = 0 at t = 0. In our approximation, Neptune, having been scattered by the other giant planets, effectively instantaneously appears and imparts a forced eccentricity of e_forced. From these initial conditions, e_free = e_forced. Then, the KBO's forced eccentricity vector remains fixed and the object's total free eccentricity vector precesses about the forced eccentricity. Thus, its total eccentricity varies sinusoidally from e = 0 to e = 2e_forced on a timescale set by g_KBO.

A hot object scattered into the classical region, in contrast, has an eccentricity e at t = 0. The magnitude of its free eccentricity is a value between max (0, e − e_forced) and e + e_forced, depending on the initial location of its periapse relative to Neptune's. Over a timescale set by g_KBO, its total eccentricity oscillates. Depending on the initial conditions, it may reach an eccentricity low enough so that its orbit no longer crosses Neptune's and/or so that it is stable under the current configuration of the giant planets.

An example of the secular evolution of cold objects "going up" and hot objects "going down" in eccentricity is shown in Figure 7, highlighting the tension between the evolution of hot objects to low eccentricities and the evolution of cold objects to high eccentricities. The cold objects (red) begin with e = 0 (see Section 3.3 for a general description of the integrations we performed). The hot objects in the integration (blue) all begin with e = 0.2 and ϖ = ϖ_N + π/3, for the purposes of illustrating secular evolution.⁴ As time progresses through three snapshots, the cold objects become excited and the hot objects reach low eccentricities. The analytical model from Equation (2) matches well except near mean motion resonances, where the secular evolution is much faster than predicted. We also overplot a more accurate analytical expression that includes a resonant correction term, which we will derive in Section 4.2.2.

Zoom In Zoom Out Reset image size

We note that throughout the paper, we perform numerical integrations of objects with initial semimajor axes out to 60 AU to give a better conceptual picture of the secular excitation. Moreover, depending on where the initial population was truncated in the solar system's planetesimal disk, it is possible that additional classical KBOs will be discovered beyond 48 AU in the future, and we would like to make testable predictions. Finally, the integration results should be interpreted as examples: since we cannot show a figure for every possible combination of parameters for Neptune, we plot out to 60 AU to let the reader imagine the results if Neptune's semimajor axis were smaller.

4.2.2. Refined Secular Expression

The secular expression in Section 4.2.1 is valid to first order in eccentricity, neglects the effects of orbital resonances, and applies in the case in which Neptune's orbit does not apsidally precess. However, these neglected effects can significantly alter a secularly evolving KBO's behavior.

• 1.  
  Neptune's high eccentricity makes terms of order e²_N (ignored in deriving Equation (2)) non-negligible. As we will show, these extra terms result in a faster secular forcing frequency of the KBO.
• 2.  
  Proximity to mean motion resonances with Neptune significantly alters the secular frequencies of KBOs (as shown for other solar system bodies in Hill 1897; Malhotra et al. 1989; Minton & Malhotra 2011). Following Malhotra et al. (1989), we incorporate resonance correction terms, described in our Appendix B.1.2. These resonance correction terms are very important for objects near resonance but not valid for objects librating in resonance (we are not considering resonant objects in this paper⁵).
• 3.  
  As Neptune's eccentricity is damped, the particle's forced eccentricity goes to zero. If the damping occurs over a timescale $\tau _{e_{\rm N}}$ shorter than the secular oscillation time, the particle's eccentricity is frozen at approximately the value it reaches at the eccentricity damping time. If the damping occurs over a longer timescale, the particle's eccentricity evolves to its initial free eccentricity.
• 4.  
  Apsidal precession of Neptune alters the forced eccentricity, keeping it low when Neptune precesses quickly. It also alters the oscillation timescale of the total eccentricity, because now both the free and forced eccentricities are precessing.
• 5.  
  Migration alters the secular frequencies and shifts the locations of the resonances.

Based on these considerations, we modify the standard expression (Equation (2)) for the secular evolution of a test particle's h and k under the influence of an eccentric Neptune to incorporate these effects:

$$\begin{eqnarray} h &=& e_{\rm free, 0} \sin (g_{\rm KBO}t + \beta _0) + \bar{e}_{\rm forced} \sin (\varpi _{{\rm N},0} + \dot{\varpi }_{\rm N}t) \nonumber \\ k &=& e_{\rm free, 0} \cos (g_{\rm KBO}t + \beta _0) + \bar{e}_{\rm forced} \cos (\varpi _{{\rm N},0}+ \dot{\varpi }_{\rm N}t) \qquad \end{eqnarray} \tag{ 4 }$$

$$\begin{equation} g_{\rm KBO}= \left(1+\frac{f_5}{f_2}e_{\rm N}^2\right)\alpha b_{3/2}^{(1)}(\alpha) \frac{m_{\rm N}}{m_\odot } \frac{n}{4} + \epsilon \delta g_{\rm KBO} \end{equation} \tag{ 5 }$$

$$\begin{equation} \delta g_{\rm KBO}= \alpha \big(C x e_{\rm N}^{x-1}\big)^2 \frac{m_{\rm N}}{m_\odot } n \end{equation} \tag{ 6 }$$

$$\begin{equation} \bar{e}_{\rm forced} = \sin ({\rm min} (g_{g_{\rm KBO}} \tau _{e_{\rm N}},\pi /2)) \frac{g_{\rm KBO}^{\prime }}{\dot{\varpi }_{\rm N}- g_{\rm KBO}} e_{\rm N}(t) \end{equation} \tag{ 7 }$$

$$\begin{equation} g_{\rm KBO}^{\prime } = {-} \left(1+\frac{f_{10}}{f_{11}}e_{\rm N}^2\right) \alpha b_{3/2}^{(2)}(\alpha) \frac{m_{\rm N}}{m_\odot } \frac{n}{4}. \end{equation} \tag{ 8 }$$

Compare Equation (4) to Equation (2). The form is the same, but Equation (4) has several key differences and new variables, which we will proceed to discuss and define throughout the remainder of this subsection. One distinction is that several quantities that were fixed in Equation (2) (i.e., ϖ_N, α, e_N) can now vary with time. Additionally, we have incorporated corrections for Neptune's high eccentricity and for the potential proximity of the KBO to mean motion resonance with Neptune. The impact of the time-varying e_N is the most complicated, so we leave it for last.

As our first correction, we allow the longitude of pericenter of Neptune, ϖ_N, to precess. Neptune's precession adds an extra term $\dot{\varpi }_{\rm N}$ to Equation (4) (where we rewrite ϖ_N as a linear function of time: $\varpi _{{\rm N},0} + \dot{\varpi }_{\rm N}t$) and to Equation (7), in analogy to the standard secular theory for the four-planet case. We note that the other giant planets impact the secular evolution of the KBOs indirectly by causing Neptune's eccentricity to precess (see Section 4.2.3 for discussion).

The precession rate $\dot{\varpi }_{\rm N}$ has the same role (and same place, in the denominator of the forced eccentricity) in the just-Neptune secular theory as in the standard four-planet secular theory (Appendix B.2) for the particle except that we are specifying Neptune's evolution via $\dot{\varpi }_{\rm N}$ instead of constructing a secular theory for the planets that produces a particular $\dot{\varpi }_{\rm N}$. When $\dot{\varpi }_{\rm N}$ is large, the effective forced eccentricity $e_{\rm forced} = |\bar{e}_{\rm forced} |$, defined below, remains low because the forced eccentricity is inversely proportional to the precession rate for $|\dot{\varpi }_{\rm N}| \gg | g_{\rm KBO}|$ (Equation (7)).

When Neptune migrates, α changes with time. In Wolff et al. (2012), we found that migration occurs in three regimes, relative to the eccentricity damping timescale $\tau _{e_{\rm N}}$: fast, comparable, and slow. When Neptune's migration timescale $\tau _{a_{\rm N}}$ is slow relative to the damping time, the secular evolution of the KBOs effectively takes place as if Neptune remains at its initial location. When Neptune's migration is fast relative to the damping time, the secular evolution of the KBOs effectively takes place as if Neptune were always at its final location. In the intermediate case, in which $\tau _{e_{\rm N}} \sim \tau _{a_{\rm N}}$, modeling the secular evolution at the location Neptune reaches after half a damping time is a fair approximation. In this work, we therefore use a fixed α, which should be chosen according to these principles for a given evolution of Neptune. Figure 8, which we will describe after discussing our treatment of Neptune's eccentricity, provides an example showing that this approach is effective.

Zoom In Zoom Out Reset image size

Before considering the impact of e_N varying with time, we discuss the correction terms for Neptune's high eccentricity and for resonances. The correction terms for Neptune's high eccentricity do not change the form of the secular evolution. We have applied a correction term, (f₅/f₂)e²_N, in our expression for g_KBO (Equation (5)) and another such factor, (f₁₀/f₁₁)e²_N, in g'_KBO (Equation (8)), which is another eigenfrequency. These terms are derived, the f factors (which are of order unity) defined, and their necessity demonstrated in Appendix B.

Proximity to resonance changes the secular frequency g_KBO (Equation (5)), as described by Malhotra et al. (1989). Orbital resonances greatly increase the secular forcing frequency because terms in the disturbing function that depend on the resonant angle can no longer be averaged over. The amplitude of the resonant correction term is defined in Appendix B and depends on how close the particle is to the location of mean motion resonance. The frequency, δg_KBO, defined in Equation (6), depends on the order of the resonance x and a constant C, of order unity, that is different for each resonance. See Figure 7 for a demonstration of the resonance correction terms.

Finally, we turn to the impact of eccentricity damping, which alters e_N. Instantaneously, the KBO has the eccentricity components $h = e \sin \varpi = e_{\rm free} \sin \phi + \bar{e}_{\rm forced} \sin \varpi _{\rm N}$ and $k = e \cos \varpi = e_{\rm free} \cos \phi + \bar{e}_{\rm forced} \cos \varpi _{\rm N}$, where ϕ = g_KBOt + β. However, now the forced eccentricity vector ($\bar{e}_{\rm forced} \cos \varpi _{\rm N}, \bar{e}_{\rm forced} \sin \varpi _{\rm N}$) is changing. We have already accounted for the apsidal precession, but now the magnitude of the forced eccentricity vector is changing as well. When e_forced does not change, e_free is a constant determined by initial conditions. This remains true when e_forced evolves slowly compared to the secular forcing time of the KBO. Otherwise, e_free changes. Instead of allowing both e_free and e_forced to change with time, we use e_{free, 0} and define an "effective" forced eccentricity, $| \bar{e}_{\rm forced}|$ (Equation (7)). The "effective" forced eccentricity changes with time proportionally to e_N. Throughout the rest of the paper, we will refer to e_{free, 0} as e_free and the "effective" forced eccentricity as e_forced. The "effective" forced eccentricity includes a factor $\sin ({\rm min} (g{_{g_{\rm KBO}}} \tau _{e_{\rm N}},\pi /2))$. When Neptune's eccentricity damping timescale is long compared to the KBO's secular oscillation period (i.e., $g_{\rm KBO}\tau _{e_{\rm N}} > \pi /2$), the particle's total eccentricity damps to its initial free eccentricity, e_{free, 0}, and this factor is unity. However, when Neptune's eccentricity damps quickly, the particle's total eccentricity damps to a value near the eccentricity it reached after one damping time. The empirical correction factor allows us to model the particle's behavior without altering the form of the secular evolution. This empirical factor provides a match to the integrations (see Figure 8).

The modified secular theory, Equation (4), matches the integrations even when damping and migration are included (Figure 8). Figure 8 shows an example of a case in which $\tau _{e_{\rm N}} = 0.3$ Myr is shorter than its migration timescale, $\tau _{a_{\rm N}} = 5 \,{\rm Myr}$. (See Section 3.3 for a description of how we implemented the damping and migration of Neptune's orbit.) In the top row, Neptune undergoes eccentricity damping from e = 0.3 at a constant semimajor axis a_N  = 28 AU. In the middle row, Neptune migrates from 28 AU to 30 AU on a timescale $\tau _{a_{\rm N}} = 5 \,{\rm Myr}$, without its eccentricity damping. Unlike any other model in this paper, we use a time-dependent α(t) for the analytical model. In the bottom row, Neptune undergoes both damping and migration, but we plot again the same analytical model as in row 1, in which Neptune undergoes only eccentricity damping (no migration). The model overplotted in row 3, even though it in no way includes the effects of migration, matches well. When the migration is slow compared to the damping, the change in the secular frequency g_KBO is negligible over the timescale during which Neptune's eccentricity is high. It is as if Neptune's eccentricity damps while Neptune remains at its initial a_N. Thus, when the migration is slow compared to the damping timescale, we can model the KBOs' secular evolution as if Neptune's eccentricity damps while Neptune remains in place.

4.2.3. Effects of Other Planets

Our approach of modeling only Neptune, undergoing a range of orbital histories that would be caused by interactions with the other giant planets and with the solar system's planetesimal disk, is sufficient because the other planets primarily affect the KBOs only indirectly, through influencing Neptune. The influence on the Kuiper Belt of a single, apsidally precessing Neptune matches the influence of multiple planets in both integrations and theory.

An illustrative case is shown in Figure 9. The initial conditions for the particles are the same as in Figure 8, and Neptune has a_N = 28 and e_N = 0.2. The top panel (four planets) and middle panel (just Neptune undergoing precession) are very similar, keeping the objects at lower eccentricities than in the bottom panel (just Neptune, no precession). Thus, precession must be included in the parameter space exploration, and including precession successfully accounts for the influence of the other giant planets. We note that in the time of the snapshot (1.4 Myr), the particles have secularly evolved to have eccentricities large enough so that their orbit, at the proper orientation, could intersect Neptune's. In cases in which Neptune precesses (top two panels), the objects are scattered by Neptune as Neptune's orbit precesses to intersect the orbits of the KBOs with eccentricities above the scattering line. The cutoff is at 45 AU because, interior this location, particles are secularly evolving quickly due to their proximity to resonance and have thus reached high eccentricities, allowing them to scatter. In the final panel, the particles are not scattered because Neptune's orbit does not precess to intersect the orbits of the particles.

Zoom In Zoom Out Reset image size

The other planets matter in that they affect Neptune, but their direct effect on the KBOs is negligible. Their main effect is to cause precession of Neptune. Interactions between Uranus and Neptune, if they closely approach the 2:1 resonance while Neptune's eccentricity is high, cause additional effects that we discuss in Section 4.5. These effects can also be modeled using only Neptune, with appropriate orbital variations.

Treating only Neptune reduces the number of free parameters, allowing a more thorough exploration of the restricted space. The constraints we will develop can be applied to more extensive models that include the other giant planets. See Appendix B.2 for a mathematical discussion of how the full four-planet secular theory reduces to the just-Neptune case.

4.2.4. Constraints from Secular Excitation of Cold Objects

It follows from the expressions in Section 4.2.2 that the excitation of the cold classicals happens on timescales of millions of years (or shorter near resonances) with an amplitude and timescale that depend on Neptune's semimajor axis, eccentricity, eccentricity damping timescale, migration rate, and precession rate. Complementarily, hot objects scattered into the classical region can evolve to lower eccentricities on similar timescales.

The observations require that the cold classicals not be excited above e > 0.1 in the region 42.5 AU < a < 45 AU, as demonstrated in Section 2. From Equations (4) and (5) derived in Section 4.2, it follows that the cold classicals will not be excited above e > 0.1 at a given location if (Constraint 1)

$$\begin{equation} \sin ({\rm min}(g_{\rm KBO}\tau _{e_{\rm N}},\pi /2)) | \frac{g_{\rm KBO}^{\prime }}{\dot{\varpi }_{\rm N}- g_{\rm KBO}} | e_{\rm N}< 0.1. \end{equation} \tag{ 9 }$$

Thus, there are three possible regimes to "preserve" an in situ population of cold classicals through Neptune's wild days.

• 1.  
  The eccentricity of Neptune is small enough that the region's forced eccentricity, proportional to e_N, is below 0.1 (i.e., $|{g_{\rm KBO}^{\prime }}/({\dot{\varpi }_{\rm N}- g_{\rm KBO}}) | e_{\rm N}< 0.1$).
• 2.  
  Neptune's periapse precesses quickly enough that the region's forced eccentricity, inversely proportional to $\dot{\varpi }_{\rm N}$, is below 0.1 (i.e., $\dot{\varpi }_{\rm N}$ is large).
• 3.  
  The eccentricity of Neptune damps quickly enough that the objects are not excited above 0.1 (i.e., τ_e is small).

Hot objects that have been scattered into the classical region from the inner disk (Section 4.1) will undergo secular evolution when they arrive in the classical region. They will reach a given semimajor axis a in the classical region with eccentricities between 1 − r_{a, N}/a and 1 − r_{p, N}/a. Not all of these eccentricities are consistent with stable orbits over 4 Gyr (Figure 3). If a particle is scattered to a high eccentricity above the stable region, under certain conditions—if $\tau _{e_{\rm N}}$ is not too fast and Neptune imparts a forced eccentricity that is large enough relative to the particle's free eccentricity—the particle can reach a region of long-term survival through secular evolution.

In Figure 10, we show two examples of KBOs that are scattered into the classical region from an integration resembling Levison et al. (2008) Run B. In this integration, Neptune begins with a_N = 28.9 AU, e_N = 0.3. Its eccentricity damps on a timescale of τ_e = 2 Myr, and it undergoes migration to 30.1 AU with τ_a = 10 Myr. Uranus begins with a = 14.5 and e = 0 and also undergoes migration, to 19.3 AU, on the same timescale. Jupiter and Saturn begin with a = 5.2 AU, e = 0.05 and a = 9.6 AU, e = 0.05, respectively, and do not undergo migration or eccentricity damping. The integration includes 24,000 test particles, half of which (following Levison et al. 2008) begin in the region from 20 to 29 AU with e = 0.2 and half of which begin in the region from 29 to 34 AU with e = 0.15. The two example particles shown were among the group of particles found in the stable classical region in this integration after 4 Gyr and exhibit typical behavior. The particles are scattered into the classical region above the region of survival (Figure 3; below dotted line in right panel) but secularly evolve down into the stable region. After 4 Gyr, the particles remain in this location. Neither of these example objects is librating in an orbital resonance.

Zoom In Zoom Out Reset image size

Therefore, Neptune's apoapse r_{p, N} must be large enough so that the hot classicals either reach the region of long-term survival immediately or can evolve there before the eccentricity of Neptune damps, after which the eccentricity of the particle is frozen.

However, the post-scattering evolution depends strongly on ϖ − ϖ_N. Particles are not actually scattered to orbits with independent, random ϖ − ϖ_N. This is because, by definition, after each scattering, the particle's new orbit fulfills the condition that at θ, the angle at which the orbit of the particle and the orbit of Neptune intersect, r = r_N:

$$\begin{equation} \frac{a(1-e^2)}{1+e\cos (\theta -\varpi)} = \frac{a_{\rm N}\big(1-e_{\rm N}^2\big)}{1+e_{\rm N}\cos (\theta -\varpi _{\rm N})}. \end{equation} \tag{ 10 }$$

This is important because once the particle is scattered, it begins to undergo secular oscillations and the initial phase of the oscillation, β₀ (Equation (4)), depends on ϖ − ϖ_N.

Since in the classical region most of the particle's orbit is outside of Neptune's orbit, the orbits will intersect close to the particle's periapse, i.e., the interior part of its orbit, so θ ≈ ϖ. When θ = ϖ, Equation (10) simplifies to

$$\begin{equation} q = a (1-e) = \frac{a_{\rm N}\big(1-e_{\rm N}^2\big)}{1+e_{\rm N}\cos (\varpi -\varpi _{\rm N})}, \end{equation} \tag{ 11 }$$

and thus ϖ maps exactly to the particle's post-scattering periapse q. For example, particles scattered to the minimum q = r_{p, N} have ϖ − ϖ_N = 0, while particles scattered to the maximum q = r_{a, N} have ϖ − ϖ_N = π. Particles scattered to an intermediate q = a_N(1 − e²_N) have ϖ − ϖ_N = ±π/2.

In Figure 11, we plot the longitude of periapse relative to that of Neptune (ϖ − ϖ_N) versus periapse q of test particles that were scattered in an integration we performed. The ϖ − ϖ_N is from the first (3000 yr) time step after the particle's scattering. The integration lasted for 1 Myr and included Neptune, with a_N = 30 AU and e_N = 0.3, and 11,600 test particles evenly spaced in semimajor axis from 29 to 34 AU, with e = 0.15. The ϖ − ϖ_N values are well matched by Equation (11).

Zoom In Zoom Out Reset image size

The secular evolution of the particle after its scattering will depend on its semimajor axis, eccentricity, and longitude of pericenter relative to Neptune's (Equation (2)). In particular, from Equation (2), it follows that the initial rate of change of the particle's eccentricity depends on ϖ − ϖ_N. By calculating the particle's total eccentricity $e = \sqrt{h^2+k^2}$, differentiating with respect to time, and evaluating the time derivative at t = 0, we find

$$\begin{equation} \dot{e} (t = 0) = - e_{\rm N}\sin (\varpi - \varpi _{\rm N}) \alpha b_{3/2}^{(2)}(\alpha) \frac{m_{\rm N}}{m_\odot } \frac{n}{4}. \end{equation} \tag{ 12 }$$

An analogous expression follows from Equation (4).

Along a scattering line q, the KBOs do not have random ϖ but ϖ close to the value dictated by Equation (11). Thus, the particles scattered to the minimum q = r_{p, N}, which have ϖ − ϖ_N = 0, and maximum q = r_{a, N}, which have ϖ − ϖ_N = π, are turning over in their secular evolution cycles $(\dot{e} = 0)$. However, particles scattered to an intermediate q = a_N(1 − e²_N), which have ϖ − ϖ_N = π/2, will be decreasing in eccentricity at the maximum rate in the cycle.

Moreover, the particle's free eccentricity also depends on ϖ − ϖ_N:

$$\begin{equation} e_{\rm free}^2 = \bar{e}_{\rm forced}^2+(e(0))^2 - 2 e(0) \bar{e}_{\rm forced} \cos (\varpi -\varpi _{\rm N}), \end{equation} \tag{ 13 }$$

where e(0) is the KBO's eccentricity at t = 0.

Thus, ϖ − ϖ_N sets not only the particle's initial phase in its secular evolution cycle but also the amplitude of its free eccentricity. The particles with the phase to achieve the lowest possible total eccentricity (e_free close to e_forced) may not initially be evolving downward in their cycle. Therefore, in order to calculate the minimum time for a particle to reach the stable region, we consider all values of ϖ − ϖ_N.

Thus, particles with ϖ − ϖ_N = 0 and ϖ − ϖ_N = π have the smallest and largest free eccentricities, respectively, while those with ϖ − ϖ_N = ±π/2 have an intermediate value. Unfortunately for particles trying to reach low eccentricities, the ones with largest free eccentricity (ϖ − ϖ_N = π) are initially going up (i.e., e is increasing).

From numerical integrations (Figure 11), it appears that the maximum deviation in ϖ − ϖ_N as a function of q from Equation (11) is ±40° for conditions relevant to Neptune and the classical region KBOs (due to the fact that θ in Equation (10) is not always exactly ϖ). In setting the initial conditions for the secular evolution of scattered particles, we employ this mapping between ϖ − ϖ_N and the particle's q, including the uncertainty.

Thus, the criterion for delivering the hot classicals developed in Section 2.2 requires the following:

Constraint 2. Neptune's apoapse is large enough (r_{a, N} > 34 AU) so that particles are immediately scattered into the stable region, or Neptune imparts an e_forced large enough relative to e_free so that it is possible for particles with semimajor axes in the range 42.5–47.5 AU to evolve to q > 34 AU in less than τ_{e, N}.

This constraint ensures that it is possible for at least some hot objects in the region from 42.5 to 47.5 AU to be delivered into the region of long-term stability.

When two bodies are near resonance, the secular eccentricity forcing happens on a much faster timescale. This effect has been recognized as the cause of Saturn and Jupiter's fast precession, which Hill (1897) attributed to their period ratios being near 5:2. More recently, Minton & Malhotra (2011) recognized that the 2:1 resonance also contributes to Jupiter and Saturn's fast precession, and Malhotra et al. (1989) identified the classical Uranian satellites' proximity to resonance as the cause of their deviation from their predicted ephemerides.

For the cold objects, the accelerated secular forcing near resonance quickly excites the eccentricities of these objects (as seen in Figure 7), disrupting the confinement of the cold population. Near resonance, the correction term, δg_KBO, to g_KBO is large (Equation (5)). Thus, the secular frequency g_KBO is very high. If resonances overlay the cold classical region at early times, Neptune's eccentricity would have to damp on unrealistically short timescales to fulfill Constraint 1 (Equation (9)). As shown in Figure 8, the objects near resonance remain dynamically disrupted even after Neptune's eccentricity damps. Since the cold objects in the region 42.5 AU < a < 45 AU are confined to low eccentricities (Section 2.2), they cannot have been excited by accelerated secular forcing near resonance while Neptune's eccentricity was high.

Constraint 3. Resonances cannot overlie the region 42.5 AU < a < 45 AU while Neptune's eccentricity is high. This constraint is a special case of Constraint 1 (Equation (9)) and is quantified in Section 5.1.2.

For the hot objects, accelerated secular forcing near resonance can drive down their eccentricities once they have been scattered into the classical region. Figure 12, inspired by Figure 3 of Levison et al. (2008), shows two example integrations of particles beginning at large eccentricities and evolving down to smaller eccentricities. The initial conditions (panel (i)) match those in Levison et al. (2008), Figure 3 (their top left panel). The snapshots in panels (ii) and (iii) are after 1.4 Myr. The integration in panel (ii) includes all four giant planets. Neptune has initial conditions a_N = 30 AU and e_N = 0.2, and the other planets have their modern orbital elements. The integration in panel (iii) is the same except without Uranus. The particles are undergoing secular evolution, as demonstrated by their vertical paths in (a, e) space (a does not change under secular evolution) and their eccentricity oscillations (panel (iv)). If Uranus is present near the 2:1 resonance with an eccentric Neptune, the evolution of the particles is chaotic (panel (ii)), as we will discuss in more detail in Section 4.5. For objects capable of reaching low eccentricities through secular evolution (Equation (13)), proximity to resonance significantly decreases the delivery time. Fast secular evolution near resonances can also assist in capturing objects into resonance, in addition to the previously identified mechanisms of chaotic capture (Levison et al. 2008) and smooth migration (Malhotra 1995).

Zoom In Zoom Out Reset image size

When Neptune's eccentricity is large, the resonances are widened and potentially overlap in what Levison et al. (2008) describe as a "chaotic sea." Levison et al. (2008) argue that this region extends to the 2:1 resonance when Neptune's eccentricity e_N  >  0.15. We have investigated the circumstances for chaos and reached several conclusions, which we will state and then justify.

• 1.  
  Even when widened by Neptune's high eccentricity, the resonances between the 5:3 and 2:1 do not overlap except for particles at high eccentricities.
• 2.  
  Variations in Neptune's semimajor axis on timescales of order a KBO libration time combine with widened resonances to cause a chaotic region.
• 3.  
  The existence and extent of the chaotic region depend on the details of Neptune's interactions with Uranus.
• 4.  
  The potential for chaos does not impose immediate additional constraints beyond those described in the previous subsection (Section 4.2.2). We will show that we can still model just Neptune, taking into account its potential evolution under the influence of the other giant planets.

4.5.1. Just Neptune: No Chaotic Sea

In addition to examining each particle individually, a qualitative way to distinguish between (1) a chaotic sea and (2) particles secularly evolving until they are scattered is to plot the individual paths of a collection of particles through (a, e) space (Figure 13). First, we consider an integration that includes only Neptune (row 1). Each starting at e = 0, the particles move straight upward vertically in (a, e) space until they reach the scattering line (solid black line, e = 1 − r_{a, N}/a), rather than moving both horizontally and vertically as they would in a chaotic sea.

Zoom In Zoom Out Reset image size

In integrations including just Neptune, with or without apsidal precession, the eccentricity of an initially cold particle grows secularly until its orbit crosses Neptune's at e > 1 − r_{a, N}/a. For a particle at 42.5 AU, when a_N = 30 AU and e_N = 0.2, this threshold is 0.15. After reaching this threshold, the particle then undergoes scattering events. Even particles near resonances evolve secularly, with the increased secular frequency defined by Equation (5). These particles appear separated from the coherent excitation of the other cold particles (Figure 7) because they undergo secular evolution so quickly and because the secular oscillation rate depends quite steeply on semimajor axis near resonance. It appears that, for particles that begin at low eccentricities, the resonances do not overlap even when e_N = 0.2. Therefore, the constraints from Section 4.2.2 hold in the just-Neptune case, even when the planet is precessing.

4.5.2. Neptune, Jupiter, and Saturn: No Chaotic Sea

4.5.3. Neptune and Uranus: Chaotic Sea

By modeling the four dynamical processes that result from Neptune's high eccentricity—scattering (Section 4.1), secular forcing (Sections 4.2–4.3), accelerated secular forcing near resonances (Section 4.4), and a chaotic sea (Section 4.5)—we have translated the conservative criteria imposed by the observed eccentricity distributions of the hot and cold classicals (Section 2.2) into the following constraints:

• 1.  
  Neptune's apoapse must be large enough to deliver hot objects to the long-term stable classical region immediately (r_{a, N} > 34 AU) or e_forced must be large enough (relative to e_free) to evolve the particle's e to <0.3 in less than Neptune's eccentricity damping time $\tau _{e_{\rm N}}$.
• 2.  
  The final value for the eccentricities of planetesimals in the region from 42.5 AU < a < 45 AU must be less than e = 0.1: $\sin ({\rm min} (g_{\rm KBO}\tau _{e_{\rm N}},\pi /2)){g_{\rm KBO}^{\prime }}/({\dot{\varpi }_{\rm N}- g_{\rm KBO}}) e_{\rm N}< 0.1$ (or the forced eccentricities must be kept below 0.1 by fast precession).
• 3.  
  Resonances cannot overlie the region 42.5 AU < a < 45 AU while Neptune's eccentricity is high.

Thus, Neptune's eccentricity must be high in order to deliver hot objects to the classical region and yet will disrupt the cold objects quickly unless Neptune's high eccentricity damps quickly or the planet's orbit apsidally precesses quickly. In all cases, mean motion resonances with Neptune cannot overlie the region 42.5 AU < a < 45 AU while Neptune's eccentricity is large.

In Section 5.1.1, we identify which regions of parameter space fulfill Constraint 1 (Equation (9)), preserving the cold classicals below e < 0.1 (as summarized in Section 4.6), without including the effects of orbital resonances or precession. In Section 5.1.2, we incorporate Constraint 3, the effects of orbital resonances. Finally, in Section 5.1.3, we consider the special case of fast precession.

5.1.1. Constraints on Neptune's Eccentricity and Damping Time

5.1.2. Constraints on Neptune's Dynamical History, Including the Effects of Resonances

In regions near orbital resonance with Neptune, KBOs undergo significantly faster secular evolution, as demonstrated in Section 4.2. Here we incorporate the resonance correction terms for the secular excitation times. In Figure 15, analogous to Figure 14, we plot contours of eccentricity damping time as a function of (a_N, e_N) that fulfill the criteria set by Constraint 1 (Equation (9)). We calculate the damping time of Neptune's eccentricity, τ_e, for which 80% of initially cold objects in the region 42.5 AU <a < 45 AU remain below e = 0.1. Note the key difference between the two figures: in Figure 14, τ_e varied smoothly through (a_N, e_N) space, but in Figure 15, there are dark regions where the eccentricity damping time is substantially reduced due to resonances overlying the classical region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The damping times are unrealistically short for parameters of Neptune for which resonances overlie the cold classical region. Neptune is unlikely to have spent substantial time with high eccentricity in these a_N ranges. Ford & Chiang (2007) find that a Neptune-mass planet at 20 AU with an eccentricity of 0.3 undergoes eccentricity damping on $\tau _{e_{\rm N}} =$ 0.6–1.6 Myr (depending on whether Neptune's orbit intersects the planetesimal disk at pericenter/apocenter or at quadrature) if the surface density σ of planetesimals is 1 g/cc. This is roughly the surface density needed to grow planetesimals in the region of Neptune (Kenyon & Luu 1998). Since $\tau _{e_{\rm N}} \propto \sigma ^{-1},$ it is unlikely to be less than 0.1 Myr (10⁵ yr), an order of magnitude faster than the $\tau _{e_{\rm N}}$ calculated by Ford & Chiang (2007). With resonances incorporated, certain regions cannot satisfy Constraint 1 (Equation (9)) in the zero precession case without unphysically low values of $\tau _{e_{\rm N}}$. Constraint 3 is a qualitative statement of this result. Thus, by applying Constraint 1, including resonances, we recover Constraint 3.

Incorporating the resonance correction terms, we plot the maximum $\tau _{e_{\rm N}}$ as a function of the KBO's semimajor axis for several combinations of a_N and e_N (Figure 16). The maximum eccentricity damping time is larger for lower eccentricities, as shown for several illustrative values of e_N in the top panel. The dips where $\tau _{e_{\rm N}}$ approaches 0 correspond to regions of the Kuiper Belt near mean motion resonances where the secular excitation time is very short. The dips are wider when Neptune's eccentricity is higher. When these dips overlie the cold classical region from 42.5 AU  <  a  <  45 AU, the damping time requirement is unphysically short. In the bottom panel, we see how the resonances shift with Neptune's semimajor axis. At a_N = 30 AU, several resonances overlie the cold classical region. When a_N = 28 AU, the resonances are shifted interior and the cold classical region is sandwiched between two resonances, the 9:5 and the 2:1. For a_N = 27.5 AU (not shown), the 2:1 resonance would be on top of the cold classicals. When a_N = 26 AU, all the major resonances are interior to the cold classical region.

Zoom In Zoom Out Reset image size

5.1.3. Special Regimes of Fast Precession

In this subsection, we explore the special case of fast precession. Batygin et al. (2011) first suggested that if Neptune were to precess sufficiently quickly, the cold classical population would remain unexcited because fast precession lowers the forced eccentricity (Equation (5)). In Figure 17, we plot the forced eccentricity as a function of the particle's semimajor axis for a range of precession rates when a_N = 30 AU, e_N = 0.3. Thus, $\tau _{e_{\rm N}}$ in the allowed region can be arbitrarily long if the precession period is sufficiently short. A precession period of 0.9 Myr keeps e_forced < 0.1 in the cold classical region for these parameters of Neptune (compare to the current g₈ precession period of 2 Myr). However, an even faster apsidal precession period for Neptune may be required if not only secular evolution but scattering and/or chaos excite the cold classicals.

Zoom In Zoom Out Reset image size

Unlike in the case of fast damping (Section 4.1.1), in the case of fast precession, scattering cold KBOs by Neptune may impose an important constraint. In the case of fast precession (and slow damping), even though the KBOs eventually evolve to e_free once Neptune's eccentricity damps, the KBO can reach a maximum value of e_max = 2e_forced while e_N is still high. Thus, the KBO is vulnerable to being scattered. For example, for a_N = 30 AU, e_N = 0.3, a KBO at 42.5 AU can be scattered if it reaches e > 0.08. Thus, the KBO needs e_forced < (1/2)0.08 = 0.04 to avoid being scattered, requiring a fast precession period of 0.4 Myr for Neptune.

Chaos (Figure 13) can still excite the cold classicals even if their forced eccentricity is low. In this case, resonances cannot overlie the cold classical region if Neptune's semimajor axis oscillates strongly due to interactions with Uranus (Section 4.5). Therefore, if the chaotic sea is present, the dark regions in Figure 15 are still forbidden.

Figure 18 demonstrates how fast precession caused by Neptune's interaction with the other planets is unable to keep the cold objects at low eccentricities in resonance regions. We performed two integrations; in each of which Neptune has a_N = 30.06 and e_N = 0.2, and 600 test particles begin with e = 0. The first integration (top panel) includes the other three giant planets, on their current orbits; they cause Neptune to undergo apsidal precession with a period of 1.2 Myr. In the second integration, Neptune is forced to precess at this rate without any other planets included (see Section 3.3 for details on the implementation). We plot the maximum eccentricity reached by each particle. In the top panel, chaos has dynamically disrupted particles that were preserved at low eccentricities by fast precession in the bottom panel (i.e., in the region interior to 45 AU).

Zoom In Zoom Out Reset image size

Unfortunately, the most obvious configuration that causes fast precession—Neptune and Uranus near their 2:1 resonance—also causes large oscillations in Neptune's semimajor axis. The larger the oscillations in Neptune's semimajor axis, the wider the chaotic sea. Although Batygin et al. (2011) suggest that a massive planetesimal disk could contribute to fast precession, this precession rate would be roughly (m_disk/m_☉)n_N, where n_N is Neptune's orbital frequency. Thus, a disk mass of roughly 60 Earth masses would produce a precession period of 0.9 Myr. Therefore, it is unlikely that the planetesimal disk alone could produce sufficiently fast precession. When Neptune is in the region from 28 to 29 AU—and the classical region is sandwiched between the 9:5 and 2:1 resonances—oscillations in Neptune's semimajor axis due to resonant interactions with Uranus could cause a chaotic sea that extends into this "sandwich" region. Thus, fast precession may fail for parameters for which fast damping has the possibility of working.

We will present a contour map for Neptune's precession rate, analogous to Figure 15, in Section 5.3.

Now we consider which parameters of Neptune will allow the transport of the hot objects from the inner disk into the classical region. We consider the full range of (a, e) into which Neptune can scatter objects from the inner disk and the resulting secular oscillations of the objects once they reach the classical region. We calculate a minimum eccentricity damping time as a function of (a_N, e_N) by requiring that particles can reach the stability region defined by q > 34 AU (Constraint 2). The "minimum time" (contoured) is the time at which it is possible for objects in 50% of semimajor axis intervals, Δa = 0.063 AU, to reach this stable region.

Figure 19, a map of constraints on Neptune's parameters that deliver the hot objects, illustrates that there are three outcomes for the transport of hot objects.

• 1.  
  White region. If r_{a, N} is sufficiently large, some objects are scattered into the stable region. This criterion can be quantified as r_{a, N} = a_N(1 + e_N) > 34 AU.
• 2.  
  Black region with red stripes. If r_{a, N} is too small, the KBOs are scattered into a region that is not stable over 4 Gyr and their forced eccentricity is too small (relative to e_free) to allow them to secularly evolve down into the stable region. This second outcome can be understood from Equation (2). If the forced eccentricity is small, the free eccentricity will be close to the particle's initial (large) eccentricity after scattering. Thus, the free and forced eccentricity vectors cannot destructively cancel to achieve a low total eccentricity.
• 3.  
  Middle region with gray-scale timescale contours. For intermediate values of r_{a, N}, the particles are scattered above the stable region (i.e., to eccentricities too large for stability) but can secularly evolve down into the stable region if Neptune's eccentricity damping timescale is long enough. The time contours refer to the minimum time (defined above) for particles to secularly evolve to the stable region. We note that, because the analytical model is only first order in the KBO's eccentricity, the timescales in the intermediate region may be in error by up to a factor of two (see Appendix 4.2.2 for discussion). However, the particles most relevant for this scenario—those rapidly declining in eccentricity—are matched with substantially smaller error. We also note that this regime (Outcome 3) is a narrow region of parameter space and that most parameters for Neptune fall within one of the two regimes above.

Zoom In Zoom Out Reset image size

The minimum values for (a_N, e_N) in the intermediate range of r_{a, N} have a messy analytical solution. The minimum eccentricity a particle can achieve is

$$\begin{equation} e_{\rm min} = \big(e_{\rm free}^2+e_{\rm forced}^2-2e_{\rm free}e_{\rm forced}\big)^{1/2}. \end{equation} \tag{ 14 }$$

Combining with Equations (11) and (13), for a given a_N one can solve for the cutoff e_N, below which particles cannot evolve into the long-term stable region.

5.2.1. Effects of Resonances and Chaos

5.2.2. Effects of Precession

5.2.3. Scattering Efficiency of the Hot Objects

We placed constraints on parameters of Neptune that preserve the confinement of cold objects to low eccentricities (Section 5.1) or allow transportation of the hot objects from the inner disk into the classical region (Section 5.2). These constraints are useful separately, and they confirm in situ formation as a feasible origin for the cold objects and transport from the inner disk as a feasible origin for the hot objects. In this section, we investigate which parameters permit the combination of these two origins, producing both a cold and a hot population. This may be possible if there is overlap between the parameter space that preserves the cold classicals (Section 5.1) and the parameter space that transports the hot classicals (Section 5.2).

In Figure 20, we combine the constraints from the hot and cold classicals. First, as a function of (a_N, e_N) we plot the contours for the maximum eccentricity damping time necessary to keep the cold classicals at e < 0.1 in the region from 42.5 to 45 AU. In the white region, Neptune's eccentricity is low enough such that the final eccentricities of the cold objects will be below 0.1, no matter how long the damping timescale. Then we overplot red, diagonal stripes in the region where Neptune cannot deliver the hot classicals, neither by direct scattering nor post-scattering secular evolution. In creating this region, we took into account the sliver of parameter space (Figure 6, contoured region) for which Neptune cannot scatter objects directly into the region q > 34 AU but for which objects can evolve into the stable region through secular evolution. We compare the time required for the hot objects to reach the long-term stable region to the maximum eccentricity damping time to preserve the cold classicals. In most cases, Neptune's eccentricity must damp before any hot objects reach the long-term stable region. From all these considerations, there is only a small region of parameter space where Neptune can deliver the hot classicals (no red, diagonal stripes), while fast secular evolution near resonances does not quickly excite the cold classicals (light regions).

Zoom In Zoom Out Reset image size

The allowed region in (a_N, e_N) space, in which the hot objects can be transported without disrupting the cold population, is bounded by a strip extending from (23 AU, 0.4) to (30 AU, 0.10) on the left and (27.5 AU, 0.4) to (30 AU, 0.2) on the right. To the top right of this strip, the series of resonances extending from the 5:2 to the 9:5, widened by Neptune's high eccentricity, are on top of the cold classical region. Within this bounding strip is a forbidden region near a_N = 27.5 AU, at which the 2:1 resonance overlies the cold classical region.

Based on the considerations above, lower limits can be placed on Neptune's migration timescale while Neptune's eccentricity is high. During Neptune's period of high eccentricity, Neptune's migration timescale must be slow enough to keep Neptune in a region that will not quickly excite the cold classicals. Neptune should not spend substantial time with high eccentricity near a_N = 27.5 AU or a_N = 30 AU, lest resonances disrupt the classical region. Thus, Neptune is constrained to migrate no more than a few AU during Neptune's eccentricity damping time. Because of the discrete ranges of consistent semimajor axes, when resonances are accounted for, a "damp first and then migrate" scenario is consistent with preserving the cold classicals, while a "migrate first and then damp" scenario is not.

In Figure 21, we plot analogous constraints for the special case of fast precession. The contours represent the precession rate of Neptune necessary to keep e_forced sufficiently low (Equation (7)). In all cases, the KBOs' forced eccentricities must be below e_forced < 0.1 (Constraint 1). We impose an additional constraint to avoid the scattering of cold particles as they reach their maximum eccentricities of e_max = 2e_forced. We require e_max to be below the scattering line at the inner edge of the cold classical region; therefore, e_forced < (1/2)(1 − q_{a, N}/42.5). For large values of Neptune's apoapse q_{a, N}, this constraint is stricter than e_forced < 0.1. This constraint may, in fact, be too strict because if a cold classical KBO is scattered by Neptune, it is more likely to be scattered out of the classical region altogether than to end up between 42.5 AU < a < 45 AU at an eccentricity too large to be consistent with the observations but small enough to survive over 4 Gyr. We leave a detailed investigation into the role of scattering in the case of fast precession for future work. The red, diagonally striped region is where hot objects cannot reach the long-term stable region.

Zoom In Zoom Out Reset image size

Next, we consider the constraints for delivering the hot classicals. For each (a_N, e_N), we calculate whether objects can be directly scattered into the region q > 34 AU. If not, we assume that Neptune imparts a forced eccentricity small enough to be consistent with the constraints above for keeping the cold classicals unexcited (i.e., e_forced = 0.1, or, if smaller, e_forced < (1/2)(1 − q_{a, N}/42.5)) and determine whether the hot object can reach the stable region through secular evolution, using Equation (14). If Neptune cannot directly scatter objects into the stable region and if objects secularly evolving cannot reach q > 34 AU in more than 50% of semimajor axis intervals, (a_N, e_N) is not consistent with delivering the hot objects. In the cases where fast precession is necessary to keep the cold classicals confined to low eccentricities, hot classicals scattered into the region have the same low forced eccentricities; therefore, they do not experience significant secular evolution down into the classical region but remain at their post-scattering eccentricities.

Finally, in blue vertical stripes, we overplot ranges of a_N for which the centers of one or more resonances of fourth order or lower lie in the cold classical region from 42.5 to 45 AU. These a_N are not necessarily forbidden, but chaos can occur if Neptune's eccentricity is high and interactions with Uranus cause oscillations in a_N.

Though we are applying constraints that are lenient and conservative in ruling out regions of parameter space, a substantial fraction of parameter space is ruled out. The following scenarios have not been ruled out and may work, depending on the details of the interactions between Uranus and Neptune. In both of these scenarios, resonances cannot overlie the region of 42.5 AU < a < 45 AU.

• 1.  
  Short τ_e and larger_{p, N}. An eccentric Neptune transports objects into the classical region until its eccentricity damps, which occurs quickly enough that the cold classicals are not excited (Figure 20).
• 2.  
  Fast $\dot{\varpi }_{\rm N}$ and larger_{p, N}. An eccentric, quickly precessing Neptune transports objects into the classical region and its quick apsidal precession keeps the forced eccentricity of the particles low, preserving the cold classicals (Figure 21).

In principle, an intermediate r_{p, N} could deliver hot objects to high eccentricities (e > 0.3) in the classical region but allow them to evolve to lower eccentricities (e < 0.3) before the cold objects are excited. In practice, we found that the timescales are not compatible: if the eccentricity damps quickly enough to preserve the cold classicals from secular excitation, there is not time for hot objects to secularly evolve down into the Belt. However, for parameters of Neptune for which hot objects are scattered directly in the stable region, secular evolution can allow the objects to reach even lower eccentricities, especially the objects undergoing fast secular evolution near resonances.

Therefore, in practice, the consistent dynamical histories are ones in which Neptune has a large enough apoapse to transport hot objects immediately into the stable region, with its eccentricity damping quickly enough or precessing quickly enough so that the cold classicals remain at low eccentricities consistent with the observations.

The goal of this paper is to determine which parameters for Neptune allow the planet to deliver the hot classicals from the inner disk into the classical region without dynamically disrupting the in situ cold classicals. Throughout its dynamical history, Neptune must satisfy the constraints presented here to avoid excessively exciting the cold classicals (Section 5.1). At some point, Neptune must also spend time in a region of parameter space where it can also deliver the hot classicals (Section 5.3). We clarify that Neptune did not necessarily form at the location where hot classical delivery takes place or arrive there after undergoing a single, instantaneous scattering. Before and after hot classical delivery, Neptune can potentially spend time in any region of parameter space as long it obeys our constraints against not excessively exciting the cold classicals. The constraints developed here, which can serve as a "road map" for Neptune's path through parameter space, hold for realistic models that include multiple scatterings and have a straightforward interpretation.

In the context of the Nice model, Neptune may have undergone a series of scatterings, spending time at a variety of spots in (a_N, e_N) space. In each of these spots, Neptune must obey the constraints we place to avoid disrupting the cold classicals. Perhaps the scattering occurs quickly compared to the excitation time for the cold classicals. If not, Neptune could pre-excite the cold classicals (but only to below the observational limit) before it reaches a region where it can deliver the hot classicals. In imposing our conservative constraints, we are assuming that the cold classicals begin with e = 0, but if they are pre-excited, the constraints will be stricter.

Another possibility is that Neptune may spend a long time at a location where it creates a very small region of stability in the classical region, clearing out most of the cold classicals. This is a potential solution to the mass efficiency problem, which we will discuss in Section 6. However, as shown in Section 3.1.1, Neptune cannot deliver the hot classicals in this regime. Therefore, after this period has ended, when Neptune is delivering the hot classicals, the constraints we will place on not exciting the cold classicals will hold.

The scattering(s) Neptune undergoes are quick changes in its orbit. After its period of high eccentricity ends—or during a temporary period of low eccentricity—Neptune can undergo slow evolution, including slow migration and slow damping of its (now small) eccentricity. The KBOs will maintain their free eccentricities throughout this slow evolution. If Neptune's eccentricity were excited very gradually, on a timescale much longer than the secular evolution times of the cold classicals, the cold classicals would keep their initial low free eccentricities. However, we expect the excitation of Neptune's eccentricity via scatterings to take place on a timescale shorter than millions of years (e.g., Thommes et al. 1999).

We have ruled out much of parameter space with the constraints that Neptune cannot excite the cold objects above e = 0.1 and must be able to deliver at least a few hot objects to q > 34 AU in the region from 42 to 47.5 AU. We note that we are "ruling out" parts of parameter space where Neptune cannot deliver the hot objects without disrupting the cold, not "ruling out" that the planet can ever spend time there (see above). In Figure 22, we look for greater consistency with the observations (Figure 3): a forced eccentricity less than 0.075 for the cold objects and q_N > 36 AU (meaning that a hot object could be scattered to an eccentricity as low as 0.24 at 47.5 AU). The parameter space shrinks where Neptune can deliver the hot objects without disrupting the cold (i.e., light gray regions with no red, diagonal constraints). Overplotted on this figure is an example (arrows) of Neptune's path through parameter space. In this example, Neptune first undergoes multiple scatterings, spending a short enough time at each (a_N, e_N) to avoid exciting the cold classicals. Then it reaches 28 AU with an eccentricity of 0.3; here it delivers the hot classicals, as its eccentricity damps quickly enough to avoid exciting the cold classicals. Then it migrates at a low eccentricity to its current location.

Zoom In Zoom Out Reset image size

On this figure, we overplot some parameters as symbols. For the triangles and circles, we will show example integrations in Section 5.5. The pluses are parameters taken from the literature. The plus at a_N = 20, e_N = 0.02 marks the initial condition of Hahn & Malhotra (2005), from which Neptune undergoes migration to its current location. Neptune remains in a region that is white (never excites the cold classicals) but also red, diagonally striped (does not deliver the hot classicals); thus, this simulation maintained the low eccentricities of the cold classicals but did not deliver the hot classicals. The plus at a_N = 28.9, e_N = 0.3 indicates the initial condition for Levison et al. (2008), Run B. In this part of parameter space, Neptune can deliver the hot classicals, which were indeed produced by their simulation. However, Neptune's eccentricity damped on a timescale of τ_e = 2 Myr, too slow to avoid exciting an in situ cold population.

Example integrations illustrating the constraints we have derived are shown in Figures 23 and 24. Each integration lasts for 1.6 Myr unless otherwise noted (see Section 3.3 for a general description of the integrations). In addition to the 600 test particles in the region from 40 to 60 AU, we add 59,600 test particles in the region from 18 to 38 AU, each with initial e = i = 0, representing the inner planetesimal disk from which the hot classicals are scattered.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The parameters for the four integrations shown in Figure 23 are plotted in the parameter space map in Figure 22 as circles and correspond to regions of parameter space where Neptune cannot both deliver the hot objects and keep the cold objects at low eccentricities. In row 1, Neptune is at a_N = 24, e_N = 0.02 and undergoes no eccentricity damping. The cold classicals remain confined to low eccentricities, but the hot classicals are delivered to high eccentricities and cannot secularly evolve down to lower eccentricities. For the same reason (i.e., that Neptune's apoapse is too small), the hot classicals are also not delivered in row 2 (a_N = 22, e_N = 0.35). The damping timescale in this integration is longer (2 Myr) than half a secular evolution time, so the eccentricities of the cold objects converge to the forced eccentricities, which are above 0.1. This integration lasts 6 Myr. In the third row (a_N = 27.5, e_N = 0.35, τ_e = 0.2 Myr) and fourth row (a_N = 30.06, e_N = 0.35, τ_e = 0.2 Myr), the hot classicals are delivered but the cold classicals are excited by fast secular evolution near resonances.

The parameters for the two integrations shown in Figure 24 are plotted in the parameter space map in Figure 22 as triangles and correspond to a set of parameters in each of the two viable regions of parameter space. In both cases, Neptune can deliver the hot classicals. Moreover, since it is in a region of parameter space where no resonances overlie the region of 42.5 AU <a < 45 AU, and since its eccentricity damps quickly enough to obey our constraints, the cold classicals remain at e < 0.1. In the top panel, we plot the observed objects from Figure 3 for comparison. In the middle panel, the 2:1 resonance is interior to the cold classical region. In the bottom panel, the cold classicals are sandwiched between two regions where the cold classicals are excited by fast secular evolution near resonance.

We now consider the free, or proper, elements of the observed KBOs, which have been computed for a subset of KBOs by Knežević & Milani (2000), Knezevic et al. (2002), and Knežević & Milani (2003). The free elements precess about the forced values, which are set by the current configuration of the giant planets, and thus provide a better window to the history of the solar system than the instantaneous orbital elements. In Figure 25, we plot the proper element eccentricities and inclinations of observed KBOs on top of the survival maps of Lykawka & Mukai (2005), which are formulated in terms of instantaneous eccentricity and inclination. Qualitatively, we see the same features as in Figures 3 and 4: the cold classicals (red squares) are confined below e < 0.1 in the region from 42.5 to 45 AU. Throughout the region, the hot objects (blue triangles) occupy the upper portion of the stability region.

Zoom In Zoom Out Reset image size

We repeat the same statistical tests (Tables 5 and 6) as above for the subset of observed KBOs that have computed proper elements. We use the proper inclinations to classify the objects with i < 2 as cold. For the second test, for which the sample is probabilistically selected, we use the inclination distribution from Volk & Malhotra (2011), which considers the inclinations with respect to the invariable plane. The results are consistent with those from the instantaneous orbital elements above. We caution that we are comparing the proper eccentricities of observed objects to a stability map of instantaneous eccentricities. This comparison should be repeated when a stability map formulated in terms of proper elements becomes available.

In Section 4.2.2, we relegated to the appendix the derivation of several additional terms in the modified secular theory. In Section B.1.1, we derive the factors proportional to e²_N that appear in the extra factors used in Equations (5) and (8) as coefficients to e²_N. In Section B.1.2, we follow Malhotra et al. (1989) to derive resonance correction terms.

B.1.1. High-order Eccentricity Terms

The basic secular theory includes only the lowest order eccentricity terms. However, when Neptune's eccentricity is high, terms containing e²_N are no longer negligible. Therefore, g_KBO and $\bar{e}_{\rm forced}$ must be modified. Here we define the extra terms and factors used in Equation (5), which come from additional terms in the disturbing function (see chap. 7 of Murray & Dermott 2000, for a standard derivation). The disturbing function has the form, up to second order, of

$$\begin{eqnarray} R &=& n^2 a^2 \frac{m_{\rm N}}{m_\odot } \big[e^2 \big(f_2 + f_5 e_{\rm N}^2 + f_6 e^2\big)\nonumber \\ && + \cos (\varpi -\varpi _{\rm N}) e e_{\rm N}\big(f_{10} + f_{11} e_{\rm N}^2 + f{12} e^2\big) \nonumber \\ &&+ \cos (2(\varpi -\varpi _{\rm N})) e^2 e_{\rm N}^2 f_{17}\big]. \end{eqnarray} \tag{ B1 }$$

Fully incorporating all e and e_N to second order would modify the functional form of the secular theory. However, if we treat e_N as a constant, we can modify the f₂ term in the secular forcing frequency g_KBO to f₂ + f₅e²_N (Equation (5)) and the f₁₀ term in the forced eccentricity to f₁₀ + f₁₁e²_N (Equation (8)). Because the form of the secular evolution Equation (2) is derived by differentiating R with respect to the particle's h and k, treating e_N as a constant does not modify the form of the secular evolution equations but simply adds extra correction factors. The f factors are defined in Appendix B of Murray & Dermott (2000). The success and necessity of these extra terms are illustrated in Figure 26.

Zoom In Zoom Out Reset image size

B.1.2. Modification Near Resonances

Proximity to resonance modifies the secular frequencies. Following Malhotra et al. (1989), we define the resonance correction terms C and in Equations (5) and (6). Malhotra et al. (1989) developed equations for a pair of moons near first-order resonance. Instead of a pair of moons, we treat a massless KBO and Neptune, and we include resonances above the first order.

Resonances add additional important terms to the disturbing function, R. The factor C in Equation (6) is proportional to the coefficient of the direct part of the disturbing function, with argument (j + x)λ − jλ_N − xϖ_N, where x is the order of the resonance. Since the coefficients of this part of the disturbing function R are proportional to e^x_N, and since $\dot{h} \propto (\frac{\partial R}{\partial h})$ and $\dot{k} \propto \frac{\partial R}{\partial k}$, a factor of xe^{x − 1}_N comes in (Equation (6)), which was not explicitly included in Malhotra et al. (1989) because they treated only the x = 1 case. The factor C for each order x is tabulated in Table 7. These coefficients were taken from the expansion of the disturbing function in Appendix B of Murray & Dermott (2000).

Table 7. Coefficients for Equation (6)

x	C
1	$\frac{1}{2} \big(-2(j+x) -\alpha \frac{d}{d\alpha }\big)b_{1/2}^{(j+x)}(\alpha)$
2	$\frac{1}{8} \big((-5(j+x)+4(j+x)^2) + (-2+4(j+x))\alpha \frac{d}{d\alpha }\,{+}\,\alpha ^2 \frac{d^2}{d\alpha ^2}\big) b_{1/2}^{(j+x)}(\alpha)$
3	$\frac{1}{48} \big((-26(j+x)+30(j+x)^2-8(j+x)^3)+ (-9+27(j+x)-12(j+x)^2)\alpha \frac{d}{d\alpha }\,{+}\,(6-6(j+x))\alpha ^2 \frac{d^2}{d\alpha ^2} - \alpha ^3 \frac{d^3}{d\alpha ^3}\big)b_{1/2}^{(j+x)}(\alpha)$
4	$\frac{1}{384} \big((-206(j+x)+283(j+x)^2-120(j+x)^3+16(j+x)^4)+ (-64+236(j+x)\,{-}\,168(j+x)^2+32(j+x)^3)\alpha \frac{d}{d\alpha }$
	$+(48-78(j+x)+24(j+x)^2)\alpha ^2 \frac{d^2}{d\alpha ^2} + (-12+8(j+x)) \alpha ^3 \frac{d^3}{d\alpha ^3}\,{+}\, \alpha ^4 \frac{d^4}{d\alpha ^4}\big)b_{1/2}^{(j+x)}(\alpha)$

Download table as:  ASCIITypeset image

The factor in Equation (5) depends on the proximity to resonance. Extending Malhotra et al. (1989) to resonances of arbitrary order, we obtain

$$\begin{eqnarray} \omega &=& j n_{\rm N} - (j+x) n \left(1+\frac{m_{\rm N}}{m_\odot }\left(1+\alpha \frac{d}{d\alpha }\right)b_{1/2}^{(0)}\right)\nonumber \\ \epsilon &=& \frac{3}{2} \frac{m_{\rm N}}{m_\odot } (j+x)^2 \frac{1+\frac{m_{\rm N}}{m_\odot } \alpha \big(1+\frac{7}{3}\alpha \frac{d}{d\alpha }+\frac{2}{3}\alpha ^2 \frac{d^2}{d\alpha ^2}\big)b_{1/2}^{(0)} }{(\omega /n)^2}, \nonumber \\ \end{eqnarray} \tag{ B2 }$$

where n_N is the mean motion of Neptune and n is the mean motion of the particle.

The n:1 resonances have indirect terms not explicitly included in Malhotra et al. (1989) (R. Malhotra 2011, private communication). However, the only relevant n:1 resonance in the region of the Kuiper Belt we are studying is the 2:1 resonance, and its indirect terms result in expressions that, when incorporated above, are directly proportional to the mass of the KBO and thus assumed to be negligible.

In the case of multiple planets, the forced eccentricity of a small body on an external orbit is given by (Murray & Dermott 2000)

$$\begin{eqnarray} h_0&=& {-} \sum _{i=1}^N \frac{\nu _i}{A-g_i} \sin (g_i t + \beta _i) \nonumber \\ k_0&=& {-} \sum _{i=1}^N \frac{\nu _i}{A-g_i} \cos (g_i t + \beta _i), \end{eqnarray} \tag{ B3 }$$

where

$$\begin{eqnarray} \nu _i &=& \sum _{j=1}^N A_j e_{ji} \nonumber \\ A_j &=& {-}n \frac{1}{4} \frac{m_j}{m_\odot } \alpha _jb_{3/2}^{(2)}(\alpha _j) \\ A &=& \sum _{j=1}^N n \frac{1}{4}\frac{m_j}{m_c} \alpha _j b_{3/2}^{(1)}(\alpha _j),\nonumber \end{eqnarray} \tag{ B4 }$$

where the particle's forced eccentricity $e_{\rm forced} = \sqrt{\vphantom{A^A}\smash{\hbox{${h_0^2 + k_0^2}$}}}$ and g_i and e_ji are the eigenfrequencies and eigenvector components of the planetary system. We compare Equation (B3) with Equation (7). For KBOs in the classical region, Neptune's orbit dominates A, the precession rate of the particle's free eccentricity. The precession rate of a particle's free eccentricity due to Neptune alone agrees with the four-planet case to within 30%. Thus, A ≈ g_KBO. Neptune dominates the ν_i term, and thus the four-planet secular theory reduces to the single-planet secular theory with an extra $g_i = \dot{\varpi }_{\rm N}$ term for Neptune's precession. Therefore, the four-planet Equation (B3) reduces to Equation (7). The quantity ν_i corresponds to g'_KBOe_N(t), g_i corresponds to $\dot{\varpi }_{\rm N}$, and A corresponds to g_KBO. The extra sin term in Equation (7) is an empirical factor to account for eccentricity damping and is discussed in detail in the main text. This conclusion is consistent with the result of Chiang & Choi (2008) that the current forced eccentricities of the KBOs are largely determined by Neptune's orbit.