Astronomical The Journal, 123:28622883, 2002 May
# 2002. The American Astronomical Society. All rights reserved. Printed in U.S.A.
THE FORMATION OF URANUS AND NEPTUNE AMONG JUPITER AND SATURN E. W. Thommes
Department of Astronomy, University of California at Berkeley, Berkeley, CA 94720; [e-mail]
M. J. Duncan
Department of Physics, Queen's University, Kingston, ON K7L 3N6, Canada; [e-mail]
and H. F. Levison
Southwest Research Institute, 1050 Walnut Street, Boulder, CO 80302; [e-mail] Received 2001 November 14; accepted 2002 January 30
ABSTRACT The outer giant planets, Uranus and Neptune, pose a challenge to theories of planet formation. They exist in a region of the solar system where long dynamical timescales and a low primordial density of material would have conspired to make the formation of such large bodies ($15 and 17 times as massive as Earth, respectively) very difficult. Previously, we proposed a model that addressed this problem: Instead of forming in the trans-Saturnian region, Uranus and Neptune underwent most of their growth among proto-Jupiter and proto-Saturn, were scattered outward when Jupiter acquired its massive gas envelope, and subsequently evolved toward their present orbits. We present the results of additional numerical simulations, which further demonstrate that the model readily produces analogs to our solar system for a wide range of initial conditions. We also find that this mechanism may partly account for the high orbital inclinations observed in the Kuiper belt. Key words: celestial mechanics -- Kuiper belt -- planets and satellites: formation -- solar system: formation
1. INTRODUCTION
The growth of Uranus and Neptune in the outer solar system is not readily accounted for by conventional models of planet formation. A low primordial density of planetesimals and weak solar gravity would have made the process of accretion slow and inefficient. In direct N-body simulations of accretion among (approximately) Earth-mass bodies beyond 10 AU, performed with three different computer codes, little accretion is found to take place over timescales of 108 yr and, by extrapolation, over the age of the solar system (Levison & Stewart 2001). Earlier simulations by Brunini & Fernandez (1999) showed accretion of the ice giants in several times 107 yr with the same initial conditions, but later simulations, performed with an improved integrator, required that bodies be enhanced in radius by a factor of at least 10 relative to bodies having the density of Uranus and Neptune in order to recover the previous result (Brunini 2000). Therefore, Uranus and Neptune are unlikely to have formed from a late stage of mergers among large protoplanets, analogous to the putative final phase of planet formation in the terrestrial zone (e.g., Wetherill 1996; Chambers & Wetherill 1998). The oligarchic growth model, in which the principal growth mode is accretion of small planetesimals by a protoplanet, also produces timescales that are too long (Kokubo & Ida 2000; Thommes 2000), unless the feeding zones of Uranus and Neptune can be replenished quickly enough with low random velocity planetesimals from elsewhere in the nebula (Bryden, Lin, & Ida 2000). Thommes, Duncan, & Levison (1999, hereafter TDL99) develop an alternative model to in situ formation for the origin of Uranus and Neptune. Beginning with four or more planetary embryos of 1015 M in the Jupiter-Saturn region, they explore through N-body simulation the evolution of the system after one of these bodies (and in one case, 2862
two at the same time) accretes a massive gas envelope to become a gas giant. They find that the remaining giant protoplanets are predominantly scattered outward. Dynamical friction with the planetesimal disk subsequently recircularizes their orbits, which leads in about half the simulations performed to a configuration quite similar to the present outer solar system, with the scattered giant protoplanets taking the roles of Uranus, Neptune, and Saturn. These results suggest that Uranus and Neptune are actually potential gas giant cores that formed in the same region as Jupiter and Saturn but lost the race to reach runaway gas accretion. Here we explore this model in more detail and perform further simulations: xx 2 and 3 motivate our choice of initial conditions for the simulations; x 4 discusses the mechanism for transporting a proto-Uranus and proto-Neptune to the outer solar system. N-body simulation results are presented in x 5. The effect on the asteroid and Kuiper belts is discussed in x 6. We summarize and discuss our findings in x 7.
2. AVAILABLE MATERIAL IN THE JUPITER-SATURN REGION
Hayashi (1981) estimated the minimum primordial surface density of solids in the outer solar system to be 3=2 a g cm2 : 1 min a 2:7 5 AU The requirement that the cores of Jupiter and Saturn be massive enough to have initiated runaway gas accretion suggests that they are $10 M in mass (Mizuno, Nakazawa, & Hayashi 1978; Pollack et al. 1996). Interior models are consistent with such a core mass but also allow a coreless Jupiter (Guillot 1999). We assume in this work that both Jupiter and Saturn began as $10 M bodies; putting gas giant cores
URANUS AND NEPTUNE FORMATION and ice giant cores on the same footing is necessary in the `` strong '' version of our model, although it is not essential to the basic mechanism; we discuss variations on the model in x 7. A surface density of 2.7 g cm2 was likely too low to form an $10 M body at 5 AU before the gas was removed from the protoplanetary disk. Lissauer (1987) finds that a surface density of 1530 g cm2 is needed to allow formation of Jupiter's core on a timescale of 5 105 5 106 yr, while the model of Pollack et al. (1996), which includes concurrent accretion of solids and gas, produces Jupiter in less than 107 yr with 10 g cm2. The formation of giant planet cores may have been triggered at least in part by the enhancement in the solids surface density beyond the `` snow line,'' where water goes from being a gaseous to a solid constituent of the protoplanetary disk. In fact, outward diffusion and subsequent freezing of water vapor from the inner solar system may have resulted in a large local density enhancement around 5 AU, perhaps yielding a surface density even higher than 30 g cm2 (Stevenson & Lunine 1988). Here we assume a power-law surface density, a a 0 : 2 5 AU The above discussion suggests $ 1030 g cm2 as a plausible surface density at 5 AU. Allowing also the exponent to vary between 1 and 2, one obtains a total mass in the Jupiter-Saturn region in excess of 40 M , and as high as 180 M . Therefore, it is likely that the region originally contained significantly more solids than ended up in the cores of Jupiter and Saturn.
3. IMPLICATIONS OF OLIGARCHIC GROWTH
2863
where a is the protoplanet semimajor axis, n is the spacing between adjacent protoplanets in units of rH, and p is the fraction of the total mass in the zone a nrH =2 < a < a nrH =2 incorporated into the protoplanet (Kokubo & Ida 2000). Using numerical simulations, Kokubo & Ida (1998) show that n $ 510. Models of giant planet formation by concurrent planetesimal and gas accretion suggest a surface density profile of / a2 (Pollack et et 1996). Using this in equation (3), one obtains a protoplanet mass independent of semimajor axis. Adopting 0 10 g cm2 and a spacing of n 7:5rH, one must then set the accreted mass fraction to about 0.8 in order to obtain a protoplanet mass of 10 M . With this spacing, between three and five such bodies fit between 5 and 10 AU. It is likely, therefore, that this region originally contained more than just the future solid cores of Jupiter and Saturn. At the same time, the more recent simulations of Thommes et al. (2002) indicate that accretion is quite inefficient in the outer solar system; even in the Jupiter-Saturn region, the above value of p is probably overly optimistic, and a correspondingly higher density of solids in this region may have been necessary to form gas giant cores. But as we find in x 5.5, a higher density disk actually tends to increase the `` success rate '' of the model.
4. SCATTERING AND CIRCULARIZATION OF THE PROTOPLANETS
Runaway growth (e.g., Wetherill & Stewart 1989; Kokubo & Ida 1996) transitions to a slower, self-limiting mode called `` oligarchic growth '' (Ida & Makino 1993; Kokubo & Ida 1998, 2000) when the largest protoplanets are still orders of magnitude less than an Earth mass everywhere in the nebula. Oligarchic growth has previously been demonstrated to take place only interior to about 3 AU (Weidenschilling & Davis 2000). Although Kokubo & Ida make estimates of protoplanet mass and growth timescales in the giant planet region, they point out that since their simulations are restricted to annuli that are narrow compared to their radii, they cannot make strong predictions about how oligarchic growth works over a wide range in the semimajor axis (Kokubo & Ida 2000). However, their most recent simulations (Ida & Kokubo 2000) span a range of 0.51.5 AU and show oligarchic growth proceeding in an outward-expanding `` wave '' over time. Also, Thommes (2000) and Thommes, Duncan, & Levison (2002) perform numerical simulations that suggest that oligarchic growth does take place in the outer solar system and proceeds on approximately the timescales predicted using the approach of Kokubo & Ida (2000). Assuming that protoplanets remain approximately evenly spaced in Hill radii, rH [1rH M=3 M 1=3 a], while growing--as is characteristic of oligarchic growth--their final mass is given by 1=2 2 3=2 p3=2 n3=2 3=2 a3 ; 3 M 2 3 M
In the nucleated instability picture of gas giant formation, there are three distinct phases of growth (Pollack et al. 1996): In phase 1 (which itself encompasses the subphases of runaway and oligarchic growth; see x 3 above), a solid core grows until it depletes most of the material in its feeding zone. Phase 2 is characterized by much slower growth, with gas accretion gradually coming to dominate over solids accretion. After a time on the order of several million years, the protoplanet contains comparable masses of gas and solids; around this time, the third phase of runaway gas accretion sets in and proceeds on a timescale of $105 yr. One of the protoplanets in the Jupiter-Saturn region must have been the first to reach this point. Shorter formation timescales at smaller heliocentric radii argue for Jupiter having formed before Saturn. The relatively long time spent by a protoplanet on the `` plateau '' of phase 2 means that several protoplanets could plausibly find themselves there by the time the first of them makes it to phase 3. Furthermore, if the winner has a broad enough margin of victory, its rivals will have only had time to accumulate a few M of gas, thus resembling the present-day ice giants in both solid and gas mass. A body increasing in mass from 10 M to Jupiter's mass, MJ 314 M , expands its Hill radius by a factor of about 3. The adjacent protoplanets will therefore have a high likelihood of being gravitationally scattered. In a three-body system--Sun, Jupiter, and a single protoplanet--a Jupitercrossing protoplanet would continue to have close encounters with its much more massive scatterer and so would remain coupled to it unless it were scattered onto an unbound orbit and left the system altogether. In reality, however, one expects that as it crosses beyond the Jupiter-Saturn region, a body will encounter a less accretionally evolved part of the protoplanetary disk, consisting predominantly of much smaller bodies. As a result, the protoplanet will experience dynamical friction. This will tend to
2864
THOMMES, DUNCAN, & LEVISON
Vol. 123
reduce the eccentricity of its orbit. If the eccentricity decays enough, the protoplanet's perihelion will be lifted away from Jupiter, thus decoupling it from its scatterer. Insofar as one can neglect interactions with other scattered protoplanets, its eccentricity will then monotonically decrease until it reaches equilibrium with the planetesimals. The eccentricity rate of change due to dynamical friction on a body of mass M in a swarm of mass m bodies can be expressed as (Weidenschilling et al. 1997) de2 M Cmhe2 i Me2 Ke ; m M dt
20 40 60 80 100
4
and a similar expression gives the time evolution of the inclination; he2 i1=2 is the rms eccentricity of the mass m bodies, m and Ke is a definite integral depending on the inclinations and eccentricities, defined in Stewart & Wetherill (1988). The coefficient C is given by C 16G2 sw1 ln v3 he2 i he2 i3=2 K 1 2 ; 5
1000
100
10 20 40 a (AU) 60 80 100
where G is the gravitational constant, vK the local Keplerian velocity, sw1 the spatial mass density of the swarm, and is approximately the ratio of the maximum encounter distance between the body M and a member of the population of the swarm, to the maximum separation that results in a physical encounter (Stewart & Wetherill 1988). An equilibrium is therefore reached when e2 M ' he2 im : M m This also means that as long as e2 M4he2 im ; M m 7 6
Fig. 1.--A (rough) approximation of the eccentricity-halving timescale due to dynamical friction, Tdf e=2 , for a 10 M body in the outer solar system having an initial eccentricity such that it is Jupiter crossing, with its pericenter at 5 AU [e0 1 5 AU=a]. Top panel: Time in years. Bottom panel: Time in units of the local orbital period. The planetesimal disk is taken to have a surface density of r 10r=5 AU2 g cm2 , an rms eccentricity em 0:05, and a thickness 2h $ 2aim $ aem ; ln is given a typical value of 3 (e.g., Kokubo & Ida 2000). In computing Tdf e=2 , we approximate the effective surface density of planetesimals encountered on an orbit with semimajor axis a as a, regardless of the orbit's eccentricity.
the eccentricity decay rate of the body M will depend only on the spatial mass density of the planetesimal disk and will be essentially independent of the individual masses of the planetesimals. This feature is exploited in the numerical simulations below. One can estimate the time for the eccentricity to be reduced by half as eM0 =2 e2 0 M : 8 Tdf e=2 $ deM0 =dt de2 0 =dt M Figure 1 plots this timescale estimate for a Jupiter-crossing 10 M body over a range in the semimajor axis. The planetesimal disk parameters used are those for the baseline simulation in x 5.1. As can be seen, the eccentricity-halving timescale reaches $2 Myr, or just over 104 orbital periods, at the outer edge of the present-day giant planet region (30 AU).
simulations to assess the effect of other planetesimal disk density profiles and of proto-Jupiter being initially not the innermost giant protoplanet. We also investigate how sensitive this scenario is to the timing of Saturn's final stage of growth, relative to that of Jupiter. Simulations are once again performed using the proven SyMBA symplectic integrator (Duncan, Levison, & Lee 1998). This integrator is able to handle close encounters among massive bodies while preserving the symplectic properties of the method of Wisdom & Holman (1991). 5.1. Initial Conditions: Set 1 For the baseline set of simulations, a planetesimal disk of surface density 2 a g cm2 9 1 10 5 AU is used, with the disk extending from 5 to 60 AU. The above density profile is steeper than those of TDL99 ( / a1 and a3=2 ). The total disk mass between 10 and 60 AU (106 M ) is about half that of the first two sets of runs in our earlier work (216 M ) and similar to that in the third run (119 M ). The truncation at 60 AU is to keep the number of bodies in the simulations tractably small; in reality the planetesimal disk may have extended for hundreds of AU, far beyond the presently visible Kuiper belt. Observational evidence does point to the possibility of a truncation of the belt at $50 AU (e.g., Chiang & Brown 1999). But if a `` Kuiper cliff '' exists, it is not necessarily primordial.
5. N-BODY SIMULATIONS
Putting the pieces together, one can now envision a scenario in which (1) a rapidly growing Jupiter scatters its smaller neighbors outward and (2) these `` failed cores '' are decoupled from Jupiter and ultimately evolve onto circular, low-inclination orbits in the outer solar system. TDL99 performed three series of numerical simulations to test this model; about half of the simulations produced a stable final configuration qualitatively similar to the present-day solar system. Here we expand on the earlier work, performing
No. 5, 2002
URANUS AND NEPTUNE FORMATION
2865
Since we truncate the disk at a location where the planetesimal surface density is already very low (0.07 g cm2 for the density profile of eq. [9]), one expects this not to have a strong effect on the evolution of any scattered protoplanets that cross the region r > 60 AU. As a test, we perform four pairs of runs comparing the evolution of a 10 M body with initial a 100 AU, e 0:90.95 in planetesimal disks with a surface density as prescribed by equation (9), in one case extending to 60 AU and in the other to 200 AU. No systematic difference in the 10 M body's semimajor axis evolution over the first few million years is apparent between the two different disk sizes. The eccentricity, however, tends to be damped more rapidly in a 200 AU disk. We expect, therefore, that the subsequent runs somewhat underestimate the effectiveness with which those scattered protoplanets that cross beyond 60 AU are circularized, if the disk was in fact larger. In particular, a planetesimal disk with a radius of hundreds of AU may allow retention of protoplanets that in our runs are scattered strongly enough to become unbound; we will explore this possibility in future work. The simulated disk is made up of equal-mass `` planetesimals,'' each having a mass of 0.2 M . At twice the mass of Mars, these bodies far exceed the actual characteristic mass of planetesimals in the early outer solar system. In reality, planetesimals are thought to have been on the order of 1 100 km in size, and thus with a mass of $1012 to 106 M (e.g., Lissauer 1987). The unrealistically large masses are chosen to keep the number of bodies manageably low, at slightly over 500. As mentioned in x 4, when the large bodies have eccentricities high enough that equation (7) is satisfied, the eccentricity decay timescale will be effectively independent of the planetesimal masses. Thus, despite the large planetesimal masses, the initial strength of eccentricity damping will be realistic. The equilibrium eccentricity condition equation (6), however, does depend on the planetesimal mass, so the equilibrium eccentricity reached by a large body among the planetesimals in the simulation will tend to be unrealistically high. Of course, a true equilibrium between protoplanets and planetesimals will not exist anyway, since mutual perturbations among the protoplanets will also have an effect. The initial planetesimal eccentricities and inclinations are given a Rayleigh distribution in e and i (Kokubo & Ida 1992), with he2 i1=2 0:05, hi2 i1=2 0:025 1=4. In the numerical integrations, the planetesimals are treated as a `` second-class,'' nonself-interacting population. Thus, they are perturbed in their Keplerian orbits only by forces from the protoplanets, not from each other. The protoplanets, on the other hand, are subject to forces from each other as well as from the planetesimals. This serves two purposes: It makes the simulations run much faster, since for N second-class bodies the computation time scales as N instead of as N 2 . Also, it prevents unrealistically strong selfstirring of the disk. Of course, not modeling planetesimal interactions means that collective planetesimal effects are not accounted for. Wave phenomena could have had an important effect on the evolution of the planetesimal disk velocity distribution (e.g., Ward & Hahn 1998), provided the disk was sufficiently massive and dynamically cold. However, it is unlikely that significant wave phenomena could persist once the initial scattering has taken place and the planetesimal disk has been stirred by eccentric 10 M bodies.
Also unmodeled is nebular gas, as a source of either aerodynamic drag (relevant for small planetesimals; e.g., Adachi, Hayashi, & Nakazawa 1976) or tidal forces (relevant for bodies e0.1 M ; e.g., Ward 1997). The former mechanism will keep the planetesimal disk more dynamically cold, but this does not affect the strength of dynamical friction on larger bodies until e2 M $ e2 m. The latter effect is thought m M to cause the inward migration of $110 M objects on timescales short compared to their formation times. This of course constitutes a potentially severe problem not only for our model but for any model of giant planet formation that requires the accumulation of large solid cores. Addressing this problem is beyond the scope of our present work, but we summarize some reasons why it may in reality have been less severe in x 7. The 10 M bodies are given densities of 0:25 1:5 g cm3 , roughly equal to those of Uranus and Neptune. This gives them radii of 2:18 104 km. Four such bodies are put in the simulation, initially on nearly circular and uninclined orbits. The orbits are spaced by 7.5 mutual Hill radii, starting from an innermost distance of 6 AU to allow for later inward migration by Jupiter. Thus, the bodies' initial semimajor axes are 6.0, 7.4, 9.0, and 11.1 AU. Between 5 and 12 AU, the disk is depleted in planetesimals so that the surface density is still given by equation (9). Since the large bodies are spaced proportionally to their semimajor axes, their distribution is consistent with a surface density proportional to a2. Integrating the surface density from 5 to 12 AU, the total mass is 51.5 M . With 40 M of this in the large bodies, this leaves 10.5 M in planetesimals in the Jupiter-Saturn region, in addition to 94.7 M between 12 AU and the outer disk edge at 60 AU. In summary, then, the initial conditions amount to a state where phase 2 of giant planet formation has been reached between 5 and 12 AU, with $80% of the planetesimals having accreted through oligarchic growth into four bodies of 10 M each, while no large bodies have yet formed beyond this region (Fig. 3). The adoption of equal-mass protoplanets is a simplification; apart from stochastic variations among the oligarchic growth end products, one can expect some intermediate bodies, with masses perhaps up to a few M . Such bodies are likely to ultimately be cleared, along with the planetesimals, from the giant planet region. However, they may end up playing a role in the dynamics of the trans-Neptunian region; see Petit, Morbidelli, & Valsecchi (1999). For computational reasons, the inner simulation radius is chosen as 1 AU; any body whose orbit penetrates this boundary is eliminated from the system. This is done because a limitation of the SyMBA integrator used here is its inability to handle close perihelion passages. Although a new version of SyMBA has since been developed that removes this restriction (Levison & Duncan 2000), the runs presented here predate this development. In any case, this limitation has little relevance for the present runs; typically, less than 10 planetesimals are lost at the inner boundary over the course of a 5 Myr simulation. The base time step is chosen as 0.05 yr, giving 20 steps per orbital period for an orbit with its semimajor axis at the inner radius and over 200 steps per orbit for Jupiter. Experimentation shows that this time step is small enough that the energy of the system is well conserved. Runs initially go to 5 Myr; in cases where the system still appears to be undergoing rapid evolution at this point, the runs are extended by another 5 Myr.
2866
0.8
THOMMES, DUNCAN, & LEVISON
0.5 0.4
Vol. 123
0.6 0.3 0.4 0.2 0.2 0.1 0 0 0 40 15 30 10 20 5 10 20 40 60 80 0 20 a (AU) 40 60
0 0 20 40 Semi-major axis (AU) 60 80
0 0 20 a (AU) 40 60
Fig. 2.--Eccentricities (top) and inclinations (bottom) in the outer solar system at the present epoch, showing the giant planets as well as all KBOs and Centaurs (objects with d30 AU) that have been observed at multiple oppositions as of 2001 September. Curve in the top panel shows the locus of orbits with perihelia at the semimajor axis of Neptune. KBO and Centaur data are taken from the Minor Planet Center site (Williams 2002).
Fig. 3.--Initial state for runs in set 1, showing eccentricity (top) and inclination (bottom) vs. semimajor axis. Larger circles denote the four 10 M protoplanets, and each of the small dots represents a 0.2 M planetesimal. Planetesimal density in the vicinity of the protoplanets is decreased to keep the density of protoplanets plus planetesimals consistent with the surface density given by eq. (9).
5.2. Set 1 Results To model gas accretion, SyMBA was modified to allow a subset of bodies to have artificially time-varying masses (in addition to any changes in mass resulting from the accretion of other bodies). For the runs in subset 1, it is assumed that the innermost protoplanet undergoes runaway gas accretion first and grows into Jupiter. The simulations start at the time when this happens, which should be a few million years into the life of the solar system (see x 4). Runaway gas accretion is simulated by increasing the body's mass over the first 105 yr of simulation time, from its original mass of 10 M to 314 M , approximately the present Jupiter mass. A linear growth in mass is used; this is deemed appropriate since the actual time evolution of mass during runaway growth is highly uncertain. Also, as the simulations show, 105 yr is roughly the response time of the system, and the system's subsequent evolution is therefore unlikely to be affected in a systematic way by the exact form of the time evolution of the runaway phase mass growth. Set 1 consists of eight alternate realizations of a run, differing only in the initial phases of the four protoplanets; for each, the angles , !, and M are randomly generated. As is seen, the stochasticity of the system ensures that this difference in phases is sufficient to bring about a very different evolution in each of the versions of the run. Figure 4 shows the evolution of one of the eight runs, denoted as run 1F, which after 5 Myr produces final protoplanet orbits that bear a particularly close resemblance to those of the giant planets in the present-day solar system. Semimajor axis, perihelion distance, and aphelion distance are plotted versus time for each of the four protoplanets; the innermost one has grown into Jupiter after the first 105 yr. By this time, the protoplanet orbits begin to mutually cross and strong scattering occurs. The protoplanet plotted in red
briefly has its semimajor axis increased to greater than 100 AU. However, dynamical friction acts to reduce eccentricities universally, decoupling the protoplanets from Jupiter and from each other. By about 1.2 Myr, none of the protoplanets are on crossing orbits anymore. After about 3 Myr, the bodies no longer undergo any changes in semimajor axis greater than a few AU on a million year timescale. At this point the orbits are well spaced and all eccentricities are d0.05, with no large fluctuations. Because of the large stochastic variations among runs, the exact timescales differ, but we find that orbits typically become noncrossing after less than 5 Myr in these and subsequent runs. This is consistent with Figure 1, which predicts timescales on the order of millions of years for the circularization of bodies scattered into the outer solar system. Subsequent semimajor axis evolution proceeds by scattering of planetesimals by protoplanets rather than scattering of the protoplanets off each other and Jupiter. As planetesimals are scattered among Jupiter and the protoplanets, the former experiences a net loss of angular momentum while the latter experience a gain. Thus, Jupiter's orbit shrinks, while those of the protoplanets expand (Fernandez & Ip 1996; Hahn & Malhotra 1999). This phase takes place over a timescale of several tens of millions of years and ends, at the very latest, when the planetesimals have been cleared from among the planets. Because the length of migration during this phase is only a few AU and in order to save time, most of the runs are stopped after 5 Myr. As an example, Figure 5 shows run 1F continued to 50 Myr; the net migration of the outer two protoplanets subsequent to 5 Myr is only $14 AU outward. The semimajor axes at which the scattered bodies end up are very noteworthy if one compares them to the present orbits of the giant planets (Fig. 2). At 5 Myr, the outer two
No. 5, 2002
URANUS AND NEPTUNE FORMATION
2867
Fig. 4.--Run 1F: Evolution of semimajor axis (bold lines), perihelion distance q (thin lines) and aphelion distance Q (dotted lines) of the four 10 M protoplanets. The protoplanet that grows to Jupiter mass (314 M ) over the first 105 yr of simulation time is shown in black.
protoplanets are at 16 and 26 AU. The innermost one is at 11 AU, while Jupiter is at 5.5 AU. This configuration of orbits is very similar to that of the present solar system, where Uranus and Neptune are at 19 and 30 AU, respectively, Saturn is at 9.6 AU, and Jupiter is at 5.2 AU. And at 5 107 yr, after some more net outward migration, the outer two protoplanets' semimajor axes are even closer to those of Uranus and Neptune (although `` proto-Saturn,'' having also moved outward, is further away from its present orbit). Eccentricities and inclinations are likewise very close to their present values. The end states of all the runs in set 1 are summarized in Figure 6. Depicted are snapshots of eccentricity versus semimajor axis at 5 Myr, except for runs 1A and 1H, which were continued to 10 Myr. These were extended because one or more of the protoplanets still had a high but decreasing eccentricity at 5 Myr. Runs 1A, 1D, 1F, 1G, and 1H result in a final ordering of orbits that is at the very least broadly consistent with the present solar system: Jupiter is the innermost body, with the
other three bodies interior to the region of the Kuiper belt and eccentricities low enough that no protoplanets cross each other or Jupiter. Out of these five runs, 1D and 1F in particular resemble the present solar system. Of course, to actually reproduce the solar system, another important event has to take place: the next protoplanet beyond Jupiter must also undergo a runaway gas accretion phase to acquire an envelope of mass e80 M and become Saturn. With the innermost core on a stable orbit in the vicinity of Saturn's present location though, as it is in runs 1D and 1F, the time delay between the runaway phases of Jupiter and Saturn appears to not be strongly constrained. A more detailed investigation of the role of Saturn's growth follows in x 5.8. In contrast, Uranus and Neptune must somehow have been prevented from later reaching runaway gas accretion. One possibility is that they simply ran out of time, still caught in the long plateau of the second giant planet growth phase (Pollack et al. 1996; see x 4) when the nebular gas was removed. Alternatively, it may be that the gas disk was truncated early on by photoevaporation between the orbit of
2868
THOMMES, DUNCAN, & LEVISON
Vol. 123
Fig. 5.--Run 1F continued to 50 Myr. Between 5 and 50 Myr, the net migrations for Jupiter and the protoplanets, going from inside to outside in semimajor axis, are 0.2, 1.5, 4, and 1.3 AU.
Saturn and the eventual orbit of Uranus (Shu, Johnstone, & Hollenbach 1993). This possibility is discussed further in x 7. The final state of the planetesimal disk differs substantially among the runs. In 1B and 1C, the planetesimal disk is largely unperturbed over most of its radial extent, while in 1F and 1H, eccentricities have been greatly increased throughout the entire disk. The rest of the runs are intermediate cases. The radial extent of the perturbation simply depends on how much of the disk is crossed by the protoplanets over the course of the run. In run 1F, for example, Figure 3 shows that one body's aphelion spends some time beyond 100 AU. On the other hand, in 1C, no protoplanet's aphelion ever goes further out than Qmax ' 18 AU. In those runs where Qmax < 60 AU (the disk radius), Qmax corresponds closely to where the planetesimal disk makes a transition from perturbed to unperturbed. The outer limit of the planetesimal disk's excitation by scattered protoplanets has been referred to as the `` fossilized scattered disk '' (TDL99). The contemporary scattered disk, by contrast, consists of
objects that were scattered by Neptune after the latter had attained its current orbit (Duncan & Levison 1997). A more detailed discussion of the fossilized scattered disk follows in x 6.2. Previous work has shown that a slow, steady outward migration of Neptune, like that which would accompany the scattering of planetesimals among the giant planets (see above), results in the capture of planetesimals into Neptune's exterior mean motion resonances (Malhotra 1995). One would expect the same thing to happen in the model presented here once the giant planets are no longer strongly perturbing each other and can migrate smoothly. None of the simulations, however, show any evidence of resonant capture of planetesimals. The problem is that the exaggerated planetesimal masses make the migration too jittery for any planetesimals to be entrained in the resonances. The simulations of Hahn & Malhotra (1999) clearly show the trend of smoother migration for smaller planetesimal masses; however, even with their smallest mass of 0.01 M , they find no resonant capture; the disk is still too discretized.
No. 5, 2002
URANUS AND NEPTUNE FORMATION
2869
1
1
B
0.5
0.5
0 0 1 C 0.5 20 40 60
0 0 1 D 0.5 20 40 60
0 0 1 E 0.5 20 40 60
0 0 1 F 0.5 20 40 60
0 0 1 G 0.5 20 40 60
0 0 1 20 40 60
0.5
0 0 20 40 60
0 0 a (AU) 20 40 60
Fig. 6.--End states of the eight set 1 runs after 5 Myr of simulation time, except for 1A and 1H, which were continued on to 10 Myr. Eccentricity is plotted vs. semimajor axis. Three different sizes of points denote planetesimals (smallest), 10 M protoplanets (medium), and Jupiter (largest). Planetesimal orbits crossing Jupiter or any of the protoplanets are generally unstable on timescales short compared to the age of the solar system; thus, the region among the protoplanets would be essentially cleared of planetesimals long before the present epoch.
However, in our model, there is the additional feature that the planetesimal disk tends to get strongly stirred during the initial scattering of the giant protoplanets. Since the probability of a body to be resonantly captured decreases with its initial eccentricity (e.g., Borderies & Goldreich 1984), the subsequent capture efficiency will be lower than if Neptune's resonances sweep through an unperturbed disk, even if Neptune's migration is perfectly smooth in both cases. Three of the runs produce systems at 5 Myr that are irreconcilably different from our own. In 1B, one of the protoplanets has merged with Jupiter. In 1C, a protoplanet has been scattered onto an orbit interior to Jupiter, in the region of the present-day asteroid belt, with its semimajor axis at 3.4 AU, its perihelion at 2.8 AU, and its aphelion at 3.9 AU. It attains a stable orbit not crossing Jupiter even though the only planetesimals available for damping interior to Jupiter
are those few that are also scattered there. This is possible because the protoplanet is scattered not just by Jupiter but also by the other protoplanets. Even assuming the protoplanet could be subsequently removed from this region, such an event would most likely have cleared much of the asteroid belt. Finally, in 1E, two of the protoplanets have merged. It should be noted that mergers--or, indeed, ejections--that reduce the number of protoplanets are not intrinsically a problem, since extra ones may have existed. Scenarios with five protoplanets (Jupiter+Saturn+3) are explored in x 5.4. 5.3. Set 2: Dependence on Initial Ordering How strongly does the final configuration of the system depend on the initial ordering of the protoplanets? The
2870
THOMMES, DUNCAN, & LEVISON
Vol. 123
next set of runs, set 2, uses the same initial conditions as set 1, within a random variation in the protoplanet phases. However, for these runs it is the second-innermost protoplanet, rather than the innermost one, that undergoes simulated runaway gas accretion. One can reasonably expect that this will favor an outcome like 1C (Fig. 6), where a protoplanet is scattered inward into the region of the asteroid belt. The end states of the runs are shown in Figure 7. In six of the eight runs, a protoplanet has indeed ended up interior to Jupiter. However, in two cases (2C and 2D), Jupiter is the innermost body. Thus, it appears that if Jupiter does not grow from the innermost protoplanet, the likelihood of ending up with a final configuration similar to our solar system declines, although such an outcome continues to be quite possible.
5.4. Set 3: Dependence on Number of Cores How sensitively does the end state depend on the initial number of core-sized bodies? In the next set of simulations, an extra 10 M protoplanet is added. All protoplanets are more tightly spaced, by 6.5 instead of by 7.5 mutual Hill radii. Starting again from 6.0 AU, the outermost protoplanet is therefore initially at 12.2 AU. The surface density of planetesimals in the region of the protoplanets is reduced to keep the average surface density unchanged at 1 10a=5 AU2 g cm2. The innermost protoplanet is, again, increased in mass to that of Jupiter over the first 105 yr. Figure 8 shows the end states of the runs. This time all runs initially have a length of 107 yr, since the larger number of bodies take longer to decouple from each other. Run 3A
1
1
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 20 40 60
0 0 a (AU) 20 40 60
Fig. 7.--End states of the eight set 2 runs (second-innermost protoplanet becomes Jupiter) after 5 Myr of simulation time, except for 2C, which was continued on to 15 Myr. Eccentricity is plotted vs. semimajor axis. Three different sizes of points denote planetesimals (smallest), 10 M protoplanets (medium), and Jupiter (largest).
No. 5, 2002
URANUS AND NEPTUNE FORMATION
2871
1
1
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 20 40 60
0 0 20 40 60
Fig. 8.--End states of the eight set 3 runs (one extra protoplanet) after 10 Myr of simulation time, except for 3A, which was continued on to 15 Myr
is continued to 1:5 107 yr because 107 after yr some of the protoplanets are still on crossing orbits. Eccentricities are uniformly low, and the protoplanet orbits are well spaced for the most part; thus, the systems have a good chance of being stable indefinitely. However, once all the planetesimals have been scattered from among the protoplanets so that the latter are no longer subject to dissipative forces, some of these systems may still become unstable. This caveat applies to all the runs presented in this work, but generally speaking, larger numbers of bodies increase the potential for instability (Levison, Lissauer, & Duncan 1998). All of the protoplanets remain in six of the eight runs, thus resulting in systems with one too many planets relative to the present solar system. However, in Runs 3B and 3F, one of the protoplanets is ejected from the solar system, leaving the right number of bodies behind. In both these runs, a protoplanet ends up with a semimajor axis within
10% of present-day Uranus and Neptune, respectively, although both `` Saturns '' are too far out. One may conclude that with one extra initial protoplanet in the JupiterSaturn region, scattered protoplanets continue to be readily circularized and the resulting systems tend to look like ours with one extra outer planet. However, a system with four giant planets remains a possible outcome. Also, the subsequent formation of Saturn from one of the protoplanets may trigger more ejections, especially in those cases where the inner protoplanets are more closely spaced, such as 3D and 3H. 5.5. Set 4: A More Massive Planetesimal Disk In this set of runs, a number of parameters are changed to simulate a system with a more massive planetesimal disk. The protoplanets are now 15 M bodies. The planetesimal disk surface density profile still varies as a2 , but it is scaled
2872
THOMMES, DUNCAN, & LEVISON
Vol. 123
up to be 15 g cm2 at 5 AU; that is, 2 a 2 15 g cm2 : 5 AU
10
Other minor differences are a legacy of chronologically earlier runs. The innermost protoplanet is initially at 5.3 AU, and successive protoplanets are spaced by only 5.8 mutual Hill radii. Thus, the outermost protoplanet is initially at 9.0 AU. The individual planetesimals have a mass of 0.24 M . The innermost protoplanet, as before, has its mass increased to 314 M over the first 105 simulation years. The end states of the runs are shown in Figure 9. All except C yield the correct number and ordering of bodies, and eccentricities are uniformly low. This `` success rate '' is higher than that of set 1, in which only five out of eight runs yield qualitatively the correct orbital configuration. This is accounted for by the more massive planetesimal disk; it pro-
vides stronger dynamical friction so that scatterings of protoplanets tend to be less violent, and subsequently, orbits tend to be circularized and mutually decoupled more quickly. Runs A, D, and H end up with protoplanet orbits that are particularly close to those of Saturn, Uranus, and Neptune. Jupiter systematically ends up at too small a heliocentric distance, indicating that the initial distance of 5.3 AU is too small. Also, the larger disk mass gives the protoplanets and Jupiter more planetesimals to scatter and thus increases the distance they travel due to angular momentum exchange (x 5.1). 5.6. Set 5: A Shallower Disk Density Profile In this set of runs, a shallower planetesimal disk surface density, varying as a3=2 , is used. The disk now begins at 10 AU, with the interior region being initially occupied solely by the protoplanets. In other words, it is assumed that in the epoch from which the runs start, all but a negligible mass of
1
1
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0 0 20 40 60
0 0 20 40 60
Fig. 9.--End states of the eight set 4 runs (a more massive planetesimal disk). First four are run to 107 yr; the last four are run to 5 106 yr.
No. 5, 2002
URANUS AND NEPTUNE FORMATION
2873
the planetesimals among the protoplanets has been swept up or scattered from the region. The surface density profile is given by 3=2 a g cm2 ; 11 3 1:8 10 AU and thus a total mass of 123 M is contained in the disk between 10 and 60 AU. Extending this planetesimal disk inward to 5 AU would yield a surface density there of only 5 g cm2 and a total mass between 5 and 10 AU of only 25 M . Despite this, we still put four 15 M protoplanets in this region; thus, it is assumed that the original planetesimal surface density profile in this region was steeper, perhaps due in part to redistribution of water vapor from other parts of the disk to the vicinity of the snow line (Stevenson & Lunine 1988). This set consists of twelve runs, each to 5 Myr. As before, the inner body's mass is increased to 314 M over the first
105 yr of simulation time. The end states are shown in Figure 10. Eight out of 12 runs possess the right number and ordering of bodies. Eccentricities of the protoplanets and Jupiter are d0.1 in five of these. This `` success rate '' is close to that of set 1 (which has a similar total disk mass); this model thus does not appear to be highly sensitive to changes in density profile alone. Various degrees of disruption of the planetesimal disk can again be seen. Most show a sharp transition between a disrupted region crossed by the scattered planetesimals and a largely undisturbed outer region. In run 5A, for example, this transition occurs at slightly below 40 AU, while in run 5F it is located between 45 and 50 AU. Runs 5B and 5L show strong disruption throughout the entire disk, indicating that all of it was crossed by one or more protoplanets. Another state can be seen in run 5G, where all eccentricities in the outer part of the belt are uniformly raised. This occurs when a protoplanet crosses the outer disk with an inclination high enough that it spends most of its orbit above or
1
1
1
0.5
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0.5
0 0 1 20 40 60
0 0 1 20 40 60
0 0 1 20 40 60
0.5
0.5
0.5
0 0 20 40 60
0 0 20 40 60
0 0 20 40 60
Fig. 10.--End states of the 12 set 5 runs (a shallower disk density profile). All are run to 5 106 yr.
2874
THOMMES, DUNCAN, & LEVISON
Vol. 123
below, rather than inside, the disk. The disk planetesimals are then excited primarily by long-range secular effects rather than by short-range scattering encounters (TDL99). Another (less dramatic) example of this effect can be seen in run 1G above. 5.7. Set 6: A Minimum Mass Disk How well does the model work in a `` minimum mass '' planetesimal disk with a surface density as given by equation (1)? Such a disk contains only about 13 M in solids between 5 and 10 AU. This is inconsistent with the assertion that the Jupiter-Saturn region is the source of all the giant planets, unless one invokes a very large density enhancement in that area. Nevertheless, it is useful to investigate how effective dynamical friction is in this extreme case. We repeat the baseline run (set 1; x 5.2) but with a planetesimal disk of the surface density given by equation (1). All eight runs are stopped after 10 Myr. The end states are shown in
Figure 11. In only two of the eight runs--6D and 6G--are all protoplanets retained. However, only in 6G are all the protoplanets noncrossing by 10 Myr. One may conclude that the circularization of scattered protoplanets is marginally effective even in a minimum mass disk. Therefore, this model places a stronger lower limit on the disk mass in the Jupiter-Saturn region than in the trans-Saturnian region. 5.8. Set 7: The Role of Saturn Thus far, the only gas giant in the simulations has been Jupiter. The innermost of the scattered protoplanets does tend to end up near the present location of Saturn. However, to reproduce the solar system, it must at some point accrete $80 M of nebular gas. We therefore investigate what effect the subsequent growth of a Saturn-mass object has on our model. A set 1 run that produced a good solar system analog (1F; Fig. 3) is used as a starting point. At about 1 Myr, the proto-
1
A
1
B
0.5
0.5
0 0 1 C 0.5 20 40 60
0 0 1 D 0.5 20 40 60
0 0 1 E 0.5 20 40 60
0 0 1 F 0.5 20 40 60
0 0 1 G 0.5 20 40 60
0 0 1 H 0.5 20 40 60
0 0 20 40 60
0 0 a (AU) 20 40 60
Fig. 11.--End states of the eight runs in set 6, with a minimum mass initial planetesimal disk. All are run to 10 Myr.
No. 5, 2002
URANUS AND NEPTUNE FORMATION
2875
planets and Jupiter are no longer on crossing orbits. At this time, the initially outermost protoplanet (plotted in blue) has become the innermost one, closest to Jupiter at $10 AU. We perform a set of five runs that branch off from this point. In these runs, the innermost protoplanet has its mass increased to that of Saturn over a 105 yr interval, starting at 1, 1.2, 1.4, 1.5, and 1.6 Myr, respectively. We choose to start only after the protoplanets are on noncrossing orbits in order to avoid ending up with an eccentric `` Saturn ''; the eccentricity evolution of an initially eccentric giant planet in a gas disk is uncertain (Lin et al. 2000). The end states are shown in Figure 12. No protoplanets have been lost from the system by the end of the runs. A protoplanet still has a high eccentricity in the 1 Myr case; this is because this protoplanet's perihelion is still very close to the innermost protoplanet's orbit at 1 Myr, and it suffers strong perturbations as the latter grows to Saturn's mass. In the
other cases, however, Saturn's formation does not cause large eccentricities in the protoplanets. This is as one would expect; at its final mass and at 10 AU Saturn's Hill radius is 0.45 AU, and by 1:2 106 yr, the closest protoplanet's perihelion is at $14 AU, almost 9rH away, and is thus unlikely to be in reach of strong scattering. We leave for future work the effect of Saturn's gas accretion on systems with more protoplanets, such as those in set 3 (x 5.4). Such systems tend to have weaker stability and thus may be more susceptible to disruption by Saturn's final growth spurt. One effect visible in Figure 12 is that the semimajor axis of Saturn at 5 Myr tends to be smaller than that of the innermost scattered protoplanet--the putative proto-Saturn--in those runs where the gas accretion of Saturn is not modeled (see, e.g., Fig. 6). When a protoplanet grows to Saturn's mass, its subsequent migration speed is much slower, since the rate of migration depends on how much mass in planet-
1 0.8 0.6 0.4 0.2 0 0 20 40 60
1 0.8 0.6 0.4 0.2 0 0 20 40 60
1 0.8 0.6 0.4 0.2 0 0 20 40 60
1 0.8 0.6 0.4 0.2 0 0 20 40 60
1 0.8 0.6 0.4 0.2 0 0 20 40 60
Fig. 12.--End states of the set 7 runs, in which `` Saturn '' commences growing after the protoplanets have largely decoupled from each other. In each case, Saturn is the second innermost of the largest points, the innermost being Jupiter. The start time of Saturn's growth, tS , is denoted on each panel.
2876
THOMMES, DUNCAN, & LEVISON
Vol. 123
esimals it scatters relative to its own mass. This counteracts the tendency of the innermost scattered protoplanet to end up at a semimajor axis larger than that of present-day Saturn in the other runs, where only Jupiter grows.
6. SCATTERED PROTOPLANETS AND THE SMALL-BODY BELTS
The simulations presented above show that scattering of giant planet coresized protoplanets is a violent event that leaves a strong dynamical signature on the surrounding planetesimal disk. The asteroid belt and the Kuiper belt are therefore the natural places to look for evidence of large scattering events in the solar system's early history. 6.1. The Asteroid Belt It is those members of the asteroid belt larger than about 50 km in diameter that are of interest in inferring properties of the early solar system; smaller bodies cannot be primordial because they could not have survived intact for the age of the solar system (e.g., Petit et al. 1999). This population displays a high degree of dynamical excitation and a severe mass depletion, containing at most a thousandth of its original mass, as extrapolated from the terrestrial region (Weidenschilling 1977). In the context of the model presented here for the early evolution of Uranus and Neptune, the initial violent scattering of protoplanets is perhaps the most obvious candidate to look to for perturbation of the asteroid belt. Indeed, in numerous runs, protoplanets spend some time interior to the orbit of Jupiter, crossing part of the asteroid belt region for up to $104 yr at a time. One might expect that such an occurrence would wreak havoc in the asteroid belt. This is analogous to the model investigated by Petit et al. (1999), in which 1 M bodies are scattered by Jupiter to cross the Kuiper belt; the difference is that the bodies here are a factor of 10 more massive. To assess the effect of scattered Uranus/Neptunesized bodies on the asteroid belt, we perform six additional runs with planetesimals added in the asteroid belt region. The individual bodies have masses of 0:024 M and are distributed with a surface density of 1 a belt 8:0 g cm2 12 1 AU between 2.5 and 4.5 AU. This shallow density profile is the same as the one used by Chambers & Wetherill (1998) in the terrestrial region, which in turn was chosen to be more consistent with the large densities required at larger heliocentric distances to form Jupiter and Saturn. The protoplanet masses in this case are 15 M . We find that eccentricities can get excited to their present values in this way, although only down to the crossing protoplanet's minimum perihelion distance, which in none of the runs performed reaches the inner edge of the belt. Inclinations fare more poorly; protoplanets seldom raise them much above 10 , which is less than the present median inclination of asteroids beyond 2.5 AU. Also, very little mass is scattered out of the belt while the protoplanets are crossing it. Figure 13 shows `` before and after '' snapshots of the run that out of the six displays the strongest disruption of the region interior to Jupiter by scattered protoplanets. As can be seen, few bodies, apart from those that are nearly Jupiter crossing, attain inclinations
above 15 . Thus, despite the much more massive Jupiterscattered bodies, our result is similar to that of Petit et al. (1999): the bodies cross too little of the asteroid belt, and for too short a time, to reproduce the observed high inclinations and mass depletion. The scattered giant protoplanet model may have figured more indirectly in the dynamical sculpting of the asteroid belt region. If proto-Saturn was very close to Jupiter when it accreted its massive gas envelope--a situation favored by an initially compact region of giant planet formation--then the subsequent migration of the gas giants may have been enough to sweep the inclination-exciting 16 secular resonance through most of the asteroid belt (e.g., Gomes 1997; Levison et al. 2001). Gomes initially places Jupiter and Saturn at 5.4 and 8.7 AU, respectively, which positions the 16 resonance at 2.7 AU. As the gas giants migrate away from each other, the resonance sweeps inward, reaching 2 AU when they reach their present orbits. If Saturn were to complete its growth even closer to Jupiter, the 16 resonance would start at a larger heliocentric distance, and more of the asteroid belt would be swept during the subsequent migration. If large (lunar to martian sized) protoplanets originally accreted in the asteroid belt, then there is little or no need for any exterior sources of the belt's excitation and mass depletion; Petit, Morbidelli, & Chambers (2001) show that in this case both can be well reproduced. In their simulations, protoplanets in the belt scatter the surrounding small bodies and are eventually themselves removed through mutual scattering and perturbation by Jupiter. In this scenario, a limited role may still have been played by Jupiterscattered bodies in the outer part of the asteroid belt, where Jupiter's influence would have impeded the formation of belt protoplanets. 6.2. The Kuiper Belt In the present solar system, a new class of Kuiper belt object (KBO) has recently been identified (Duncan & Levison 1997; Luu et al. 1997). These objects have semimajor axes and eccentricities such that they lie near the locus of Neptune-encountering objects shown in Figure 5. They are thought to be part of a population referred to collectively as the scattered disk: formerly low-eccentricity KBOs that have had their orbits changed by close encounters with Neptune. Many of the simulations in x 5 show an analogous class of planetesimals in their `` Kuiper belt '' regions. However, these fall on the locus of orbits crossing not the final semimajor axis of the outermost protoplanet, but the furthest aphelion distance of any of the protoplanets during their initial high-eccentricity phase. Since these orbits are no longer being crossed by a protoplanet, they will be stable over long times. One can refer to these structures as `` fossilized '' scattered disks (TDL99) because they preserve part of the dynamical history of the planetesimal disk. Such structures appear only in runs where the initial scattering was strong enough that one or more protoplanets had their aphelia increased to well beyond the final semimajor axis of the (ultimately) outermost protoplanet. Observations of our solar system's Kuiper belt do indeed reveal an anomalously high degree of excitation (e.g., Petit et al. 1999; Malhotra, Duncan, & Levison 2000). The eccentricities and, to a lesser degree, inclinations of bodies in mean motion resonances with Neptune, particularly the 2 : 3
No. 5, 2002
1 0.8 0.6 0.4 0.2 0 0 1 2
URANUS AND NEPTUNE FORMATION
1 0.8 0.6 0.4 0.2 0 3 a (AU) 40 40 4 5 0 1 2 3 a (AU) 4 5
2877
30
30
20
20
10
10
0 0 1 2 3 a (AU) 4 5
0 0 1 2 3 a (AU) 4 5
Fig. 13.--Eccentricities and inclinations of planetesimals in the asteroid belt region, interior to Jupiter (large dot), at 3 104 yr (top panel) and 2 105 yr (bottom panel ). These times are, respectively, just before and just after the period during which a protoplanet repeatedly crossed the region interior to Jupiter. Jupiter in this run has moved inward to $4.8 AU, 0.4 AU less than its present semimajor axis. Curve marks the locus of Jupiter-crossing orbits
resonance at 39.5 AU, can be explained by resonance sweeping during Neptune's migration, as can the paucity of objects on nonresonant orbits interior to 39 AU (Malhotra 1995). However, the high inclinations found beyond $41 AU in what is commonly called the `` classical '' Kuiper belt, cannot be explained in this way. Petit et al. (1999) propose large (up to 1 M ) Neptune-scattered planetesimals as the mechanism that stirred and cleared the belt. But even when such bodies remain in the belt for 100 Myr, the inclinations they raise are almost always less than 20 . Can the excitation of the Kuiper belt be better accounted for if we are thus far only seeing the part of it that is interior to the locus of a fossilized scattered disk? The inclinations raised by a protoplanet can be directly obtained from the simulations of x 5. The runs of set 1, set 3, and set 5 are used for comparison. Their inclinations are shown in Figures 14, 15, and 16, respectively. Inclinations beyond the outermost protoplanet are excited up to a maximum of about 30 , similar to those observed in the classical Kuiper belt today (see Fig. 5). However, observations of the Kuiper belt are biased against high-inclination objects. Brown (2001) derives a debiased inclination distribution function for the classical
belt: i2 i2 ft i sin i a exp 2 1 a exp 2 ; 21 22 13
with a 0:93 0:02, 1 2:20:2 , and 2 17 3. One 0:6 can define a parameter i0 cos 1 cos i to give a measure of the characteristic inclination of a population. For the debiased distribution above, i0 21 . However, in the runs presented here, the largest i0 in the region corresponding to the classical belt (between the outermost planet's 2 : 3 and 1 : 2 mean motion resonance) is only 15 . Thus, although higher inclinations are produced here than in the large Neptune-scattered planetesimals model of Petit et al. (1999), the inferred full velocity distribution of the classical Kuiper belt still cannot be accounted for. It is possible to estimate with a simple numerical experiment if a planet as large as Uranus or Neptune can in principal excite the Kuiper belt to observed values. This experiment consists of a single Uranus-mass planet on an orbit with a 45 AU, e 0:2, and i 25 . The planet is embedded in a swarm of 400 massless test particles infor-
2878
THOMMES, DUNCAN, & LEVISON
Vol. 123
30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 20 40 60 20 40 60 20 40 60 20 40 60
30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 20 40 60 20 40 60 20 40 60 20 40 60
Fig. 14.--Counterpart to Fig. 6, showing inclination vs. semimajor axis for the end states of the runs in set 1 (baseline)
mally spread from 35 to 55 AU, with initial e 0:01 and i 1 . Figure 17 plots the i0 of the particles as a result of scattering off of the planet. It shows that a planet the mass of the ice giants can indeed excite the Kuiper belt to i0 20 in a million years. However, in none of the runs we performed was such a high inclination imparted to a protoplanet. Alternatively, a less inclined protoplanet may be able to reproduce the Kuiper belt inclinations if it remains on a belt-crossing orbit significantly longer than 1 Myr (and thus longer than in our runs). Right at the outset, though, there is a caveat: If one wants to appeal to a truncation in the gas disk beyond Saturn as the mechanism that forestalls the evolution of Uranus and Neptune into gas giants, then a circularization time of several million years may cause their pericenters to spend a perilously long time in the gas. In any case, one way to increase the dynamical friction timescale is by decreasing the mass of the trans-Saturnian planetesimal disk. However, this also increases the chances of protoplanets being lost
altogether during scattering, and makes it more difficult to (plausibly) fit enough solids for forming all the giant planets into the Jupiter-Saturn zone (see x 5.7). A longer circularization time will also result if a protoplanet can be decoupled from Jupiter and the other protoplanets at a larger semimajor axis. For the latter situation to have a better chance of occurring in simulations, the planetesimal disk needs to be extended to larger heliocentric distance. We will investigate this issue further in future work.
7. DISCUSSION
The conventional picture of the formation of Uranus and Neptune, whereby the ice giants accrete near their present heliocentric distances, has grave problems. Numerical simulations have not been able to produce $10 M objects in the trans-Saturnian region in the lifetime of the solar system without significantly increasing protoplanet radii or invoking dissipational forces and planetesimal disk densities too
No. 5, 2002
URANUS AND NEPTUNE FORMATION
2879
30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 20 40 60 20 40 60 20 40 60 20 40 60
30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 30 20 10 0 0 20 40 60 20 40 60 20 40 60 20 40 60
Fig. 15.--Counterpart to Fig. 8, showing inclina...