CHAOTIC CAPTURE OF NEPTUNE TROJANS

Numerical studies showed that Neptune Trojans (NTs) can survive the standard planetary migration (Gomes 1998; Kortenkamp et al. 2004) and be orbitally stable over 4 × 10⁹ yr (Weissman & Levison 1997; Nesvorný & Dones 2002). This led Nesvorný & Dones (2002) to make detailed predictions about the sky density of NTs near Neptune's L₄ and L₅ Lagrangian points. Various orbit distributions were used to generate these maps, each corresponding to a different assumption on the NT orbit distribution. This study showed that the cloud of NTs, as viewed from the Earth, can be widely dispersed on the sky about the positions of Neptune's Lagrangian points, mainly because the NT orbits can be stable up to inclinations i ≈ 35°.

The first NT, 2001 QR₃₂₂, was discovered by the Deep Ecliptic Survey (Chiang et al. 2003). The orbit of this object has a small libration amplitude around the L₄ point (A ≈ 70°, where A is the full longitude range experienced by the particle as it orbits about a Lagrange point; see below) and very small inclination (i = 18). Based on these orbital characteristics, Chiang & Lithwick (2005) proposed that "the largest NTs represent the unadulterated outcome of dispersion-dominated oligarchic accretion" (a process that presumably took place in situ) and that "the large NTs may today have a dispersion in orbital inclination of less then ∼10°, despite the niches of stability at higher inclinations." These predictions have not been confirmed.

As of 2008 November we know of two NTs (out of six NTs in total; Table 1), 2005 TN₅₃ and 2007 VL₃₀₅, that have very high inclinations (i = 25° and i = 28°, respectively). Both of these objects have been observed during several oppositions (three and four, respectively); their orbit determination is therefore reasonably precise. They must represent the "tip of an iceberg" because the detection of 2005 TN₅₃ by the Magellan–Baade low-ecliptic-latitude survey implies that the population of NTs with highly inclined orbits is numerous (Sheppard & Trujillo 2006). This shows that the NTs occupy a thick disk, which is indicative of capture instead of the in situ formation proposed by Chiang & Lithwick (2005).

Table 1. Known Neptune Trojans as of 2009 January

ID	H	Opps.	Ln	A (deg)	$\bar{i}$ (deg)	$\bar{e}$
2001 QR322	7.8	5	L4	71	1.8	0.030
2004 UP10	8.9	3	L4	37	2.7	0.027
2005 TN53	9.1	3	L4	22	25.3	0.058
2005 TO74	8.5	3	L4	24	4.8	0.053
2006 RJ103	7.4	3	L4	22	6.7	0.023
2007 VL305	8.0	4	L4	39	27.9	0.064

Notes. The columns are: (1) object's ID (ID); (2) absolute magnitude (H); (3) number of oppositions (Opps.); (4) Lagrangian point about which the orbit librates (L_n); (5) full libration amplitude (A); (6) mean inclination $(\bar{i})$; and (7) mean eccentricity $(\bar{e})$. All known NTs librate about L₄ because L₅ was not looked at as it is currently in the Galactic plane. The libration amplitude, mean inclination and mean eccentricity were determined from a 17-Myr-long numerical integration of orbits. The uncertainties of these values are ≈5°, 005, and 0.001, respectively. The uncertainties were estimated by tracking the evolution of clones of each object with the initial orbit elements distributed within the 1σ orbit-uncertainty region around the nominal orbit. The source of initial orbit elements was AstDyS (http://hamilton.dm.unipi.it/astdys/).

Download table as:  ASCIITypeset image

Here we investigate the capture of Neptune's Trojans in the framework of the Nice model (Tsiganis et al. 2005) because this is the most complete model currently available for the early evolution of the outer solar system. In the Nice model, the giant planets are assumed to have formed in a compact configuration (all were located at 5–18 AU). Slow migration was induced in these planets by gravitational interaction with planetesimals leaking out of a massive primordial transplanetary disk. After a long period of time, most likely some 700 Myr after formation of the giant planets (Gomes et al. 2005), planets crossed a major mean-motion resonance (MMR; e.g., the 2:1 resonance between Jupiter and Saturn; Tsiganis et al. 2005). This event triggered a global instability that led to a violent reorganization of the outer solar system. Uranus and Neptune penetrated the transplanetary disk, scattering its inhabitants throughout the solar system. Finally, the interaction between the ice giants and the planetesimals damped the orbits of these planets, leading them to evolve onto nearly circular orbits at their current locations.

The Nice model is compelling because it can explain many of the characteristics of the outer solar system, such as the current orbits of four giant planets (Tsiganis et al. 2005), orbital distribution of the irregular satellites (Nesvorný et al. 2007), and the main properties of the Kuiper Belt (Levison et al. 2008). The late instability between planets in the Nice model may have produced the Late Heavy Bombardment (LHB) of the Moon (Gomes et al. 2005). Moreover, the model can explain the large inclinations of the Jupiter Trojan (JT) population as these objects became trapped from a planetesimal disk that was stirred by planets (Morbidelli et al. 2005).

The capture of JTs in the Nice model occurs during discrete episodes when the co-orbital region is swept by secondary resonances associated with the mean-motion commensurabilities between Jupiter and Saturn (Morbidelli et al. 2005; Marzari & Scholl 2007; Robutel & Bodossian 2009). Each time one of these resonances sweeps over Jupiter's Lagrangian points, the local dynamics become chaotic and objects can freely flow in (and out). They can become permanent JTs when the resonances move away and the planetary migration ceases. This mechanism, very different from the one that produced populations of KBOs in the exterior MMRs with Neptune (such as Plutinos; see Malhotra 1995), is known as the chaotic capture (Morbidelli et al. 2005).

Tsiganis et al. (2005) noticed in their simulations that several particles were trapped on long-lived orbits characteristic of Neptune's Trojans (with 0.2% trapping efficiency, on average, and greater than 80 Myr dynamical lifetimes). These particles were eventually removed from the NT region, but this could have been an artifact of the graininess of Neptune's migration as only 1000 to 5000 particles were used to model the planetesimal disk. The capture was thought to have been produced as the eccentric orbit of Neptune shrank due to dynamical friction from the planetesimal disk (H. F. Levison 2008, private communication; see also Li et al. 2007). We investigate these issues in more detail here.

The structure of this paper is as follows. In Section 2, we describe our model. We used ≈27,000 particles to represent the planetesimal disk and performed 90 simulations of the Nice model (with selected runs extending to time t = t₀ + 10⁹ yr; see below for the definition of t₀) to test the variability of the results on different assumptions. As it turned out that the efficiency and orbit distribution of captured NTs are both sensitive to these assumptions, we selected only a few Nice-model migration runs where the architecture of the final planetary systems was very similar to that of our own. The NT populations produced in these runs were analyzed in detail.

The results are described in Section 3. We found that NTs are produced by chaotic capture with changes of Neptune's eccentricity having only a minor effect. The capture efficiency that we estimate from our model indicates that large NTs may be (at least) ∼10 times more populous than large JTs, in agreement with the results of Chiang et al. (2003) and Sheppard & Trujillo (2006, 2008). Broader implications of our work, such as constraints on the migration histories of the outer planets, are discussed in Section 4.

We used the results of Gomes et al. (2005) to record the orbits of the outer planets and the state of the planetesimal disk at time t₀. We chose t₀ so that the state of the system represents the situation shortly before the orbits of Jupiter and Saturn migrate through the 2:1 MMR, become dynamically excited, and initiate the epoch of planetary encounters. We used one of the published simulations of Gomes et al. (2005) and t₀ = 868 Myr. We do not need to study the epoch prior to t₀ because a population of primordial NTs captured at t < t₀ would not survive the epoch of planetary encounters in the Nice model when Neptune's orbit was violently changing.

In the original simulation of Gomes et al. (2005), the pre-LHB planetesimal disk was represented by several hundred objects. We cloned each original particle so that the resulting disk was represented by 27,028 equal-mass objects. This increased resolution was necessary because the capture of Trojans is a low-probability process. By increasing the resolution we also wanted to avoid artificial graininess of the migration from an underresolved disk.

The cloning was done as a small perturbation to the positions and velocities of the original particles. The masses of the clones were adjusted so that the total mass of the disk (≈35 M_E, where M_E = 6 × 10²⁷ g is the Earth mass) and its distribution were preserved (the original planetesimal disk was located between 21 and 35 AU). In total, we created 60 distinct initial states of the original disk. Below we identify these states by the seed value that was used to initialize the random number generator.

The positions and velocities of four outer planets (Jupiter to Neptune) at t₀ = 868 Myr were preserved. The starting semimajor axis of Jupiter, Saturn, and the two ice giants at t₀ were 5.4, 8.4, 12.3, and 18.0 AU, respectively. The orbits of disk planetesimals at t₀ had low eccentricities and inclinations (e ≲ 0.05 and i ≲ 1°) except for the inner part of the disk where bodies were scattered by Neptune and passing resonances.

In addition to the nominal migration runs described above, we also performed experiments with the planetesimal disk truncated at ≈30 AU. This adjustment was motivated by a frequent outcome of our nominal migration runs in which Neptune migrated beyond 30 AU and Uranus/Neptune passed through their mutual 2:1 MMR. This produced problems because the NT population was removed in these simulations as a result of late crossing of this resonance. In reality, Uranus and Neptune stopped just before reaching the 2:1 MMR. By setting the disk edge at ≈30 AU, we thus hoped to stabilize Neptune near 30 AU in our simulations and obtain realistic final configurations of the planets.

These new migration simulations were divided into three sets, each having a different mass of the planetesimal disk. We used disks with 25 M_E, 35 M_E, and 45 M_E. In each case, we performed 10 simulations with slightly distinct initial states of the truncated disk. As in our 60 nominal simulations, we created these initial states by a small perturbation to the positions and velocities of the original particles and identified them by the seed value that was used to initialize the random number generator. To distinguish the nominal and truncated disk simulations, we use the following notation: "n" (nominal), "tl" (truncated, low mass), "tm" (truncated, medium mass), and "th" (truncated, high mass). Thus, for example, the nominal simulation with seed value 22 is denoted by "n22," while the truncated disk simulation with a 35-M_E-mass disk and seed value 6 is denoted by "tm6."

In each simulation, we numerically tracked the orbits of planets and planetesimals for at least 100 Myr. We used the symplectic integrator known as SyMBA (Duncan et al. 1998). SyMBA is a highly efficient symplectic N-body integrator similar to the Wisdom–Holman map (Wisdom & Holman 1991) which is able to handle close encounters between massive bodies. In the code, planets gravitationally interact with each other and also act on other bodies in the simulation. The planetesimals do not interact with each other but affect the orbits of planets. (This approximation may lead to a slightly faster migration rate than that expected for planet interactions with the real self-gravitating planetesimal disk (Moore et al. 2009) but it should not matter otherwise.)

In total, we performed 90 numerical integrations with individual runs that started with slightly different disk states (see above) and produced different migration outcomes. Such a large number of experiments is needed because the planet evolution in the Nice model is chaotic and difficult to predict. We must therefore perform a large number of runs to produce at least a few in which all important details of the current outer solar system architecture are just right. Since each simulation required two to four months on a single CPU, the total computer use in this work was 270 CPU-months. The results of these runs are summarized in Section 3.

We selected cases out of our 90 migration simulations in which the migrating planets at t₀ +  100 Myr resembled the current architecture of the outer solar system. Specifically, we required that none of the planets be ejected from the solar system and that all four planets were near their current semimajor axis locations with low eccentricities and low inclinations. We paid particular attention to the ice giants because, for the reasons described above, their orbits needed to be tuned with care. The selected migration runs were then continued to t₀ + 1 Gy to check on the migrating planets's behavior in a progressively more depleted disk and to filter out the captured NTs that were not stable in the long run. We found that the migration rate of planetary orbits and escape/capture rate of NTs at t₀ + 1 Gy was minimal and thus represented a good proxy for the final state of the system.

For each of these endstates we determined the mean values of eccentricity, $\bar{e}$, and inclination, $\bar{i}$, of objects captured near Neptune's Lagrangian points, and their libration amplitudes, A. Specifically, we calculated Δλ = λ − λ_N, where λ and λ_N were the object's and Neptune's mean longitudes, respectively, and defined the full amplitude A = Δλ_max − Δλ_min, where Δλ_max and Δλ_min are the maximum and minimum values attained by Δλ during the object's libration. Using numerical integrations we also determined $\bar{e}, \bar{i}$, and A for the known NTs and reported them in Table 1.

The populations of NTs produced in our successful simulations were compared with the real population of NTs predicted from observations (Chiang et al. 2003; Sheppard & Trujillo 2006, 2008). Specifically, we compared (1) the distributions of libration amplitudes, eccentricities, and inclinations and (2) the NT population size expected from our simulations and the one inferred from the observations. We also analyzed the orbital histories of captured objects to determine (3) the mechanism of capture, and (4) the orbital evolution of objects before and after the instant of capture.

3.1. Planetary Migration

A vast majority of our nominal migration jobs (35 M_E disk truncated at 35 AU) had to be discarded because they did not produce a good analog for the outer solar system. In eight cases out of 60 (13%), Jupiter and Saturn have not reached the 2:1 resonance within 100 Myr after t₀. In 20 cases (33%), the simulation ended with less then four planets (with at least one having been ejected) or produced planetary orbits with large e and i, much larger than those in the present solar system. While these migration results show the true range of possible evolution outcomes for a planetary system in the Nice model, they are not useful here.

The main problem with most of the remaining nominal jobs was that Uranus and Neptune migrated beyond their mutual 2:1 MMR (26 cases; 43%), which has not happened in the solar system where currently P_N/P_U ≈ 1.96 with P_U and P_N being the orbital periods of Uranus and Neptune, respectively. As we discussed in Section 2, these cases cannot be considered here because the 2:1 MMR crossing by Uranus and Neptune would unrealistically deplete the NT population.

Discarding several additional migration simulations with minor inconsistencies (e.g., Jupiter and Saturn ending too close to the 5:2 resonance or Uranus and Neptune in the 7:4 resonance; no NTs were produced in these simulations because the NT orbits were unstable in the long run), we found that only one case, n22, has all required characteristics (Figure 1). This shows that it is rather difficult to obtain precisely the correct planetary orbits with our nominal Nice-model setup (one good case out of 60; 1.7%).

Zoom In Zoom Out Reset image size

Next, we discuss the results of simulations in which the disk was truncated at 30 AU (instead of the original 35 AU) to limit the total radial migration span.

None of 10 finished tl simulations has been successful: six ended with Uranus and Neptune migrating beyond the 2:1 resonance and four with large planetary eccentricities and inclinations. In the tm simulations, however, the success rate was better with two jobs out of 10 (20%) ending with good planetary orbits and P_N/P_U < 2.0. We discuss these cases, tm6 and tm7, in more detail below. In addition, one simulation with the massive disk, th10, has also been successful. It therefore seems that truncating the disk at 30 AU helps in achieving a better success rate as far as the disk is sufficiently massive (≳35 M_E).

The following analysis is based on simulations n22, tm6, tm7 and th10 (Table 2). We also include simulations n12 and n25 in which Uranus and Neptune barely crossed the 2:1 resonance. To make these two runs useful, we stopped them when P_N/P_U = 1.96. The masses of disk planetesimals were then set to zero and the simulations were continued to t₀ + 1 Gy. We denote these new runs by n12b and n25b. In four successful cases listed at the top of Table 2, n22, tm6, tm7, and th10, massive planetesimals were used for the full length of the simulation.

Table 2. Selected Migration Simulations

ID	Mass (ME)	Edge (AU)	PN/PU	No. of NTs	No. of Stable	Prob. (10−4)
n22	35	35	1.89	6	1	0.4
tm6	35	30	1.87	6	2	0.7
tm7	35	30	1.86	24	9	3.3
th10	45	30	1.92	22	6	2.2
						*1.6
n12b	35	35	1.96	18
n25b	35	35	1.96	21

Notes. The columns are: (1) job's ID (ID); (2) planetesimal disk's mass (Mass); (3) outer disk's edge (Edge); (4) final orbital period ratio of Uranus and Neptune (P_N/P_U); (5) total number of captured NTs; (6) number of stable NTs with A < 70° and e < 0.08; and (7) probability of capture per one disk particle (Prob.; asterisk denotes the average from n22, tm6, tm7, and th10 simulations). The capture probability estimate is based on column 6, which is the number of NTs stable over 1 Gy. Simulations n12 and n25 were stopped at t₀ + 40 Myr and t₀ +  54 Myr, respectively, when P_N/P_U = 1.96. The masses of disk planetesimals were set to zero at this point and the simulations were then continued to t₀ + 1 Gy. We denote these new runs by n12b and n25b. In four successful cases listed at the top of the table, n22, tm6, tm7 and th10, massive planetesimals were used for the full length of the simulation.

Download table as:  ASCIITypeset image

The basic information about selected runs is listed in Table 2. Figures 1, 2, and 3 show the evolution of planets and disk planetesimals in n22. Figure 4 shows P_N/P_U for n22.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In n22, the instability due to the 2:1 resonance crossing between Jupiter and Saturn occurred at ≈t₀ + 8 Myr. The resonant crossing and subsequent encounters of planets produced evolution that is characteristic for the Nice model: planetary eccentricities and inclinations increased, Neptune moved into the outer planetesimal disk and the planets started migrating. The dynamical friction from the planetesimal disk then damped planetary e and i so that the final stages of the planet evolution resemble that of standard planetary migration (Malhotra 1995).

Interestingly, our successful simulations share some similarities. For example, in all these runs Neptune was ejected by planetary encounters to a ≈ 25 AU from where it smoothly migrated to 30 AU by interaction with the planetesimal disk. The migration of Neptune from 25 to 30 AU may be required to explain the observed characteristics of the Kuiper belt (Levison & Morbidelli 2003; see also Malhotra 1995).

Several slight discrepancies were noted when we compared the final results of our calculations with the real planetary orbits. First, in all four successful runs Saturn ended below its real orbit with a ≈ 9 AU (the real mean semimajor axis is 9.55 AU). The same issue has been noted in the simulations of Tsiganis et al. (2005) and Nesvorný et al. (2007). This suggests that the initial orbits of the planets prior to LHB were slightly different from those tested so far and/or that the planetesimal disk was slightly more massive (≳50 M_E), thus producing a larger overall drift of the Saturn's orbit.

Second, we noted that the final orbits of Uranus and Neptune obtained in the simulations were dynamically colder than the present orbits of ice giants. For example, the final mean eccentricities of Uranus and Neptune were typically ∼0.005 while the real values are $\bar{e}_{\rm U}=0.042$ and $\bar{e}_{\rm N}=0.01$. This difference stems from the strong dynamical friction that these planets experience in our simulations from the dynamically cold planetesimal disk (Tsiganis et al. 2005).

Interestingly, in the original simulations of Tsiganis et al. (2005), Uranus and Neptune ended on dynamically more excited orbits than those we obtained here. We believe that this difference stems from the graininess of the planet–disk interaction in the simulations of Tsiganis et al. (2005), in which the disk was resolved by 1000–5000 particles only, each having mass ∼3-14 M_P, where M_P ≈ 1.3 × 10²⁵ g is Pluto's mass. In our simulations, in which the disk was resolved by ≈27,000 particles with ∼0.5 M_P, the graininess was reduced. Also, the results of Nesvorný et al. (2007), where ≈7000 disk particles were used, were intermediate between the results obtained here and those of Tsiganis et al. (2005).

These results could indicate that the real migration was more grainy than the one assumed in this work. For example, Charnoz & Morbidelli (2007) suggested that the planetesimal disk may have contained ∼1000 planetesimals with ∼1 M_P; interactions with these objects could have caused grainy migration. A simple interpolation between the migration results discussed above shows, however, that it is unlikely that ∼1000 Pluto-mass objects would have caused the desired effect. Instead, interactions of Uranus and Neptune with more massive, planetary-sized bodies would be needed.

Alternatively, and perhaps more likely, massive disk objects may have stirred the planetesimal disk prior to the LHB thus reducing the disk's ability to damp planetary eccentricities and inclinations by dynamical friction effects. We did not account for these effects because our integration scheme neglects interaction between planetesimals. The study of planetary migration into a planetesimal disk containing massive objects is left for future work.

3.2. Chaotic Capture

Table 3 lists six disk objects that were captured in our simulation n22 in permanent orbits around Neptune's L₄ and L₅ points. Figures 5 and 6 show the orbit histories of two of these objects. These plots illustrate that captures occur when Uranus and Neptune cross their mutual MMRs, such as the 5:3, 7:4, 9:5, and 11:6 MMRs. This happens as follows. Each MMR is surrounded by a number of secondary resonances where the NT's libration period is commensurable with the Uranus/Neptune MMR period. These secondary resonances, when they occur during the planetary migration, produce chaos in Neptune's co-orbital region and have destabilizing effects on the NT population (Figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. The Orbits of Six NTs Obtained in the n22 Simulation

Particle		A (deg)	${\bar{i}}$ (deg)	${\bar{e}}$	Res.
2428	L5	93.2	18.5	0.015	5:3
3570	L5	26.3	11.5	0.085	5:3–7:4
3699	L5	91.5	22.5	0.090	11:6
3895	L5	15.7	8.5	0.065	5:3
2981	L4	97.9	5.0	0.060	3:2–8:5
4178	L4	81.2	6.0	0.105	7:4–9:5

Notes. The columns are: (1) particle's ID; (2) Lagrangian point about which the orbit librates (L_n); (3) libration amplitude (A); (4) mean inclination $(\bar{i})$; (5) mean eccentricity $(\bar{e})$; and (6) MMRs between Uranus and Neptune that contributed to the capture (Res.). If range is given in column 6 (e.g., 5:3–7:4 for particle 3570) this indicates that weaker resonances in the given range contributed to the capture. The libration amplitude, mean inclination and mean eccentricity were determined at t₀ + 1 Gy.

Download table as:  ASCIITypeset image

The destabilizing effects of chaos produced by secondary resonances were studied by Kortenkamp et al. (2004). They showed that the secondary resonances associated with the above listed MMRs are very powerful and would severely deplete any original NT population during planetary migration. What their study missed, however, is that the effects of chaos can also produce NTs because "if there is a way out, there is a way in." As shown by Morbidelli et al. (2005) for JTs, the chaotically evolving orbits can be stabilized when resonances move away and the planetary migration ceases. Here we find that NTs have been captured by a similar mechanism.

In the n22 simulation, Uranus and Neptune crossed resonances at ≈t₀ +  25 Myr (5:3 MMR), ≈t₀ +  70 Myr (7:4), ≈t₀ +  100 Myr (9:5), and ∼t₀ + 140 Myr (11:6) (see Figure 4). For a comparison, particle no. 3895 (Figure 5) was captured into libration around L₅ at t₀ + 23 Myr. This nicely corresponds to the time of crossing of secondary resonances between the NT's libration frequency and Uranus/Neptune 5:3 resonant frequency. The capture of particle no. 2428 as NT was aided by the same resonances (Table 3).

Figure 6 shows the orbit history of particle no. 3699. Particle no. 3699 was captured on an orbit co-orbital with Neptune at t₀ + 142 Myr. This corresponds to crossings of secondary resonances associated to the 11:6 MMR between Uranus and Neptune (Figure 4). Later on, the orbit history of the particle showed two episodes, at t₀ +  250 Myr and t₀ +  680 Myr, where the libration amplitude increased and the object temporarily resided on a horseshoe trajectory. During the t₀ + 250-Myr episode it even temporarily librated around the L₄ point before moving back to L₅.

Figure 7 shows that NTs can be removed (or captured) not only at the time of crossing of the main MMRs between Uranus and Neptune but also during the interim migration periods. In fact, a large fraction of NT orbits with small libration amplitudes are unstable when Uranus and Neptune migrate between the 3:2 and 7:4 resonances (1.5 < P_N/P_U < 1.75; Figure 7). This instability is probably generated by secondary resonances associated with weaker MMRs between Uranus and Neptune. This shows that NTs can be captured at any time during the early phase of planetary migration. Many orbits reached NT libration in our migration model during this period (Table 3).

3.3. Orbital Distribution of Captured NTs

Several particles captured as NT in our simulations had very large inclinations. For example, particles no. 2428 and 3699 in n22 ended up with $\bar{i}=18\mbox{$.\!\!^\circ $}5$ and $\bar{i}=22\mbox{$.\!\!^\circ $}5$, respectively (Table 3). These values are comparable to those of known high-inclination NTs (Table 1). This result is encouraging. It shows that chaotic capture is capable of producing the observed high-i orbits of NTs, which would not be expected from other NT formation mechanisms proposed so far (e.g., Chiang & Lithwick 2005). Based on this result, and the efficiency of chaotic capture discussed in Section 3.4, we believe that observed NTs are captured objects rather than having formed in situ.

The large inclinations of captured NTs reflect the similarly large inclinations of disk planetesimals that have become scattered by planets. Figure 3 shows that the inclination distribution of planetesimals at t₀ + 100 Myr is nearly uniform between 0 and 35 deg. Therefore, late captures are expected to produce very large inclinations. Indeed, particle no. 3699, which was captured at t₀ + 142 Myr, ends up with i = 225 (Figure 6). On the other hand, particles that were captured relatively early are expected to have lower inclinations because it takes time for the planets to stir up the planetesimal disk. For example, the mean inclination of planetesimals on Neptune-crossing orbits at t₀ + 11 Myr is 102 while it is 152 at t₀ + 100 Myr.

It is difficult to compare the model and observed orbits of NTs in detail because: (1) only six NTs are known and their orbit distribution is strongly biased toward low i; and (2) only a small number of disk objects were captured as NTs in our simulations. One important discrepancy, however, is already apparent from the existing results.

Sheppard & Trujillo (2006) estimated from the discovery rate of NTs with the Magellan–Baade telescope that the high-inclination NTs (i ∼ 25°) should be ∼4 times more numerous than the low-inclination NTs (i ∼ 5°). We are not obtaining such an extreme inclination distribution in our simulations. In the n22 simulation, the numbers of low-i and high-i NT orbits are roughly similar. In other successful simulations listed in Table 2, the fraction high-i orbits is even smaller. For example, 15 and seven orbits ended with i < 10° and i > 10°, respectively, in the th10 simulation (Figure 8). We discuss possible explanations for this problem in Section 4.

Zoom In Zoom Out Reset image size

Figure 8 compares the observed orbits of known NTs with those obtained in the th10 simulation. The model and observed distributions of eccentricities and libration amplitudes are broadly similar. The captured objects have, however, typically slightly higher eccentricities and libration amplitudes than the real NTs. Specifically, the eccentricities of captured objects range up to ≈0.1–0.13, while the real NTs have e < 0.07. Also, their libration amplitudes range up to ≈100°, while the real NTs have A ≲ 70°.

These differences probably stem from the fact that the final orbit of Neptune obtained in our simulations has lower e and i values than those of the present Neptune's orbit ($\bar{e}_N=0.004$ and $\bar{i}_N=0\mbox{$.\!\!^\circ $}12$ as opposed to present values $\bar{e}_N=0.01$ and $\bar{i}_N=0\mbox{$.\!\!^\circ $}66$; inclinations given with respect to the invariant plane). We discussed this issue in Section 3.1. With the present larger values, the NT orbits with e > 0.08 and/or A > 70° are not stable (Nesvorný & Dones 2002), while these orbits survive in our simulations. To resolve this problem, the model needs to be improved so that the final orbits of the outer planets produced by simulations become more similar to the current ones (see Section 4).

3.4. Capture Probability and Mass of NTs

We performed 90 simulations of the Nice model to study the origin of the recently discovered population of NTs. The initial orbits of planets and state of the planetesimal disk were taken from Gomes et al. (2005), but we cloned the original particles to obtain a better resolution of the disk. The disk was set to be dynamically cold with inclinations less than 1°. We found that, if the disk has the outer edge at ≈30 AU (rather than at 35 AU), there is a better chance that Uranus and Neptune stop migrating before reaching the 2:1 resonance, as it happened in reality. The disk's edge at ≈30 AU is therefore favored by our results.

NTs were captured in our migration model by a process similar to chaotic capture of JTs proposed by Morbidelli et al. (2005). We showed that chaotic capture can produce orbital inclinations of NTs in excess of 20° and thus explains the high-inclination orbits of 2005 TN₅₃ and 2007 VL₃₀₅. The original orbital radii of captured bodies range nearly over the full radial extent of the planetesimal disk. This can explain why the observed colors of known NTs are similar to those of JTs, irregular satellites and gray objects in the classical Kuiper belt (Sheppard & Trujillo 2006), because all these populations may have formed in the same source region (Morbidelli et al. 2005; Nesvorný et al. 2007; Levison et al. 2008).

We found evidence that the planetesimal disk may have been dynamically excited before the Nice-model planetary instability occurred because: (1) capture from a dynamically cold disk cannot explain the 4:1 ratio between the high- and low-inclination NTs (Sheppard & Trujillo 2006); (2) capture probability from a dynamically cold disk (∼10⁻⁴ per disk particle) is too large and inconsistent with observations; and (3) final orbits of Uranus and Neptune have incorrectly low e and i values in our simulations, because these planets suffer strong dynamical friction from the dynamically cold planetesimal disk that we used in this work.

We believe that these issues can be resolved if the planetesimal disk became dynamically excited prior to the planetary instability in the Nice-model framework. The required stirring of the planetesimal disk could have been produced by Pluto-sized and larger objects that may have formed in the disk. Since the escape velocity from Pluto (≈1.2 km s⁻¹) is about one fifth of the orbital speed at 25 AU (≈6 km s⁻¹), numerous Pluto-sized objects in the disk (Charnoz & Morbidelli 2007) could have easily produced orbits of planetesimals with e ≳ 0.2 and i ≳ 12°. Chaotic capture from this population should have lower efficiency than the one discussed for the cold disk in Section 3.4 because e < 0.08 is required for the NT stability. Large inclinations would also help to resolve (1). We plan to perform these new simulations in the near future.

This article is based on work supported by NSF. The work of D.V. was supported by research grant 205/08/0064 of the Czech Grant Agency. We thank Alessandro Morbidelli, Hal Levison, Luke Dones, and Bill Bottke for their insightful comments on this work. A. Morbidelli has provided initial conditions for our migration simulations.