NEPTUNE ON TIPTOES: DYNAMICAL HISTORIES THAT PRESERVE THE COLD CLASSICAL KUIPER BELT

The solar system is often used as a case study for the formation of planetary systems from proto-planetary disks. The current configuration of Kuiper Belt objects (KBOs) provides a map for how the dynamical evolution of the giant planets in our solar system sculpted the disk of planetesimals. Therefore, it is possible to use the orbital properties of the Kuiper Belt to constrain how the orbits of the giant planets evolved in the early solar system, particularly for Neptune, the primary sculptor of the Kuiper Belt. Models employing N-body integrations to trace the effects of the giant planets' migration and orbital eccentricity evolution on the planetesimal disk have enjoyed substantial success in reproducing the dynamical populations of KBOs observed today (e.g., Malhotra 1993, 1995; Hahn & Malhotra 1999, 2005; Gomes 2003; Levison et al. 2008; Morbidelli et al. 2008). These populations include objects near orbital resonance with Neptune, objects scattering off Neptune, and "classical" objects decoupled from Neptune. (See Gladman et al. 2008 for definitions of the dynamical classes.) Yet substantial discrepancies still exist between the simulations and observations, particularly for the classical population. The bias-corrected inclination distribution of observed classical KBOs is bimodal (Brown 2001; Gulbis et al. 2010; Volk & Malhotra 2011). The low-inclination, dynamically "cold" component and the high-inclination, dynamically "hot" component also have distinct physical properties, including colors (Tegler & Romanishin 2000; Trujillo & Brown 2003; Peixinho et al. 2008), sizes (Levison & Stern 2001; Fraser et al. 2010), albedos (Brucker et al. 2009), and binary fractions (Stephens & Noll 2006; Noll et al. 2008). To date no simulations have been able to produce both the high- and low-inclination classical objects while qualitatively matching their observed eccentricity distribution.

Several theories of the dynamical history of the giant planets in our solar system have been proposed, inspired by the dynamical populations within the Kuiper Belt. One of the most widely accepted models, the Nice Model, stems from a postulated large-scale instability in the early solar system (e.g., Thommes et al. 1999; Gomes 2003; Morbidelli et al. 2008). Inspired by the Nice Model, Levison et al. (2008) propose a scenario in which Neptune is scattered outward onto an eccentric orbit to near its current location from a formation location closer to the Sun (e.g., Thommes et al. 1999, 2002). In this scenario, Neptune's eccentricity subsequently damps due to dynamical friction with a remnant disk of planetesimals. This model reproduces the observed resonant population, scattered population, and the hot classical population, but, we argue, does not satisfactorily match the low eccentricities of the observed population of dynamically cold objects in the classical region. Furthermore, it does not produce a sufficient number of high-inclination objects. Other incarnations of the Nice Model (e.g., Gomes 2003; Morbidelli et al. 2008) include an in situ population of cold objects, but, over the course of Neptune's evolution, these objects become excited to higher eccentricities. Whether or not the particular history described by the Nice Model is correct, the prevalence of large eccentricities among extrasolar giant planets suggests that large-scale orbital excitation is common during the formation of planetary systems.

In contrast to the upheaval of the Nice Model, another model (Malhotra 1993) proposes a period of extensive, smooth migration of the giant planets. Planetesimal-driven migration (Fernandez & Ip 1984) is likely important in shaping the architecture of many planetary systems. Malhotra (1993, 1995) and Hahn & Malhotra (1999) demonstrated that the migration of the giant planets in our solar system on flat, circular orbits would have perturbed the disk of planetesimals, scattering some to more eccentric orbits and capturing others into resonance. They also found that migration results in a large number of planetesimals in the location known today as the "Scattered Disk," which encompasses those objects that have had gravitational interactions with one or more of the giant planets and have high orbital eccentricities and inclinations. The key inconsistency between this model and current observations is that it cannot produce, without additional processes such as stochasticity (e.g., Levison & Morbidelli 2003; Murray-Clay & Chiang 2006), both the cold and hot classical populations from a single set of initial disk conditions or account for the differences in their physical properties.

Because detailed N-body simulations are computationally expensive, previous works have investigated a limited number of solar system planetary histories. Currently, no comprehensive parameter study has been done exploring the evolution of Neptune's orbit given constraints from the sculpting of the Kuiper Belt. In addition, although planet–planet scattering often produces large mutual inclinations (Chatterjee et al. 2011), as observed in the Upsilon Andromedae system (McArthur et al. 2010), no serious, detailed treatment has included the possibility that Neptune underwent a period of high orbital inclination, a conceivable outcome of the orbital instability in the early solar system. Here we consider a general model that can encompass the two detailed models described above—the Nice-Model-inspired scenario in which Neptune is scattered to a high eccentricity (Levison et al. 2008) and the scenario of extensive migration of Neptune on a low-eccentricity, low-inclination orbit (e.g., Malhotra 1995)—as well as potential scenarios in which Neptune undergoes a period of high inclination. In this generalized model, Neptune undergoes some combination of migration and/or evolution of its eccentricity and/or inclination. Here we take a step toward understanding the qualitative differences in the dynamics of the Kuiper Belt generated by a wide range of planetary histories. In a series of three papers, we perform a parameter study of the effects on a disk of planetesimals of the migration and the eccentricity and inclination evolution of Neptune in the early solar system. No component of this parameterization is new; rather, our goal is to comprehensively consider all possible parameters for Neptune's history.

In this first paper, we introduce our method for computationally and analytically modeling Neptune's orbital evolution and its effects on the planetesimal population. Then we demonstrate several concepts that will allow us to thoroughly explore the parameter space of Neptune's dynamical history:

• 1.  
  If Neptune undergoes a period of elevated orbital eccentricity, it will secularly excite the eccentricities and inclinations of an in situ planetesimal population. We can analytically model this secular excitation using a simple expression.
• 2.  
  As Neptune's eccentricity and inclination damp, the planetesimals evolve to final eccentricities and inclinations that depend on Neptune's initial eccentricity and inclination and damping timescales.
• 3.  
  The effects of Neptune's eccentricity and inclination evolution on the planetesimals can be treated separately to first order.
• 4.  
  The migration rate sets Neptune's effective location for the secular evolution of the planetesimals.

Having demonstrated these points, we can place robust constraints on Neptune's semimajor axis, eccentricity, and inclination during its late evolution (after any period of planet–planet scattering), and on its migration, eccentricity damping, and inclination damping rates, identifying which parameters are consistent with maintaining the low eccentricities and inclinations of the cold classicals. In the second paper (Dawson & Murray-Clay 2012), we place even stronger constraints by incorporating additional effects, including the effects of other planets on Neptune's orbital evolution and the effects of proximity to mean-motion resonance with Neptune on the secular excitation of the planetesimals, as well as additional constraints from the hot classical population. A third paper (R. I. Dawson & R. A. Murray-Clay, in preparation) will focus on the inclinations of the classicals.

In Section 2, we present our parameterization of Neptune's orbital evolution. In Section 3, we discuss the observational constraints that today's cold classical KBOs place on past sculpting of the planetesimal disk and describe our computational and analytical models of the excitation of the planetesimal disk by an inclined and eccentric Neptune. We present the results of the parameter study in Section 4. In Section 5, we present our conclusions and describe how the companion papers will expand upon this work.

As a first step toward a comprehensive study of the impact of a wide range of dynamical histories of the outer solar system on the Kuiper Belt, we thoroughly explore the parameter space of a general model for Neptune's dynamical history. This model encompasses specific, previously proposed solar system history models (e.g., Malhotra 1995; Levison et al. 2008). In our parameterization, Neptune is instantaneously scattered to an initial location with an initial eccentricity and inclination. Subsequently, it undergoes planetesimal-driven migration and its eccentricity and inclination are damped by dynamical friction from the disk of planetesimals. Resonant relaxation may also contribute significantly to the planet's eccentricity and inclination damping if the Toomre parameter of the planetesimal disk is large or if the random velocities of the planetesimals are large (Tremaine 1998). The parameters of this model are Neptune's "initial" (defined below) semimajor axis, eccentricity and inclination, the planet's migration rate, and the timescales for Neptune's eccentricity and inclination damping. In this model, we include the effects of only one planet (Neptune), an approach we justify briefly in Section 3.3 and more thoroughly in Dawson & Murray-Clay (2012).

2.1. Parameters

2.2. Computational Model

Direct computational modeling of the effect of planetesimals on Neptune's orbit would be computationally expensive, so instead we apply fictitious forces (Appendix) to evolve Neptune's semimajor axis a_N, eccentricity e_N, and inclination i_N, with any specified functional form. Following Malhotra (1993) and Levison et al. (2008), we use the functional forms:

$$\begin{eqnarray} e_N &=& e_0 \exp {(-t/\tau _e)} i_N = i_0 \exp {(-t/\tau _i)} \nonumber \\ a_N &=& a_f + (a_0 - a_f) \exp {(-t/\tau _a),} \end{eqnarray} \tag{ 1 }$$

where a₀ is the initial semimajor axis of Neptune, a_f = 30 AU is the final semimajor axis, and τ_e, τ_i, and τ_a are the eccentricity damping timescale, inclination damping timescale, and migration timescale, respectively. Our results do not depend on the specific form of Equation (1). As we will demonstrate in Section 4, sometimes the instantaneous rate of change of the variables (${\dot{a}}/{a}, {\dot{e}}/{e}, {\dot{i}}/{i}$) is most relevant, while in other cases the total evolution matters most. We have verified these statements with integrations (not shown) using an alternative migration form $({\dot{a}}/{a}) \propto ({\dot{e}}/{e}) \propto ({\dot{i}}/{i}) \equiv {\rm constant}$.

We model an initially unexcited disk of planetesimals that becomes today's cold classical population. In Section 3.1, we present the observational constraints on the excitation of this population. In Section 3.2, we present an analytical model for the evolution of this planetesimal disk under the influence of Neptune, which we use to predict and interpret the results of numerical simulations. In Section 3.3, we justify directly modeling only Neptune instead of all four giant planets.

3.1. Constraints from the Observations of Cold Classical Objects

The cold classicals are a class of dynamically "cold" objects on low-eccentricity, low-inclination orbits, with positions starting at 42.5 AU, the region interior to which is unstable, and falling off quickly beyond 45 AU (Kavelaars et al. 2009). We assume that today's cold classical KBOs are remnant planetesimals that formed in situ, and we use these terms interchangeably. Strong constraints can be placed on the dynamical history of the solar system by requiring that Neptune does not disrupt these objects as it migrates outward on an inclined and/or eccentric orbit. In this section, we discuss our specific criteria for "preserving" the cold classicals, which we will use to place constraints on the set of parameters (Section 2.1) defining Neptune's orbital history.

There is evidence that the classical KBOs have a bimodal inclination distribution (Brown 2001; Gulbis et al. 2010; Volk & Malhotra 2011). The cold classicals are defined as the class of objects with a distribution of inclinations i centered on a low inclination with a small width in inclination. The functional form of this distribution is typically modeled as a Gaussian multiplied by sin i. One (Gulbis et al. 2010) of the three proposed models for the de-biased inclination distribution differs substantially from the other two (Brown 2001; Volk & Malhotra 2011) with respect to the relative populations of the cold and hot classicals and the width of the hot (high-i) component (Figure 1, top panel). However, all three distributions are similar for the cold classicals (Figure 1, bottom panel). The number of cold classicals per inclination bin falls off almost entirely by i = 6° for the models of Brown (2001) and Gulbis et al. (2010) and by i = 4° for the model of Volk & Malhotra (2011). Therefore, we require that Neptune's dynamical history should not excite the cold classicals above an inclination of 6°.

Zoom In Zoom Out Reset image size

The disruption criterion for the cold classical KBO eccentricities is more subtle. The cold population is defined by its inclination distribution, not its eccentricity distribution. We could imagine a cold population of objects which have inclinations below 6° but a uniform distribution of eccentricities. If this were true (and we will demonstrate that it is not), the initial planetesimal disk could be excited to arbitrarily large eccentricities. Moreover, the eccentricity distribution could be shaped entirely by the long-term stability of the KBO orbits. In this case, objects excited in eccentricity during Neptune's high-eccentricity period would be ejected from the system over billions of years.

To test whether the eccentricities of the cold classicals are sculpted solely by stability, we compared observed cold classical objects to a stability map created by Lykawka & Mukai (2005) (Figure 2). Lykawka & Mukai (2005) do not use proper elements, so we use the instantaneous orbital elements of the observed objects. We have confirmed that the features of the distributions we identify below are qualitatively the same using proper elements (Dawson & Murray-Clay 2012, Appendix A.1). Because the hot and cold populations overlap (Figure 1), we cannot definitively determine to which distribution any particular object belongs. However, for all three model distributions, less than 10% of objects with inclinations i < 2° are hot. We find that the eccentricities of i < 2° objects in the region from 42.5 to 45 AU are confined well below the survival limit. Consequently, we can conservatively constrain the dynamical history of Neptune: Neptune cannot excite the cold classical objects in this region above e = 0.1.

Zoom In Zoom Out Reset image size

Thus we can impose two conservative criteria for preserving the cold classicals:

• 1.  
  In the region from 42.5 to 47.5 AU, the inclinations of the cold classicals must not be excited above i < 6°.
• 2.  
  In the region from 42.5 to 45 AU, the eccentricities must not be excited above e < 0.1.

3.2. Secular Evolution Model for the Planetesimals

Excitation of the in situ planetesimal population by a perturbing planet on an inclined and/or eccentric orbit occurs through secular evolution. Here we present simple analytical expressions that we will use to predict and interpret the results of our integrations in Section 4.

3.2.1. Dynamics of Secular Evolution

Due to forcing from Neptune,⁴ both the eccentricity and inclination of a planetesimal undergo secular oscillations on timescales of order a million years. The planetesimal's total eccentricity is the vector sum of its forced eccentricity e_forced, imparted by Neptune, and its free eccentricity e_free, set by initial conditions (Figure 3). The e_free vector precesses about the e_forced vector at the angular frequency g_KBO. The secular evolution of the planetesimal's inclination is analogous to—yet, to lowest order, separable from—the eccentricity evolution. The planetesimal's total inclination is the sum of its forced and free inclinations, and the free inclination precesses about the forced inclination. In the case of inclination, we can think of the particle's orbit being inclined by i_free with respect to the "forced" plane and precessing about the forced plane at the rate g_KBO. See Murray & Dermott (2000), Chapter 7, for a pedagogical presentation of secular evolution.

Zoom In Zoom Out Reset image size

The vector components of the planetesimal's eccentricity are h = esin ϖ and k = ecos ϖ, where ϖ is the planetesimal'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 _{N}) \nonumber \\\ms k &=& e_{\rm free}\cos (g_{\rm KBO}t + \beta) + e_{\rm forced}\cos (\varpi _{N}), \end{eqnarray} \tag{ 2 }$$

where

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

The constants e_free and β are determined from the initial conditions. 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 planetesimal, a, all of which are assumed to be constant. The functions b are standard Laplace coefficients. The secular frequency of the KBO is g_KBO, m_N is the mass of Neptune, m_sun is the mass of the Sun, and n = (Gm_sun/a³)^1/2 is the planetesimal's mean motion.

Similarly, the vector components of the planetesimal's inclination are p = isin Ω and q = icos Ω, where Ω is the planetesimal's longitude of ascending node. Secular forcing by Neptune causes p and q to evolve as (to first order in i and i_N)

$$\begin{eqnarray} q &=& i_{\rm free}\sin (-g_{\rm KBO}t + \gamma) + i_{\rm forced}\sin (\Omega _{N}) \nonumber \\\ms p &=& i_{\rm free}\cos (-g_{\rm KBO}t + \gamma) + i_{\rm forced}\cos (\Omega _{N}), \end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray} &&i_{\rm forced}= i_N. \end{eqnarray} \tag{ 5 }$$

The constants i_free and γ are determined from the initial conditions. Here, Ω_N is the longitude of the ascending node of Neptune and i_N is the inclination of Neptune. While e_forced depends on α (Equation (3)), to first order i_forced is equal to i_N, regardless of the particle's semimajor axis.

Three quantities are plotted in Figure 4 as a function of a (for two different a_N): (1) e_forced/e_N, (2) i_forced/i_N, and (3) the secular period 2π/g_KBO. These values describe the amplitude and timescale of the secular excitation of the in situ planetesimal population. For example, when Neptune is at 20 AU, a planetesimal at 45 AU has a forced eccentricity of 0.55e_N and a forced inclination of i_N. If the planetesimal begins with e = i = 0, it will reach its maximum eccentricity e = 2 × 0.55e_N = 1.1e_N and maximum inclination i = 2i_N on a timescale of (1/2)(2π/g_KBO) = 13 Myr (i.e., half a secular oscillation period).

Zoom In Zoom Out Reset image size

3.2.2. Secular Excitation of the Planetesimal Disk

The secular excitation of a planetesimal disk can be modeled using the expressions above. Here we consider that Neptune is temporarily on an inclined and/or eccentric orbit after undergoing scattering on an effectively instantaneous timescale (i.e., much less than the secular timescale). Then Neptune imparts a forced eccentricity e_forced and inclination i_forced on the initially cold planetesimal disk. Thus each planetesimal has e_free ∼ e_forced and i_free ∼ i_forced. The planetesimal's total eccentricity and inclination, each a vector sum of the free and forced, now oscillates from the initial values of e(0) ∼ 0 and i(0) ∼ 0, reaching maximum values of e = e_free + e_forced ∼ 2e_forced and i = i_free + i_forced ∼ 2i_forced, on a timescale set by the secular evolution rate g_KBO. Thus an initially cold planetesimal disk will become excited.

3.3. The Case for Modeling Only Neptune's Effects on the Planetesimals

We limit our parameter study to include only the planet Neptune. We have several reasons for choosing this approach. First, unlike previous approaches, we are not attempting to create a single model that reproduces the entire Kuiper Belt in detail, but to constrain which histories of Neptune are consistent with a major qualitative feature: the unexcited orbits of the cold classicals. Our constraints will feed into more detailed models. Second, restricting the parameter study to just Neptune drastically reduces the number of parameters, allowing us to thoroughly explore the remaining parameter space. Third, we find (Dawson & Murray-Clay 2012) that the primary effects of the other planets on the Kuiper Belt are indirect: they alter the orbit of Neptune, which in turn affects the Kuiper Belt. Our modifications to the Mercury 6.2 integrator (Appendix) allow us to model any orbital evolution of Neptune without needing to include the other planets. Figure 5 provides one example pair of integrations, demonstrating that cold planetesimals evolve similarly in the presence of only Neptune and in the presence of all four giant planets.

Zoom In Zoom Out Reset image size

One main effect of other planets is to cause Neptune's longitude of periapse to precess. However, the disruption of the cold classicals is not significantly affected if Neptune's precession period is comparable to or longer than the secular excitation time (Figure 5). If Neptune precesses very quickly, the cold classicals could be preserved outside our constraints, as proposed by Batygin et al. (2011). In addition, strong interactions between Neptune and Uranus can cause Neptune's semimajor axis to oscillate, leading to orbital chaos in the classical region, which can cause additional excitation. We explore these additional complications in Dawson & Murray-Clay (2012).

Here we do not explore the g₈ and ν₈ secular resonances, which over time remove low-eccentricity and low-inclination objects, respectively. The location of these secular resonances will impose additional constraints on preserving the cold population. However, since they remove unexcited objects, leaving excited objects behind, these resonances will not allow for parameters of Neptune that we rule out.

3.4. Considerations for a Massive Planetesimal Disk

Throughout the paper, we treat the planetesimals as massless test particles. We parametrically model the orbital evolution of Neptune caused by the planetesimals—Neptune's migration and the damping of its eccentricity and inclination—but do not explicitly consider the effects of planetesimal self-gravity. In this section, we consider some of these effects and how they might impact our conclusions.

One important consideration is whether the transfer of angular momentum between Neptune and the disk is sufficient to excite a massive disk to the extent we assume for massless test particles. Since the total angular momentum of the system is conserved, the angular momentum deficit (AMD) due to the secular excitation of the cold classicals cannot exceed Neptune's initial AMD (see Laskar 1997 for a description of the concept of AMD). Using the expansion of the AMD by Hahn (2003, Equations (26a) and (26b)), the ratio of the total AMD of the cold objects to the AMD of Neptune is

$$\begin{eqnarray} &&\frac{\rm AMD_{\rm cold}}{\rm AMD_{\rm Neptune}} = \frac{ m_{\rm cold} }{m_N} \sqrt{\frac{a_{\rm cold}}{a_N}} \left(\frac{e_{\rm cold}^2+i_{\rm cold}^2}{e_N^2+i_N^2}\right). \end{eqnarray} \tag{ 6 }$$

Since the factor $\sqrt{{a_{\rm cold}}/{a_N}}$ is roughly unity, if the mass in the cold classicals exceeds Neptune's mass, the AMD of Neptune would be too small to excite the cold classicals to e_N and i_N. Thus if the mass in cold classicals were large enough, the objects could remain at low eccentricities and inclinations even if Neptune's eccentricity and inclination were large.

The outer solar system likely contained several tens of Earth masses at early times, comparable to Neptune's mass of 17 M_⊕. Such a massive disk is required to form large KBOs such as Pluto by coagulation (Stern & Colwell 1997; Kenyon & Luu 1998) and to drive substantial migration (Fernandez & Ip 1984; Hahn & Malhotra 1999) and/or eccentricity and inclination damping of Neptune's orbit. In contrast, Fuentes & Holman (2008) combined the results of several surveys to estimate that the mass of the current classical Kuiper Belt is only 0.008 ± 0.001 M_⊕, 40% of which is the cold component (Kavelaars et al. 2009).

We assume throughout the rest of this paper that the total mass—and thus AMD—of the cold classicals is smaller than that of Neptune (i.e., closer to today's mass). This choice implicitly assumes that the primordial disk was partially truncated or that the cold classicals were depleted before the era we treat in this paper. If the cold classical belt began with a low mass—compared to the inner disk, responsible for Neptune's migration and for dynamical friction—truncation of the planetesimal disk may explain why Neptune's migration halted at 30 AU (Levison & Morbidelli 2003). Dynamical depletion would tend to excite the cold classical population to higher eccentricities than observed, but could be consistent if all excited objects were subsequently scattered out of the region from 42.5 to 45 AU. Thus this dynamical depletion would need to occur before the delivery of the observed hot classical population, which has higher eccentricities. During the era we consider in this paper, Neptune must deliver the hot objects without disrupting the cold ones.

The constraints presented here will require revision if Neptune's primary sculpting of the Kuiper Belt occurred while the disk was massive. For example, this could be true if the cold classical belt was depleted through collisions following Neptune's orbital evolution. Collisional depletion occurs through collisional grinding, followed by ejection of small grains by radiation forces, and Kenyon & Luu (1998) and Kenyon & Bromley (2001) argue for this scenario in the context of coagulation models for the growth at KBOs. A break in the size distribution of the cold classicals at R = 25–50 km has been attributed to collisional grinding (Pan & Sari 2005). However, Nesvorný et al. (2011) argue that, given this interpretation, a large fraction of binaries with small components (R < 50 km) would have been disrupted. Such binaries are observed, implying that collisional grinding did not generate the break. Hence collisional grinding could only have substantially depleted the cold Kuiper Belt if the majority of the Belt's mass was initially sequestered in much smaller objects.

Another potentially important effect is the propagation of spiral density waves and torsion waves in a massive disk (Ward & Hahn 1998; Hahn 2003), excited at secular resonances and at the edge of the disk. In particular, spiral density waves tend to "smear out" the excitation of the eccentricities and inclinations of the planetesimals near secular resonances. Spiral density waves are potentially important in sculpting the dynamical structure of the Kuiper Belt—Ward & Hahn (1998) use them to place constraints on the mass of the Kuiper Belt beyond 50 AU. The exact behavior of the density waves depends on a variety of disk parameters—including the size distribution of planetesimals, the surface density profile in the disk, and the thickness of the disk—so we refer the reader to Hahn (2003) for a detailed exploration. (For ease of comparison, note that assuming a surface density profile Σ∝a^−1.5, the current mass of the cold Kuiper Belt corresponds to a Σ roughly a factor of three smaller than the lowest surface density simulated by Hahn 2003.) Without the influence of density waves, very large forced eccentricities and inclinations (e → ∞, i → ∞) occur only close to secular resonances. The density waves spread the elevated eccentricity and inclination, confined to the secular resonance in the massless disk case, across the Kuiper Belt from 40 to 50 AU. However, our constraints do not include the extra excitation caused by secular resonances, which excite the cold classicals even more than we consider. Therefore, the constraints we will place on ruling out parameters of Neptune that cause excessive excitation will still hold and, as future work, the additional excitation caused by secular resonances plus density waves may place additional constraints.

In the absence of secular resonances, waves launched at the disk edge could similarly smear eccentricities and inclinations. We assume that the eccentricities and inclinations of planetesimals are not self-damped by spiral density waves launched at the disk edge. The eccentricities and inclinations of planetesimals could also be damped by dynamical friction caused by smaller, unobserved bodies whose random velocities were damped by collisions (Goldreich et al. 2004) or by resonant relaxation (Tremaine 1998). We assume that neither of these effects contribute significantly during the era of interest.

Finally, the eccentricities and inclinations of the planetesimals may not be excited in the first place if the planetesimal disk is massive enough to cause Neptune's eccentricity to quickly precess, lowering the forced eccentricity of the cold classicals. In order for the disk to significantly affect the precession rate of the cold classicals, its mass must be at least that of Neptune. A scenario involving a quickly precessing Neptune is explored in Batygin et al. (2011) and discussed in our companion paper, Dawson & Murray-Clay (2012).

We present here the results of a set of integrations containing Neptune and a collection of test particles designed to represent the in situ planetesimals that become today's cold classical KBOs. The orbital evolution of Neptune is modeled as described in Section 2.2 and the Appendix. In Section 4.1, we describe how we computationally model the planetesimals. In Section 4.2, we characterize the secular excitation of the planetesimals under the influence of an inclined and eccentric Neptune and place limits on the eccentricity and inclination of Neptune. In Section 4.3, we determine the effect of damping of Neptune's eccentricity and inclination on the excitation of the planetesimals. We demonstrate that there are two regimes for damping Neptune's orbit—in the slow-damping regime, the planetesimals evolve to their initial free eccentricity and inclination, set by Neptune's initial eccentricity and inclination, and in the fast-damping regime, their eccentricity and inclination are frozen at the value reached at the damping time. We also show that the evolution of the inclination and eccentricity can be treated separately. In Section 4.4, we consider the effects of Neptune's migration. Finally, we present two example scenarios in Section 4.5 demonstrating the principles we have derived: a pathological scenario in which the evolution of Neptune's semimajor axis, eccentricity, and inclination each occur on a different timescale, and a realistic scenario consistent with preserving the cold population.

4.1. Computational Model of Planetesimals

4.2. Secular Excitation of Planetesimals Under the Influence of an Inclined and Eccentric Neptune

In scenarios of the early solar system in which Neptune is scattered onto an eccentric and/or inclined orbit, Neptune can potentially disrupt an in situ population of planetesimals above the confined eccentricities and inclinations we observe (Section 3.1) in the cold classical Kuiper Belt today. Neptune forces the planetesimals to undergo secular evolution (Section 3.2), in which they oscillate through high eccentricities and inclinations. The rate of secular evolution depends on the location of the planetesimal relative to Neptune (Figure 4). Figure 6 shows snapshots of the secular evolution of test particles under the influence of Neptune when the planet is stationary at 20 AU (top), migrates from 20 to 30 AU (middle), and is stationary at 30 AU (bottom). In their secular evolution, the planetesimals reach maxima e = 2e_forced and i = 2i_forced. In the middle panel, the evolution transitions from matching the top panel (early snapshots, left) to matching the bottom panel (late snapshots, right). Thus, in the case of migration, we can estimate the secular evolution using a constant a_N. When e_N or i_N has not yet damped, we can estimate the secular evolution as if a_N were constant at the location where Neptune spent most time. Because we have chosen a migration law for which ${\dot{a}_N}/{a_N}$ decreases with time, planetesimals behave as if a_N has always had the value a_N(t). We will discuss migration in detail in Section 4.4.

Zoom In Zoom Out Reset image size

If Neptune's eccentricity and inclination are small (Figure 7), the test particles remain below the observational limit of e < 0.1, i < 6° indefinitely. Since, to first order, i_forced is independent of α, the planetesimal's position relative to Neptune, the planetesimals will remain below i < 6° if i_N < 3°. The forced eccentricity does depend on α (Figure 4), decreasing with the planetesimal's distance from Neptune. Thus to satisfy our conservative criterion established in Section 3.1, which requires planetesimals with 42.5AU < a < 45 AU to remain at e < 0.1, Neptune's eccentricity must stay below 0.09 at 20 AU or 0.06 at 30 AU. However, if e_N and i_N damp due to dynamical friction with the planetesimals, the initial values of Neptune's eccentricity and inclination can be higher. In the next subsection, we discuss the effects of damping.

Zoom In Zoom Out Reset image size

4.3. Damping of Neptune's Eccentricity and Inclination

The constraints we placed previously were on how high Neptune's eccentricity and inclination can remain indefinitely. During secular evolution, planetesimals reach eccentricities up to 2e_forced and inclinations up to 2e_forced. However, the planetesimals evolve to final values less than these maximum values, when Neptune's eccentricity and inclination damp. Here we refine the limits on Neptune's eccentricity and inclination in light of damping, considering a range of timescales for the dynamical-friction-driven damping rates of e_N and i_N. The actual damping timescales depend on the surface density and size distribution of the planetesimals.

4.3.1. Eccentricity and Inclination can be Treated Separately

To first order, the effects on the planetesimals of Neptune's inclination and eccentricity, including damping, can be treated independently. In first-order secular theory (Section 3.2), the evolution of a planetesimal's e is independent of Neptune's inclination i_N and of a planetesimal's i is independent of Neptune's eccentricity e_N. Thus we can place constraints on the evolution of Neptune's eccentricity without taking into account its inclination or vice versa. Figure 8 illustrates the independence of the eccentricity and inclination parameters.

Zoom In Zoom Out Reset image size

4.3.2. Two Regimes for Damping: Slow and Fast

The damping of Neptune's eccentricity e_N or inclination i_N affects the secular excitation of the planetesimals in two regimes: slow and fast. We summarize the constraints on Neptune's eccentricity and inclination in Table 1.

Table 1. Constraints on Neptune's Eccentricity and Inclination

	No Damping	Slow Damping	Fast Damping
eN	$2 \frac{b_{3/2}^{(2)}(\alpha)}{ b_{3/2}^{(1)}(\alpha) }e_{N} < 0.1$	$\frac{b_{3/2}^{(2)}(\alpha)}{ b_{3/2}^{(1)}(\alpha) }e_{N} < 0.1$	$\frac{b_{3/2}^{(2)}(\alpha)}{ b_{3/2}^{(1)}(\alpha) }e_{N} \sin (g_{\rm KBO}\tau _e) < 0.1$
iN	2iN < 6°	iN < 6°	iNsin (gKBOτi) < 6°

Download table as:  ASCIITypeset image

In the slow regime, e_N or i_N damps on a timescale τ_e or τ_i, respectively, longer than the time ((1/2)2π/g_KBO) for e or i to reach its maximum value. Because each planetesimal has an initial i(t = 0) = e(t = 0) = 0, its free inclination i_free and eccentricity e_free are set by Neptune's initial e_N(t = 0) and i_N(t = 0). In this slow regime, the planetesimal's i_free and e_free are conserved. As e_N and i_N damp, the planetesimal's forced eccentricity e_forced and inclination i_forced decrease, and the total eccentricity e and inclination i of the planetesimal approach e_free and i_free.

Thus, in the slow regime, the planetesimals will evolve to a value below i < 6° if i_N < 6°. To satisfy our conservative criterion established in Section 3.1, which requires planetesimals with 42.5 AU < a < 45 AU to remain at e < 0.1, Neptune's eccentricity must stay below 0.18 at 20 AU and 0.12 at 30 AU. Thus in the slow-damping case, there is a strong constraint (Table 1) on the maximum eccentricity and inclination at which Neptune can remain over long timescales. Note that these values are twice the values given in Section 4.2. This is because, in the slow-damping regime, the planetesimals will evolve to e_free = e_forced(t = 0) and i_free = i_forced(t = 0), whereas if e_N and i_N remain high indefinitely, the planetesimals' e and i continue to oscillate, reaching a maximum of 2e_forced and 2i_forced.

If e_N and i_N damp in the fast regime, on a timescale shorter than the secular excitation time, i_free and e_free are not conserved. The planetesimal's total e and i are frozen at the values they reach after approximately one damping time. Figure 9 illustrates the behavior of the planetesimals in this regime in integrations in which Neptune's eccentricity or inclination damps.

Zoom In Zoom Out Reset image size

Neptune's orbit could have been even more eccentric and inclined than in the slow regime if e_N and i_N damped quickly (Table 1, right column). If the planetesimal has not yet reached the maximum of its secular cycle after one damping time, instead of converging to the e_free and i_free set by the initial conditions, the planetesimal evolves to e_final = e_freesin (g_KBOτ_e) and i_final = i_freesin (g_KBOτ_i) (see Dawson & Murray-Clay 2012 for a detailed explanation and justification). Consider again a planetesimal at 42.5 AU. If Neptune's eccentricity and inclination damp on a timescale of τ = 0.32 Myr, the final eccentricity and inclination of the planetesimal are reduced by a factor of sin (g_KBOτ) = 0.5 compared to the slow-damping case. Thus the planetesimals will evolve to a value below i < 6° if i_N < 6°/0.5 = 12°. To satisfy our conservative criterion established in Section 3.1, which requires planetesimals with 42.5 AU < a < 45 AU to remain at e < 0.1, e_N must stay below 0.18/0.5 = 0.36 at 20 AU and 0.12/0.5 = 0.24 at 30 AU.

We note that because "slow" and "fast" damping are defined relative to the secular evolution time, and because the secular evolution time increases with the planetesimal's semimajor axis, it may be that the damping is "fast" for particles in the outer disk yet "slow" for particles in the inner disk.

4.4. Migration of Neptune

In the process of secular excitation of the cold classicals, migration alters a planetesimal's secular evolution timescale and forced eccentricity, which depend on α, the ratio of the planetesimal's semimajor axis to Neptune's (the forced inclination is independent of α). For a planetesimal at 42.5 AU, the forced eccentricity is 40% larger when Neptune is at 30 AU versus 20 AU and the secular evolution period is five times shorter. Figure 10 shows the effect of migration on the secular evolution of the planetesimals. The excitation at a given time is well modeled by the secular theory (Equations (2) and (4)) using Neptune at its snapshot location. Recall that we are using a migration law that slows with time, so that Neptune spends increasing amounts of time at each, subsequent location.

Zoom In Zoom Out Reset image size

For retaining the cold classicals, what matters is the ratio of the migration timescale to the damping timescale. When this ratio is large ("slow" migration), Neptune effectively damps at its initial position and we can model the secular evolution at this location (Figure 11, top two rows). When this ratio is small ("fast" migration), Neptune effectively damps at its final position and we can model the secular evolution there (Figure 11, bottom two rows). When the ratio is of order unity, modeling the secular evolution at the location Neptune reaches after half a damping time is a decent approximation (Figure 11, middle row). However, we have not explored the regime in which τ_a ∼ τ_e or τ_a ∼ τ_i in detail. Typically, we expect the damping and migration to occur in either the fast or slow regime, for reasons we will now state. The models in which Neptune is scattered from its location of formation to close to its current location (e.g., Levison et al. 2008) can be considered fast migration. For extensive migration in a planetesimal disk, we expect the migration timescale to be significantly longer than the damping timescale because the "random momentum" is a fraction of order e or i of the Keplerian momentum. (Though we note that Neptune only has to migrate some fraction of its semimajor axis.)

Zoom In Zoom Out Reset image size

4.5. Two Example Scenarios

Consider two example dynamical histories for Neptune. In each, the planet starts at an initial location, eccentricity, and inclination and undergoes migration, eccentricity damping, and inclination damping to evolve to its current orbit. First, we consider a pathological scenario in which the evolution of Neptune's semimajor axis, eccentricity, and inclination each occur on a different timescale (Figure 12). Neptune begins at a_N = 20 AU, e_N = 0.2, and i_N = 10°. On a timescale of τ_a = 3 Myr, it migrates to its current location at 30 AU. The eccentricity damps on a longer timescale of τ_e =30 Myr and the inclination on a shorter timescale of τ_i = 0.3 Myr. For the inclination damping, the migration occurs in the slow regime (τ_a/τ_i > 1), and we model the damping as taking place at 20 AU. With Neptune at 20 AU, the secular timescale of a planetesimal at 45 AU is 2π/g_KBO = 24 Myr and the initial forced inclination of the planetesimal is i_forced = i_N = 10°. The value of the inclination to which the planetesimal evolves is reduced by a factor of sin (g_KBOτ_i) = 0.07. Thus the planetesimal at 45 AU evolves to a final value of 07. For the eccentricity, the migration occurs in the fast regime. Because the eccentricity damping timescale is much longer than the migration timescale (fast migration), we model the eccentricity damping as taking place at 30 AU. With Neptune at 30 AU, the secular timescale of a planetesimal at 45 AU is 2π/g_KBO = 6 Myr and its initial forced eccentricity is e_forced = 0.78e_N = 0.16. The eccentricity damping is in the slow-damping regime, so the planetesimal's final eccentricity evolves to its initial forced eccentricity, e = 0.16. This predicted value and those for planetesimals at other semimajor axes (orange curve, Figure 12) match the integrations well, validating our method.

Zoom In Zoom Out Reset image size

Finally, we consider a more realistic scenario consistent with preserving the cold population (Figure 13). Neptune begins at 26 AU, e_N = 0.35, and i_N = 14°. On a timescale of τ_a = 10 Myr, it migrates to its current location at 30 AU. Its eccentricity and inclination damp on a timescale of τ_e = τ_i = 0.3 Myr. In this slow migration regime, we model the damping taking place with Neptune at 26 AU. Thus the secular timescale of a planetesimal at 45 AU is 2π/g_KBO  ∼ 11 Myr, its initial forced eccentricity is 0.7e_N = 0.24 and its initial forced inclination is i_N = 14°. The values of eccentricity and inclination to which the planetesimal evolves are reduced by a factor of sin (g_KBOτ_e) = sin (g_KBOτ_i) = 0.16. Thus the planetesimal at 45 AU evolves to a final eccentricity of 0.04 and inclination of 24, values satisfying the observational constraints established in Section 3.1.

Zoom In Zoom Out Reset image size

As a first step in a comprehensive study of the impact on the Kuiper Belt of a wide range of possible dynamical histories of the outer solar system, we have performed a suite of numerical integrations probing the impact of the orbital evolution of a single planet on a disk of planetesimals. We have presented the observational evidence for a population of dynamically cold objects in the Kuiper Belt in the region from 42.5 to 45 AU that are confined to e < 0.1 and i < 6°. We argued that recent models of Kuiper Belt sculpting—which explain many of the observed dynamical and physical properties of the Kuiper Belt—do not generate or preserve sufficiently low eccentricities for the cold classicals (e.g., Gomes 2003; Levison et al. 2008; Morbidelli et al. 2008). Our results have revealed several principles key to constraining which orbital histories of Neptune are consistent with preserving the in situ planetesimal population at the low eccentricities and inclinations required by the observations:

• 1.  
  If Neptune is scattered onto an eccentric and/or inclined orbit, it will secularly excite the eccentricities and inclinations of an in situ planetesimal population.
• 2.  
  Planetesimals starting with e = i = 0 reach eccentricities and inclinations up to twice their forced eccentricity and inclination on timescales given by Equations (2) and (4).
• 3.  
  As Neptune's eccentricity and inclination damp, the planetesimals evolve to their final eccentricities and inclinations. If the damping timescales τ_e and τ_i of Neptune's orbit are slow compared to the secular evolution time, a planetesimal evolves to its initial free eccentricity and inclination, which are set by Neptune's initial eccentricity and inclination. If the damping is fast compared to the secular evolution time, the planetesimal's e and i effectively freeze at the values they reach after one damping time, reaching final values of e_freesin (g_KBOτ_e) and e_forcedsin (g_KBOτ_i), respectively. See Table 1 for constraints on Neptune's eccentricity and inclination in the two damping regimes.
• 4.  
  The effects of Neptune's (1) eccentricity evolution and (2) inclination evolution on the planetesimals can be treated separately to first order.
• 5.  
  At a given location in the planetesimal disk, the secular excitation timescales and forced eccentricity (but not the forced inclination) depend on Neptune's location, which is altered by Neptune's migration. When Neptune's migration is slow relative to the damping time, the secular evolution effectively takes place at Neptune's initial location. When Neptune's migration is fast relative to the damping time, the secular evolution effectively takes place at Neptune's final location.

From these principles, it is evident that the three models described in Levison et al. (2008) would not be able to retain a cold classical population. Neptune begins with an eccentricity of 0.3 and a semimajor axis of 27.5 AU (runs A and C) or 28.9 AU (run B) and damps on a timescale of 1 Myr (runs A and B) or 3 Myr (run C) in the slow migration regime. When Neptune is at 27.5 AU, a planetesimal at 42.5 AU has a forced eccentricity of e = 0.76e_N = 0.23 and a secular evolution timescale of 6 Myr. The final eccentricity the planetesimal evolves to is reduced by a factor of sin (g_KBOτ_e) = 0.85 for a damping timescale of 1 Myr (fast damping) and is not reduced for a damping timescale of 3 Myr (slow damping). Thus the final eccentricity of the planetesimal is 0.19 (run A) or 0.23 (run C), well above the observational limit. When Neptune is at 28.9 AU (run B), a planetesimal at 42.5 AU has a forced eccentricity of e = 0.79e_N = 0.24 and a secular evolution timescale of 5 Myr. The final eccentricity the planetesimal evolves to is reduced by a factor of sin (g_KBOτ_e) = 0.97 for a damping timescale of 1 Myr. Thus the final eccentricity of the planetesimal is 0.23, well above the observational limit. According to the constraints established in our paper, any of these three initial conditions for Neptune could retain the cold classicals if Neptune's eccentricity were to damp more quickly.

Based on these principles, we can place robust constraints on Neptune's dynamical history. For example, in the regime of slow damping and slow migration, if Neptune's initial semimajor axis is 20 AU, then its eccentricity must stay below 0.18 and inclination below 6°. If Neptune's initial semimajor axis is 30 AU—or if it migrates quickly, relative to the damping time, from its initial location to 30 AU—then its eccentricity must stay below 0.12. In the case of fast damping—on a timescale shorter than a planetesimal's secular excitation time—the initial eccentricity and inclination of Neptune can be even higher.

Having established these principles, we complete our parameter study in two companion papers: Dawson & Murray-Clay (2012) and R. I. Dawson & R. A. Murray-Clay (in preparation). In the first companion paper, we consider more generally the constraints on Neptune's dynamical history from the eccentricity distribution of the classical KBOs and incorporate several other important effects in the constraints from the cold classicals:

• 1.  
  A more accurate model for secular excitation that includes higher-order terms.
• 2.  
  The greatly increased secular frequency near Neptune's mean-motion resonances, which places strong constraints on Neptune's semimajor axis when its eccentricity is high.
• 3.  
  The effects of the other giant planets, including precession of Neptune, which can reduce a planetesimal's forced eccentricity (Batygin et al. 2011), and oscillations in Neptune's semimajor axis, which can create a chaotic sea in the Kuiper Belt region.

Then we combine these constraints for retaining the cold classicals with constraints for creating the hot classical population and identify which regions of parameter space of Neptune's dynamical history can produce the hot classical population without disrupting the cold population. In the second companion paper, we place constraints on Neptune's inclination and inclination damping time.

Our parameterization of Neptune's orbital history allows us to constrain which orbital histories are consistent with maintaining the cold classical population. By combining the observational constraints (Section 3.1) with the principles derived in this paper, we can immediately check whether a particular set of initial conditions for Neptune's semimajor axis, eccentricity, and inclination, and the rates of its migration and eccentricity and inclination damping, are consistent with maintaining the observed dynamically cold population, without performing computationally expensive integrations. A picture, both qualitative and quantitative, is emerging of what dynamical histories of Neptune allow a promising general scenario—Neptune's delivery of the hot classicals from the inner disk to classical region, where the cold population has formed in situ—to be consistent with observed low eccentricities and inclinations of the cold classical population. Barring fast precession of Neptune's orbit, which we do not consider in this work, the existence of the cold classical population implies that Neptune could have spent only a limited time at high eccentricity and/or inclination during its dynamical history.

We thank an anonymous referee for helpful comments, including recommending the addition of Section 3.4. S.W. gratefully acknowledges funding by the National Science Foundation Research Experiences for Undergraduates (REU) and the Department of Defense Awards to Stimulate and Support Undergraduate Research Experiences (ASSURE) programs under grant no. 0754568. R.I.D. is supported by a National Science Foundation Graduate Research Fellowship under grant DGE-1144152 and R.M.C. by the Smithsonian Institution. The computations in this paper were run on the Odyssey cluster supported by the Harvard FAS Sciences Division Research Computing Group.

In this Appendix, we follow Lee & Peale (2002) to derive modifications to the equations of motion of an orbiting body to allow any form of evolution of the body's orbital elements. We use these modified equations to model the early solar system evolution of Neptune's orbit caused by interactions with the other giant planets and with the planetesimal disk. Directly modeling Neptune, the other giant planets, and tens of thousands of massive planetesimals would be computationally prohibitive. Instead, we model Neptune alone, with its orbit evolving as it would under the influence of the other planets and the planetesimal disk, including migration, damping of its eccentricity, and damping of its inclination. However, the equations we derive can be used to implement any type of orbital evolution, not just migration and damping. For example, in a companion paper (Dawson & Murray-Clay 2012), we use the modifications to cause Neptune's semimajor axis to oscillate, as it would due to resonant interactions with Uranus.

We modify the Mercury 6.2 (Chambers 1999) N-body integration code to allow for arbitrary orbital evolution by adding extra terms to the equations of motion at each time step. We add these terms via a user-defined force, already supported by Mercury 6.2, and an additional velocity modification described below. We note that Lee & Peale (2002) use their derived modifications for a non-symplectic code and use a symplectic implementation when employing a symplectic algorithm. However, we do not add the extra terms in a symplectic manner, even for a symplectic integration algorithm; thus these modifications must represent small perturbations on the planet. We determined that a time step of 200 days was sufficiently small for our models, yielding the same results as smaller step sizes.

We make these modifications as follows:

• 1.  
  At each time step, Mercury 6.2 calls the user-defined force routine. We have modified this routine to return not only a user-defined acceleration but a user-defined velocity. We also modify the routine to only apply these corrections for Neptune, not for the planetesimals.
• 2.  
  The user-defined force routine calculates the user-defined velocity additional terms from Equations (A9)–(A11) below and the user-defined acceleration additional terms from Equations (A12)–(A14) below.
• 3.  
  The acceleration additional terms are used to advance the velocity of Neptune and the velocity additional terms to advance the position of Neptune, in addition to the usual gravitational forces.

To derive these additional terms, we follow Lee & Peale (2002) but instead of holding the orbital inclination i constant, we allow it to vary with time.

Equation (A1) gives the position vector components (x, y, and z) as a function of the orbital elements (compare to Equation (A5) of Lee & Peale). The variable r is the distance of the planet from the central body, Ω is the longitude of the ascending node in the xy plane, ω is the argument of periapse, f is the true anomaly, and i is the inclination:

$$\begin{eqnarray} x &= r \cos {\Omega } \cos {(\omega + f)} - r \cos {i} \sin {\Omega } \sin {(\omega + f)} \nonumber \\ y &= r \sin {\Omega } \cos {(\omega + f)} + r \cos {i} \cos {\Omega } \sin {(\omega +f)}\nonumber \\ z &= r \sin {i} \sin {(\omega + f).} \end{eqnarray} \tag{ A1 }$$

The velocity vector components are given in Equation (A2) (compare to Equation (A6) of Lee & Peale). A dot over a variable indicates its derivative with respect to time:

$$\begin{eqnarray} \dot{x} &=& \cos {\Omega } \,[\dot{r} \cos {(\omega + f)} - r \dot{f} \sin {(\omega + f)}]\nonumber\\ && -\, \sin {\Omega } \,[\dot{r} \cos {i} \sin {(\omega + f)} + r \dot{f} \cos {i} \cos {(\omega + f)} -z \dot{i}] \nonumber \\ \dot{y} &=& \sin {\Omega } \,[\dot{r} \cos {(\omega + f)} - r \dot{f} \sin {(\omega + f)}]\nonumber\\ && +\, \cos {\Omega } \,[\dot{r} \cos {i} \sin {(\omega + f)} + r \dot{f} \cos {i} \cos {(\omega + f)} - z \dot{i}]\nonumber \\ \dot{z} &=& \dot{r} \sin {i} \sin {(\omega + f)} + r \dot{f} \sin {i} \cos {(\omega + f)}]\nonumber\\ &&+r \dot{i} \cos {i} \sin {(\omega + f)}. \end{eqnarray} \tag{ A2 }$$

Below are the x-components of the "velocity" and "acceleration" additional terms used to update the body's position and velocity, respectively, at each time step (compare to Equations (A3) and (A4) of Lee & Peale). The variable a is the semimajor axis and e is the orbital eccentricity:

$$\begin{equation} \frac{dx}{dt} \bigg |_{\dot{a}} + \frac{dx}{dt} \bigg |_{\dot{e}} + \frac{dx}{dt} \bigg |_{\dot{i}} = \frac{\partial x}{\partial a} \dot{a} + \frac{\partial x}{\partial e} \dot{e} + \frac{\partial x}{\partial i} \dot{i} \end{equation} \tag{ A3 }$$

$$\begin{equation} \frac{d \dot{x}}{dt} \bigg |_{\dot{a}} + \frac{d \dot{x}}{dt} \bigg |_{\dot{e}} + \frac{d \dot{x}}{dt} \bigg |_{\dot{i}} = \frac{\partial \dot{x}}{\partial a} \dot{a} + \frac{\partial \dot{x}}{\partial e} \dot{e} + \frac{\partial \dot{x}}{\partial i} \dot{i}. \end{equation} \tag{ A4 }$$

Similar expressions can be derived for the other coordinates. It should be noted that these are simply the additional terms to the equations of motion resulting from the orbital evolution—$\dot{a}$, $\dot{e}$, and/or $\dot{i}$—and do not describe the overall motion of the system. All coordinates, velocities, and accelerations must be calculated with respect to the central body. In order to compute these partial derivatives, we must first define the variables r, $\dot{r}$, and $r \dot{f}$ following Murray & Dermott (2000):

$$\begin{eqnarray} &&r = \frac{a (1 - e^{2})}{1 + e \cos {f}}, \quad \dot{r} = \frac{n a}{\sqrt{1 - e^{2}}} e \sin {f}, \nonumber\\ &&\qquad r \dot{f} = \frac{n a}{\sqrt{1 - e^{2}}} (1 + e \cos {f}). \end{eqnarray} \tag{ A5 }$$

Next, we calculate the partial derives of $r, \dot{r},$ and $r\dot{f}$ with respect to the orbital elements a, e, and i (compare to Equation (A7) of Lee & Peale):

$$\begin{equation} \frac{\partial r}{\partial a} = \frac{r}{a}, \quad \frac{\partial r}{\partial e} = \bigg [ - \frac{2 e r}{1 - e^{2}} - \frac{r^{2} \cos {f}}{a (1 - e^{2})}\bigg ], \quad \frac{\partial r}{\partial i} = 0 \end{equation} \tag{ A6 }$$

$$\begin{equation} \frac{\partial \dot{r}}{\partial a} = {-} \frac{\dot{r}}{2 a}, \quad \frac{\partial \dot{r}}{\partial e} = \frac{\dot{r}}{e (1 - e^{2})}, \quad \frac{\partial \dot{r}}{\partial i} = 0 \end{equation} \tag{ A7 }$$

$$\begin{eqnarray} &&\frac{\partial (r \dot{f})}{\partial a} = {-} \frac{r \dot{f}}{2 a}, \quad \frac{\partial (r \dot{f})}{\partial e} = \frac{r \dot{f} (e + \cos {f})}{(1 - e^{2}) (1 + e \cos {f})}, \nonumber\\ &&\qquad\qquad\qquad\frac{\partial (r \dot{f})}{\partial i} = 0. \end{eqnarray} \tag{ A8 }$$

Combining Equations (A6) through (A8) with Equation (A3) yields an expression for change in position (compare to Equation (A8) of Lee & Peale):

$$\begin{equation} \frac{d x}{d t} \bigg |_{\dot{a}} + \frac{d x}{d t} \bigg |_{\dot{e}} + \frac{d x}{d t} \bigg |_{\dot{i}} = \frac{x}{a} \dot{a}+ \bigg [ \frac{r}{a (1 - e^{2})} - \frac{1 + e^{2}}{1 - e^{2}} \bigg ] \frac{x}{e} \dot{e} \,{+}\, (z \, \sin {\Omega }) \, \dot{i}, \end{equation} \tag{ A9 }$$

$$\begin{eqnarray} \frac{d y}{d t} \bigg |_{\dot{a}} + \frac{d y}{d t} \bigg |_{\dot{e}} + \frac{d y}{d t} \bigg |_{\dot{i}} &=& \frac{y}{a} \dot{a} + \bigg [ \frac{r}{a (1 - e^{2})} - \frac{1 + e^{2}}{1 - e^{2}} \bigg ] \frac{y}{e} \dot{e}\nonumber\\ && - (z \, \cos {\Omega }) \, \dot{i}, \end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray} \frac{d z}{d t} \bigg |_{\dot{a}} + \frac{d z}{d t} \bigg |_{\dot{e}} + \frac{d z}{d t} \bigg |_{\dot{i}} &=& \frac{z}{a} \dot{a} + \bigg [ \frac{r}{a (1 - e^{2})} - \frac{1 + e^{2}}{1 - e^{2}} \bigg ] \frac{z}{e} \dot{e}\nonumber\\ && +\,(- \sin {\Omega } x + \cos {\Omega } y) \, \dot{i}. \end{eqnarray} \tag{ A11 }$$

We can now use these additional terms to update the position of a body undergoing any arbitrary orbital evolution in a, e, and/or i, including migration, eccentricity damping, and inclination damping. Below we combine Equations (A6) through (A8) with Equation (A4) to obtain the acceleration terms used to update the body's velocity (compare to Equation (A9) of Lee & Peale):

$$\begin{eqnarray} &&\frac{d \dot{x}}{dt} \bigg |_{\dot{a}} + \frac{d \dot{x}}{dt} \bigg |_{\dot{e}} + \frac{d \dot{x}}{dt} \bigg |_{\dot{i}} = \frac{- \dot{x} + 3 z \dot{i} \sin {\Omega }}{2 a} \dot{a}\nonumber\\\ms &&\quad + \bigg (\cos {\Omega } \, \bigg [ \frac{\partial \dot{r}}{\partial e} \cos {(\omega + f)} - \frac{\partial (r \dot{f})}{\partial e} \sin {(\omega + f)} \bigg ] \nonumber \\\ms &&\quad -\, \sin {\Omega } \bigg [ \frac{\partial \dot{r}}{\partial e} \cos {i} \sin {(\omega + f)} + \frac{\partial (r \dot{f})}{\partial e} \cos {i} \cos {(\omega + f)}\nonumber\\\ms &&\quad - \dot{i} \frac{\partial r}{\partial e}\frac{z}{r} \bigg ] \bigg) \dot{e} + \sin {\Omega } \bigg [ \dot{z} + \frac{\dot{i}}{i} z \bigg ] \dot{i}, \end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray} &&\frac{d \dot{y}}{dt} \bigg |_{\dot{a}} + \frac{d \dot{y}}{dt} \bigg |_{\dot{e}} + \frac{d \dot{y}}{dt} \bigg |_{\dot{i}} =\, \frac{- \dot{y} - 3 z \dot{i} \cos {\Omega }}{2 a} \dot{a}\nonumber \\\ms &&\quad + \bigg (\sin {\Omega } \, \bigg [ \frac{\partial \dot{r}}{\partial e} \cos {(\omega + f)} - \frac{\partial (r \dot{f})}{\partial e} \sin {(\omega + f)} \bigg ] \nonumber \\\ms &&\quad+ \cos {\Omega } \bigg [ \frac{\partial \dot{r}}{\partial e} \cos {i} \sin {(\omega + f)} + \frac{\partial (r \dot{f})}{\partial e} \cos {i} \cos {(\omega + f)}\nonumber \\\ms &&\quad - \dot{i} \frac{\partial r}{\partial e}\frac{z}{r} \bigg ] \bigg) \dot{e} - \cos {\Omega } \bigg [ \dot{z} + \frac{\dot{i}}{i} z \bigg ] \dot{i}, \end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray} &&\frac{d \dot{z}}{dt} \bigg |_{\dot{a}} + \frac{d \dot{z}}{dt} \bigg |_{\dot{e}} + \frac{d \dot{z}}{dt} \bigg |_{\dot{i}}\nonumber \\\ms &&\quad = \bigg [ \frac{- \dot{r}}{2 a} \sin {i} \sin {(\omega + f)} + \frac{ - r \dot{f}}{2 a} \sin {i} \cos {(\omega + f)}\nonumber \\\ms &&\qquad + \frac{r}{a} \dot{i} \cos {i} \sin {(\omega + f)} \bigg ] \dot{a} + \bigg [ \frac{\partial \dot{r}}{\partial e} \sin {i} \sin {(\omega + f)}\nonumber \\\ms &&\qquad + \frac{\partial (r \dot{f})}{\partial e} \sin {i} \cos {(\omega + f)} + \frac{\partial r}{\partial e} \dot{i} \cos {i} \sin {(\omega + f)} \bigg ] \dot{e}\nonumber \\\ms &&\qquad + \bigg [ \dot{r} \cos {i} \sin {(\omega + f)} + r \dot{f} \cos {i} \cos {(\omega + f)}\nonumber \\ &&\qquad + r \frac{\dot{i}}{i} \cos {i} \sin {(\omega + f)} - \dot{i} z \bigg ] \dot{i}. \end{eqnarray} \tag{ A14 }$$

Arbitrary functions for the evolution of the orbital elements can be assigned to the $\dot{a}/a$, $\dot{e}/e$, and $\dot{i}/i$ terms. With the additional terms derived above (Equations (A9)–(A14)), the full motion of a body undergoing arbitrary evolution in semimajor axis, eccentricity, and/or inclination can be implemented. When the eccentricity e reaches a small value, machine round-off error can cause the eccentricity damping terms to damp the eccentricity below 0, leading to failure of the integrator. Therefore, we did not apply the eccentricity damping terms if e < 10⁻⁴. If the orbit of the body precesses, ω is no longer constant (apsidal precession) and/or Ω is no longer constant (nodal precession). Since in this paper we do not consider precession (see Batygin et al. 2011; Dawson & Murray-Clay 2012), we leave the additional modifications to the equations of motion to account for precession for future work.