RETENTION OF A PRIMORDIAL COLD CLASSICAL KUIPER BELT IN AN INSTABILITY-DRIVEN MODEL OF SOLAR SYSTEM FORMATION

The quest to understand the origins of the solar system dates back centuries. The last two decades, however, have seen a renewed interested in the problem, as the discovery of the Kuiper Belt (Jewitt & Luu 1993) has provided important new clues to the physical processes that took place during the early stages of our solar system's evolution. The continued acquisition of new information gave rise to a multitude of new formation models (see Morbidelli et al. 2008 for a comprehensive review). Among the newly proposed scenarios, an instability model, termed the "Nice" model (Tsiganis et al. 2005; Gomes et al. 2005; Morbidelli et al. 2005), has been particularly successful in reproducing the observed properties of planetary orbits and the Kuiper Belt (Levison et al. 2008).

Within the context of the narrative told by the Nice model, planets start out in a multi-resonant configuration (Morbidelli et al. 2007; Batygin & Brown 2010), and, driven by planetesimal scattering, begin migrating divergently (Fernandez & Ip 1984). Eventually, the planets encounter a low-order mean-motion resonance (MMR), which results in a transient period of instability. During this period, the ice giants scatter outward and settle roughly onto their current semimajor axes but with high eccentricities (Tsiganis et al. 2005; Thommes et al. 2008). Neptune's excited eccentricity gives rise to a chaotic sea between its exterior 3:2 and 2:1 MMRs, allowing planetesimals to random walk into the "classical" region (Levison et al. 2008). Subsequently, as the planets circularize due to dynamical friction (Stewart & Wetherill 1988), the scattered and resonant populations of the Kuiper Belt are sculpted.

An outstanding problem within the Nice model lies in the formation of the cold classical population of the Kuiper Belt, which is the central theme of this study. The cold population is distinctive from the rest of the Kuiper Belt in a number of ways. First and foremost, as its name suggests, the orbital distribution is dynamically unexcited. When Neptune scatters planetesimals, it tends to pump up their inclinations to tens of degrees. Yet the cold population resides on nearly co-planar orbits, with inclinations not exceeding ${\sim}5 \deg$ (Brown 2001; Gladman et al. 2008). The eccentricities of the cold population, on average, also tend to be diminished in comparison with the hot population, but the division there is not as apparent. Figure 1 shows the eccentricities of the current aggregate of observed Kuiper Belt objects (KBOs) between 30 and 60 AU. Cold classical objects, whose inclinations are below $5 \deg$, are plotted as black dots, while all other objects with inclinations above $5 \deg$ are plotted as blue dots. Note that the cold population's eccentricity distribution is not monotonic in semimajor axes. Between 42 AU and 45 AU, planar KBOs have roughly isentropic eccentricities. However, low-eccentricity objects progressively disappear beyond 45 AU. We refer to this feature of the Kuiper Belt as the "wedge" (see Figure 1).

Zoom In Zoom Out Reset image size

A second distinction is the colors of cold classical KBOs. In general, the Kuiper Belt exhibits a vast diversity of colors, from neutral gray to deep red. Within this range, cold classical KBOs readily stand out as clumps of exclusively red material (Trujillo & Brown 2002; Lykawka & Mukai 2005). In a similar manner, the size distribution of the cold population differs significantly from that of the hot classical population (Fraser et al. 2010). Finally, the fraction of binaries present in the cold population is uniquely large (Stephens & Noll 2006). Moreover, it has been shown that the wide binaries of the cold population in particular would have been disrupted by encounters with Neptune (Parker & Kavelaars 2010), and thus must have never been scattered. While it is difficult to interpret each of these observational facts as conclusive evidence for a particular history, their coherence suggests that the cold classicals are a unique population whose dynamical similarity has been maintained through the dramatic evolution of the outer solar system (Morbidelli & Brown 2004).

A number of formation mechanisms for the cold classical population have been suggested. Within the context of a smooth migration scenario (Malhotra 1995; Murray-Clay & Chiang 2005; Hahn & Malhotra 2005), a primordially cold population can in principle escape dynamical excitation. However, other drawbacks of the smooth migration scenario, such as the inability to reproduce secular architecture of the planets and difficulties in forming the hot classical belt, render it unlikely (Morbidelli et al. 2009; Brasser et al. 2009). Levison et al. (2008) advocated a similar emplacement history for the cold classicals as the hot classicals (i.e., via MMR overlap). The cold population that is produced in such simulations, however, is not cold enough and not physically distinct from the hot population. Subsequently, Morbidelli et al. (2008) showed that if a local, cold population is implemented into the orbital solution of Levison et al. (2008), it will have the same orbital distribution as the implanted population after the instability, so the problem remains. Thus, no coherent picture of the formation of the cold population exists.

In this work, we show that in situ formation of the cold population is consistent with an instability model and all observed dynamical properties of the population, including the low inclinations and the wedge (shown in Figure 1), can be formed. The aim here is not to replicate the Kuiper Belt and the orbits of the planets in a complex N-body simulation, but rather to identify the dynamical processes responsible for the sculpting of the region. The plan of our paper is as follows: in Section 2, we construct an analytical model for secular excitation of a primordially unexcited belt, and thus derive the conditions for retention of dynamically cold orbits. Moreover, we show that the now-fossil wedge is a result of a temporary slowdown in orbital precession. In Section 3, we perform self-consistent N-body simulations that confirm our analytical results. Numerical simulations show that while the cold population can remain undisrupted, similar objects immediately interior to the 3:2 MMR get scattered all over the Kuiper Belt. We conclude and discuss our results in Section 4.

Here, we seek to develop a simplified analytical model that describes the long-term interactions between Neptune and an initially dynamically cold population of KBOs, residing between its exterior 3:2 and 2:1 MMRs, during a transient phase of high eccentricity. Prior to the instability, the planets sit in a compact configuration on near-circular orbits. As long as the orbital separation between the planets and KBOs remains large, their mutual interactions are extremely weak, and so the KBOs maintain their dynamically cold orbits. Consequently, this period is unimportant to the problem at hand.

When planet–planet resonance crossing (or some other mechanism) causes the instability, the gain in semimajor axes and acquisition of high eccentricities and inclination of the planets takes place on a timescale that is considerably shorter than Neptune's apsidal precession period (i.e., less than a million years or so). As a result, it can be viewed as instantaneous within the context of a secular approximation. Thus, in an orbit-averaged sense, it is as if Neptune suddenly appears at 30 AU with a high e and i and begins interacting with the KBOs. Since we seek to show that, after the transient phase of high eccentricity, the KBOs can end up on dynamically cold orbits, we must restrict Neptune from penetrating the region beyond 40 AU. This places the maximum eccentricity attainable by Neptune below e_max < (4/3 − 1) = 1/3. This is, however, a weak constraint, since an eccentric, inclined Neptune can still cause large modulations in the eccentricities and inclinations of the KBOs on a secular timescale (Murray & Dermott 1999). Let us now develop a mathematical model for these secular interactions.

We begin by modeling Neptune's evolution. In our model, we take the mass of the cold KBOs to be negligible, so they have no effect on Neptune's orbit (this is not necessarily true, at all times, for other Kuiper Belt populations). The lack of mass in the primordial cold belt is a requirement for our model that brings up concerns about its formation. We shall discuss this in some detail in Section 4. Since we seek to retain the majority of the local population, and we know that the mass of the current cold classical population is much less than that of the Earth, this is a reasonable assumption. The other planets, as well as the massive component of the Kuiper Belt, will cause apsidal and nodal precession of Neptune's orbit, which we write as $g = \left<\dot{\varpi }_N \right>$ and $f= \langle\dot{\Omega }_N\rangle$, respectively. Note that we are only accounting for the average precessions. We express dynamical friction as exponential decay of e and i with constant timescales τ_e and τ_i. These timescales are different, and their numerical values in N-body simulations tend to be of order ∼10⁷ years (Levison et al. 2008). We neglect the modulation of Neptune's e and i by the other planets. In other words, we only retain the free elements.

In terms of complex Poincaré variables (x = eexp (ıϖ), y = iexp (ıΩ)), we can formulate the first-order Lagrange's equations for Neptune as follows:

$$\begin{equation} \frac{d x_n}{d t} = \imath g x_n - \frac{x_n}{\tau _e} \ \ \ \ \frac{d y_n}{d t} = \imath f y_n - \frac{y_n}{\tau _i}. \end{equation} \tag{ 1 }$$

It is trivial to show that these equations admit the solutions

$$\begin{equation} x_n = e_n^0 \exp [(\imath g - 1/\tau _e)t ] \ \ \ \ y_n = i_n^0 \exp [(\imath f - 1/\tau _i)t ], \end{equation} \tag{ 2 }$$

where e⁰_n and i⁰_n are the initial (maximum) eccentricity and inclination of Neptune, respectively, and $\imath = \sqrt{-1}$. Here, we take e⁰_n = 0.25 and $i_n^0=10 \deg$, in accord with results of numerical simulations (Tsiganis et al. 2005; Levison et al. 2008; Batygin & Brown 2010). In our simple model, the secular evolution of the KBOs is dictated entirely by Neptune's evolution. In the spirit of Laplace–Lagrange secular theory, we only retain terms up to second order in eccentricity and inclination in the disturbing function of the KBOs to ensure a decoupled, analytical solution. The resulting first-order Lagrange's equations read (Wu & Goldreich 2002)

$$\begin{equation} \frac{d x_{{\rm KBO}}}{d t} = \imath A x_{{\rm KBO}} + \imath A_{n} x_n \ \ \ \ \frac{d y_{{\rm KBO}}}{d t} = \imath B y_{{\rm KBO}} + \imath B_{n} y_n, \end{equation} \tag{ 3 }$$

where A, A_n, B, and B_n are constants that depend only on the planetary masses and semimajor axes ratios of Neptune to KBOs (e.g., chap. 7 of Murray & Dermott 1999). Note that in the free precession terms, (A, B), the presence of other planets can also be accounted for with ease.

From here, let us focus only on the eccentricity evolution, since the derivation of the inclination evolution follows an identical procedure. Setting the initial orbital state vector of the KBO to zero ($\left[ x, y \right] = \vec{0}$), the solution to the above equation reads

$$\begin{equation} x_{{\rm KBO}}=\frac{e_n^0 \tau _e A_{n}(\exp [\imath A t] - \exp [(\imath g - 1/\tau _e) t]) }{A \tau _e -g \tau _e -\imath }. \end{equation} \tag{ 4 }$$

The controlling parameter in this solution is Neptune's precession, g. Four solutions, for a KBO at a = 45 AU, with different g's are presented in Figure 2. A natural unit of g is the g₈ eigenfrequency of the Laplace–Lagrange secular solution for the solar system, which physically corresponds to Neptune's average precession rate in the current solar system (see Murray & Dermott 1999). Incidentally, the same unit can be used for the nodal recession rate in the inclination solution, since quantitatively g₈ ≈ −f₈ ≈ 065 yr. As can be seen in Figure 2, varying g leads to dramatically different results. In particular, if low eccentricities are to be retained, g must significantly exceed g₈.

Zoom In Zoom Out Reset image size

After a sufficient amount of time, when Neptune's eccentricity has decayed away (i.e., t ≫ τ_e), the second exponential in the numerator of Equation (4) can be neglected. Such a solution represents a precessing KBO with a constant eccentricity. Accordingly, the time dependence of the solution only governs its angular part. Since we are solely interested in the final orbits of the KBOs, we must extract only the radial part of the solution. Let us write the t ≫ τ_e solution as an exponential of an arbitrary number, ξ:

$$\begin{equation} \exp (\xi)=\frac{e_n^0 \tau _e A_{n} \exp [\imath A t] }{A \tau -g \tau -\imath }. \end{equation} \tag{ 5 }$$

Solving for ξ, and complex expanding the logarithm, we have

$$\begin{eqnarray} \xi &=&\ln \bigg(\frac{e_n^0 \tau _e A_{n} }{\sqrt{1+\tau _e^2 (g - A)^2 }}\bigg) \nonumber \\ &&+\imath \arg \bigg(-\frac{e_n^0 \tau _e A_{n} \exp [\imath A t] }{\imath - \tau _e (A-g)}\bigg). \end{eqnarray} \tag{ 6 }$$

The argument of the logarithm in the above equation is the radial part of the complex solution, which corresponds to the final eccentricity of the KBO, with an equivalent expression for the inclination:

$$\begin{eqnarray} e_{{\rm KBO}}^{{\rm final}}&=&\frac{e_n^0 \tau _e A_{n} }{\sqrt{1+\tau _e^2 (g - A)^2 }}, \nonumber \\ i_{{\rm KBO}}^{{\rm final}}&=&\frac{i_n^0 \tau _i B_{n} }{\sqrt{1+\tau _i^2 (f - B)^2 }}. \end{eqnarray} \tag{ 7 }$$

In principle, we could have arrived at the same answer by complex expanding the solution and taking the square root of the sum of the squares of the real and imaginary parts, although the intermediate expressions would have been considerably more messy.

The above equations can be simplified even further by considering their limiting regimes. If the decay timescale is much longer than the beat frequency (g − A, f − B), we can Taylor expand the equations to first order in (1/τ) around zero. The answer then becomes independent of τ:

$$\begin{eqnarray} &&e_{{\rm KBO}}^{{\rm final}} \simeq e_n^0 \frac{ A_{n} }{ g - A } \ \ \ \ i_{{\rm KBO}}^{{\rm final}} \simeq i_n^0 \frac{B_{n} }{f - B }. \end{eqnarray} \tag{ 8 }$$

This procedure is equivalent to assuming that τ²_e(g − A)² ≫ 1 or τ²_i(f − B)² ≫ 1 and throwing away the 1 under the square root in the denominator.¹ It is clear from Figure 2, where the approximate solutions are plotted as big dots, that quantitative agreement with the "full" solution (Equation (4)) is excellent in the parameter regime of interest.

Figures 3 and 4 show the secular excitation of initially cold KBOs' eccentricities and inclinations between the 3:2 and 2:1 MMRs with g = 0, g₈, 2g₈, and 3g₈ as solid lines. These solutions suggest that if one is to retain an eccentricity below e < 0.1 and inclination below $i < 5 \deg$, Neptune's average orbital precession and nodal recession rates must have exceeded ∼3g₈ during the eccentric phase. The enhanced precession is primarily a consequence of Uranus. When Neptune scatters, it arrives somewhat closer to the Sun than its current orbit and migrates to ∼30 AU by scattering KBOs (here, we have implicitly omitted this effect by stating that the coefficients A, A_n, B, and B_n are constant). Thus, at the time of scattering, the semimajor axis ratio of Neptune to Uranus may be lower, leading to an enhanced precession. Additionally, the mass contained in the Kuiper Belt may also play a role in inducing secular precession of Neptune.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The solution described above gives eccentricities that monotonically decrease with semimajor axes. However, as already discussed above, the observed cold population exhibits a somewhat different behavior, with low-eccentricity objects progressively disappearing in the vicinity of the 2:1 MMR. This dynamically unique structure (i.e., the wedge—see Figure 1) is an essential feature to any proposed formation mechanism for the cold classicals.

A wedge-like structure cannot be reproduced by a sweeping 2:1 MMR in an instability-driven formation model. Unlike the smooth migration scenario, where resonant capture is possible (Henrard & Lamaitre 1983; Malhotra 1995; Ketchum et al. 2011), when Neptune is eccentric, the chaotic motion that arises from resonant splitting (Wisdom 1980) leads to an effective randomization of the eccentricities (Quillen 2006). In other words, the KBOs that are temporarily captured do not form a coherent structure such as the wedge. An alternative scenario for formation of the wedge is one where the local population ends at 45 AU, and the wedge is a result of an extended scattered disk with q ∼ 40 AU (Gladman et al. 2008). It is unlikely, however, that, in the extended scattered disk scenario, the low inclinations of scattered objects could be preserved.

Here, we propose the formation of the wedge to be a consequence of secular perturbations. Thus, we seek to modify the above secular solution such that it yields eccentricities that are not monotonically decreasing with semimajor axes in the region of interest. As already described above, the controlling parameter in the secular solution is g. So far, we have kept g constant. However, since Neptune scatters numerous KBOs during its circularization, and the orbits of other planets (particularly Uranus) are changing as well, one would expect Neptune's precession to vary considerably, in a chaotic manner.

It is difficult to predict the exact nature of this variation without a detailed calculation, so here we consider an extreme case as a proof of concept. Namely, we set g = 4g₈ at all times, except τ < t < 1.1τ, where we set g = 0. Note that the precession of Neptune need not necessarily stop. We are choosing g = 0, rather than a diminished precession rate (such as, say g = g₈) merely for the sake of argument. An analytical solution is attainable in a similar fashion as above, by breaking up the integration into three separate time intervals. If g is not held constant throughout Neptune's circularization, the final eccentricity and inclination take on a different character. Qualitatively, this can be understood as follows: when Neptune stops precessing, it starts to induce considerable oscillations in eccentricities of KBOs; however, once the precession becomes rapid again, the modulation stops and the eccentricities become frozen-in. These solutions are plotted as dashed curves in Figures 3 and 4. The details of the non-monotonic solution depend on when and for how long Neptune's precession is halted, and change further if the precession is merely slowed down, rather than stopped. Furthermore, the dashed curves in Figures 3 and 4 shift to larger semimajor axes if the free precessions of the KBOs (A, B) are enhanced. While it is understood that these calculations do not reproduce the cold classical population in detail, they do show that primordially unexcited objects can retain cold orbits in face of dynamical excitation, and coherent structure can be formed in the context of a purely secular solution.

Having motivated in situ formation of the cold classical population with analytical arguments, we now turn to numerical N-body simulations for confirmation of the above results and inclusion of omitted physics (such as close encounters, MMRs, and higher-order secular terms in the disturbing function). In this study, the integrations were performed using the mercury6 integration software package (Chambers 1999) utilizing the "hybrid" algorithm. The disk was composed of two components. The massive planetesimal swarm, containing 3000 particles, resided between the immediate stability boundary of the initial multi-resonant configuration and ∼35 AU. This was followed by a disk of another 3000 massless particles that extended to 60 AU. Thus, we are assuming that a significant density gradient exists, in the vicinity of Neptune's final orbit, such that the mass in the outer disk is insufficient to drive Neptune's migration. However, the numbers of particles were chosen due to considerations of computational cost and are not intended to be representative of the relative fraction of bodies in the planetesimal disk in any way. The initial conditions were drawn from the eight multi-resonant states that were identified by Batygin & Brown (2010) as being compatible with an instability formation model. The planetesimals were initialized on near-coplanar, near-circular orbits (e ∼ sin i ∼ 10⁻³). The self-gravity of the planetesimal swarm was neglected to reduce the computational cost of the experiments, as 30 permutations of each initial condition were integrated.

Batygin & Brown (2010) used the presence of scattering events between an ice giant and a gas giant, followed by a transient phase of high eccentricity, as a proxy for whether successful formation of the classical Kuiper Belt can occur. Further constraints on the initial conditions can be placed by considering the reproduction of the outer solar system's secular eigenmodes. Particular difficulty has been found in ensuring that the amplitude of Jupiter's g₅ mode is larger than that of the g₆ mode (Morbidelli et al. 2009). Having completed all of the integrations, we checked the relative amplitudes of the g₅ and the g₆ modes in all solutions. Surprisingly, we found that despite a transient period of instability and gas giant/ice giant scattering, the (3:2 3:2 4:3)² and the (5:3 4:3 4:3) initial conditions did not reproduce the secular architecture of the planets, in neither this set nor in the set of integrations of Batygin & Brown (2010). If Jupiter and Saturn were indeed initially locked in the 3:2 MMR, as hydrodynamic simulations suggest (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Pierens & Nelson 2008), only the (3:2 3:2 5:4) and (3:2 4:3 4:3) initial conditions are left as viable options for the starting state of the solar system.

As already discussed in Section 2, interactions between the cold outer disk and the outermost ice giant are largely independent of the starting condition, since scattering in a successful simulation always sets the planets onto orbits that are close to that of the current solar system, but with moderate eccentricities. Consequently, we did not restrict our analysis to any particular initial condition. Out of our set of 180 integrations, in eight cases primordially cold objects were able to retain unexcited orbits in addition to the gas-giant eigenmodes being reproduced correctly. Here, we focus on two representative integrations: one starting from the (2:1 4:3 4:3) initial condition (Figure 5) and another starting from the (5:3 4:3 3:2) initial condition (Figure 6). In both cases, the cold classical population is produced, but the wedge is only formed in the simulation that starts from the (2:1 4:3 4:3) initial condition (although it is somewhat smaller than its observed counterpart). Note that the formation of the wedge has little to do with the initial condition—rather, its production is a random process. Similarly, the exact degree of excitation of the cold population's inclinations is sensitively dependent on the details of Neptune's evolution, which is chaotic. Thus, the fact that the wedge is reproduced in one simulation and the degree of excitation of the inclinations is reproduced in another are unrelated results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A vast majority (>90%) of the objects in the cold classical region (i.e., a ≳ 42 AU) are retained in our simulations on stable orbits. On the contrary, only about a few thousandths of the particles in the inner disk are emplaced onto stable orbits in the Kuiper Belt region. This implies that in order to self-consistently study the formation of the Kuiper Belt, N ≫ 3000 is needed. Unfortunately, the required resolution is not computationally feasible. However, the problem can still be addressed by the use of "tracer" simulations, an approach already utilized in the context of Kuiper Belt formation by Levison et al. (2008).

In a tracer simulation the planets and planetesimals are not self-consistently evolved in time. Rather, the evolution of the planets is pre-loaded from a master simulation and the planetesimals, which are treated as test particles, are evolved subject only to gravitational interactions with the planets. At the beginning of a tracer simulation, the tracer disk is initialized to have the same distribution as the massive component of the planetesimal disk. Consequently, at all times during the integration, the tracer particles also have an identical orbital distribution to that of the massive planetesimals. Each simulation was seeded with 100 test particles and integrated on Caltech's PANGU supercomputer. We employed the Bulirsch–Stoer algorithm (Press et al. 1992) in our tracer integrations.

We performed 200 tracer simulations for each of the evolutions presented in Figures 5 and 6. This amounts to evolving a primordial disk of ∼26, 000 particles, including the outer belt. After the ∼30 Myr simulations were completed, approximately 7% of the particles that originated interior to 35 AU had semimajor axes in the range 35 AU <a < 60 AU, shown as pale blue dots in Figures 7–10. We further cloned the populations³ of tracer particles in the Kuiper Belt region to effectively increase the number of implanted hot classical, scattered, and resonant particles by another factor of 6. The resulting Kuiper Belt, including the test particles that originate beyond 35 AU, was then evolved for an additional 500 Myr to ensure that all unstable particles have time to eject. At the end of the 500 Myr, only ∼5% of the implanted objects that were present at the end of the 30 Myr simulations ended up on stable orbits. Consequently, the cumulative fraction of objects that are implanted into the Kuiper Belt from the inner disk is ∼0.3%. The stable objects are shown as dark blue dots in Figures 7–10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that, in our simulations, the resonant populations are considerably diminished in number. This is largely a cost of performing self-consistent simulations with planetesimals that are unrealistically massive. Every time Neptune scatters a KBO, its resonances jump unrealistically far, disturbing the resonant KBOs, leading to their eventual ejection (Murray-Clay & Chiang 2006). In a suite of customized simulations where the instability still occurs, but planets are analytically guided to their final orbits and gravity is softened (Levison et al. 2008), Neptune's MMRs end up overpopulated. This leads one to believe that the true parameter regime of Neptune's migration resided somewhere between what is presented in this work and that of Levison et al. (2008, A. Morbidelli 2011, private communication).

Although both of the integrations presented here produce a cold classical belt, it is immediately apparent that the wedge is only produced in the integration that starts from the (2:1 4:3 4:3) initial condition, although again the process has little do with the choice of initial condition. Furthermore, from Figure 7, it can be readily inferred that the production of the wedge must be a secular effect since the structure in this simulation extends beyond the 2:1 MMR, i.e., the unswept region. Note that owing to the enhanced free precession of the KBOs (due to the presence of a massive Kuiper Belt), the wedge structure is shifted to the right compared with the analytical estimates presented in the previous section.

In the context of these integrations, we are further able to confirm that the formation of the wedge is due to a considerable slowdown in Neptune's precession. During circularization in the integration that starts from a (5:3 4:3 3:2) initial condition, Neptune's precession is always roughly g ≈ 4.7g₈ while it is eccentric. On the contrary, in the integration that starts from the (2:1 4:3 4:3) initial condition, Neptune's precession rate varies considerably (1.6g₈ ≲ g ≲ 3.7g₈) between 23 Myr and 25 Myr (Figure 11). The presence of a mechanism for the successful formation of the wedge from a local population is an important argument for confirmation of the in situ formation of the cold classical population in the context of an instability model.

Zoom In Zoom Out Reset image size

It is noteworthy that, in the results of the simulation, the wedge appears much less coherent at semimajor axes interior to the 2:1 MMR. This is a consequence of eccentric resonant sweeping. Because of Neptune's considerable eccentricity, the KBO multiplet and the Neptune multiplet of the resonance overlap even for small KBO eccentricities. This allows the KBO to randomly explore the phase space occupied by both sections of the resonance. However, as Neptune's eccentricity is monotonically decreasing, so is the phase space volume occupied by Neptune's multiplet of the 2:1 MMR, making capture impossible (Quillen 2006). Moreover, because of different precession rates, the nominal location of Neptune's multiplet of the resonance lags (i.e., smaller semimajor axis) that of the KBO. Thus, if a KBO exits the resonance shortly after it enters, it tends to get transported closer to the Sun, since it enters at the KBO multiplet and exits at the Neptune multiplet. The change in semimajor axes, however, is only the resonant splitting width, so it is rather small (δa < 0.1 AU). This randomization of the orbital elements causes the inner part of the wedge to appear less coherent in Figure 7.

Finally, it is worth noting that although KBOs that become the cold classical population are able to roughly retain their primordial orbital distribution, the objects between 35 AU and 40 AU inevitably get scattered by Neptune during the instability. Indeed, in both simulations presented here, the scattered cold classicals (shown as yellow triangles in Figures 7–10) join the scattered disk as well as the hot classical population, while some particles get temporarily trapped in resonances during their evolution.

The fact that these lifted objects mostly get emplaced onto stable orbits is suggestive that the results of intrusion of inclined populations by cold-classical-like objects, that took place during the instability, should still be observable today. In other words, the in situ formation scenario for cold classicals presented here predicts that a class of objects, occupying the same unique color region as the cold classicals, should be present in the excited populations.

In this paper, we present a self-consistent dynamical model for the evolution of a primordial cold classical population of the Kuiper Belt, in the context of an instability-driven formation scenario for the solar system. We show, from simple analytical considerations, that the cold belt can survive the transient period of dynamical instability, inherent to the planets. In order for a primordially cold population of KBOs to maintain an unexcited state, the average apsidal precession and nodal recession rates of Neptune during the transient phase of instability must have been considerably faster than what is observed in today's solar system. Simultaneously, successful formation of the wedge (see Figure 1) requires that the apsidal precession rate drops by a factor of a few for a short period of time. Numerical integrations presented in this work confirm the results of the analytical calculations and reveal a particular result: the formed cold population and the wedge closely resemble their observed counterparts. The dynamical evolution of cold classicals we propose here is in close agreement with the uniqueness of cold classicals' physical characteristics.

In situ formation of cold classicals brings to light the issue of truncation of the classical belt near the 2:1 MMR. In the chaotic capture mechanism proposed by Levison et al. (2008), the outer edge comes about naturally as the 2:1 MMR sculpts the belt. In our solution, however, a cold belt that extends further out is surely possible. Thus, we are forced to attribute the proximity of the edge and the 2:1 MMR to a mere coincidence. Another question of interest is the fate of primordially cold binaries in the 35–40 AU region. It is likely that many of these binaries will get disrupted by close encounters with Neptune, although the exact fraction will depend on the details of Neptune's evolution. Consequently, an in-depth analysis of the evolution of the scattered cold KBOs may open up an avenue toward further constraining the orbital history of Neptune.

Although in situ formation of cold classicals resolves a pressing dynamical problem within the Nice model, it gives rise to a new issue that requires attention. Namely, the outstanding question of importance is planetesimal formation beyond ∼35 AU, given the steep size distribution of the cold classicals. In other words, how is the formation of planetesimals up to ∼200 km in size accomplished in such a low-density environment?

Although the answer to this question is by no means trivial, one possible solution to this problem lies within the context of streaming instabilities (Youdin & Goodman 2005). Streaming instabilities have already been suggested as the dominant formation process in the cold classical population, as gravitational collapse has been shown to yield wide binaries (Nesvorný et al. 2010). Importantly, in the proposed picture, planetesimal formation is a threshold process that "turns on" only when gas drag accumulates a critical amount of dust in a given location within the solar nebula. Thus, one can in principle envision a system where most of the dust gets carried inward of ∼35 AU by gas drag, but, infrequently, the dust surface density reaches a critical value in the outer nebula, causing a few, but sizable planetesimals to be born. Such a scenario would likely result in a very sharply decreasing surface density profile in the outer nebula, at the epoch of disappearance of the gas. As a result, this picture would imply the existence of a steep density gradient in the primordial planetesimal disk, such as the one we require in our model, consistent with preventing Neptune's extended migration.

Another possibility for the formation process is hierarchical coagulation, where planetesimal growth is accomplished by collisions among smaller objects in a quiescent environment (Kenyon 2002). In particular, it has been suggested that, if aided by turbulent concentration, hierarchical coagulation could yield the desired mass of the cold classical population (Cuzzi et al. 2010). In fact, even if the original mass of the cold belt exceeded its current value, erosion by collisional grinding could in principle be invoked to reduce the overall mass. However, this process may prove problematic in reproducing the observed wide binary fraction of the cold belt (Nesvorný et al. 2011).

Whatever the formation process for the cold classicals is, the results presented here have considerable implications. First and foremost, the successful retention of the cold classical population in the context of an instability-driven model fixes the most significant drawback of the Nice model. Second, our scenario suggests that the cold classical Kuiper Belt is the only population of objects in the outer solar system that has not been transported away from its formation site. Furthermore, assuming that collisional grinding has played a negligible role in the cold population's evolution (as suggested by the observed binary fraction; Nesvorný et al. 2011), the cold classical population essentially yields the surface density of the solar nebula at a ∼ 45 AU, since the majority of the KBOs are retained in place. This potentially makes the cold classicals a unique laboratory for the study of surface processes as well as the chemistry of the primordial solar nebula. Third, based upon the results of the numerical simulations, we expect that objects that are physically similar to the cold classicals should be scattered throughout the Kuiper Belt and as a result may explain the spectral similarity between cold classicals and some objects at higher inclinations.

In conclusion, it appears that quantitative evaluation of planetesimal formation beyond ∼35 AU is required to draw a complete picture of the in situ formation and evolution scenario for the cold classicals. However, the considerable improvement of the model for early dynamical evolution of the Kuiper Belt presented here supports the overall validity of the hypothesis.

We thank Alessandro Morbidelli, Hal Levison, Darin Ragozzine, and Peter Goldreich for useful conversations.