NEPTUNE'S ORBITAL MIGRATION WAS GRAINY, NOT SMOOTH

Our numerical integrations consist of tracking the orbits of four giant planets (Jupiter to Neptune) and a large number of test particles representing the outer planetesimal disk. To set up an integration, Jupiter and Saturn were placed on their current orbits. Uranus and Neptune were placed on inside of their current orbits and were migrated outward. The initial semimajor axis ${a}_{{\rm{N,0}}}$, eccentricity ${e}_{{\rm{N,0}}}$, and inclination ${i}_{{\rm{N,0}}}$ define Neptune's orbit before the main stage of migration/instability. In most of our simulations we used ${a}_{{\rm{N,0}}}=24\;{\rm{au}}$, because the wide inclination distribution of HCs and resonant populations requires that Neptune's migration was long range (a_N,0 ≲ 25 au; Nesvorný 2015a). We also set ${e}_{{\rm{N,0}}}=0$ and ${i}_{{\rm{N,0}}}=0$. All inclination values reported in this article are referred to the invariant plane of the solar system.

The swift_rmvs4 code (Levison & Duncan 1994) was used to follow the evolution of planets (and massless disk particles; see Section 2.2). The code was modified to include fictitious forces that mimic the radial migration and damping of planetary orbits. These forces were parametrized by the exponential e-folding timescales, ${\tau }_{a}$, ${\tau }_{e}$, and ${\tau }_{i}$, where ${\tau }_{a}$ controls the radial migration rate, and ${\tau }_{e}$ and ${\tau }_{i}$ control the damping rate of e and i. Specifically, the semimajor axis of Neptune changes from ${a}_{{\rm{N,0}}}$ to its current average of ${a}_{{\rm{N,c}}}=30.11\;{\rm{au}}$ as

$$\begin{eqnarray}&&{a}_{{\rm{N}}}(t)={a}_{{\rm{N,c}}}+({a}_{{\rm{N,0}}}-{a}_{{\rm{N,c}}})\mathrm{exp}(-t/{\tau }_{a}),\end{eqnarray} \tag{ 1 }$$

and the eccentricity and inclination of Neptune evolve according to

$$\begin{eqnarray}&&{e}_{{\rm{N}}}(t)={e}_{{\rm{N,0}}}\mathrm{exp}(-t/{\tau }_{e})\ {\rm{and}}\ {i}_{{\rm{N}}}(t)={i}_{{\rm{N,0}}}\mathrm{exp}(-t/{\tau }_{i}).\end{eqnarray} \tag{ 2 }$$

The expressions for ${e}_{{\rm{N}}}(t)$ and ${i}_{{\rm{N}}}(t)$ differ from those used in Morbidelli et al. (2014), where the damping rate $({{de}}_{{\rm{N}}}/{dt})/{e}_{{\rm{N}}}$ was chosen to be proportional to $\mathrm{exp}(-t/{\tau }_{i})$. Here we set ${\tau }_{a}\sim {\tau }_{e}\sim {\tau }_{i}$ $(={\tau }_{1})$, because such roughly comparable timescales were suggested by previous work.

The numerical integrations of the first migration stage were stopped when Neptune reached ${a}_{{\rm{N,1}}}\simeq 28\;{\rm{au}}$. Then, to approximate the effect of planetary encounters during the instability (NM12, Nesvorný 2015b), we applied a discontinuous change of Neptune's semimajor axis and eccentricity, ${\rm{\Delta }}{a}_{{\rm{N}}}$ and ${\rm{\Delta }}{e}_{{\rm{N}}}$. Motivated by the NM12 results (see Figure 3), we set ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ or 0.5 au, and ${\rm{\Delta }}{e}_{{\rm{N}}}=0$, 0.05 or 0.1. The purpose of ${\rm{\Delta }}{a}_{{\rm{N}}}={\rm{\Delta }}{e}_{{\rm{N}}}=0$ is to have a reference case for comparison purposes. We use ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$, because Nesvorný (2015b) showed that this jump size would be needed to explain the Kuiper Belt kernel. Note that the resonant objects are released from resonances with ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$, because the typical resonance width is just smaller than the jump size. No change was applied to the orbital inclination of Neptune, because a small inclination change should not critically affect the processes studied here.

The second migration stage starts with Neptune having the semimajor axis ${a}_{{\rm{N,2}}}={a}_{{\rm{N,1}}}+{\rm{\Delta }}{a}_{{\rm{N}}}$. We apply the swift_rmvs4 code, and migrate the semimajor axis (and damp the eccentricity) on an e-folding timescale ${\tau }_{2}$. The migration amplitude was adjusted such that the final semimajor axis of Neptune ended to be within 0.05 au of its current mean ${a}_{{\rm{N,c}}}=30.11\;{\rm{au}}$, and the orbital period ratio, ${P}_{{\rm{N}}}/{P}_{{\rm{U}}}$, where ${P}_{{\rm{N}}}$ and ${P}_{{\rm{U}}}$ are the orbital periods of Neptune and Uranus, ended within 0.5% of its current value (${P}_{{\rm{N}}}/{P}_{{\rm{U}}}=1.96$). A strict control over the final orbits of planets is important, because it guarantees that the mean motion and secular resonances reach their present positions.

As for the specific values of ${\tau }_{1}$ and ${\tau }_{2}$ used in our model, we found from NM12 that the orbital behavior of Neptune can be approximated by ${\tau }_{1}\simeq 10\;{\rm{Myr}}$ and ${\tau }_{2}\simeq 30\;{\rm{Myr}}$ for a disk mass ${M}_{{\rm{disk}}}=20$ ${M}_{{\rm{Earth}}}$, and ${\tau }_{1}\simeq 20\;{\rm{Myr}}$ and ${\tau }_{2}\simeq 50\;{\rm{Myr}}$ for ${M}_{{\rm{disk}}}=15$ ${M}_{{\rm{Earth}}}$. Slower migration rates are possible for lower disk masses. Moreover, we found from NM12 that the real migration is not precisely exponential with the effective ${\tau }_{2}$ being longer than the values quoted above during the very late migration stages (τ₂ ≳ 100 Myr). This is consistent with constraints from Saturn's obliquity, which was presumably exited by late capture in a spin–orbit resonance (Vokrouhlický & Nesvorný 2015; see Hamilton & Ward 2004 and Ward & Hamilton 2004 for the original work that proposed the spin–orbit resonance as the means of exciting Saturn's obliquity). Much shorter migration timescales than those quoted above do not apply, because they would violate constraints from the wide inclination distribution of HCs and resonant populations (Nesvorný 2015a). We therefore used ${\tau }_{1}=10\;{\rm{Myr}}$ or 30 Myr, and ${\tau }_{2}=30$ or 100 Myr. These cases should bracket the range of possible migration timescales.

The planetesimal disk was divided into two parts. The inner part of the disk, from just outside Neptune's initial orbit to ${r}_{{\rm{edge}}}$, was assumed to be massive. We used ${r}_{{\rm{edge}}}=28\;{\rm{au}}$ or 30 au, because our previous simulations in NM12 showed that the massive disk's edge must be at 28–30 au for Neptune to stop at $\simeq 30\;{\rm{au}}$. If the edge of the massive disk were at >30 au, Neptune would continue migrating past 30 au (Gomes et al. 2004). The solar nebula could have become truncated, for example, by photoevaporation from the UV and FUV irradiation by background stars in a cluster (e.g., Adams 2010; see also discussion in Petit et al. 2011). In fact, a recent study of the dynamics of planetesimals embedded in a gas disk suggested that the solar nebula was truncated (or else it would act to produce very high orbital inclinations, $i\gt 40^\circ$, in the Kuiper Belt; Kretke et al. 2012). The estimated mass of the planetesimal disk below 30 au is ${M}_{{\rm{disk}}}\simeq 20$ ${M}_{{\rm{Earth}}}$ (NM12). As shown in Levison et al. (2008), the massive disk is the main source of HCs, Plutinos, and other resonant populations. It therefore has a crucial importance for the resonance overpopulation problem considered here.

The planetesimal disk probably had a low mass extension reaching from 30 au to at least $\simeq 45\;{\rm{au}}$. The low mass extension of the disk beyond 30 au is presumably the source of CCs (Batygin et al. 2011; Wolff et al. 2012; Nesvorný 2015b). It is needed to explain why the CCs have several unique physical and orbital properties (see Section 3.4). The disk extension should not substantially contribute to the present populations of the hot and resonant KBOs (Nesvorný 2015b)⁴ , because the orbital inclinations of bodies native to $a\gt 40\;{\rm{au}}$ remain small during Neptune's migration. Here we therefore do not initially consider the disk extension, and return to it only in Section 3.4, where we test whether a grainy migration is consistent with the formation of the Kuiper Belt kernel.

Each of our simulations included one million disk particles distributed from outside Neptune's initial orbit to ${r}_{{\rm{edge}}}$. The radial profile was set such that the disk surface density ${\rm{\Sigma }}\propto 1/r$, where r is the heliocentric distance. A large number of disk particles was needed because the implantation probability in different parts of the Kuiper Belt is expected to be ∼10⁻³–10⁻⁴. With 10⁶ disk particles initially, this yields ∼100–1000 implanted particles, and allows us to perform a detailed comparison of the model results with observations. The disk particles were assumed to be massless such that their gravity does not interfere with the migration/damping routines. This means that the precession frequencies of planets are not affected by the disk in our simulations, while in reality they were. This is an important approximation (Batygin et al. 2011). The direct gravitational effects of the fifth planet on the disk planetesimals were ignored (see the discussion at the end of Section 4). These effects could be especially important for the CCs (Batygin et al. 2012).

An additional uncertain parameter concerns the dynamical structure of the planetesimal disk. It is typically assumed that the disk was dynamically cold with orbital eccentricities $e\lesssim 0.1$ and orbital inclinations $i\lesssim 10^\circ$. Some dynamical excitation could have been supplied by scattering off of Pluto-sized and larger objects that presumably formed in the disk (Stern & Colwell 1997; Kenyon et al. 2008). The magnitude of the initial excitation is uncertain, because it depends on several unknown parameters (e.g., the number of large objects in the disk). The initial eccentricities and initial inclinations of disk particles in our simulations were distributed according to the Rayleigh distribution with ${\sigma }_{e}=0.1$ and ${\sigma }_{i}=0.05$, where σ is the usual scale parameter of the Rayleigh distribution (the mean of the Rayleigh distribution is equal to $\sqrt{\pi /2}\sigma$).

We developed an analytic method to represent the jitter that Neptune's semimajor axis experiences due to close encounters with massive planetesimals. The method has the flexibility to use any smooth migration history of Neptune on the input, include a certain number of the massive planetesimals in the original disk, and generate a new migration history where the random element due to massive planetesimal encounters is included. This approach is useful, because we can easily control how grainy the migration is, while preserving the global orbital evolution of planets from the smooth simulations. We then proceed to test how the simulation results depend on various parameters, such as the number and mass of the massive planetesimals in the original disk.

We start with a specific migration run in which Neptune's semimajor axis evolves smoothly, except for a possible jump by ${\rm{\Delta }}{a}_{{\rm{N}}}$ due to an encounter with another planet. The migration parameters, namely ${a}_{{\rm{N,0}}}$, ${a}_{{\rm{N,1}}}$, ${\rm{\Delta }}{a}_{{\rm{N}}}$, ${\tau }_{1}$, and ${\tau }_{2}$ (Section 2.1), are specified at this point. As mentioned above, each run also includes 10⁶ test particles that represent the disk planetesimals. We first scan through the simulation output in small steps ${\rm{\Delta }}t$, and extract the orbit of Neptune and the orbital distribution of disk planetesimals at each step. We then apply the Öpik-type collision probability code (Bottke et al. 1994; see also Greenberg 1982) to calculate how many encounters between Neptune and planetesimals happen for encounter distances $r\lt R$, where R is some threshold. The gravitational focusing by Neptune is taken into account in this calculation.

The goal is to find how the number of encounters with $r\lt R$ depends on R. We find that for small values of R this dependence is linear (while a quadratic dependence would be expected without gravitational focusing). This can be understood from the following expression for the impact parameter:

$$\begin{eqnarray}&&{b}_{{\rm{max}}}^{2}={R}^{2}\left[1+\displaystyle \frac{{R}_{{\rm{N}}}}{R}{\left(\displaystyle \frac{{v}_{{\rm{esc}}}}{{v}_{\infty }}\right)}^{2}\right]\end{eqnarray} \tag{ 3 }$$

(e.g., Bertotti et al. 2003), where ${R}_{{\rm{N}}}$ = 24,622 km is Neptune's mean radius, v_esc ≃ 23.5 km s⁻¹ is the escape velocity from Neptune's cloudtops, and ${v}_{\infty }$ is the encounter speed "at infinity." Parameter ${b}_{{\rm{max}}}(R)$ is the maximal impact parameter value for which the minimal encounter distance is lower than specified R, when ${v}_{\infty }$ is fixed. Since ${v}_{\infty }$ ≃ 1–2 km s⁻¹, and thus ${v}_{\infty }\ll {v}_{{\rm{esc}}}$, the second term in Equation (3) is greater than 1 for all encounters with $R\lt {R}^{* }$, where ${R}^{* }={R}_{{\rm{N}}}{({v}_{{\rm{esc}}}/{v}_{\infty })}^{2}$. The number of encounters with $r\lt R$ in this regime is therefore proportional to R. For $R\gt {R}^{* }$, on the other hand, the first term in Equation (3) prevails, and the number of encounters with $r\lt R$ becomes proportional to R². In practice, we find it satisfactory to neglect the effect of distant encounters, because the distant encounters do not (individually) induce any significant changes of Neptune's semimajor axis. We therefore only consider encounters with $r\lt {R}^{* }$, where the scaling is linear. Note that ${R}^{* }\gt 140\ {R}_{{\rm{N}}}$ for ${v}_{\infty }\lt 2$ km s⁻¹.

The Öpik code gives us the number of planetesimals having encounters with Neptune, $n(R,t;{\rm{\Delta }}t)$, as a function of the distance R, time t, and time interval ${\rm{\Delta }}t$. Obviously, $n(R,t;{\rm{\Delta }}t)\propto {\rm{\Delta }}t$ for the small intervals used here (Δt = 10³ years). The time profile of the number of encounters depends on the timescale of the planetesimal disk dispersal, which in turn is related to Neptune's migration speed. We find from our simulations that $n(R,t;{\rm{\Delta }}t)\propto R\;{t}^{-\alpha }$ with some exponent α. For example, in the simulation with ${\tau }_{1}=30$ Myr and τ₂ = 100 Myr, we have α ≃ 1.1–1.2. Then, if we approximate $n(R,t;{\rm{\Delta }}t)=C\;R\;{t}^{-\alpha }$, the constant C will depend on the initial number of massive planetesimals in the disk (${N}_{{\rm{mp}}}$). Here we consider ${N}_{{\rm{mp}}}=1000$, 2000, and 4000, and rescale the original $n(R,t;{\rm{\Delta }}t)$ obtained with 10⁶ disk particles to ${N}_{{\rm{mp}}}$. This gives us an approximate record of the history of close encounters between the massive planetesimals and Neptune.

Equipped with the calibrated $n(R,t;{\rm{\Delta }}t)$ function, we sequentially consider each ${\rm{\Delta }}t$ interval. With ${\rm{\Delta }}t={10}^{3}\;{\rm{years}}$, $n(R,t;{\rm{\Delta }}t)$ is always much smaller than unity. In each step, we generate a random variate X with a uniform distribution between 0 and 1. We then set $n(R,t;{\rm{\Delta }}t)=X$ and solve for $R=R(t,X;{\rm{\Delta }}t)$. If ${R}_{{\rm{N}}}\lt R\lt {R}_{{\rm{max}}}$, where ${R}_{{\rm{max}}}\lt {R}^{* }$ is the maximum encounter distance considered here (in practice we set ${R}_{{\rm{max}}}={R}_{{\rm{Hill}}}/10$, where ${R}_{{\rm{Hill}}}=$ is the Hill radius), the encounter is recorded, and we proceed with the next timestep ${\rm{\Delta }}t$. In the end, the method allows us to generate an encounter sequence that closely approximates the reality. Note that for the typical encounter speeds considered here it only takes up to several months to cross ${R}_{{\rm{max}}}$.

We proceed by computing the effect of individual encounters on the semimajor axis of Neptune. This is done using a hyperbolic approximation of the planetesimal's trajectory relative to Neptune. The hyperbolic approximation is adequate for deep encounters considered here, and in a regime when the encounter duration is much shorter than the orbital period of Neptune. The deflection angle θ of the asymptotes of the hyperbola describing the planetesimal encounter trajectory can be obtained from R and ${v}_{\infty }$. Specifically, introducing the half-angle ${\theta }_{1/2}=\tfrac{1}{2}\;\theta$, we find (e.g., Bertotti et al. 2003)

$$\begin{eqnarray}&&\mathrm{sin}{\theta }_{1/2}={\left[1+2\displaystyle \frac{R}{{R}_{{\rm{N}}}}{\left(\displaystyle \frac{{v}_{\infty }}{{v}_{{\rm{esc}}}}\right)}^{2}\right]}^{-1}.\end{eqnarray} \tag{ 4 }$$

Expressed in the inertial frame, the change of Neptune's velocity vector is

$$\begin{eqnarray}&&\delta {\boldsymbol{V}}\simeq 2\;\displaystyle \frac{m}{{M}_{{\rm{N}}}}\;{v}_{\infty }\mathrm{sin}{\theta }_{1/2}[({{\boldsymbol{e}}}_{1}\times {{\boldsymbol{e}}}_{2})\mathrm{cos}{\theta }_{1/2}-{{\boldsymbol{e}}}_{2}\;\mathrm{sin}{\theta }_{1/2}],\end{eqnarray} \tag{ 5 }$$

where ${M}_{{\rm{N}}}$ and m are Neptune's and planetesimal masses, ${{\boldsymbol{e}}}_{1}$ is a unit vector along the angular momentum vector of planetesimal's planetocentric trajectory, and ${{\boldsymbol{e}}}_{2}$ is a unit vector along the incoming trajectory of the planetesimal (relative to the planet). Since we do not propagate the information about ${{\boldsymbol{e}}}_{1}$ and ${{\boldsymbol{e}}}_{2}$ from the original simulation, we assume that ${{\boldsymbol{e}}}_{1}$ and ${{\boldsymbol{e}}}_{2}$ are randomly (isotropically) oriented in space (obviously, these vectors are perpendicular to each other). This should be a reasonable approximation in a situation when the disk of planetesimals is dispersed by Neptune. On the other hand, the actual distribution of encounters is likely to be, to some degree, anisotropic, because Neptune migrates outward due to these encounters. A study of the dynamical effects of anisotropic encounters on the Kuiper Belt goes beyond the scope of this paper.

Finally, we use the Gauss equations to compute the orbital effect of $\delta {\boldsymbol{V}}$ on Neptune's orbit. The eccentricity and inclination changes are neglected. In the approximation of a near-circular orbit of Neptune, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\delta {a}_{{\rm{N}}}}{{a}_{{\rm{N}}}}\simeq 2\;\displaystyle \frac{\delta {V}_{\parallel }}{{V}_{{\rm{N}}}},\end{eqnarray} \tag{ 6 }$$

where $\delta {a}_{{\rm{N}}}$ is the change of Neptune's semimajor axis, $\delta {V}_{\parallel }$ is a projection of $\delta {\boldsymbol{V}}$ onto the direction of Neptune's orbital motion, and ${V}_{{\rm{N}}}$ is Neptune's orbital speed. To compute $\delta {V}_{\parallel }$ from $\delta {\boldsymbol{V}}$, we assume that $\delta {\boldsymbol{V}}$ has a random orientation. The change $\delta {a}_{{\rm{N}}}$ is computed for all encounters with massive planetesimals and recorded in a file. In a simulation, we then use the modified swift_rmvs4 code with smooth migration (Section 2.1), and apply $\delta {a}_{{\rm{N}}}$ every time that is recorded in the encounter file.

Figure 4 shows an example of the sequence of $\delta {a}_{{\rm{N}}}$ in the case with ${\tau }_{1}=30\;{\rm{Myr}}$, ${\tau }_{2}=100\;{\rm{Myr}}$, and ${N}_{{\rm{mp}}}=1000$, where each massive planetesimal was assumed to have Pluto's mass. The total number of encounters recorded in this case is ∼10⁴, implying that a massive planetesimal had on average ∼10 encounters with $r\lt {R}_{{\rm{max}}}={R}_{{\rm{Hill}}}/10$. As expected, most encounters happen during the initial migration stages when the planetesimal disk is still massive. The slight preference for a somewhat larger range of the $\delta {a}_{{\rm{N}}}$ values at late epochs reflects the outward migration of Neptune's orbit. The root mean square of $\delta {a}_{{\rm{N}}}/{a}_{{\rm{N}}}$ is of the order of $\sim (m/M)({v}_{\infty }/{V}_{{\rm{N}}})\sim {10}^{-4}$, as expected from Equations (5) and (6). By design, the $\delta {a}_{{\rm{N}}}$ distribution has a zero mean, while in reality the distribution should be skewed toward positive values, because the massive planetesimals contribute to Neptune's outward migration.

Zoom In Zoom Out Reset image size

Figure 5 shows an example of grainy migration of Neptune produced by the method described in this section. The effects of encounters are clearly visible in panel (b), where the ratio of orbital periods, ${P}_{{\rm{N}}}/{P}_{{\rm{U}}}$, shows an irregular pattern. During the late stages, when the smooth migration nearly stalls, the random effect of encounters with massive planetesimals can cause Neptune to move slightly inward during some time intervals. This may be important for Neptune Trojans, because their stability sensitively depends on the maximal ${P}_{{\rm{N}}}/{P}_{{\rm{U}}}$ reached during the planetary migration (Gomes & Nesvorný 2016).

Zoom In Zoom Out Reset image size

We used the CFEPS detection simulator (Kavelaars et al. 2009) to compare the orbital distributions obtained in our simulations with observations. CFEPS is one of the largest Kuiper Belt surveys with published characterization (currently 169 objects; Petit et al. 2011). The simulator was developed by the CFEPS team to aid the interpretation of their observations. Given intrinsic orbital and magnitude distributions, the CFEPS simulator returns a sample of objects that would have been detected by the survey, accounting for flux biases, pointing history, rate cuts, and object leakage (Kavelaars et al. 2009). In the present work, we input our model populations in the simulator to compute the detection statistics. We then compare the orbital distribution of the detected objects with the actual CFEPS detections using the Kolmogorov—Smirnov (K–S) test (Press et al. 1992).

This is done as follows. The CFEPS simulator takes as an input: (1) the orbital element distribution from our numerical model, and (2) an assumed absolute magnitude (H) distribution. As for (1), the input orbital distribution was produced by a short integration starting from the final model state of the Kuiper Belt. The orbital elements of each object were recorded at 100 yr intervals during this integration until the total number of the recorded data points reached ≃10⁵. Each data point was then treated as an independent observational target. We rotated the reference system such that the orbital phase of Neptune in each time frame corresponded to its ecliptic coordinates at the epoch of CFEPS observations. This procedure guaranteed that the sky positions of bodies in Neptune's resonances were correctly distributed relative to the pointing direction of the CFEPS frames.

The magnitude distribution was taken from Fraser et al. (2014). It was assumed to be described by a broken power law with $N(H){dH}\quad =\quad {10}^{{\alpha }_{1}(H-{H}_{0})}\;{dH}$ for $H\lt {H}_{{\rm{B}}}$ and $N(H){dH}\quad =\quad {10}^{{\alpha }_{2}(H-{H}_{0})+({\alpha }_{1}-{\alpha }_{2})({H}_{{\rm{B}}}-{H}_{0})}\;{dH}$ for $H\gt {H}_{{\rm{B}}}$, where ${\alpha }_{1}$ and ${\alpha }_{2}$ are the power-law slopes for objects brighter and fainter than the transition, or break magnitude ${H}_{{\rm{B}}}$, and H₀ is a normalization constant. Fraser et al. (2014) found that ${\alpha }_{1}\simeq 0.9$, ${\alpha }_{2}\simeq 0.2$, and ${H}_{{\rm{B}}}\simeq 8.0$ (r band) for the HCs. In the context of a model where the HCs formed at <30 au, and were implanted into the Kuiper Belt by a size-independent process (our integrations do not have any size-dependent component), the HC magnitude distribution should be shared by all populations that originated from <30 au (Morbidelli et al. 2009; Fraser et al. 2014). We varied the parameters of the input magnitude distribution to understand the sensitivity of the results to various assumptions. We found that small variations of ${\alpha }_{1}$, ${\alpha }_{2}$, and ${H}_{{\rm{B}}}$ within the uncertainties given in Fraser et al. (2014) have essentially no effect. Note that because we compare our results with the CFEPS survey, we work with the absolute magnitudes in the g band (${H}_{{\rm{g}}}\simeq {H}_{{\rm{r}}}+0.6$).

We first discuss whether Neptune's jump, as suggested in Nesvorný (2015a), can help to resolve the resonance overpopulation problem. Figure 6 shows the final distribution of orbits implanted into the Kuiper Belt in case 1 with smooth migration. For ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ (i.e., no jump) we obtain ${N}_{{\rm{HC}}}/{N}_{3:2}=0.14$, while for ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$ we find ${N}_{{\rm{HC}}}/{N}_{3:2}=0.35$ (Table 1). Thus, the ratio increased by a factor of 2.5 when Neptune's jump was accounted for in the model. The main difference between these two cases is that the probability of capture on a stable orbit in the 3:2 resonance, ${P}_{3:2}$, is ${P}_{3:2}=2.0\times {10}^{-3}$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ and ${P}_{3:2}=6.8\times {10}^{-4}$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$. (The probabilities are normalized here to one particle in the original disk below 30 au.)

Zoom In Zoom Out Reset image size

Table 1.  The Capture Statistics of the HCs and Plutinos Obtained in Different Dynamical Models

	${P}_{{\rm{HC}}}$	${P}_{3:2}$	${N}_{{\rm{HC}}}/{N}_{3:2}$
	(×10−4)	(×10−4)
Smooth Migration
T30	1.9	5.3	0.36
C1-0.0	2.8	20	0.14
C1-0.5	2.4	6.8	0.35
C2-0.0	5.0	11	0.45
C2-0.5	6.6	10	0.66
Grainy Migration
C1-0.0-1000P	4.6	8.9	0.52
C1-0.5-1000P	4.2	3.2	1.3
C1-0.0-2000P	5.2	3.9	1.3
C1-0.5-2000P	5.7	3.2	1.8
C1-0.0-4000P	6.6	2.1	3.1
C1-0.5-4000P	6.2	1.9	3.3
C1-0.0-1000P2	5.5	1.3	4.2
C1-0.5-1000P2	6.0	1.6	3.8
C2-0.0-1000P	8.5	4.7	1.8
C2-0.5-1000P	9.2	3.7	2.5
C2-0.0-2000P	11	1.9	5.8
C2-0.5-2000P	13	1.9	6.8
Observations
CFEPS/DES	$\simeq 5$	$\simeq 1.5$	≃2–4

Note. T30 is a case from Nesvorný (2015a), where Neptune started at ${a}_{{\rm{N,0}}}=24\;{\rm{au}}$ and smoothly migrated to 30 au with an e-folding timescale of $\tau =30\;{\rm{Myr}}$. C1 stands for Case 1 with ${\tau }_{1}=30\;{\rm{Myr}}$ and ${\tau }_{2}=100\;{\rm{Myr}}$, C2 stands for case 2 with ${\tau }_{1}=10\;{\rm{Myr}}$ and ${\tau }_{2}=30\;{\rm{Myr}}$. Labels 0.0 and 0.5 denote the cases with ${\rm{\Delta }}{a}_{{\rm{N}}}=0.0$ and ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$. The simulations with 1000, 2000, and 4000 Plutos are labeled by 1000P, 2000P and 4000P, Respectively. The Case with 1000 Twoplutos is denoted by 1000P2. The Columns give the probability of capture as HC (${P}_{{\rm{HC}}}$) and in the 3:2 resonance (${P}_{3:2}$), and the ratio between the two populations (${N}_{{\rm{HC}}}/{N}_{3:2}$). The last row lists observational contraints. The ${P}_{{\rm{HC}}}$ value reported in this row was computed from the estimated mass of the HCs, ∼0.01 ${M}_{{\rm{Earth}}}$ according to Fraser et al. (2014). Assuming a ${M}_{{\rm{disk}}}=20$ ${M}_{{\rm{Earth}}}$ disk from NM12, this gives ${P}_{{\rm{HC}}}\simeq 5\times {10}^{-4}$. According to the CFEPS and DES the population of Plutinos in the 3:2 resonance is ≃2–4 smaller than the HCs. Thus, ${P}_{{\rm{HC}}}\simeq 1.5\times {10}^{-4}$.

Download table as:  ASCIITypeset image

We looked into this issue in detail and found that the change of ${P}_{3:2}$ was mainly contributed by Neptune's jump. As a consequence of Neptune's jump, the 3:2 resonant objects captured during the previous stage were released from the resonance. Interestingly, however, many bodies were captured into the 3:2 resonance from the scattered disk immediately after Neptune's jump, when the 3:2 resonance suddenly moved into a new orbital location, and during the subsequent slow migration of Neptune. This explains why the ${N}_{{\rm{HC}}}/{N}_{3:2}$ ratio did not change more substantially. The probability of capture in the main belt, ${P}_{{\rm{HC}}}$, also changed, but the change was minor (${P}_{{\rm{HC}}}=2.8\times {10}^{-4}$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ and ${P}_{{\rm{HC}}}=2.4\times {10}^{-4}$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}};$ Table 1).

A similar result holds for case 2 with smooth migration, where ${N}_{{\rm{HC}}}/{N}_{3:2}=0.45$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ and ${N}_{{\rm{HC}}}/{N}_{3:2}=0.66$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$. These values are somewhat higher that those obtained in case 1 perhaps suggesting that ${N}_{{\rm{HC}}}/{N}_{3:2}$ might increase further if even shorter migration timescales were used. The short migration timescales are not plausible, however, because they do not satisfy the inclination constraint (Nesvorný 2015a). Given these results, we conclude that the effect of Neptune's jump can help, but it is insufficient in itself to resolve the resonance overpopulation problem. This is because even in the smooth migration cases with ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$, the 3:2 resonance is still strongly overpopulated, by a factor of ∼5–10 relative to HCs. Other resonances, such as the 2:1 or 5:2, are overpopulated by a significant factor as well. We therefore proceed by considering the cases with grainy migration.

Figure 7 shows some of our best results for the grainy migration. These results were obtained for case 1 (${\tau }_{1}=30\;{\rm{Myr}}$ and τ₂ = 100 Myr) and 1000 Twoplutos. Here, ${N}_{{\rm{HC}}}/{N}_{3:2}\simeq 4$ for both ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ and ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$, in a close match to CFEPS observations. This is encouraging. Neptune's jump does not seem to have much to do with this result, because ${N}_{{\rm{HC}}}/{N}_{3:2}\simeq 4.2$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ and ${N}_{{\rm{HC}}}/{N}_{3:2}\simeq 3.8$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$. These two values are similar and the small difference between them probably reflects some minor difference in the planetary migration histories. Neptune's jump has only a small effect here, because the case with grainy migration encourages late captures with many bodies being captured after the instability.

Zoom In Zoom Out Reset image size

The results for the case-1 grainy migration with 4000 Plutos are similar to those reported above for 1000 Twoplutos. The ones obtained for 1000 and 2000 Plutos are intermediate, showing a clear progression from the smooth migration case to cases with an increased migration graininess (Table 1). The main trend is that ${P}_{3:2}$ drops when more and more Plutos are included. For case 1 with 1000 Plutos, ${P}_{3:2}=8.9\times {10}^{-4}$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ and ${P}_{3:2}=3.2\times {10}^{-4}$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$. With 2000 or 4000 Plutos, or 1000 Twoplutos, P_3:2 ≃ (1.5–4) × 10⁻⁴, which represents only a small fraction of ${P}_{{\rm{3:2}}}$ obtained with the smooth migration. This change is attributed to the jittery motion of the 3:2 resonance that accompanies Neptune's grainy migration (Murray-Clay & Chiang 2006). Due to this jitter, many bodies captured into the resonance during the earlier stages were later released from the resonance, because their libration amplitude increased beyond the limits required for the orbital stability (${A}_{\sigma }\simeq 130^\circ ;$ Nesvorný & Roig 2001). The distribution of the libration amplitudes in the 3:2 resonance is discussed in Section 3.3.

The HC capture probability also changes when the migration graininess is increased in the model. With case 1 and 1000 Plutos, ${P}_{{\rm{HC}}}$ ≃ 4–5 × 10⁻⁴, about a factor of 2 higher than in the smooth migration case. The capture probability increases further, to ${P}_{{\rm{HC}}}\simeq 6\times {10}^{-4}$, when 2000 or 4000 Plutos are included. Together, the increasing ${P}_{{\rm{HC}}}$ and decreasing ${P}_{3:2}$ lead to ${N}_{{\rm{HC}}}/{N}_{3:2}$ values that are more in line with observations. The best fit to observations, with N_HC/N_3:2 ≃ 2–4, occurs for the case-1 simulations with 2000 and 4000 Plutos (Table 1). It is difficult to infer the precise number of Pluto-mass bodies with more confidence. On one hand, additional observations are needed to better constrain the ${N}_{{\rm{HC}}}/{N}_{3:2}$ ratio in the present Kuiper Belt. On the other hand, the massive planetesimals in the original disk must have had a range of masses, while here we represent their size distribution by a delta function. A more realistic modeling with a continuous mass distribution of massive planetesimals is left for future work. See Gladman & Chan (2006) for a modeling work of the effect of very massive planetesimals.

For an outer disk with mass M_disk = 15–20 ${M}_{{\rm{Earth}}}$, ${P}_{{\rm{HC}}}\simeq 6\times {10}^{-4}$ in our reference case with grainy migration implies the HC mass M_HC = 0.008–0.012 ${M}_{{\rm{Earth}}}$, while the mass inferred from observations is ${M}_{{\rm{HC}}}\simeq 0.01$ ${M}_{{\rm{Earth}}}$ (Fraser et al. 2014). This is an excellent agreement. Since the model population of Plutinos in the 3:2 resonance has the right proportion relative to the HCs, as discussed above, this implies that the Plutino mass obtained in the simulation is also approximately correct. (The inner classicals below the 3:2 resonance and Neptune Trojans provide additional constraints that are not considered here.)

The case 2 with ${\tau }_{1}=10\;{\rm{Myr}}$ and ${\tau }_{2}=30\;{\rm{Myr}}$ shows trends in many ways similar to those discussed for case 1 above. The principal effect of faster migration rates is to produce the ${P}_{{\rm{HC}}}$ values that are a factor of ∼2 larger than in case 1 (for the same level of graininess), and ${P}_{3:2}$ values that are somewhat smaller. Together, these trends make it easier to obtain the observed N_HC/N_3:2 ≃ 2–4 with a lower level of graininess. This is illustrated in Figure 8, where we show the orbital distributions obtained with 2000 Plutos. We find that ${N}_{{\rm{HC}}}/{N}_{3:2}\simeq 1.5$ in case 1 and ${N}_{{\rm{HC}}}/{N}_{3:2}\simeq 6$ in case 2. Thus, while the case 1 with 2000 Plutos produces a ratio that is slightly lower than the one indicated by observations, the case 2 with 2000 Plutos overshoots it. The best results in case 2 were obtained with 1000 Plutos in the original disk. In this case, ${N}_{{\rm{HC}}}/{N}_{3:2}=1.8$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0$ and ${N}_{{\rm{HC}}}/{N}_{3:2}=2.5$ for ${\rm{\Delta }}{a}_{{\rm{N}}}=0.5\;{\rm{au}}$ (Table 1).

Zoom In Zoom Out Reset image size

We conclude that the statistics inferred from observations of the resonant and non-resonant populations in the Kuiper Belt implies that the massive planetesimal disk below 30 au contained 1000–4000 Plutos. The combined probability that a planetesimal from the original disk below 30 au evolves on a Kuiper Belt orbit is ∼10⁻³. With 1000–4000 Plutos in the original disk, we would therefore expect that ∼1–4 Pluto-class objects should exist in the Kuiper Belt today, while two such objects are known (Pluto and Eris). This is a reasonable agreement, but note that neither Pluto or Eris is a member of the HC population, while we would expect from our model that the Pluto-size objects are preferentially deposited in the HC population. The expectations would be slightly different if a continuous mass distribution of massive planetesimals were considered. For example, it is plausible that the needed migration graininess was produced by the combined effect of 1000 Plutos and 500 Twoplutos. This would yield ∼1 Pluto and ∼0.5 Twoplutos in the Kuiper Belt today. Obviously, these considerations are subject to the small number statistics. Their main point is to show that the needed graininess from the ${N}_{{\rm{HC}}}/{N}_{3:2}$ constraint is not contradictory to having two Pluto-class objects in the Kuiper Belt today.

The implantation of the disk planetesimals into the Kuiper Belt is a multi-step dynamical process that was first pointed out in Gomes (2003), and is hereafter called the Gomes mechanism. In the Gomes mechanism, disk planetesimals are first scattered by Neptune to >30 au, where they can evolve onto orbits with large libration amplitudes in mean motion resonances. The secular dynamics inside the mean motion resonances, including the Kozai cycles (Kozai 1962), can subsequently act to raise the perihelion distance and decouple the orbits from Neptune. Finally, if Neptune is still migrating, bodies can be released from the resonances onto stable non-resonant orbits. While the Gomes mechanism can operate for a wide range of migration parameters, Nesvorný (2015a) found that the inclination constrain requires that Neptune is given sufficient time to act on the scattered bodies and increase their orbital inclinations before bodies decouple from Neptune (the Kozai cycles inside mean motion resonances also contribute to increasing the orbital inclinations, but they are not the principal factor). Hence it is required that Neptune's migration was slow.

Nesvorný (2015a) used a slightly different migration setup from the one utilized here. The migration in their case was smooth (i.e., no massive bodies in the disk) and characterized by a single migration phase with the starting position of Neptune ${a}_{{\rm{N,0}}}$ and e-folding migration timescale τ. They found that the inclination constraint implies that ${a}_{{\rm{N,0}}}\lesssim 25\;{\rm{au}}$ and $\tau \gtrsim 10\;{\rm{Myr}}$. The migration recipe used in this work was described in Section 2. Here we have two migration stages with a slow migration during the first stage, and even slower migration during the second stage. The migration timescales used for the two phases satisfy the inclination constraint because ${\tau }_{\mathrm{1,2}}\geqslant 10\;{\rm{Myr}}$. At the end of the first phase, we assumed that the orbit of Neptune may have changed discontinuously. And, in addition, our preferred migration mode is grainy. These differences could affect the inclination distribution. Here we therefore test whether the orbital distribution obtained with our favored migration parameters satisfies the inclination constraint.

Figure 9 illustrates how the orbital distribution of bodies obtained in our case-1 simulation (${\tau }_{1}=30\;{\rm{Myr}}$, ${\tau }_{2}=100\;{\rm{Myr}}$, Δa_N = 0.5 au) compares with observations. Two cases are shown: (1) a grainy migration case with 1000 Twoplutos, and (2) a smooth migration case for a reference. We used the CFEPS detection simulator, as described in Section 2.4, and compared the simulated orbits with the actual CFEPS detections. The agreement is satisfactory in the grainy case, where the distribution of the non-resonant orbits roughly follows the lines of constant perihelion distance such that a larger value of the semimajor axis implies larger eccentricity. This trend is a characteristic property of the HC population (see Nesvorný 2015a for a more detailed comparison of our model with the CFEPS observations). The 3:2 resonance population obtained in the grainy migration case has a correct distribution of orbital eccentricities. Moreover, as we discussed in Section 3.1, the HCs and resonant populations appear in the right proportion. For comparison, the smooth migration case shown in Figures 9(a) and 9(b) leads to an excessive number of detections in resonances, which is clearly incorrect.

Zoom In Zoom Out Reset image size

The inclination distribution obtained in the model is wide and roughly comparable to the one inferred from observations. A more careful comparison of the inclination distributions is presented in Figure 10, where the model distributions are shown to follow very closely the CFEPS distributions. This is especially true for the inclination distribution of Plutinos (Figure 10(a)), for which the K–S test gives a 84% probability that the simulated and observed distributions are being derived from the same underlying distribution. The agreement is somewhat less satisfactory for HCs, where the model distribution is slightly wider than the observed one and the K–S test gives a 23% probability. Still, this is a satisfactory match. Also, note that the inclination distribution is sensitive to the migration timescale, and slightly shorter migration timescales should lead to a better agreement (Nesvorný 2015a).

Zoom In Zoom Out Reset image size

Indeed, our case 2 with ${\tau }_{1}=10\;{\rm{Myr}}$ and ${\tau }_{2}=30\;{\rm{Myr}}$ yields a narrower inclination distributions of the HCs (Figure 11(b)). In this case, the K–S test gives the 57% probability for Plutinos and 80% probability for the HCs. The main difference with respect to case 1 is that the model distribution of HCs now represents a reasonable match to observations all the way down to $i\simeq 5^\circ$, while in case 1 we were able to produce a satifactory fit only for $i\gt 10^\circ$ (Nesvorný 2015a). This may indicate that the real migration timescales were closer to case 2 than to case 1.

Zoom In Zoom Out Reset image size

In case 1 with 2000 and 4000 Plutos the inclination distributions are similar to the one shown in Figure 10. All other cases studied here show narrower inclination distributions for Plutinos. Cases 1 and 2 with a smooth migration also show narrower inclination distributions of HCs, which nicely fit the observed distribution for $5^\circ \lt i\lt 10^\circ$, where the CFEPS inclination distribution steeply raises, but fail to match the wide distribution for $i\gt 10^\circ$. This shows that there is some trade-off between the level of graininess and the migration timescale. Any future attempt to closely match the inclination distribution will thus need to explore both these parameters with more resolution. Here we content ourselves with showing that the general results published in Nesvorný (2015a) are valid even if Neptune's migration was grainy.

The distribution of libration amplitudes in the 3:2 resonance was characterized by CFEPS (Gladman et al. 2012). According to Figure 3 in Gladman et al. (2012), the cumulative distribution of libration amplitudes, A_σ, appears to be steadily raising from ${A}_{\sigma }\simeq 20^\circ$ to ${A}_{\sigma }\simeq 100^\circ$, and tails off for ${A}_{\sigma }\gt 100^\circ$. In terms of a differential distribution, Gladman et al. (2012) suggested an asymmetric triangle model, where the number of orbits in a ${\rm{\Delta }}{A}_{\sigma }$ interval linearly increases from zero at ${A}_{\sigma }=20^\circ$ to a maximum at ${A}_{\sigma }\simeq 90^\circ$, and then linearly drops to zero at ${A}_{\sigma }\simeq 130^\circ$, because the orbits with ${A}_{\sigma }\gt 130^\circ$ are unstable.

Gladman et al. (2012) also pointed out that the distribution of A_σ inferred from CFEPS observations is very similar to the theoretical distribution reported in Nesvorný & Roig (2001), where the 3:2 resonance was randomly populated, and the orbital distribution was dynamically evolved over 4 Gyr (with Neptune on its current orbit). The number of resonant orbits increases with A_σ, because the resonant orbits with larger libration amplitudes represent a larger phase-space volume than the orbits with smaller libration amplitudes, and are therefore more populated to start with. There are fewer orbits with ${A}_{\sigma }\gt 100^\circ$, because this is already close to the stability limit, and the original population was depleted when the orbits with ${A}_{\sigma }\gt 100^\circ$ evolved out of the resonance.

Figure 12 shows the distributions of libration amplitudes in the 3:2 resonance obtained in the models with smooth and grainy migrations. The grainy migration case matches the observed distribution much better that the smooth case. The K–S probabilities obtained for the grainy case are 32% for ${\rm{\Delta }}a=0 \%$ and 17% for ${\rm{\Delta }}a=0.5\;{\rm{au}}$, while the probability in the smooth case is ∼10⁻³. This result, in itself, could be used rule out the smooth migration case, where the libration amplitudes are significantly larger than the ones found by CFEPS (Gladman et al. 2012). This is a consequence of the capture mechanism in the 3:2 resonance, which tends to produce large amplitudes if the migration is smooth. The amplitudes in the grainy case are a bit larger than what would be ideal for ${\rm{\Delta }}a=0.5\;{\rm{au}}$, and the distribution is somewhat shallower for ${\rm{\Delta }}a=0$, but we do not consider these slight differences being significant.⁶

Zoom In Zoom Out Reset image size

The distributions of A_σ obtained for the smooth and grainy migration cases are significantly different in all cases studied here. The amplitude distribution obtained for the smooth-migration case 2 is similar to the one shown in Figure 12 for the smooth-migration case 1 (blue line). This is independent of whether or not Neptune jumped during the instability. The smooth migration cases therefore produce, in general, the amplitude distributions that do not agree with observations. With the grainy migration corresponding to 4000 Plutos, on the other hand, the model distributions of A_σ become slightly narrower than what is inferred from observations. From this we conclude that 1000–4000 Plutos give the best fit to observations. It is encouraging to see that while this argument is independent of the one based on the population statistics (Section 3.1), it leads to a similar inference about the number of Pluto-sized objects. In any case, these results represent a significant improvement over those shown in Figure 7 of Gladman et al. (2008), where the amplitude distribution obtained in the model of Levison et al. (2008; ${a}_{{\rm{N,0}}}=28\;{\rm{au}}$, smooth migration with τ = 1 Myr) was shown to be strongly discordant with the CFEPS observations.

The CCs have low orbital inclinations ($i\lt 5^\circ$) and several physical properties (red colors, large binary fraction, steep size distribution of large objects, relatively high albedos) that distinguish them from all other KBO populations.⁷ The most straightforward interpretation of the unique physical and orbital properties is that CCs formed and/or dynamically evolved by different processes than other trans-Neptunian populations. Here we consider a possibility that the CCs formed at >40 au and survived Neptune's early "wild days" relatively unharmed (e.g., Kavelaars et al. 2009; Batygin et al. 2011; Wolff et al. 2012). This requires that the massive planetesimal disk at <30 au had a low-mass extension beyond 30 au, as already discussed in Section 2.2. Nesvorný (2015b) studied this model and found that the original disk at 42–47 au only contained the mass ∼6 × 10⁻³ M_Earth. The surface density of solids in this region, Σ_s ∼ 2 × 10⁻⁵ g cm⁻², was probably therefore some ∼3000 times lower than in the massive part of the disk below 30 au.⁸ This implies that the CCs must have formed by an efficient accretion mechanism that was capable of building ∼100 km planetesimals in a low-mass environment (e.g., Johansen et al. 2009).⁹

According to Petit et al. (2011), the CC population can be divided into the "stirred" and "kernel" components. The stirred orbits have the semimajor axes $42.4\lt a\lt 47\;{\rm{au}}$, inclinations $i\lt 5^\circ$, and small eccentricities with an upper limit that raises from $e\simeq 0.05$ for a = 42 au to $e\simeq 0.2$ for a = 47 au. The kernel is a narrow concentration of low-inclination orbits with $a\simeq 44\;{\rm{au}}$, $e\simeq 0.05$, and a ≃0.5–1 au width in the semimajor axis. Figure 13 illustrates a model of the orbital distribution inferred from the CFEPS observations.

Zoom In Zoom Out Reset image size

Nesvorný (2015b) suggested that the Kuiper Belt kernel can be explained if Neptune's otherwise smooth migration was interrupted by a discontinuous change of Neptune's semimajor axis when Neptune reached ≃28 au. Before the discontinuity happened, planetesimals located at ∼40 au were swept into Neptune's 2:1 resonance, and were carried with the migrating resonance outwards. The 2:1 resonance was at $\simeq 44\;{\rm{au}}$ when Neptune reached ≃28 au. If Neptune's semimajor axis changed by fraction of au at this point, perhaps because Neptune was scattered off of another planet (see Figure 3), the 2:1 population would have been released at $\simeq 44\;{\rm{au}}$, and would remain there to this day. The orbital distribution produced in this model provides a good match to the orbital properties of the kernel.

Nesvorný (2015b) model assumptions and migration parameters were the same as in this work, except that (1) they considered a low mass extension of the planetesimal disk at 30–50 au, and (2) their migration was ideally smooth, while here we showed that the population statistics inferred from observations requires that the migration was grainy. We therefore repeat the simulations of Nesvorný (2015b) to test whether the kernel can form even if the migration was grainy.

Each simulation included 5000 test particles distributed from 30 to 50 au. Their radial profile was set such that the disk surface density ${\rm{\Sigma }}\propto 1/r$. There is therefore an equal number of particles (250) in each radial au. A larger resolution is not needed, because a significant fraction of particles in the CC region survive. The disk extension was assumed to be dynamically cold with low orbital eccentricities and low orbital inclinations. The initial inclinations were set to be similar to those inferred for the present population of CCs. Specifically, we used $N(i){di}=\mathrm{sin}i\mathrm{exp}(-{i}^{2}/2{\sigma }_{i}^{2}){di}$, with ${\sigma }_{i}=2^\circ$ (Brown 2001; Gulbis et al. 2010). The initial eccentricities were set according to the Rayleigh distribution with ${\sigma }_{e}=0.01$, 0.05, or 0.1.

Figure 14 shows the orbital distribution of particles obtained in a model with ${\sigma }_{e}=0.01$, and the case-1 migration parameters with the graininess corresponding to 1000 Twoplutos. The result is similar to those published in Nesvorný (2015b) for the smooth migration. The concentration of orbits near 44 au has orbital properties comparable to those of the CFEPS kernel. This shows that the grainy migration required to explain the population statistics of resonant and non-resonant populations also allows for the formation of the Kuiper Belt kernel. These results are therefore consistent with each other. The concentration of orbits obtained in the model near 44 au becomes slightly more fuzzy for ${\sigma }_{e}=0.05$ or 0.1, following the trends described in Nesvorný (2015b). The case-2 parameters with grainy migration also lead to the formation of the kernel. We therefore conclude that the model of kernel formation described in Nesvorný (2015b) does not require that the migration was smooth. Instead, it works even if the migration was grainy.

Zoom In Zoom Out Reset image size

Adams et al. (2014) found from the Deep Ecliptic Survey (DES) that ${N}_{3:2}/{N}_{2:1}\simeq {N}_{3:2}/{N}_{5:2}\simeq 2$, while Gladman et al. (2012) suggested from the CFEPS that N_3:2/N_2:1 ≃ 3–4 and ${N}_{3:2}/{N}_{5:2}\simeq 1$. Part of these differences between the DES and CFEPS may stem from differences in observational strategies and/or debiasing approach. While there is obviously a significant uncertainty in these estimates, it is probably fair to say that observations suggest roughly comparable populations in the 2:1 and 5:2 resonances (to within a factor of ∼2 or so), both of which are ≃1–4 times smaller than Plutinos in the 3:2 resonance.

Here we find that the smooth migration cases produce ${N}_{3:2}/{N}_{2:1}\sim 15$ and ${N}_{3:2}/{N}_{5:2}\sim 10$. The 2:1 and 5:2 resonances are therefore clearly underpopulated, relative to the 3:2 resonance, if the migration is assumed to be smooth. Much better results were obtained for the grainy migration. In case 2 and 2000 Plutos, where the lowest resonant ratios were found, N_3:2/N_2:1 ≃ 2–2.5 and ${N}_{3:2}/{N}_{5:2}\simeq 1$. It is therefore plausible that the population of bodies in the 5:2 resonance can be as large as Plutinos (Volk et al. 2016). Other results obtained here are intermediate. For example, case 1 gives ${N}_{3:2}/{N}_{2:1}\sim 10$ and ${N}_{3:2}/{N}_{5:2}\simeq 6$ for 2000 Plutos, and N_3:2/N_2:1 ≃ 3–5 and ${N}_{3:2}/{N}_{5:2}\simeq 2$ for 4000 Plutos.

The increased level of migration graininess therefore leads to lower values of ${N}_{3:2}/{N}_{2:1}$ and ${N}_{3:2}/{N}_{5:2}$, which are in better agreement with observations (Gladman et al. 2012; Adams et al. 2014; Volk et al. 2016). This trend is mainly contributed by the lower value of ${P}_{3:2}$ when the migration is assumed to be grainy. On the other hand, if different migration timescales are considered for the same level of graininess, then the cases with faster migration rates tend to produce larger populations in the 2:1 (a ≃ 47.8 au) and 5:2 (a ≃ 55.5 au) resonances than the cases with slower migration rates. This trend is clearly visible in Figure 8, where the 2:1 and 5:2 resonances are more populated in case 2 (${\tau }_{1}=10\;{\rm{Myr}}$ and τ₂ = 30 Myr) than in case 1 (${\tau }_{1}=30\;{\rm{Myr}}$ and τ₂ = 100 Myr).