Numerical simulation of the final stages of terrestrial planet formation

ICARUS 4 4 , 7 0 6 - 7 2 1

(1980)

Numerical Simulation of the Final Stages of Terrestrial Planet Formation L A R R Y P. C O X AND J O H N S. L E W I S Department of Earth and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139

Received October 17, 1979; revised May 5, 1980 Three representative numerical simulations of the growth of the terrestrial planets by accretion of large protoplanets are presented. The mass and relative-velocity distributions of the bodies in these simulations are free to evolve simultaneously in response to close gravitational encounters and occasional collisions between bodies. The collisions between bodies, therefore, arise in a natural way and the assumption of expressions for the relative velocity distribution and the gravitational collision cross section is unnecessary. These simulations indicate that the growth of bodies with final masses approaching those of Venus and the Earth is possible, at least for the case of a two-dimensional system. Simulations assuming an initial uniform distribution of orbital eccentricities on the interval from 0 to emax are found to produce final states containing too many bodies with masses which are too small when emax< 0.10, while simulations with emax> 0.20 result in too many catastrophic collisions between bodies thus preventing rapid accretion of planetarysize bodies. The emax 0.15 simulation ends with a state surprisingly similar to that of the present terrestrial planets and, therefore, provides a rough estimate of the range of radial sampling to be expected for the terrestrial planets. =

I. INTRODUCTION In v i e w o f t h e effort b e i n g e x p e n d e d in t h e s t u d y o f p l a n e t e s i m a l a c c r e t i o n in t h e early stages of terrestrial planet formation ( G o l d r e i c h a n d W a r d , 1973; S a f r o n o v , 1969; W e i d e n s c h i l l i n g , 1974; G r e e n b e r g et al., 1978), it s e e m e d a p p r o p r i a t e to t r y to determine whether terrestrial-like planetary s y s t e m s c o u l d be p r o d u c e d as a r e s u l t o f the gravitationally induced collisions and a c c r e t i o n s o c c u r r i n g in a s w a r m o f p r o t o p l a n e t s a s s u m e d to h a v e a l r e a d y r e a c h e d a p p r o x i m a t e l y l u n a r m a s s . In a d d i t i o n , it s e e m e d d e s i r a b l e to e x p l o r e a v a r i e t y o f initial c o n d i t i o n s for t h e p r o t o p l a n e t a r y s w a r m in o r d e r to d e t e r m i n e w h i c h o f t h e s e c o n d i t i o n s r e s u l t in t e r r e s t r i a l - l i k e s y s t e m s o f final b o d i e s . A comparison of predictions of the commonly employed two-body gravitational model with integrations of the equations of m o t i o n w h i c h i n c l u d e t h e i n f l u e n c e o f the S u n i n d i c a t e t h a t the t w o - b o d y m o d e l is i n a d e q u a t e for t h e t r e a t m e n t o f e n c o u n t e r s 0019-1035/80/120706-16502.00/0 Copyright © 1980 by Academic Press, Inc. All rights of reproduction in any form reserved.

b e t w e e n large p r o t o p l a n e t s , p a r t i c u l a r l y when they move along low eccentricity o r b i t s r e s u l t i n g in e n c o u n t e r s w i t h l o w relative v e l o c i t i e s . I n t e g r a t i o n s p e r f o r m e d b y D o l e (1962) a n d Giuli (1968) in a t t e m p t i n g to a c c o u n t for t h e r o t a t i o n o f t h e E a r t h b y the gravitational accretion of small bodies h a d e a r l i e r r e v e a l e d nine d i s t i n c t b a n d s o f h e l i o c e n t r i c o r b i t s r e s u l t i n g in i m p a c t i n g t r a j e c t o r i e s for t h e c a s e o f b o d i e s a p p r o a c h ing a l o n g c i r c u l a r o r b i t s w h e r e o n e o f t h e b o d i e s w a s a s s u m e d to b e m a s s l e s s . Because of the inadequacy of previous twob o d y g r a v i t a t i o n a l m o d e l s in t r e a t i n g t h e c l o s e e n c o u n t e r s a r i s i n g in t h e p r e s e n t simu l a t i o n s , a n e w m o d e l for single c l o s e encounters between two small masses, ml a n d m2, o r b i t i n g a m u c h l a r g e r m a s s , M, h a s b e e n p r o p o s e d ( C o x et al., 1978, h e r e a f t e r r e f e r r e d to as P a p e r I). T h e p r e d i c t i o n s o f t h e m o d e l w e r e c o m p a r e d w i t h integrations of the three-body equations of m o t i o n in o r d e r to d e t e r m i n e for w h i c h e n c o u n t e r s t h e m o d e l is a n a d e q u a t e approximation. 706

TERRESTRIAL PLANET FORMATION II. THE BASIC ASSUMPTIONS OF THE SIMULATIONS

707

the sum o f the masses o f the terrestrial planets, with uniform spacing of the protoplanets in the ecliptic plane between the A. The Treatment o f Encounters heliocentric distances 0.5 and 1.5 AU, the In these simulations the protoplanetary typical spacing is - 1 7 0 R~. The separation bodies are a s s u m e d to m o v e along unpero f 50 R~ was chosen as a c o m p r o m i s e turbed coplanar orbits until they pass b e t w e e n the a b o v e two requirements. It within a sphere of influence distance, Rs, of should be noted that only the effects of another body. This radius is given by R~ those encounters e x p e c t e d on the basis of Rl[ml(ml + m~)/MZ] 1f~, (Battin, 1964), unperturbed orbits to result in separations where ml is taken to be the larger mass less than the sphere of influence are inwhen ml and m2 are unequal and where it is cluded in this simulation. a s s u m e d that m~,m~ ~ M. The radius R~ One can ask what fraction of the changes denotes the distance from M to m~. I f a in the orbital elements of the protoplanets close a p p r o a c h to another b o d y occurs, might be e x p e c t e d to result from close new orbital elements are c o m p u t e d for each encounters (Ar)min < Rs, as o p p o s e d to body. When the p a r a m e t e r s of the encounmore distant encounters, (Ar)min > Rs, for a ter fall within the region o f validity of the two-dimensional s y s t e m of protoplanets model described in Paper I for close ensuch as the one assumed for the initial state counters, the new orbital elements are com- of our simulations. A p p r o x i m a t e estimates puted using the model. T h o s e encounters of the relaxation times due to close and having p a r a m e t e r s outside the range o f va- more distant encounters were made (Cox, lidity o f the model and those predicted by 1978) where the relaxation time for close the model to have separations at closest encounters was taken as the time for a a p p r o a c h less than the sum of the radii of typical body to experience an encounter the two protoplanets are treated b y integra- with a deflection angle, 8 => rr/2, and where tion of the regularized equations o f motion the relaxation time for more distant en(Bettis and Szebehely, 1971). The motion of counters was taken as the time within the bodies during these latter encounters is which the root mean square velocity integrated starting at an initial separation o f change due to encounters was the same 50 R~, and continuing through close ap- order of magnitude as the typical relative proach to a final separation of 50 Rs or to velocity between bodies. impact, where new orbital elements are The encounters between bodies in these c o m p u t e d for each body. estimates were treated as t w o - b o d y interacT w o considerations influence the choice tions neglecting the presence of the Sun. o f 50 Rs as the separation distance at which Including the influence of the Sun on the integrations are initiated and terminated. encounters by means of the model deClearly, the separation must be large scribed in Paper I would have made the enough to include in the integration as estimation of relaxation times impractically m u c h as possible of the effect of the en- difficult, where an order-of-magnitude esticounter on the orbital elements o f the two mate is sufficient. These estimates indicate protoplanets. At the same time, the separa- that the cumulative effect of long-range tion should be sufficiently smaller than the encounters is of no more than the same typical spacing b e t w e e n protoplanets that order as that of encounters with (Ar)min < the effect of the nearest neighboring proto- R~. When it is also noted that the relative planets on the encounter between the two velocities are greater for long-range enunder consideration can generally be ne- counters due to differential Keplerian moglected. For the s y s t e m a s s u m e d here, 100 tion o f the bodies in the system, it can protoplanets having a total m a s s equal to reasonably be argued that encounters with

708

COX AND LEWIS

make the major contribution to the relaxation of the system and that our results should not be significantly altered by neglecting the effect of long-range encounters. (Ar)mi n < R s

B. The Treatment o f Collisions between Bodies In a collision model in which the immediate formation of a solid body is assumed, the primary mechanisms for dissipating the impact kinetic energy are heating of the resultant body and escape of the remaining material. For masses as large as those considered here, m >- Me~50, the impact kinetic energy, if converted entirely into heat, would be sufficient to melt all of the material in both bodies. The heat capacity and melting temperature of the material in the bodies were taken to be those of ordinary chondrites, 0.166-0.182 cal/g and 1180-1350°C (Alexeyeva, 1958), respectively. Further, for bodies of this size the comminution energy is far less than the gravitational binding energy. It seems much more likely that the two impacting bodies are shattered into smaller fragments. From the computer simulations it is observed that the fraction of collisions with sufficiently high impact kinetic energy to disperse the material of both bodies completely to large distances is quite small provided the orbital eccentricities are not too large (Tables I III). We then conclude that although a few high-velocity fragments may escape carrying considerable kinetic energy off with them, most of the fragments, having smaller velocities, will in general form a gravitationally bound approximately spherical swarm of fragments. Since the size distribution of fragments for shattered bodies of this size, where gravitational binding forces dominate over material binding forces is unknown, we make the plausible assumption that the mass of a typical fragment is equal to that mass for which comminution energy equals gravitational binding energy. Since the fragment clusters are relatively small ( - 1 0 -3

AU) and the fragments are, on the basis of the above assumption, small (a few tens of kilometers in radius) and, therefore, numerous, collisions between fragments will be frequent. These collisions will eventually dissipate the kinetic energy of the cluster as heat, additional comminution energy, etc. If the cluster is not disrupted in the meantime by collisions or close encounters with other clusters or large bodies, it will shrink primarily by flattening, at least at first, as its kinetic energy is dissipated, eventually forming a solid body of mass approaching that of the cluster. Order-of-magnitude estimates assuming that one-tenth of the relative kinetic energy of colliding fragments is dissipated per collision show that these clusters will shrink to solid densities in a time which is two to three orders of magnitude less than the characteristic collision time between protoplanetary bodies during even the early stages of our simulations of planar systems. Considerably greater mean collision times between bodies are expected for more realistic nonplanar systems. A more detailed treatment of the collisions between protoplanets and the subsequent evolution of their collision remnants would of course be highly desirable, but is certainly too ambitious for the present computer simulations. Matsui and Mizutani (1978) have investigated this settling process for planar protoplanetary clusters neglecting the presence of the Sun and assuming that bodies accrete upon collision by numerically integrating the gravitational N-body equations of motion. For all of their models they find accretion times corresponding to a few days for both the Earth and the Moon systems. Collisions between bodies are therefore treated in the simulations as follows. When a close approach with separation less than the sum of the two radii is found by integration, the two colliding protoplanets are assumed to accrete and form a body with mass equal to the sum of the masses of the two incoming bodies. The orbital elements of the newly formed body are determined

TERRESTRIAL PLANET FORMATION f r o m the p o s i t i o n a n d v e l o c i t y o f the c e n t e r o f m a s s o f the t w o colliding b o d i e s at the instant of impact. The relative velocity of the b o d i e s u p o n i m p a c t is r e c o r d e d a n d u s e d to c a l c u l a t e a v a l u e for (T~ + U 1 ~ ) / ( - ~ ) , w h e r e Tj is the i n c o m i n g k i n e t i c e n e r g y at i m p a c t , U12 is the m u t u a l p o t e n tial e n e r g y o f the b o d i e s , a n d f l is the s u m o f the g r a v i t a t i o n a l b i n d i n g e n e r g i e s o f the i n c o m i n g b o d i e s . T h e c o n d i t i o n (T~ + U ~ ) / ( - f ~ ) < 1 is e q u i v a l e n t to Ei < 0, w h e r e Ei is the total e n e r g y o f the t w o - b o d y syst e m at i m p a c t , a n d is t h e r e f o r e u s e d to d e c i d e w h e t h e r the a s s u m p t i o n o f a c c r e t i o n w a s a valid one. V a l u e s o f the ratio (T~ + U ~ ) / ( - f ~ ) are t a b u l a t e d for the c o l l i s i o n s o c c u r r i n g in t h e t h r e e s i m u l a t i o n s pres e n t e d here in T a b l e s I - I I I . A d e t a i l e d s t e p - b y - s t e p d e s c r i p t i o n o f the a l g o r i t h m e m p l o y e d in the p r e s e n t s i m u l a TABLE I COLLISION STATISTICS(emax = 0.05) mI

mi

(Ms/100)

(M@/100)

2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 6 6 6 6 8 14

2 4 6 8 10 12 14 22 26 4 6 8 10 12 14 16 22 6 8 10 16 8 16

Aa

D~

Ab

(A only) 34 14 4 3 2 3 2 1 1 5 3 2 I 1 1 1 1 1 I 1 2 1 1

i

0.158 0.073 0.065 0.065 0.012 0.030 0.030 0.02L 0.039 0.084 0.063 0.046 0.004 0.012 0.031 0.022 0.693 0.169 0.033 0.020 0.031 0.037 0.027

aA = accretive; D = destructive collisions. b Simulation average value of A for the accretive collisions = 0.103.

709

TABLE II COLLISION STATISTICS(e.a~ = 0.10) mj (M~/100)

mi (M~/100)

A

D

A~ (A only)

2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 6 6 6 6 8 8 8 10 10 12

2 4 6 8 10 12 18 20 22 24 30 4 6 8 12 14 16 22 26 6 8 I0 12 12 14 22 12 20 20

27 18 6 2 3

5

0.325 0.153 0.102 0.256 0.096 0.005 0.024 0.009 0.002 0.009 0.I01 0.168 0.073 0.061 0.012 0.001 0.002 0.040 0.021 0.094 0.021 0.009 0.043 0.072 0.027 0.033 0.044 0.005 0.281

1 1 1 1 1 l

2 4

a Simulation average value of A for the accretive collisions = 0.171. t i o n s is g i v e n in (Cox, 1978, pp. 103-109) for t h o s e r e a d e r s i n t e r e s t e d in f u r t h e r details. C. T h e I n i t i a l C o n d i t i o n s

First, s t a r t i n g v a l u e s are c h o s e n for the m a s s e s a n d o r b i t a l e l e m e n t s o f the b o d i e s in the s y s t e m w h o s e e v o l u t i o n is to b e s i m u lated. All s i m u l a t i o n s to b e p r e s e n t e d here b e g i n with 100 i d e n t i c a l b o d i e s with m a s s e q u a l to 0.02 E a r t h m a s s e s . This artificial c h o i c e for the initial m a s s d i s t r i b u t i o n o f b o d i e s reflects the p r a c t i c a l l i m i t a t i o n t h a t o u r p r o g r a m , at its p r e s e n t level o f efficiency, c a n be u s e d to e v o l v e o n l y sys-

710

COX AND LEWIS T A B L E III COLLISION STATISTICS (emax = 0.15)

mj (M~/IO0)

mi (M~/IO0)

A

2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 6 6 6 6 8 8 8 10 10 10 18 26 24

2 4 6 8 10 12 18 20 24 32 58 4 6 10 14 20 8 18 20 24 10 12 24 12 14 30 40 34 60

20 15 9 6 3

2 3 2 2 1 1 1 1 1 1 2

D

Aa (A only) 0.438 0.253 0.266 0.292 0.052 0.208 0.065 0.006 0.022 0.039 0.030 0.163 0.165 0.094 0.051 0.050 0.025 0.000 0.115 0.068 0.015 0.255 0.192 0.025 0.013 0.089 0.235 0.104 0.012

a Simulation average value of A for the accretive collisions = 0.238.

tems containing no more than 100 or at most 200 bodies. The total mass of the s y s t e m is then very nearly equal to the total mass of the present terrestrial planets. The semimajor axes of the orbits of these bodies are required to fall within the range 0,5 to 1.5 AU, corresponding approximately to the range of semimajor axes of the terrestrial planets. The surface density in all o f these simulations was assumed to have an R -3/2 d e p e n d e n c e on heliocentric distance between 0.5 and 1.5 A U and to be zero elsewhere. The initial eccentricities of the orbits were randomly chosen from a uniform distribution on the interval from 0 to

some m a x i m u m eccentricity, emax- The angles of perihelion for the orbits are given b y a sequence of r a n d o m deviates from a uniform distribution on the interval (0, 27r). Finally, the time since perihelion passage for each b o d y is t a k e n from a uniform distribution on the interval (0, P ) where P is the orbital period of the given body.

D. Advantages of the Present Simulations First, the mass and relative velocity distributions of the protoplanetary bodies are free to evolve simultaneously. Once the initial conditions are chosen and the algorithm for handling encounters is established, the distributions of the orbital elements of the bodies evolve in response to the gravitational close encounters and occasional collisions occurring between bodies. E v e n t s involving larger bodies and m o d e r a t e mass ratios are included and are treated in the same way as events between pairs o f bodies with large mass ratios. Collisions b e t w e e n bodies arise in a natural way in these simulations and a s s u m e d expressions for the relative velocity distribution and the gravitational collision cross section are not needed. This is a fortunate feature, since integrations of close encounters between small bodies initially in circular orbits show characteristics that are not at all explicable by t w o - b o d y gravitational models. Integrations of close encounters b e t w e e n bodies initially in small-eccentricity orbits, e <~ 0.02, say, also indicate more complicated behavior than can be accounted for by two-body models. III. R E S U L T S OF T H R E E R E P R E S E N T A T I V E NUMERICAL SIMULATIONS

Figures 1 through 5 present results of three representative numerical simulations of the growth of the terrestrial planets by accretion of large protoplanets. These simulations primarily differ in the values chosen for emax. The time evolutions of the orbital eccentricities and semimajor axes are given for the bodies in each simulation in Figs. 1A, 2A, and 3A. In Figs. 1B, 2B,

TERRESTRIAL PLANET FORMATION

0.201

711

ema = = 0 . 0 5

)I-- 0.15 n~ I'Z (~ O.lO LIJ

0.05

1.81.71.61.5--I 1.4-

1.3~ 1.2i~ I . | '~ 1.0 IE ~ 0.9 OJ °3 0.8 0.7 0.6 0.~ 2.48 2.82 3.06 3,25 3,37 3.47 3.57 3.68 3.76 3,94 4.08 4.16 4.Z5 4.34 4~56 4.81

A

LOG ( T I M E I Y E A R S I )

FIG. IA. The time evolutions of the orbital eccentricities and semimajor axes are plotted for the bodies in the emax 0.05 simulation. For clarity in plotting an equal abscissa increment was allotted to each gravitational encounter, whether resulting in a collision or not. The associated event times are given as the log of the time in years. Collisions and subsequent accretions of bodies are denoted by a vertical line connecting the values of eccentricity or semimajor axis of the two colliding bodies just prior to the current encounter and by a small box locating the orbital eccentricity or semimajor axis of the newly accreted body. =

a n d 3B t h e s e m i m a j o r a x e s , e c c e n t r i c i t i e s , and masses of the system bodies are s h o w n for nine d i f f e r e n t s t a g e s o f e a c h s i m u l a t i o n i n c l u d i n g t h e initial a n d final o n e s . S e v e r a l f e a t u r e s a r e a p p a r e n t in t h e s e p lo ts o f th e simulations.

F i r s t , in t h e e c c e n t r i c i t y e v o l u t i o n p l o t s we note a redistribution of eccentricities for all o f t h e s i m u l a t i o n s : t h e n u m b e r o f b o d i e s w i t h v e r y sm al l v a l u e s o f e c c e n t r i c i t y decreases and a " t a i l " extending toward higher values of eccentricity develops.

712

COX AND LEWIS 0,2( [-

N = tO0

N=80

N=70

emax:O.05

o. lg

-~

O.ZO

N=60

>,I--

°o~oo ~ ~°o oO~o o~

--

N=50

N=40

0,1~ I'Z

0.10 Oo

o

O,Oe

'Io 0.20 - -

N=30

"N=20

~

j::L o ' ,

N=8

O.IS

0.10 oD

[]

0.05

[j~ -I i os

F.i~[ cutI.O

•

. ~J . i.t,

.

.

i. . . . o.~

,

/ Lo

,~

i

'

1.5

o.s

i/

~!o

i!

I

1.5

SEMIMAJOR AXIS (A.U)

FIG. lB. T h e semimajor axes, eccentricities, and m a s s e s o f the s y s t e m bodies are s h o w n for nine different stages in the emax = 0.05 simulation, including the initial and final ones. T h e n u m b e r o f s y s t e m bodies for the nine stages is indicated at the top of e a c h frame. T h e boundaries o f the initial s e m i m a j o r axis and eccentricity distributions are indicated by the rectangles in each frame. T h e areas o f the small b o x e s marking the location o f each body are proportional to the m a s s e s o f the bodies.

Also, an initial increase in eccentricity is apparent for the emax = 0.05 simulation. Where the lines in these figures are not too dense, one can o b s e r v e the range of typical eccentricity changes occurring during close encounters. Later, quantitative m e a s u r e s of these eccentricity changes are tabulated and discussed for the early (N > 90) evolution of each simulation. Also note that in the late stages o f evolution o f each simulation the smallest bodies are c o m m o n l y perturbed into orbits having relatively high eccentricity values. In the plots o f semimajor axis evolution one o b s e r v e s that " e d g e effects" are fairly m o d e s t in that bodies are only infrequently scattered out of their initial range of heliocentric distance, 0.5 to 1.5 AU. Further, bodies are more frequently scattered outward b e y o n d 1.5 A U than inward of 0.5 AU. One can also m a k e qualitative estimates f r o m these plots o f the range of semimajor axis changes typically occurring

during close encounters. A detailed examination of the semimajor axis plots reveals that pairs of bodies in orbits having either nearly identical or very dissimilar values of their semimajor axes rarely collide and accrete, e v e n when the orbits a p p r o a c h each other closely enough for collisions to take place. Most collisions then o c c u r for pairs of bodies with only m o d e r a t e l y dissimilar values of their semimajor axes. M a n y instances exist in the simulations where bodies in orbits having nearly identical values for their semimajor axes repeatedly interact without collision until an event happens to scatter one or the other of the bodies into an orbit with a m o d e r a t e l y dissimilar s e m i m a j o r axis value; for example, the pair o f bodies in the emax = 0.05 simulation with a ~ 1.2 A U during the period 4.38 - log t -< 4.60. In the plots of s y s t e m states at selected times (Figs. 1B, 2B, and 3B) one can m o r e easily follow the progress of accretion in

TERRESTRIAL PLANET FORMATION

713

n,," I-Z LLI C'= IJJ

4 X e~

0 -3

m,l

A

2.17 2.51 2.75 2.93 3.03

3J4 3.25 3.33 3,42 3.50 3.56 3.66 3,79 3 . 9 9 4.20 4.38

LOG (TIMEIYEARS)) FIG. 2 A . S a m e a s F i g . 1A, b u t f o r t h e emax = 0 . 1 0 s i m u l a t i o n .

each simulation, since the masses of the bodies, represented b y the b o x areas, are shown along with the orbital eccentricities and semimajor axes of the bodies. As expected, accretion is m o r e rapid for smaller heliocentric distances since the surface density of bodies is greater in that region. The evolution of the eccentricity and semimajor axis distributions is again a p p a r e n t in these plots. Also as expected, the m o r e

massive bodies tend to have smaller orbital eccentricities since the accretion process averages the eccentricity vectors o f the orbits o f the accreting bodies. The final states of the emax = 0.05 and 0.10 simulations contain too m a n y bodies with masses which are too small relative to those of Venus and Earth. On the other hand, the emax = 0.15 simulation ends with a final state which is very suggestive o f the terres-

714

COX AND LEWIS N= I 0 0

N=70

N:80

emox=O.[O ,r

=

"

•

~

~

~

=

=~

=

= • ~

~= =

==

,

~

,~=1

f

!

e= I

,~,

,

,=1

N= 40

N= 50

N=60

).. O.tS I-I1~ o J ¢ I.-

= ~c

w ( ~ o.o~ LU

7 7-

~

~=

=~

N=20

N;30

N= I0

7 r~

r

;7

ols

~

l/

i

7'

o',

I

I,s

,!o

,/,

SEMIMAJOR AXIS (A.U.)

FIG. 2B. S a m e as Fig. 1B, but for the emax = 0.10 simulation.

than one, a collision is expected to result in the complete disruption and dispersal of the material of both bodies; however, when A is less than one, a collision is expected to result in accretion. Values of A for collisions occurring in the three simulations are tabulated in Tables I-III. The numbers of accretive and destructive collisions, A and D , are also tabulated for each of the rn~ - mj values for which collisions occurred. There are a number of interesting observations which can be made from these statistics. First, the average A ratio for all of the collisions occurring in a given simulation increases as emax increases. Further, the average value of A for the ml = m2 = 0.02 M e collisions exceeds the average value of A when all of the collisions for a given simulation are included. It is no surprise then that the rnl = ms = 0.02 M e A. The Collision Statistics and ml = 0.04 M e , m2 = 0.02 M e colliThe condition A = (T~ + /-/12)/(-~) < 1 sions account for all of the destructive is equivalent to E~ < 0 and therefore pro- collisions. Perhaps the most significant vides a convenient measure which may be observation to be made from these statisused in deciding whether a given collision is tics is that the number of destructive colaccretive or destructive. When A is greater lisions rapidly increases as emax increases trial system of planets, with two of these bodies approaching the Earth and Venus in mass. The final states of those simulations which we have presented in Figs. 1-3 are compared with the terrestrial system of planets in Fig. 5. The heliocentric gaps between orbits of adjacent bodies are expressed as fractions of A -- 2.4 R [(rnl + rnz)/Mo] lfa, where R is the heliocentric distance to the center of the gap between ml and rn~. Birn (1973) has carried out numerical integrations and presented theoretical arguments based on the Jacobian integral which purport to show that bodies in orbits separated by gaps greater than A are stable in the sense that, except for oscillations, the semimajor axes of all orbits remain constant over long time scales.

TERRESTRIAL PLANET FORMATION

715

0'35 1 ernox =0.15

0.30 0.25"

>... I"('~ 0:: 0 . 2 0 I'Z LtJ (.~ 0 . 1 5 m,I

0.10 0 . 0 ~=

.,.-: :D (/) X

0

W

2 . 2 2 2.5 .4 2.7'5 2 , 9 0 3.03 3.13, 3.24 3.39 3.52 3.64 3.78 3.89 4.01 4 . 0 9 4.19 4 . 3 5

A

LOG (TIME(YEARS)) FIG. 3 A . S a m e as Fig. I A , but for the emax = 0.15 simulation.

b e y o n d 0.10, so that for simulations with initial values o f emax > 0.15, one is clearly no longer justified in assuming that all collisions result in accretions.

B. Eccentricity Evolution in the Simulations In Table IV statistics are tabulated pertaining to the eccentricity changes occurring during the early stages o f four repre-

sentative simulations including the three presented here. The quantity denoted in the table by Ae is the average value of miAei + mjAej over those events occurring b e t w e e n the simulation starting time and the 10th accretion event, where Aei and Aej are the eccentricity changes o f the t w o bodies in the given event and the masses mi and mj are in units o f the initial mass, 0.02 M s . In these statistics w e refer to those gravita-

716

COX AND LEWIS TABLE IV

encounter)/(lAel/encounter) quantifies the rate of the eccentricity increase due to gravitational encounters. OCCURRING FOR ENCOUNTERS AND COLLISIONS F r o m the simulation starting time to the emax 10th accretion event, the m a s s average eccentricity increases quite rapidly for the 0.025 0.05 0.10 0.15 emax = 0.025 and emax = 0.05 simulations Ae/encounter +0.0027 +0.0018 +0.0007 +0.0007 and decreases for the em~x = 0.10 and 0.15 Ae/collision -0.0133 0 . 0 1 4 3 -0.0507 -0.0846 simulations. It is e x p e c t e d that an emax IAel/encounter 0.0100 0.0104 0.0071 0.0057 Ae/encounter value s o m e w h e r e b e t w e e n 0.05 and 0.10 0.27 0.17 0.10 0.12 I-~el/encouoter would on the average result in a balance XAe(encounters) +2.96 + 1.59 +0.20 +0.41 b e t w e e n the eccentricity d e v e l o p e d in grav[EAe(collisions)J Collision-encounter itational encounters and that lost in collibalance ratio 5.0 8.0 68.7 128.2 sions for the coplanar s y s t e m of bodies Number of encounters during the first 10 chosen here to initiate all the simulations. collision events 147 140 137 527 One simulation was run with the same initial conditions as those o f the three described here, but with the difference that tional interactions not resulting in collisions encounters with close a p p r o a c h distances as encounters. N o t e that in all of the simu- small enough to result in collisions were lations, at least for the early stage consid- ignored. This simulation allowed us to obered here, the gravitational encounters on serve the eccentricity increases due to the average result in eccentricity increases, close gravitational interactions and the rebut as ernax is increased the mass average suiting distribution of orbital eccentricities eccentricity increase per e n c o u n t e r de- in a s y s t e m of bodies where no collisions creases. Collisions, on the other hand, al- were allowed. The initial eccentricities ways result in eccentricity decreases, and were randomly chosen from a uniform disthe eccentricity decreases due to collisions tribution on the interval from 0 to em~x = increase as em~x is increased. The mass 0.005. A steady increase in the average average eccentricity of the s y s t e m of proto- eccentricity was o b s e r v e d until the simulaplanetary bodies is then determined by the tion was terminated at 9800 years at which balance between the eccentricity increases time the average eccentricity of the 100 resulting from gravitational encounters bodies had reached 0.054. and the eccentricity decreases resulting In a s y s t e m where nonzero orbital inclifrom collisions. The ratio Y.Ae nations are allowed, fewer close gravita(encounters)/IEA e (collisions)] provides tional interactions are e x p e c t e d to result in a quantitative m e a s u r e o f this balance; collisions. The n u m b e r of encounters per when this ratio is greater than one, the collision for a three-dimensional s y s t e m of m a s s average eccentricity is increasing bodies with the same m e m b e r s as the sysand, conversely, when it is less than tem used to initiate the two-dimensional one, the mass average eccentricity is de- simulations is e x p e c t e d to be at least one creasing. The collision-encounter balance order o f magnitude greater than the n u m b e r ratio gives the n u m b e r of encounters re- in the two-dimensional case. The eccentricquired on the average to offset the eccen- ity balance referred to a b o v e is then extricity decreases due to collisions. A mea- pected to o c c u r for higher values o f emax in sure of the average absolute eccentricity the three-dimensional case than in the twochange per e n c o u n t e r is provided by dimensional one. The necessary eccenIAel/encounter = ([miAeil + ]mjAejl)/ en- tricity to produce terrestrial-like final counter and finally, the ratio ( A e / configurations of bodies such as those in STATISTICS OF THE ECCENTRICITY C H A N G E S

TERRESTRIAL PLANET FORMATION the emax = 0.15 simulation might then be built up by the additional gravitational encounters occurring in the three-dimensional case.

lows. F o r a given b o d y , ~ 7 5 % of the material comes from the zone bounded by the orbits o f the next inner, i.e., toward the Sun, and the next outer bodies, 8% comes from the zone bounded by the orbits of the next inner and second inner bodies, and the remaining 17% comes from the zone bounded by the orbits of the next outer and second outer bodies. H o w e v e r , we must caution that this radial sampling function is based on the statistics of a small n u m b e r of c o m p o n e n t bodies. Furthermore, it undoubtedly reflects the edge effects resulting from our choice o f the initial radial distribution of material. The similarity of the final state of the emax = 0.15 simulation to the system of terrestrial planets does permit us to form the above working hypothesis for the radial sampling to be expected for at least the Earth and Venus. Much more complex " e d g e effects" have no doubt influenced the radial sampling of the material comprising Mercury and Mars. It must be pointed out that although the radial sampling observed for bodies in the simulations represented only the sampling occur-

C. Heliocentric Distance Sampling of Material It is informative to examine the range of heliocentric distance sampling of material for the bodies in the simulations at several stages o f their evolution. In Fig. 4 the original semimajor axis distributions of c o m p o n e n t bodies are shown for those bodies with mass greater than 0.1 M e at six stages, including the final one, of the emax = 0.15 simulation. As expected, the range of heliocentric distance sampling in bodies with smaller mass tends to be narrower than the range of sampling in bodies with larger mass, and the range o f radial sampling was found to be greater in the emax 0.15 simulation than in the ema x = 0.05 simulation for bodies o f a given mass. The radial sampling function observed for the bodies in the final state of the em~x = 0.15 simulation may be described as fol=

N=IO0

N=80

0.~0

emax=

o.io

N=70

0.15

"°*

o

....

717

: o!,oo ° N= 6 0

, ~,°~

N= 5 0

....

°l,,

N= 40

> " 0.=0 I'o,i~

.

.

.

.

_

° o

,..

o.o~

°o ~ o°o :ii

#

° O~o~

,;:o,t

o

....o.,of~ ,

o.io

~

N=6

,,=3o

] <1

o o.~

o -

,

~

/

J

° E] o

,

7. °°, ..

o []

o.os

B

_o

[]

~o

0

I O

i!o . . . .

o

I J " i.s

0 []

"

o.s

. . . .

oo

o iI o

~.o . . . .

SEMIMAdOR

.

'

'

o.~

to

A X I S (A.U.)

FIG. 3 B . S a m e a s F i g . 1B, b u t f o r t h e emax = 0.15 s i m u l a t i o n .

1.5

718

COX AND LEWIS z° I I0

o

N= 6

N=60

'°i-

°I

N=So

9C

emo x = O. 15

aQ

."

2 .IL:.:?

,...

O3

5: 3 0 |~

* ~

.

~°I

N=40

20 " F " "- "

eo

,~ ........... --,~

io

O.5

LO

1.15 '

,

,

-

.N=20

. . . . . •

0.5

Io

I0

1.5

SEMIMAJOR

05

I

I0

,

I

15

A X I S (A.U.)

FIG. 4. The original semimajor-axis distributions of component bodies are shown for those bodies with mass greater than 0.1 M e at six stages, including the final one, of the emax = 0.15 simulation. The original semimajor axes of the component bodies are indicated by a " + " and the current semimajor axis and mass of each body are located by an " ' × " .

ring during the last stage of accretion, the range o f radial sampling is expected to be much narrower for bodies with small values of their mass. Therefore, it is apparent that the range o f radial sampling exhibited in the final bodies is primarily determined by the last few accretion events involving the given body. Hartmann (1976) has used Wetherill's (1975) calculations of the gravitational dispersal of planetesimals from different orbits about the Sun to obtain data on the sources o f different late-accreted materials added to the terrestrial-planet surface layers. Wetherill tabulated the final impact destinations of statistical samples of particles starting in nine different orbits in the terrestrial region o f the solar system. In his Table III (p. 555), Hartmann has converted the relative numbers of impacts obtained by Wetherill to obtain the percentage of mass striking target planets from each source orbit. These results indicate that large fractions o f the mass added to the terrestrial planets in their later stages o f growth probably formed nearer to other planet's orbits, in agreement with conclusions drawn from the present simulations. Although our radial sampling function pertains to a much earlier stage of accretion than that given by

Hartmann, the two functions seem to be qualitatively similar.

D. Numerical Checks of the Simulations One numerical simulation was run to check the validity of the encounters model for those close encounters typically occurring in the present simulations. Comparison of values for the orbital element changes calculated using the model with values determined by integration o f the regularized equations o f motion revealed that the model, in agreement with the conclusions stated in Paper I, does provide a useful approximation for the orbital element changes occurring in those close encounters for which e ~> 4 and a ~< 0.80. A further check on the accuracy of the simulations was made by calculating the total angular momentum o f the system of protoplanetary bodies for the initial and final states of the emax = 0.15 simulation. The total angular momentum is closely conserved (to the fourth digit), as it should be, while the total system energy decreases by 1.23% between the initial and final states of the emax = 0.15 simulation. Finally, the numerical integration of the individual encounters was checked by integrating backwards in time to make sure

TERRESTRIAL PLANET FORMATION

mass average eccentricity values in the range 0.025 < ( e ) < 0.05 can be sustained by close gravitational encounters offsetting the eccentricity lost in collisions, at least for the early stages of evolution o f the coplanar systems chosen to begin these simulations. These values of ( e ) are somewhat too low to produce terrestrial-like systems, but higher values of ( e ) may result from the higher ratio of close gravitational encounters to collisions holding for three-dimensional systems. An estimate of the time interval required for this stage of accretion in a reasonably thin hypothetical three-dimensional protoplanetary system was made using the ratio of the t w o - b o d y collision probability for a three-dimensional system to that for a twodimensional system along with the evolution times obtained in the two-dimensional c o m p u t e r simulations. The collision probability per unit time for one body may be written as .7--1 = Uno', where U is the mean relative speed of the bodies, n is the number of bodies per unit surface area or volume, and (r is the collision cross section. The collision cross section in two dimensions was taken to be 2 ~ and in three

that the initial conditions could be reproduced. IV. CONCLUSIONS DRAWN FROM THE SIMULATIONS REGARDING THE FORMATION OF THE TERRESTRIAL PLANETS Perhaps the primary conclusion is that terrestrial-like planetary systems can be formed as a result of the gravitationallyinduced collisions and accretions occurring between bodies in coplanar systems of protoplanets such as those chosen here to initiate the simulations, provided the initial average orbital eccentricity o f the system of bodies is in the range 0.05 -< ( e ) --- 0.10. Simulations with initial values of ( e ) smaller than 0.05 are found to result in finalstate configurations having many bodies o f nearly equal mass, none o f which approach the Earth or Venus in mass, Fig. 5. On the other hand, initial values of ( e ) greater than 0.10 result in substantial numbers of catastrophically disruptive collisions between bodies and it is expected that the accretion process for this case will at least be dramatically slowed, if not reversed. The numerical simulations indicate that

emox:O.05

12 18 5 2 N~I~L5F,q 2.6

3,6 24 I-,-I 3.3 ~

6.4

2 22

emox= Ol I 5

I'-~H

60

zo

H

9.6

1.3

58 I'~

I

I

0,3

I

I 0,4

12.9

L

0.5

84

6

I-,-II

•

6 •

54 3.5

4.4

4.4

26 I-,-12.71

L~/L~ s

2 •

I

26

I

I

5.5

0.7

5

Terrest riol

Plonels

4.1

6 52 50 18 $ 2 8 30 16 I,-~1391.~.13.81-~1 I.,.=.-I I,.~1-=,-I__4.1 I-~ 3.9 ~ 1.9 1.4 0.3

emax:O, lO

719

I

0.6

81 H

L

0.7

I00 I~1

6.0

I

0.8

I

0.9

I

LO

8.5

[

I.I

I

I.;~

II •

I

I

1.3

I

1.4

I

1.5

I

1.6

I

1.7

1__

1.8

A.U.

FIG. 5. A comparison of the final states of the three representative simulations with the present terrestrial system of planets. The bars indicate the heliocentric range of each body. The numbers above the bars give the mass of the bodies in hundredths of an Earth mass, and the heliocentric distance gaps between orbits are expressed in multiples of 2~ = 2.4R [m1 + m=)/Mo] 113,where A~is a gap argued to be stable by Birn (1973).

720

COX AND LEWIS

dimensions to be ~r~ 2, where ~ is a multiple /3f times the physical radius o f the bodies, r. The factor/3f, which provides a measure of the increase of the collision cross section o v e r the geometric cross section due to "gravitational focusing," was estimated from the integrations of encounters between pairs o f bodies initially in circular orbits to be a p p r o x i m a t e l y 28 for the twodimensional case with rnl = rn2 = 0 . 0 2 M , . The surface area available to the 100 bodies in the simulation is S = rr(1.52 - 0.52) -- 2rr A U 2. The corresponding volume a s s u m e d to be available to the 100 bodies in a threedimensional system was taken to be equal to that o f the solid of revolution formed by rotating a regular trapezoid, b o u n d e d a b o v e b y ~b = / m a x , below b y ~b = - / m a x , and b y R1 = 0.5 A U and R2 = 1.5 A U in heliocentric distance, one full revolution. For a set of test bodies having a dimensionless relative velocity, U = u/Vk, where V z = G M o / a , with respect to a field b o d y in a circular orbit of radius a , there exist expressions for emax and /max in t e r m s of U (Opik, 1951). The m a x i m u m possible eccentricity occurs when the orbits are coplanar and when the perihelion of the test b o d y ' s orbit is at the field b o d y ' s orbit and m a y be written in terms of U as emax = U z + 2 U. For any U, the m a x i m u m inclination (relative to the orbital plane of the field body) is, according to 0 p i k , /max = 2 sin-l(U/2). For small values o f U, e m a x = 2 U and imax = U are good approximations. With the values er~ax = 0.15 (/max = 0.075 rad), ml = m2 = 0.02 M s , and/3f = 28, the ratio r3-D/r2-D is determined to be approximately 280. This ratio with the emax = 0.15 simulation time of 61,400 years yields an estimate of 2 × 107 years for the time required for the stage of accretion considered in our simulations. This time falls well within the limits of a few times 108 years for the formation of the terrestrial planets, based on R b - S r model ages for returned samples o f lunar soil (Nyquist, 1977) along with the K4°-A 4° and U, T h - H e 4 gas-retention ages for various meteorites.

V. LIMITATIONS OF THE PRESENT SIMULATIONS AND POSSIBLE FUTURE IMPROVEMENTS Many of the desirable features of these simulations were described earlier (Section II-D). In this section we will turn our attention to some of their limitations. First, actual systems of protoplanets, even if they are quite thin, are nevertheless three-dimensional, and this fact presents a whole new set of ditticulties in attempting to simulate such systems. It was argued earlier that the n u m b e r of encounters per collision for a three-dimensional system o f bodies having the same bodies as those in the s y s t e m used to initiate the two-dimensional simulations is e x p e c t e d to be at least one order of magnitude greater than the n u m b e r for the two-dimensional case. Therefore, significantly more c o m p u t e r time would be required to simulate the three-dimensional systems. Further, eccentricity balance is e x p e c t e d to o c c u r at higher values o f ( e ) in the three-dimensional case than in the twodimensional one, but since the character of the encounters p r e s u m a b l y also changes in three dimensions, one cannot be sure that this expectation is justified in fact. Several simulations which were carried out have not been discussed here, because the orbital spacings b e t w e e n bodies were not stable (Birn, 1973), when no further close approaches, ( A r ) m i n < Rs, were predicted and the simulations terminated. In the two-dimensional case this situation is m o s t often the result of an alignment of some n u m b e r of the orbital eccentricity vectors of the s y s t e m bodies. Inclusion of the effects o f distant encounters in some o f these " f i n a l " states would give rise to s o m e additional collisions and subsequent evolution of these systems. " A v o i d a n c e " mechanisms such as this are e x p e c t e d to be an e v e n greater problem in three-dimensional simulations which include only the effects of close encounters. Another limitation of the simulations in their present form is the small n u m b e r (N = 100) of bodies which can currently be han-

TERRESTRIAL PLANET FORMATION

721

d i e d b y t h e p r o g r a m s . T h i s m a k e s it difficult sions, and Molly Rauber for typing the manuscript. to a s s e s s t h e i n f l u e n c e o f o u r c h o i c e o f a n Computing time for this research was provided by NASA Grant NGL-22-009-521. initial s y s t e m o f o n e h u n d r e d b o d i e s o f e q u a l m a s s , rn = 0.02 M s , on t h e final c o n f i g u r a t i o n s p r o d u c e d in t h e s i m u l a t i o n s . REFERENCES Two-dimensional simulations beginning ALEXEYEVA, K. N. (1958). Physical properties of with protoplanetary systems having more stony meteorites and their interpretation in the light bodies and a variety of mass distributions of hypotheses regarding the origin of meteorites. Meteroritika 16, 67-77; in Russian. s h o u l d e l u c i d a t e t h e i n f l u e n c e o f t h e initial BATTIN, R. H. (1964). Astronautical Guidance. Mcc o n d i t i o n s o n t h e final c o n f i g u r a t i o n s p r o Graw-Hill, New York. duced. Such simulations seem feasible BETTlS, D. G., AND SZEBEHELV,V. (1971). Treatment upon some improvement and/or modof close approaches in the numerical integration of ification o f t h e w a y in w h i c h e n c o u n t e r s the gravitational problem of N bodies. Astrophys. a r e h a n d l e d in t h e s i m u l a t i o n p r o g r a m s . Space Sci. 14, 133-150. BraN, J. (1973). On the stability of the planetary W e h a v e a t t e m p t e d h e r e to s i m u l a t e t h e system. Astron. Astrophys. 24, 283-293. last stage of terrestrial planet formation Cox, L. P. (1978). Ph.D. dissertation, Department of n e g l e c t i n g a n y o u t s i d e i n f l u e n c e s . If, for Earth and Planetary Sciences, M.I.T. example, Jupiter has already attained a Cox, L. P., LEWIS, J. S., AND LECAR, M. (1978). A model for close encounters in the planetary probs u b s t a n t i a l p a r t o f its m a s s b e f o r e t h e s t a g e lem. Icarus 34, 415-427. of accretion considered here, the influence DOLE, S. H. (1962). The gravitational concentration of of bombarding material perturbed by Jupiparticles in space near the Earth. Planetary Space t e r into o r b i t s c r o s s i n g the t e r r e s t r i a l r e g i o n Sci. 9, 541. m u s t b e i n c l u d e d in a r e a l i s t i c s c e n a r i o ( s e e GIULI, R. T. (1968). On the rotation of the Earth produced by gravitational accretion of particles. W e i d e n s c h i l l i n g , 1975). Icarus 8, 301-323. F i n a l l y , s o m e c o l l i s i o n s do o c c u r in t h e s e GOLDREICH,P., AND WARD, W. R. (1973). The formasimulations which, based on their expected tion of planetesimais. Astrophys. J. 183, 1051-1061. r e l a t i v e v e l o c i t i e s at i m p a c t , s h o u l d b e GREENBERG, R., WACKER,J. F., HARTMANN, W. K., completely disruptive, and yet the simulaAND CHAPMAN, C. R. (1978). Planetesimals to planets: Numerical simulation of collisional evolut i o n s a s s u m e t h a t all c o l l i s i o n s r e s u l t in tion. Icarus 35, 1-26. a c c r e t i o n . T h i s a s s u m p t i o n a p p e a r s to b e HARTMANN,W. K. (1976). Planet formation: Compow e l l j u s t i f i e d for c o l l i s i o n s o c c u r r i n g in s y s sitional mixing and lunar compositional anomalies. tems with mass average eccentricities (e) Icarus 27, 553-559. ~< 0.05. It c a n n o t b e r u l e d o u t , h o w e v e r , MATSUI, T., AND MIZUTANI, H. (1978). Gravitational that bodies were accreted and subsequently N-body problem on the accretion process of terrestrial planets. Icarus 34, 146-172. d i s p e r s e d b y c o l l i s i o n s s e v e r a l t i m e s in t h e NYQUIST, L. E. (1977). Lunar Rb-Sr chronology. process of forming the terrestrial planets. Phys. Chem. Earth. 10, I03-142. In o r d e r to s i m u l a t e p r o p e r l y this f o r m a t i o n 0PIg, E. J. (1951). Collision probabilities with the p r o c e s s , it w o u l d b e n e c e s s a r y to c o n s t r u c t planets and the distribution of interplanetary matter. Proc. Roy. Irish Acad. Sect. A 54, 165. models of the disruptive collisions comSAFRONOV, V. S. (1969). Evolution o f the Protoplanem o n l y o c c u r r i n g in t h o s e s y s t e m s w i t h ( e ) tary Cloud and Formation o f the Earth and Planets. > 0.05. F u r t h e r m o r e , a c e r t a i n n u m b e r o f Israel Program for Scientific Translations. NASA v e r y c l o s e e n c o u n t e r s in all o f t h e s i m u l a TT F-677. t i o n s p r e s u m a b l y l e a d t o t i d a l d i s r u p t i o n o f WEIDENSCHILLING, S. J. (1974). A model for accretion of the terrestrial planets. Icarus 22, 426-435. one or both bodies. WEIDENSCHILLING,S. J. (1975). Mass loss from the region of Mars and the asteroid belt. Icarus 26, 361366. ACKNOWLEDGMENTS WETHERILL, G. W. (1975). Late heavy bombardment The authors thank Drs. Charles C. Counselman, of the Moon and terrestrial planets. Proc. Lunar Myron Lecar, and William Ward for helpful discusPlanet Sci. Conf. 6th.