Terrestrial Planet and Asteroid Formation in the Presence of Giant Planets

Icarus 143, 60–73 (2000) doi:10.1006/icar.1999.6241, available online at http://www.idealibrary.com on

Terrestrial Planet and Asteroid Formation in the Presence of Giant Planets I. Relative Velocities of Planetesimals Subject to Jupiter and Saturn Perturbations Stephen J. Kortenkamp and George W. Wetherill Department of Terrestrial Magnetism, Carnegie Institution of Washington, 5241 Broad Branch Road, Washington, DC 20015 E-mail: [email protected] Received August 24, 1998; revised August 24, 1999

1. INTRODUCTION We investigate the orbital evolution of 1013 - to 1025 -g planetesimals near 1 AU and in the asteroid belt (near 2.6 AU) prior to the stage of evolution when the mutual perturbations between the planetesimals become important. We include nebular gas drag and the effects of Jupiter and Saturn at their present masses and in their present orbits. Gas drag introduces a size-dependent phasing of the secular perturbations, which leads to a pronounced dip in encounter velocities (V enc ) between bodies of similar mass. Planetesimals of identical mass have V enc ∼1 and ∼10 m s−1 (near 1 and 2.6 AU, respectively) while bodies differing by ∼10 in mass have V enc ∼10 and ∼100 m s−1 (near 1 and 2.6 AU, respectively). Under these conditions, growth, rather than erosion, will occur only by collisions of bodies of nearly the same mass. There will be essentially no gravitational focusing between bodies less than 1022 to 1025 g, allowing growth of planetary embryos in the terrestrial planet region to proceed in a slower nonrunaway fashion. The environment in the asteroid belt will be even more forbidding and it is uncertain whether even the severely depleted present asteroid belt could form under these conditions. The perturbations of Jupiter and Saturn are quite sensitive to their semi-major axes and decrease when the planets’ heliocentric distances are increased to allow for protoplanet migration. It is possible, though not clearly demonstrated, that this could produce a depleted asteroid belt but permit formation of a system of terrestrial planet embryos on a ∼106 -year timescale, initially by nonrunaway growth and transitioning to runaway growth after ∼105 years. The calculations reported here are valid under the condition that the relative velocities of the bodies are determined only by Jupiter and Saturn perturbations and by gas drag, with no mutual perturbations between planetesimals. If, while subject to these conditions, the bodies become large enough for their mutual perturbations to influence their velocity and size evolution significantly, the problem becomes much more complex. This problem is under investigation. °c 2000 Academic Press Key Words: planet formation; asteroids; giant planets; gas drag.

In this paper we address a variation in the “standard model” for the formation of the Solar System. Specifically, we begin an investigation to determine the consequences of the early formation of Jupiter and Saturn with regard to terrestrial planet and asteroid formation. The conventional “core-accretion” model for the formation of the giant planets (see review by Lissauer 1993) postulates that protoplanetary cores initially form by collisional accretion of planetesimals. These cores eventually reach a critical mass of ∼10 M⊕ whereby they can begin accreting massive gaseous mantles in a runaway phase of growth. One long-standing problem with the core-accretion model is that the timescale for a single 10-M⊕ core to form at Jupiter’s heliocentric distance, even under very optimistic initial conditions, will be on the order of 106 years with accretion of the ≈300-M⊕ gaseous mantle requiring ∼107 years (Lissauer 1987, Pollack et al. 1996). This may be close to the time at which the solar nebula is believed to have begun dissipating (Strom et al. 1993). A further problem with the core-accretion model has been pointed out by Weidenschilling (1997), who demonstrated that protoplanetary core formation may cease at masses near 2 to 3 M⊕ , well short of the 10 M⊕ thought to be needed for runaway gas accretion. Lissauer et al. (1995) illustrate other problems associated with the formation of giant planets by the core-accretion model. Despite these serious quantitative problems the core-accretion model has enjoyed rather broad acceptance and led to a general estimate of about 107 years for the time required for Jupiter and Saturn to attain their present masses. Wetherill and Stewart (1993) have shown that in 107 years accretion in the inner Solar System is well under way, with perhaps a few dozen Mercuryto Mars-sized planetary embryos distributed between 0.5 and 3 AU. Models pertaining to the evolution of the inner Solar System typically evolve these planetesimals for 5 to 10 million

60 0019-1035/00 $35.00 c 2000 by Academic Press Copyright ° All rights of reproduction in any form reserved.

61

PLANET AND ASTEROID FORMATION

years before introducing Jupiter and Saturn (see Chambers and Wetherill 1998). At the present time the evidence supporting the core-accretion model is too weak to exclude alternative hypotheses. One such alternative is the “disk-instability” model originally proposed by Kuiper (1951) and subsequently studied by Cameron (1978). The disk-instability model postulates that gravitational instabilities in the gaseous solar nebula lead to clumping of the nebula and eventual collapse of fragments into giant gaseous protoplanets (GGPPs). These GGPPs then evolve into giant planets like Jupiter and Saturn. The disk-instability model fell out of favor primarily because of two obstacles. First, the model had difficulty explaining how Jupiter and Saturn could obtain their massive cores of ice and rock. Earlier estimates of these cores were in the range 10 to 30 M⊕ for Jupiter and 15 to 25 M⊕ for Saturn (Stevenson 1982), masses which seemed at the time to corroborate the core–accretion formation theory. Second, the observed (atmospheric) metal abundances of the giant planets are substantially enhanced relative to solar abundances, which would not be expected if the giant planets formed by diskinstability directly from the solar nebula. Each of these obstacles to GGPPs may now have become less serious due to refined models of the interiors of Jupiter and Saturn (Chabrier et al. 1992, Guillot et al. 1994). These models predict smaller core masses of 3 to 10 M⊕ for Jupiter and 1 to 13 M⊕ for Saturn. The lower end of these mass ranges may be too low to trigger the runaway phase of gas accretion vital to the core-accretion model. These same interior models also suggest the presence of a radiative zone which prevents global mixing, such that the observed metal abundances once thought to be global may actually only reflect atmospheric abundances. (See Boss (1999) for a detailed comparison of the core–accretion and disk-instability theories.) Recent discoveries of extra-solar giant planets (Mayor and Queloz 1995, Marcy and Butler 1996, Butler and Marcy 1996, Butler et al. 1997) and the refined interior models of Chabrier et al. (1992) and Guillot et al. (1994) have led to a revival of sorts for the GGPP theory. Boss (1996) has shown that low mass protoplanetary disks accreting material from a spherical halo have hot inner disks (∼1200 K inside of 4 to 5 AU) and lower temperatures in the outer disk (∼100 K outside of 4 to 5 AU). He showed that these conditions are sensitive to the disk mass accretion rate but rather insensitive to the disk mass and central protostar mass. Boss (1997) later found that under these conditions the outer disk is unstable toward GGPP formation regardless of the thermodynamics assumed for the disk (locally isothermal or adiabatic). He has since shown (Boss 1998) that in an intermediate mass disk (∼0.14 M¯ inside 10 AU) as a single GGPP forms near 6 AU, if the location of the center of mass of the protostar/disk system is preserved (i.e., if the protostar is allowed to wobble as the first GGPP forms) then a second GGPP can be triggered near 10 AU.

Boss’s (1996) disk models led him to entertain the notion of a “best of both worlds” scenario for planet formation, whereby Jupiter and Saturn formed from GGPPs in the cool outer disk while the terrestrial planets formed in the hot inner disk by the more conventional method of collisional accumulation of small planetesimals. The two GGPPs in Boss’s (1998) model appear in less than 103 years, forming toward the end of the accretion of the Sun itself, and quite possibly prior to a quiescent stage of disk evolution during which planetesimal formation is conventionally considered to occur. If this was the case then the terrestrial planets and asteroids may have begun their growth after the formation of the gas giants and would therefore be subject to long-range perturbations throughout their entire growth. This paper is the first stage of an investigation of the consequences of an early formation of Jupiter and Saturn on the formation of planetary embryos in the terrestrial planet region and the asteroid belt. Relative velocities between bodies as a function of mass are calculated for both the terrestrial planet and the asteroid regions. We also address the question of whether growth will be prevented by destructive collisions between the planetesimals. In a subsequent paper, for conditions for which growth is possible, the evolution in mass and velocity of the bodies will be considered. 2. GAS DRAG AND THE SWIFT SYMPLECTIC INTEGRATOR

Wisdom and Holman (1991) developed a symplectic mapping technique for studying the conservative gravitational N -body problem. The basic principle underlying the symplectic mapping method involves a sequence of steps alternating between pure Keplerian motion of the individual bodies along their instantaneous elliptical orbits (a “drift”) and periodic perturbations to the bodies derived from their mutual gravitational interaction (a “kick”). The drift portion of the sequence is carried out by advancing the position and velocity vectors x, v of each body from an initial time t0 to a latter time t using Gauss’ f and g functions via the initial value problem x(t) = f (t)x(t0 ) + g(t)v(t0 ),

(1)

˙ v(t) = f˙(t)x(t0 ) + g(t)v(t 0 ).

(2)

Equations for determining f and g and their time derivatives can be found in Danby (1988). Malhotra (1994) has modified the symplectic mapping method in order to study nonconservative N -body systems. The general theory developed by Malhotra is applicable to systems in which a relatively weak nongravitational force is superimposed on the nearly Keplerian motion of the bodies. Examples include perturbations of cometary orbits due to out-gassing and the influence of radiation on small micrometer-sized dust particles. Malhotra applies the theory to the study of dissipation due to gas drag. In the nonconservative

62

KORTENKAMP AND WETHERILL

method the kick portion of the “drift–kick–drift . . .” mapping sequence remains unchanged, arising only due to the mutual gravitational interaction of the N bodies. To account for weak dissipation the drift portion of the sequence is altered by introducing modified functions f˜ and g˜ and a new function h˜ such that the positions and velocities of the bodies are advanced according to ˜ ˜ x(t) = f˜(t)x(t0 ) + g(t)v(t 0 ) + h(t)[x(t0 ) × v(t0 )], ˙˜ ˙˜ v(t) = f˙˜(t)x(t0 ) + g(t)v(t 0 ) + h(t)[x(t0 ) × v(t0 )]. Following the work of Adachi et al. (1976) we describe the gas drag force at time t0 as µ ¶ 3CD ρg F0 = − (3) ku0 ku0 , 8rp ρp where CD is a nondimensional drag coefficient, ρg is the density of the gas, rp and ρp are the radius and density of the solid body, and u0 is the velocity vector of the solid body with respect to the gas at time t0 . With a time step τ , Malhotra’s expressions for f˜, ˜ and h˜ are given by g, · ¸  τ 2 xˆ 0 · F0 − (ˆx0 · vˆ 0 )ˆv0 · F0  f˜ = f + ,    2x0 1 − (ˆx0 · vˆ 0 )2    · ¸  2 vˆ 0 · F0 − (ˆx0 · vˆ 0 )ˆx0 · F0 τ + O(τ 3 ), (4) , g˜ = g +  2v0 1 − (ˆx0 · vˆ 0 )2    ¸  2 ·  (x × v ) · F τ  0 0 0  ˜h =  , 2 |x0 × v0 |2

TABLE I Model Parameters Quantity

Symbol

Value

First appearance

Planetesimal density Nebular gas density Keplerian velocity Gas velocity Ratio (Vk − Vg )/Vk Gas drag coefficient Impact efficiency Mass escape coefficient Impact energy available Critical specific energy Crushing energy Nominal Weak

ρp ρg Vk Vg η CD k1 k3 Qf Q? Qc

3 g cm−3 1.18 × 10−9 g cm−3

Eq. (3) Eq. (3) Eq. (6)

5 × 10−3 0.5 0.5 3 × 106 (cm s−1 )2.25

Eq. (5) Eq. (3) Eq. (16) Eq. (18) Eq. (16) Eq. (21)

108 erg g−1 107 erg g−1

Eq. (17) Fig. 6

on η in the expression for de/dt is from Kary et al. (1993).) Figure 1 shows the results for planetesimals with radius rp = 3 and 5 km using the model parameters listed in Table I. The points are the numerical results and the solid lines are from Eqs. (5). The difference between the analytical and numerical results corresponds to a constant factor of about 5/3 in the value of τo . The effect of this is that Adachi’s analytical decay rates for

˙˜ and h˙˜ are where x0 = x(t0 ), v0 = v(t0 ). The time derivatives f˙˜, g, taken only to first order in τ . We have introduced these modifications into the SWIFT symplectic integrator of Levison and Duncan (1994). As a test of the new code we compared the numerical results with the predicted average orbital decay rates from Adachi et al. (1976), given by ¿ À µ ¶  τo da 9.158 2   = −2η e+ I +η ,    a dt 3π π   ¿ À µ ¶   higher order terms 2.422 2 3 τo de + =− e+ I + η , in e, I, η,  e dt π π 2   ¿ À µ ¶    1 2.422 8 τo d I   =− e+ I + η , I dt 2 π 3π (5) where η is as described in Table I and τo is a characteristic timescale given by µ ¶ 8rp ρp 1 τo = , (6) 3CD ρg Vk with Vk the Keplerian circular velocity at heliocentric distance a. The numerical values in the coefficients on e are from Adachi’s evaluation of certain elliptic integrals. (The (3/2) coefficient

FIG. 1. The decay in semi-major axis (a), eccentricity (e), and inclination (I ) due to gas drag acting on planetesimals of radius 3 and 5 km. The points show the numerical results from the modified SWIFT integrator, and the solid lines are calculated from the analytical decay rates derived by Adachi et al. (1976) (Eq. (5)). The difference between the analytical and numerical results corresponds to a constant factor of about 5/3 in the radius of the body.

PLANET AND ASTEROID FORMATION

a 5-km body match our numerical results for a 3-km body—the numerical decay rates are slower than the analytical. We found this 5/3 factor to be present over a wide range of time-step sizes and for all sizes of planetesimals. To ensure that the modified SWIFT code was not responsible for this discrepancy we performed a test using Everhart’s (1985) 15th-order RADAU integrator suitably modified to account for gas drag. The RADAU results match the SWIFT results. A third numerical technique using a HERMITE integrator (Makino and Aarseth 1992) yielded a similar discrepancy to Adachi’s equations (Inaba 1999). Inaba (1999) showed empirically that the discrepancy could be resolved by modifying Eqs. (5) as µ ¶2 ¿ À ·µ ¶ ¸ 12 2 9.158 2 τo da = −2η e + I + η2 , a dt 3π π ¿ À ·µ ¶2 µ ¶2 µ ¶2 ¸ 12 2.422 2 3 τo de =− e + I + η , e dt π π 2 ¿ À ·µ ¶ ¶2 µ ¸ 12 1 2.422 2 8 τo d I 2 =− e + I +η . dt dt 2 π 3π

(7)

Figure 2 shows the orbital decay of a 5-km radius planetesimal using all four techniques. The three numerical methods and Adachi’s analytical derivations all begin with identical formulations of the gas drag force, given by Eq. (3). Careful review of

63

Adachi’s paper revealed a footnote which states that the approximations used in their derivation of Eqs. (5) would not produce errors larger than 50%, even if e, sin I , and η are comparable, provided e, sin I < 0.3. The 5/3 discrepancy established by the numerical calculations suggests that the approximate theory may be somewhat more in error. Equations (7) are qualitatively similar to Eqs. (4.21) in Adachi et al. (1976), which they estimate to be accurate to within 10%. Our tests show that Eqs. (7) provide a better fit to the numerical calculations than do Adachi’s Eqs. (4.21). Workers who have traditionally used gas drag equations similar to Eqs. (5) (e.g., Kary et al. 1993, Wetherill and Stewart 1993, Weidenschilling et al. 1997) may wish to modify their models accordingly. The allure of a more accurate analytical treatment, derived from first principles, may entice somebody into revision of the gas drag theory. While such a project might be a worthwhile cause, it is not needed in order to pursue the present work. Discrepancies even as high as the factor of 5/3 mentioned above are rather insignificant in light of the larger uncertainties in our knowledge of the conditions that existed in the early solar nebula. An alternative to Malhotra’s (1994) method was subsequently published by Cordeiro et al. (1997). (Yet another alternative is described by Mikkola (1998).) Cordeiro incorporates the weak dissipative force into the kick portion of the algorithm in a method that may be more computationally efficient than Malhotra’s method while providing at least equally accurate results. Our preferred method has little to do with the simplicity of the available techniques and is based primarily on the fact that most of our experiments had been conducted prior to the publication of—and our becoming aware of—the alternative methods. 3. ORBITAL EVOLUTION OF PLANETESIMALS

The “Standard Model”

FIG. 2. Orbital decay of a 5-km radius planetesimal calculated using SWIFT (filled circles), RADAU (open circles), results from Inaba (1999, dashed line) (Eqs. (7)) and Adachi et al. (1976, solid line) (Eqs. (5)). All four techniques begin with identical formulations of the gas drag force, given by Eq. (3).

Conventional theories put forth to explain terrestrial planet formation surmise that small micrometer-sized grains in the early solar nebula (the disk of gas and dust surrounding the proto-sun) coagulated into larger aggregates (Weidenschilling and Cuzzi 1993). These macroscopic particles further accumulated into planetesimals in the size range of tens of meters to a few kilometers. These planetesimals are largely decoupled from the gaseous nebula and orbit the Sun with Keplerian motion, leading to a relative velocity of about 60 m s−1 with respect to the partially pressure-supported gaseous nebula (Whipple 1972, Adachi et al. 1976, Weidenschilling 1977). The resulting gas drag dissipates energy and angular momentum—circularizing the planetesimal orbits and reducing their mutual inclinations while slowly causing them to decay toward the Sun. By this process the planetesimals evolve on very nearly circular coplanar orbits perturbed only by their mutual gravitational interactions. Typical orbital elements of such a population are shown in Table II. Relative encounter velocities for bodies on such orbits are ∼10 cm s−1 . In the absence of other strong gravitational forces it has been shown that after ∼105 years at 1 AU collisional

64

KORTENKAMP AND WETHERILL

TABLE II Planetesimal Initial Orbital Elements Element

Symbol

Initial value

Semi-major axis Eccentricity Inclination Arg. pericenter Long. ascending node

a e I ω Ä

Variable, see Tables III & IV 10−5 sin I = e 0◦ to 360◦ (random) 0◦ to 360◦ (random)

and orbital evolution results in the runaway growth of the largest bodies (Wetherill and Stewart 1993). These runaways become Mercury- to Mars-sized planetary embryos, separated by about ∼0.05 AU (Kokubo and Ida 1996, Weidenschilling et al. 1997). This process results in a few dozen planetary embryos distributed between 0.5 and 3 AU. Models pertaining to the final stages of terrestrial planet formation typically permit these planetesimals and embryos to evolve for 5 to 10 million years before introducing Jupiter and Saturn and dissipating the gaseous nebula (Wetherill 1996, Chambers and Wetherill 1998). The resulting simulated planetary systems have a number of features in common with the terrestrial planets of our Solar System, but also demonstrate unsolved problems. Consideration of GGPPs The distribution of encounter velocities is one of the key factors which controls the rate of planetesimal growth. The distribution of orbital elements of bodies on intersecting orbits controls the encounter velocities of these bodies. If Jupiter and Saturn formed via the GGPP scenario then they are likely to have been present during the formation of planetesimals and thus influenced terrestrial planet and asteroid formation from the “very beginning.” The gravitational influence of giant planets is important not just for orbits near the chaotic regions of resonance. Smooth secular perturbations act on orbits at all semi-major axes. In the absence of gas drag, secular gravitational perturbations from Jupiter and Saturn are dependent only on time and the semi-major axis of the orbit being perturbed (to low order in e and I ). This means that there is virtually no difference, from an orbital dynamics point of view, between the behavior of a kilometer-sized planetesimal and a full-fledged planet. Bodies of all sizes on orbits distributed over a small range in semimajor axis (1a ' 0.1 AU near 1 AU) will experience nearly identical forced variations in their elements e, I , ω, and Ä—that is, the orbits will evolve in phase. Again, this is in the absence of gas drag. Such is not the case if gas drag is acting. In this case secular gravitational perturbations from Jupiter and Saturn are dependent on time, the semi-major axes of the planetesimal orbits, and the rate of decay of the orbits. These perturbations act to increase the orbital eccentricities and inclinations of planetesimals, counteracting the effects of gas drag. For a swarm of planetesimals of a single size and initially on nearly circular

low inclination orbits distributed over a small range of semimajor axis these perturbations act identically such that the individual orbits evolve in phase; that is, they evolve as coplanar concentric elliptical orbits. Encounter velocities between these single-size bodies may remain low enough to encourage collisional growth of the bodies. However, the orbital decay rates associated with gas drag (Eqs. (5)) are dependent on the physical properties of the planetesimals (density and radius). For a given planetesimal density, a, e, and I of smaller bodies decay faster than those of larger bodies. Therefore, when a range of planetesimal sizes is studied, we find that as the orbits of smaller bodies decay past the orbits of larger bodies the orbits of the two sets are not in phase. The highly exaggerated orbital schematics of Fig. 3 provide examples of how the phasing of ω and Ä can influence the encounter velocities. The size-dependent amplitudes of e and I contribute to the encounter velocities as well. Model Parameters and Initial Conditions We include, from time t = 0, the planets Jupiter and Saturn at their present masses and in their present orbits. The model parameters are given in Table I and the initial orbital elements e, I , ω, and Ä are shown in Table II. All inclinations are referred to the invariable plane of the Jupiter–Saturn system. A number of planetesimals of various sizes were distributed over a limited range 1a = amax − amin near 1 and near 2.6 AU (see Tables III and IV). Planetesimals are treated as test particles (no mutual perturbations) and their designated masses are only used to determine the magnitude of the gas drag force. Encounter velocities of bodies on these initial orbits are in the range 20 to 40 cm s−1 . All simulations were run for a period of 25,000 years. It is not our intention in this paper to study close encounters and collisions between individual bodies, but rather, to trace the smooth secular changes in the distribution of orbital elements as a function of planetesimal size and semi-major axis. In this regard the relatively small number of bodies (25 or 50 of each size) proved more than adequate. TABLE III Initial Planetesimal Distribution near 1 AU Radius (km)

Mass (g)

Number

amin (AU)

amax (AU)

1000.0 500.0 200.0 100.0 50.0 20.0 10.0 5.0 2.0 1.0 0.5 0.2 0.1

∼1025 ∼1024 ∼1023 ∼1022 ∼1021 ∼1020 ∼1019 ∼1018 ∼1017 ∼1016 ∼1015 ∼1014 ∼1013

25 25 25 25 25 25 50 50 50 50 50 50 50

1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

1.02 1.03 1.04 1.05 1.10 1.15 1.15 1.20 1.25 1.30 1.50 2.40 2.40

PLANET AND ASTEROID FORMATION

65

FIG. 3. The top two diagrams show the orientation of two orbits with identical values of a and e but with arguments of pericenter, ω, that are nearly in phase (right) and 90◦ out of phase (left). The bottom two diagrams show the orientation of two orbits with identical inclinations but with longitudes of ascending node, Ä, that are in phase (right) and 180◦ out of phase (left). If e and I are relatively large then the phasing of ω and Ä can lead to either high or low encounter velocities, Venc .

Results of the Simulations It would not be useful to illustrate graphically the orbital evolution of each of our 13 different size classes in the two different semi-major axis regions. Therefore we have chosen to display, as representative examples, the orbital evolution of the 1- and 2-km bodies. The evolution of these bodies near 1 AU is shown in Fig. 4. For clarity only a small subset of the 100 modeled orbits are shown in the top panel. The solid lines are the 2-km bodies and the more quickly decaying dashed lines are the 1-km bodies. Of the 50 orbits in each set, only those orbits that are intersecting orbits in the other set are used to determine the mean values shown in Fig. 4. After 3000 years the eccentricities have TABLE IV Initial Planetesimal Distribution near 2.6 AU Radius (km)

Mass (g)

Number

amin (AU)

amax (AU)

1000.0 500.0 200.0 100.0 50.0 20.0 10.0 5.0 2.0 1.0 0.5 0.2 0.1

∼1025 ∼1024 ∼1023 ∼1022 ∼1021 ∼1020 ∼1019 ∼1018 ∼1017 ∼1016 ∼1015 ∼1014 ∼1013

25 25 25 25 25 25 50 50 50 50 50 50 50

2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6 2.6

2.62 2.63 2.64 2.65 2.70 2.75 2.75 2.80 3.00 3.00 3.50 4.00 4.00

risen from the initial value of 10−5 to ∼10−3 and the inclinations have risen by a similar factor. The ω and Ä panels of Fig. 4 show that the orbits of each set slowly drift out of phase after their initial orientation to the preferred ω and Ä values forced by the secular perturbations. The bottom panel of Fig. 4 shows the relative encounter velocities of the bodies on intersecting orbits. The velocities were calculated using the procedures outlined in the Appendix. The larger body is taken to be the target and the smaller body is the projectile. After 5000 years the relative encounter velocities between these 1- and 2-km bodies grows by about two orders of magnitude, to ∼10 m s−1 . Orbital evolution in the asteroid belt, between 2.55 and 2.75 AU, is shown in Fig. 5. Here, the stronger perturbations increase the eccentricities by about an order of magnitude higher than near 1 AU. In addition, orbital decay through the numerous mean-motion resonances tends to introduce more chaotic behavior. Inclinations in the asteroid belt also reach values higher than near 1 AU. Owing largely to the sharp increase in eccentricities the relative encounter velocities of these 1- and 2-km bodies in the asteroid belt grow by about three orders of magnitude, to ∼100 m s−1 . The Effects of Impacts For sufficiently energetic encounters, a collision will reduce the mass of the target. We have studied the threshold at which this occurs for bodies of sizes ranging from 0.1 to 1000 km (1013 to 1025 g). Collisional loss of mass was not found for bodies larger than 50 km in radius (1021 g). Figure 6, representing the region near 1 AU, shows the relative encounter velocities

66

KORTENKAMP AND WETHERILL

FIG. 4. The top five panels show the orbital evolution of a, e¯ , I¯ , ω, ¯ and ¯ for two different sizes of planetesimals with orbits near 1 AU. For clarity Ä only a subset of the 100 modeled orbits are shown in the top panel. The dashed lines correspond to orbits of 1-km radius bodies and the solid lines correspond to orbits of 2-km radius bodies. The mean elements are calculated only for those 1- and 2-km orbits which are intersecting and in the range 0.8 < a < 1.2. ¯ panels also show the relative phasing of the intersecting orbits, The ω¯ and Ä ¯ = |Ä ¯ 1 −Ä ¯ 2 |. The bottom panel shows the mean with 1ω¯ = |ω¯ 1 − ω¯ 2 | and 1Ä encounter velocity (±1σ ) between bodies on these intersecting orbits.

FIG. 5. Same as described in the legend to Fig. 4 except for orbits with 2.55 < a < 2.75. The stronger secular perturbations acting in the asteroid belt lead to higher e and I values and to faster precession rates of ω and Ä. Chaotic perturbations acting on the bodies as their orbits decay through the numerous jovian mean-motion resonances primarily influence the evolution of e¯ and ω. ¯

between bodies (targets) of six selected sizes, all at or below this 50-km threshold, with small bodies (projectiles) ranging from 0.1 km up to the radius of the target. The two horizontal solid lines represent the growth limits for crushing strengths of 108 ergs g−1 (nominal—upper line) and 107 ergs g−1 (weak— lower line). The horizontal dotted line is the growth limit as

PLANET AND ASTEROID FORMATION

67

FIG. 6. This figure shows encounter velocities near 1 AU for the model with the semi-major axes of Jupiter and Saturn at aJ = 5.2 AU and aS = 9.5 AU, respectively. The six panels show the cases for target planetesimals which experience impacts that straddle the growth limits imposed by the cratering criteria given in Appendix B. Shown are target planetesimals with radii of Rt = 50, 20, 10, 5, 2, 1 km and masses of Mt ≈ 1021 , 1020 , 1019 , 1018 , 1017 , 1016 g (assuming ρp = 3 g cm−3 ). The data points (±1σ ) indicate the mean encounter velocity between the target and projectiles of various sizes, averaged over a time span of 25,000 years. The top left panel labels the nominal growth limit (for Q c = 108 erg g−1 ), the growth limit for bodies composed of weaker material with a crushing strength of Q c = 107 erg g−1 , the growth limit given by Housen et al. (1983), and the surface escape velocity for the combined target–projectile. The top right panel labels the gravity regime disruption limit given by Love and Ahrens (1996).

described by Housen et al. (1983). The bold dashed diagonal line above these erosion limits represents the catastrophic disruption limit of Love and Ahrens (1996), where >50% of the target mass is lost due to the impact. The derivations and criteria used to establish these erosion and disruption limits are given in the Appendix. The horizontal dot–dash line marks the surface escape velocity of the combined target–projectile. Encounter velocities above this escape velocity preclude runaway growth of the target body.

In all cases where the relative velocity for encounters with bodies of similar mass is small a larger body can be formed. For smaller projectiles, the encounter velocities will be higher and may be above one or all of the thresholds for mass loss rather than growth. If the evolving mass distribution is such that a large fraction of the mass remains in bodies of very similar mass, then a planetesimal assemblage in which the mass of the largest bodies continues to grow may be formed, even though such a large fraction of the mass of the system may migrate

68

KORTENKAMP AND WETHERILL

FIG. 7. Same as described in the legend to Fig. 6, except for orbits near 2.6 AU.

sunward by gas drag that the planetesimal growth process may be a very inefficient one. The results of a calculation analogous to that of Fig. 6 are shown in Fig. 7, corresponding to an asteroidal semi-major axis of about 2.6 AU. This heliocentric distance was chosen to be near the region of the asteroid belt that supplies Earth-approaching asteroids and meteorites, but not within Jupiter’s 3:1 commensurability resonance, which would provide very high relative velocities that would be of no relevance to the problem under discussion. For the same reason we avoid the secular resonances. In the case shown in Fig. 7 the relative velocities are much higher than those found near 1 AU, and the inhibitions to planetary growth even more severe. The situation could be different if Jupiter and Saturn formed at somewhat greater heliocentric distances and then slowly mi-

grated inward. Results of analogous calculations in which the semi-major axes of Jupiter and Saturn were 1 AU larger are shown in Figs. 8 and 9. At least at 1 AU, the amplitudes of the perturbations become more modest, and therefore the chances of overcoming the difficulties described here would be greater. Because the original positions and migration rates of the giant planets are unknown the true locations of the mean-motion resonances cannot be determined. Likewise, the positions of secular resonances depend on unknown parameters of the gas disk such as mass, structure, and the dispersal rate. Rather than make assumptions regarding all of these factors we prefer to avoid the resonances. So while the effects of all mean-motion and secular resonances are implicitly included in our numerical integrations, we consider the size-dependent behavior described in this paper a general characteristic of resonance-free regions.

69

PLANET AND ASTEROID FORMATION

FIG. 8. Same as described in the legend to Fig. 6, except for aJ = 6.2 AU and aS = 10.5 AU.

Under the conditions we have modeled it seems very unlikely that runaway growth could occur before the escape velocities of the larger bodies of the assemblage exceed the 102 to 103 m s−1 velocities imposed by the long-range secular perturbations of Jupiter and Saturn. A thorough calculation of the actual evolution of a system of this kind is difficult, but in view of the very high relative velocities imposed by the giant planets, the situation may be sufficiently extreme, at least in the asteroid region, that simple approximations to the real problem are warranted. At present there is no statistical analytical theory of planetesimal interactions for the case in which the angular elements (ω and Ä) have not been randomized at the outset, as is assumed by Stewart and Wetherill (1988) and Stewart and Ida (1999). We are addressing this problem by parameterization of our numerical in-

tegration calculations, but cannot say anything conclusive about the results at this time. Very rough calculations using a modified form of the gas dynamic model of Wetherill and Stewart (1993) indicate that runaway growth at 1 AU would be significantly delayed, by ∼106 years, but not necessarily prevented. However, the conditions in the asteroid belt may not allow runaway growth to occur at all. 4. CONCLUSIONS

The relative encounter velocities of planetesimals at 1 AU in the presence of early formed Jupiter and Saturn and the gas of the protoplanetary disk are considerably lower than they would be if it were not for the direct effects of gas drag and the tendency of

70

KORTENKAMP AND WETHERILL

FIG. 9. Same as described in the legend to Fig. 6, except for orbits near 2.6 AU and with aJ = 6.2 AU and aS = 10.5 AU.

the perturbations of the orbital elements to remain in phase when the semi-major axes are similar. At the same time, the presence of gas reduces greatly the tendency for the perturbations to remain in phase, except for bodies very nearly equal in mass. For bodies of other masses, relative velocities of ∼100 m s−1 are found for bodies near 1 AU and several times larger near 2.6 AU. These high relative encounter velocities certainly are a serious handicap to growth of planetary embryos near 1 AU and could eliminate this possibility entirely in the asteroid belt. The Jupiter–Saturn perturbations become considerably smaller if Jupiter and Saturn had larger semi-major axes at the time of planet formation, and drifted inward to their present positions as a result of radial migrations commonly believed to be caused by gravitational interactions between the massive planets and the even more massive gas disk (Goldreich and Tremaine 1980,

Ward and Hourigan 1989). Trilling et al. (1998) describe one calculation which shows that if Jupiter formed at 5.2 AU it would drift into 0.03 AU in 107 years. Clearly this result is dependent on the assumed properties and lifetime of the solar nebula. However, Trilling’s work indicates that some radial migration of Jupiter should be expected. The migration of Jupiter has not yet been calculated for the parameters in Boss’ model, so appropriate initial distances for Jupiter and Saturn are not known at this time. It also remains an open question as to what degree of migration is compatible with the formation of terrestrial planets, and at the same time the depletion of the asteroid belt. As can be seen from Figs. 8 and 9, not much migration is required to reduce the effects of Jupiter and Saturn perturbations to modest levels in the terrestrial planet region while they remain quite high in the asteroidal region.

71

PLANET AND ASTEROID FORMATION

All the uncertainty that is inevitable in such difficult theoretical problems, the even greater uncertainty in our understanding of conditions in the presolar disk, the recent contributions by Boss, all those working on giant planet migration, and the observers of extra-solar giant planets can provide fuel for rampant speculation. We will suggest just one scenario. If giant planets form early as GGPPs then the stunted planetary system in the asteroid belt and the small mass of Mars—in contrast to the fine array of planets in the terrestrial planet region—may be a consequence of Jupiter forming slightly farther out in the solar nebula and slowly migrating to its rather accidental final position and thereby truncating the growth of planetary embryos not far beyond 1 AU. Further investigation of this problem will require understanding the evolution of the system at the stage (∼1023 g) where the mutual perturbations between the planetesimals become strong enough to significantly randomize the phases of their orbits. This is a much more difficult problem that we are attempting to address by parameterization of the numerical integration calculations. APPENDIX

A. Encounter Velocities Recall that because of the size-dependent nature of the gas drag force [Eq. (3)] the projectile orbit is decaying more quickly than the target orbit. As a consequence of this the projectile orbit will exactly intersect the target orbit at two instants in time. However, it is extremely unlikely that these intersection times will coincide with the time of output from the integrator. Therefore, given a target orbit with elements at , et , It , ωt , Ät and a target-crossing projectile orbit with ep , Ip , ωp , Äp which are output at the time step immediately after intersection, we wish to solve for the value of ap which results in the intersection of the two orbits and thus allows for collision of the bodies. From spherical trigonometry and following some of the notation of Greenberg (1982) the mutual inclination Im between the target and projectile orbits is found to be cos Im = cos Ip cos It + sin Ip sin It cos 1Ä,

(8)

where 1Ä ≡ Äp − Ät . The arguments of the mutual nodes (u t and u p ), measured from Ät and Äp respectively, are given by

Positioning the target body at the primary mutual node λt = u t − ωt gives it a heliocentric distance ¡ ¢ at 1 − et2 1 + et cos λt

(10)

and similarly, with the projectile body at λp = u p − ωp ,

d=

ap s1 − ep2 d 1 + ep cos λp

.

1 − et2 1 + ep cos λp · 1 − ep2 1 + et cos λt

# .

(12)

Following Wetherill (1967) the relative encounter velocity of the target and projectile at d is then given by

2 Venc =

" (2Ap − 1) (2At − 1) GM ap + at d2 A2p A2t # − 2f

s

ap at 1 − ep2

ds

1 − et2

dg

1/2

· (cos Im + cot αp cot αt ) ,

(13)

where G is the universal gravitational constant, M is the central mass, Ap = ap /d, At = at /d, and cot2 αp = cot2 αt =

A2p ep2 − (Ap − 1)2 A2p s1 − ep2 d A2t et2 − (At − 1)2 ¡ ¢ . A2t 1 − et2

(14)

The secondary mutual node, corresponding to u t + 180◦ , requires a slightly different value of ap and therefore results in a slightly different relative encounter velocity. Given a target and projectile of mass and radius Mt , Rt and Mp , Rp respectively the surface escape velocity of the pair is 2 Vesc =

2G(Mt + Mp ) . (Rt + Rp )

(15)

2 = Conservation of energy then gives an impact velocity of Vimp 2 2 Venc + Vesc .

B. Erosion and Disruption Conditions We follow the erosion conditions used by Wetherill and Stewart (1993). The center of mass impact energy available for fragmentation is given by Qf =

k1 2 V Mt Mp /(Mt + Mp ), 2 imp

(9)

sin u p = sin u t sin It / sin Ip .

d=

" a p = at

(16)

where the impact efficiency k1 = 0.5 is the fraction of the impact energy not lost to heating. The mass of material fragmented by the impact is

u t = arctan[−sin 1Ä/(cotIp sin It − cos 1Ä cos It )], cos u p = cos u t cos 1Ä + sin u t sin 1Ä cos It ,

The exact value of ap which results in the intersection of the target and projectile orbits can then be determined by equating these two expressions for d

Mf = Q f /Q c ,

where Q c = 108 erg g−1 is the assumed crushing strength. The mass of material assumed to escape is given by −2.25 Me = k3 Mf Vesc

(18)

(Gault et al. 1963, Greenberg et al. 1978), where k3 = 3 × 106 (cm s−1 )2.25 . For a nominal erosion limit we use Q c = 108 erg g−1 and for a weak erosion limit we use Q c = 107 erg g−1 . Considering that this may be a rather dated description for the mass of escaping material we also use the condition given by Housen et al. (1983). This relation is µ

(11)

(17)

Me = πρp

D 2

¶1.7 ¶2 µ √ g(D/2) α , Vesc

(19)

72

KORTENKAMP AND WETHERILL

where D is the diameter of the crater, given by Schmidt and Housen (1987) as à D = 2.32

2 Vimp

!0.22

2g Rp

Rp ,

(20)

g is the surface acceleration of gravity for the target, and α = 0.2 is the assumed crater depth to diameter ratio. If Me ≥ Mp there is no growth of the target. There are numerous models in use that describe how the critical specific energy (Q ? , the energy per unit target mass required for catastrophic disruption of the target) scales with the size of the target. In a review of seven recent scaling laws Durda et al. (1998) showed that the transition from the strength-dominated to gravity-dominated regimes falls between target radii of tens of meters to tens of kilometers. Clearly such a large uncertainty in the critical specific energy weakens arguments against any particular scaling law, no matter how dated the work. Nevertheless, the new models should not be ignored. For our gravity regime disruption criteria we choose to use the scaling law of Love and Ahrens (1996) primarily because it compares favorably with models used in the past by Wetherill and Stewart (1993). Love and Ahrens’ scaling law (in erg g−1 ) is given by Q ? = 24.2[Rt (cm)]1.13 .

(21)

For an impact with Q f > Q ? more than 50% of the target mass will be disrupted and dispersed (catastrophic disruption). The transition between this Q ? and most strength regime critical specific energies falls near or below the smallest target size used in our studies.

ACKNOWLEDGMENTS We thank Alan Boss and David Trilling for very useful discussions and Jack Lissauer and an anonymous reviewer for their helpful comments. This research is funded by grants from the NASA programs on Origins, Exobiology, and Planetary Geology and Geophysics.

REFERENCES Adachi, I., C. Hayashi, and K. Nakazawa 1976. The gas drag effect on the elliptic motion of a solid body in the primordial solar nebula. Prog. Theoret. Phys. 56, 1756–1771. Boss, A. P. 1996. Evolution of the solar nebula III: Protoplanetary disks undergoing mass accretion. Astrophys. J. 469, 906–920. Boss, A. P. 1997. Giant planet formation by gravitational instability. Science 276, 1836–1839. Boss, A. P. 1998. Evolution of the solar nebula IV: Giant gaseous protoplanet formation. Astrophys. J. 503, 923–937. Boss, A. P. 1999. Formation of extrasolar giant planets: Core accretion or disk instability? Earth, Moon, Planets, in press. Butler, R. P., and G. W. Marcy 1996. A planet orbiting 47 Ursae Majoris. Astrophys. J. 464, L153–156. Butler, R. P., G. W. Marcy, E. Williams, H. Hauser, and P. Shirts 1997. Three new “51 Pegasi-Type” planets. Astrophys. J. 474, L115–118. Cameron, A. G. W. 1978. Physics of the primitive solar accretion disk. Moon Planets 18, 5–40. Chabrier G., D. Saumon, W. B. Hubbard, and J. I. Lunine 1992. The molecular– metallic transition of hydrogen and the structure of Jupiter and Saturn. Astrophys. J. 391, 817–826. Chambers, J. E., and G. W. Wetherill 1998. Making the terrestrial planets: Nbody integrations of planetary embryos in three dimensions. Icarus 136, 304– 327.

Cordeiro, R. R., R. S. Gomes, and R. V. Martins 1997. A mapping for nonconservative systems. Celes. Mech. Dynamic. Astron. 65, 407–419. Danby, J. M. A. 1988. In Fundamentals of Celestial Mechanics, pp. 162–164. Willmann-Bell, Richmond. Durda, D. D., R. Greenberg, and R. Jedicke 1998. Collisional models and scaling lows: A new interpretation of the shape of the main-belt asteroid size distribution. Icarus 135, 431–440. Everhart, E. 1985. An efficient integrator that uses Gauss-Radau spacings. In Dynamics of Comets: Their Origin and Evolution (A. Carusi and G. B. Valsecchi, Eds.), pp. 185–202. D. Reidel, Dordrecht. Gault, D. E., E. M. Shoemaker, and H. J. Moore 1963. Spray ejected from the lunar surface by meteroid impact. NASA TN 1767. Goldreich, P., and S. Tremaine 1980. Disk-satellite interactions. Astrophys. J. 241, 425–441. Greenberg, R. 1982. Orbital interactions: A new geometrical formalism. Astron. J. 87, 184–195. Greenberg, R., J. Wacker, C. R. Chapman, and W. K. Hartmann 1978. Planetesimals to planets: Numerical simulation of collisional evolution. Icarus 35, 1–26. Guillot, T., G. Chabrier, P. Morel, and D. Gautier 1994. Nonadiabatic models of Jupiter and Saturn. Icarus 112, 354–367. Housen, K. R., R. M. Schmidt, and K. A. Holsapple 1983. Crater ejecta scaling laws: Fundamental forms based on dimensional analysis. J. Geophys. Res. 88, 2485–2499. Inaba, S. 1999. Modeling of the Formation of Protoplanets without Nebular Gas by Solving the Statistical Accumulation Equation. Ph.D. dissertation, Tokyo Institute of Technology, Tokyo, Japan. Kary, D. M., J. J. Lissauer, and Y. Greenzweig 1993. Nebular gas drag and planetary accretion. Icarus 106, 288–307. Kokubo, E., and S. Ida 1996. On runaway growth of planetesimals. Icarus 123, 180–191. Kuiper, G. P. 1951. On the origin of the Solar System. Proc. Natl. Acad. Sci. 37, 1–14. Levison, H. F., and M. J. Duncan 1994. The long-term dynamical behavior of short-period comets. Icarus 108, 18–36. Lissauer, J. J. 1987. Timescales for planetary accretion and the structure of the protoplanetary disk. Icarus 69, 249–265. Lissauer, J. J. 1993. Planet formation. Annu. Rev. Astron. Astrophys. 31, 129– 174. Lissauer, J. J., J. B. Pollack, G. W. Wetherill, and D. J. Stevenson 1995. Formation of the Neptune system. In Neptune and Triton (D. P. Cruikshank, Ed.), pp. 37– 108. Univ. of Arizona Press, Tucson. Love, S. G., and T. J. Ahrens 1996. Catastrophic impacts on gravity dominated asteroids. Icarus 124, 141–155. Makino, J., and S. J. Aarseth 1992. On a Hermite integrator with Ahmad–Cohen scheme for gravitational many-body problems. Publ. Astron. Soc. Jpn. 44, 141–151. Malhotra, R. 1994. A mapping method for the gravitational few-body problem with dissipation. Celes. Mech. Dynamic. Astron. 60, 373–385. Marcy, G. W., and R. P. Butler 1996. A planetary companion to 70 Virginis. Astrophys. J. 464, L147–151. Mayor, M., and D. Queloz 1995. A Jupiter-mass companion to a solar-type star. Nature 378, 355–359. Mikkola, S. 1998. Non-canonical perturbations in symplectic integration. Celes. Mech. Dynamic. Astron. 68, 249–255. Pollack, J. B., O. Hubickyj, and Y. Greenzweig 1996. Formation of the giant planets by concurrent accretion of solids and gas. Icarus 124, 62–85. Schmidt, R. M., and K. R. Housen 1987. Some recent advances in the scaling of impact and explosion cratering. Intl. J. Impact Eng. 5, 543–560.

PLANET AND ASTEROID FORMATION Stevenson, D. J. 1982. Formation of the giant planets. Planet. Space Sci. 30, 755–764. Stewart, G. R., and S. Ida 2000. Velocity evolution of planetesimals: Unified analytical formulas and comparisons with N-body simulations. Icarus 143, 28–44. Stewart, G. R., and G. W. Wetherill 1988. Evolution of planetesimal velocities. Icarus 74, 542–553. Strom, S. E., S. Edwards, and M. F. Skrutskie 1993. Evolutionary time scales for circumstellar disks associated with intermediate- and solar-type stars. In Protostars and Planets III (E. H. Levy and J. I. Lunine, Eds.), pp. 837–866. Univ. of Arizona Press, Tucson. Trilling, D. E., W. Benz, T. Guillot, J. I. Lunine, W. B. Hubbard, and A. Burrows 1998. Orbital evolution and migration of giant planets: Modeling extrasolar planets. Astrophys. J. 500, 428–439. Ward, W. R., and K. Hourigan 1989. Orbital migration of protoplanets: The inertial limit. Astrophys. J. 347, 490–495. Weidenschilling, S. J. 1977. Aerodynamics of solid bodies in the solar nebula. Mon. Not. R. Astron. Soc. 180, 57–70. Weidenschilling, S. J. 1997. Growing Jupiter’s core by runaway accretion. Proc. Lunar Planet. Sci. Conf. 28th, 1513–1514.

73

Weidenschilling, S. J., and J. N. Cuzzi 1993. Formation of planetesimals in the solar nebula. In Protostars and Planets III (E. H. Levy and J. I. Lunine, Eds.), pp. 1031–1060. Univ. of Arizona Press, Tucson. Weidenschilling, S. J., D. Spaute, and K. Ohtsuki 1997. Accretional evolution of a planetesimal swarm II. The terrestrial zone. Icarus 128, 429–455. Wetherill, G. W. 1967. Collisions in the asteroid belt. J. Geophys. Res. 72, 2429– 2444. Wetherill, G. W. 1980. Formation of the terrestrial planets. Ann. Rev. Astron. Astrophys. 18, 77–113. Wetherill, G. W. 1992. An alternative model for the formation of the asteroids. Icarus 100, 307–325. Wetherill, G. W. 1996. The formation and habitability of extra-solar planets. Icarus 119, 219–238. Wetherill, G. W., and G. R. Stewart 1993. Formation of planetary embryos: Effects of fragmentation, low relative velocity, and independent variation of eccentricity and inclination. Icarus 106, 190–209. Whipple, F. L. 1972. On certain aerodynamic processes for asteroids and comets. In From Plasma to Planet (A. Elvius, Ed.), pp. 211–232. Wiley, New York. Wisdom, J., and M. Holman 1991. Symplectic maps for the N-body problem. Astron. J. 102, 1528–1538.