The Effect of Tidal Interaction with a Gas Disk on Formation of Terrestrial Planets

Icarus 157, 43–56 (2002) doi:10.1006/icar.2001.6811, available online at http://www.idealibrary.com on

The Effect of Tidal Interaction with a Gas Disk on Formation of Terrestrial Planets Junko Kominami and Shigeru Ida Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Ookayama, Meguro-ku, Tokyo, 152-8551, Japan E-mail: [email protected] Received June 12, 2001; revised November 5, 2001

these runs, we used exponential decay model with e-folding time of 3 × 106 years. The orbits of protoplanets are stablized by the eccentricity damping in the early time. When disk surface density decays to
1% of the minimum mass disk model, the damping force is no longer strong enough to inhibit the increase of the eccentricity by distant perturbations among protoplanets so that the orbital crossing starts. In this disk decay model, a gas disk with 10−4 –10−3 times the minimum mass model still remains after the orbital crossing and accretional events, which is enough to damp the eccentricities of the Earth-sized planets to the order of 0.01. Using these results, we discuss a possible scenario for the last stage of terrestrial planet formation. c 2002 Elsevier Science (USA) Key Words: terrestrial planets; planetary formation; solar nebula; Earth; planetary dynamics.

We have performed N-body simulation on final accretion stage of terrestrial planets, including the effect of damping of eccentricity and inclination caused by tidal interaction with a remnant gas disk. As a result of runway and oligarchic accretion, about 20 Marssized protoplanets would be formed in nearly circular orbits with orbital separation of several to ten Hill radius. The orbits of the protoplanets would be eventually destabilized by long-term mutual gravity and/or secular resonance of giant gaseous planets. The protoplanets would coalesce with each other to form terrestrial planets through the orbital crossing. Previous N-body simulations, however, showed that the final eccentricities of planets are around 0.1, which are about 10 times higher than the present eccentricities of Earth and Venus. The obtained high eccentricities are the remnant of orbital crossing. We included the effect of eccentricity damping caused by gravitational interaction with disk gas as a drag force (“gravitational drag”) and carried out N-body simulation of accretion of protoplanets. We start with 15 protoplanets with 0.2M⊕ and integrate the orbits for 107 years, which is consistent with the observationally inferred disk lifetime (in some runs, we start with 30 protoplanets with 0.1M⊕ ). In most runs, the damping time scale, which is equivalent to the strength of the drag force, is kept constant throughout each run in order to clarify the effects of the damping. We found that the planets’ final mass, spatial distribution, and eccentricities depend on the damping time scale. If the damping time scale for a 0.2M⊕ mass planet at 1 AU is longer than 108 years, planets grow to Earth’s size, but the final eccentricities are too high as in gas-free cases. If it is shorter than 106 years, the eccentricities of the protoplanets cannot be pumped up, resulting in not enough orbital crossing to make Earth-sized planets. Small planets with low eccentricities are formed with small orbital separation. On the other hand, if it is between 106 and 108 years, which may correspond to a mostly depleted disk (0.01–0.1% of surface density of the minimum mass model), some protoplanets can grow to about the size of Earth and Venus, and the eccentricities of such surviving planets can be diminished within the disk lifetime. Furthermore, in innermost and outermost regions in the same system, we often find planets with smaller size and larger eccentricities too, which could be analogous to Mars and Mercury. This is partly because the gravitational drag is less effective for smaller mass planets, and partly due to the “edge effect,” which means the innermost and outermost planets tend to remain without collision. We also carried out several runs with time-dependent drag force according to depletion of a gas disk. In

1. INTRODUCTION

Planets accrete from planetesimals. Runaway growth of planetesimals results in rapid formation of protoplanets (Greenberg et al. 1978, Wetherill and Stewart 1989, 1993, Spaute et al. 1991, Barge and Pellat 1991, Kokubo and Ida 1996, Inaba et al. 2001). A small number of planetesimals run away, while most planetesimals remain small. Oligarchic growth follows the runaway growth, resulting in similar-sized protoplanets with several to 10 Hill radius orbital separations (Kokubo and Ida 1995, 1998, 2000, Weidenschilling et al. 1997). Distant perturbations among the protoplanets and dynamical friction from the small planetesimals are responsible for the configuration of the protoplanets. The dynamical friction damps the orbital eccentricities and inclinations of the protoplanets, which means their orbits are almost circular and coplanar. Kokubo and Ida (1998, 2000) showed that about 20 Marssized protoplanets are finally formed in terrestrial planet region as a result of oligarchic growth, in the case of the minimum mass disk model (Weidenschilling 1977, Hayashi 1981). Since the mass of the protoplanets is significantly smaller than the mass of Earth or Venus, further coagulation between the protoplanets would be required. The time scale of oligarchic growth is 105 –106 years in the terrestrial planet region (Kokubo and Ida 2000). The orbits of 43 0019-1035/02 $35.00 c 2002 Elsevier Science (USA)  All rights reserved.

44

KOMINAMI AND IDA

the protoplanets would not be stable on a longer time scale. Mutual gravity among the protoplanets (Chambers et al. 1996, Yoshinaga et al. 1999) and seculer resonance from giant planets (Ito and Tanikawa 1999, Nagasawa et al. 2000) may pump up the eccentricities of the protoplanets. Then, orbital crossing begins and accretion of the protoplanets starts. Chambers and Wetherill (1998) and Agnor et al. (1999) calculated accretion of protoplanets during orbital crossing by N-body simulation, starting from initial conditions consistent with the result of oligarchic growth. In many cases, some of resulting planets are as massive as Earth or Venus. However, the eccentricities of the planets are rather high (∼0.1). They are about 10 times the present eccentricities of Earth and Venus. These high eccentricities obtained by N-body simulations are the remnant of orbital crossing of the protoplanets. The time-averaged eccentricities of Earth and Venus are ∼0.03. The eccentricities of the Earth and Venus are oscillating between 0.01–0.06 in a period ∼105 years. Since the forced eccentricity by Jupiter and Saturn is ∼0.03, free eccentricities of Earth and Venus may be ∼0.03. Since we do not include the perturbation from the giant planets in the present calculations, accretion simulations should produce the eccentricity
0.03 for Earth-sized planets. This value of eccentricity is still smaller than the results in most of the previous simulations. The inclinations of Earth and Venus are also oscillating between 0–0.07 radian with forced inclinations ∼0.03 radian. Hence, inclinations
0.04 radian must be produced by accretion simulations. Since inclinations (in radian) are usually pumped up less than eccentricities, we hereafter deal with the eccentricities, although we carry out three-dimensional calculations. To make a planetary system that has planets with orbits and masses similar to those of the present terrestrial planets, the high eccentricities of the large planets have to be diminished. Collisional damping is not strong enough to diminish the eccentricities. Further drag force is needed to damp the excess eccentricities after the collisions of the planets. As will be discussed in Section 2, the damping due to tidal interaction with a protoplanetary gas disk (Ward 1989, 1993, Artymowicz 1993, Agnor and Ward 2002), which we will refer to as “gravitational gas drag,” is a good candidate; aerodynamic gas drag (e.g., Adachi et al. 1976) is too weak to diminish the eccentricities of the large planets in the final stage. Dynamical friction from leftover planetesimals can also be effective (Chambers 2001; also see Sections 2 and 4.3). However, including such drag force might prevent orbital crossing of protoplanets in the first place (Iwasaki et al. 2001, 2002) or lead to isolation between protoplanets before sufficient coagulation to form Earth-sized planets. Is there drag force that is weak enough not to prevent the growth of planets to the Earth’s size, but strong enough to damp the eccentricities of the Earth-sized planets before removal of the protoplanetary disk? The calculation we have carried out shows that such strength of the drag is possible from the gravitational gas drag due to a mostly depleted disk (with gas surface density of ∼10−4 –

10−3 times that in the minimum mass model). The full amount of disk gas in the minimum mass model prevents the onset of orbital crossing (Iwasaki et al. 2001, 2002). The decrease in the disk surface density would trigger orbital crossing and allow enough coagulation among the protoplanets. As will be shown in Section 4.1, the orbital crossing stage lasts for a few times 106 years to form Earth-sized planets and their eccentricities of the Earthsized planets are diminished on a time scale ∼107 years after the orbital crossing stage. The observationally inferred lifetime of protoplanetary disks is ∼107 years (e.g., Strom et al. 1993, Zuckerman et al. 1995). In the calculations by Chambers and Wetherill (1998) and Agnor et al. (1999), the growth time scale in the orbital crossing stage is 107 –108 years. Monte Carlo simulations (e.g., Wetherill 1992) show comparable or rather longer time scales. However, planets as large as Earth are formed on a time scale of ∼106 –107 years. The remaining 107 –108 years are required for small bodies with high eccentricities to collide with the Earth-sized planets. If the small amount of gas exists, it damps the high eccentricities of the small bodies to isolate the bodies and/or to accelerate the accretion by the Earth-size bodies. Cosmochemical evidences suggest that the Earth–Moon impact and Earth’s core formation occurred 2.5–6.5 × 107 years after the start of the Solar system (e.g., Halliday et al. 2000). Our models suggest that orbital crossing among protoplanets starts when disk gas is depleted to 10−2 –10−3 times the gas of the minimum mass model. The observationally inferred disk depletion time scale ∼107 years would correspond to decay time scale by several factors. Since the decay by a factor 102 –103 takes several e-folding times if exponential decay is assumed, our model suggests that orbital crossing starts several times 107 years after formation of the disk and it persists for
107 years, which may be consistent with the cosmochemical evidences. Chambers (2001) did some calculations with 150–160 leftover planetesimals without disk gas. It is reported that compared to the cases without planetesimals, the resulting planets have smaller eccentricities, although they are still larger than the eccentricities of Earth and Venus. As shown in Section 2, dynamical friction by planetesimals has a similar effect on planetary orbits as gravitational gas drag. We use the strength of the drag force as a parameter, and investigated the effect of the drag force extensively and systematically. The gravitational interaction between a disk and a protoplanet diminishes not only the eccentricity of the protoplanet but also its semimajor axis, which means migration toward the sun (e.g., Goldreich and Tremaine 1980, Ward 1986, 1989, 1997). However, such effect on the semimajor axis is neglected in the particular calculation in this paper for the following reason. Since the migration is caused by the strength difference of the torques of inner and outer part of the disk, while the eccentricity damping is caused by the torques themselves, the migration time scale is longer than the eccentricity damping time scale by a factor ∼(vkep /cs )2 where vkep is Keplarian velocity and cs is sound velocity of the disk gas (e.g., Ward 1993). In the minimum

45

TERRESTRIAL PLANET FORMATION

mass disk model, (vkep /cs )2 ∼ 102 –103 . Furthermore, asymmetric dissipation of density waves between inner and outer parts may make the migration time scale even longer (Tanaka et al. 2001). In conclusion, the migration time scale is much longer than the eccentricity damping time scale. In most of our calculations with a significantly depleted disk, the eccentricity damping time scale is 106 –107 years, so that the migration time scale is 108 –1010 years. Hence, we neglect the migration in the calculation here. However, the migration may be fast enough to affect orbital configuration of the protoplanets before the orbital crossing, which is the initial condition of our calculation, when the disk gas is not significantly depleted. The effect of the migration on formation of protoplanets is a future issue. In this paper, we carried out N-body simulation to show that it is possible to make planets that are about the size of Earth and have nearly circular orbits, by including the gravitational drag force due to a mostly depleted disk. The drag force is explained in the next section. In Section 3, the calculation method is explained. We show numerical results and analytical arguments to explain the results in Section 4. Some discussions are made to present a possible consistent scenario to make terrestrial planets in Section 5.

where M⊕ is the Earth’s mass, vrel is relative velocity between a body and disk gas, which is ∼evkep (vkep : Keplerian velocity), M is the mass of the body, rp is its physical radius, and ρgas is gas density. Since τaero ∝ M 1/3 , aerodynamical gas drag is more effective for a smaller body. On the other hand, gravitational gas drag is more effective for a larger body. The damping time scale of gravitational gas drag is  τgrav 

One can think of several possible drag forces for protoplanets. One of such forces is dynamical friction from leftover planetesimals (e.g., Stewart and Wetherill 1988, Ida and Makino 1992b). Dynamical friction plays an important role in runaway growth (e.g., Wetherill and Stewart 1989, Kokubo and Ida 1996). However, the result from the N-body calculation done by Kokubo and Ida (1998, 2000) shows that most planetesimals accrete to form protoplanets in the final stage of oligarchic growth; the total mass of the leftover planetesimals might not be enough for effective dynamical friction. Chambers (2001) performed N-body simulation of accretion of protoplanets including the leftover planetesimals with their total mass equal to the total mass of the protoplanets. Since the effect of dynamical friction is incorporated into gravitational gas drag as mentioned below, dynamical friction from planetesimals is not particularly considered here. Other forces are aerodynamical gas drag (e.g., Adachi et al. 1976) and gravitational gas drag (Ward 1989, 1993, Artymowicz 1993). Comparison between the damping time scales of two forces shows that aerodynamical gas drag can be neglected in this case. Damping time scale of aerodynamical gas drag is given by (e.g., Adachi et al. 1976, Tanaka and Ida 1999) τaero 

  M 1/3 Mvrel 7 −1  2 × 10 e 2 M⊕ πrp2 ρgas vrel −1   1/2 ρgas r × years, (2.1) 2 × 10−9 gcm−3 1 AU



M gasr 2



cs vkep

4

−1 kep ,

(2.2)

where cs is sound velocity of disk gas and kep is Keplerian frequency (Ward 1989, 1993, Artymowitz 1993). Supposing the minimum mass disk model with gas surface density given by min gas = 1700(r/1 AU)−3/2 gcm−2 , 

τaero

2. GRAVITATIONAL GAS DRAG

2.1. Effective Drag Force

M M

τgrav

   e −1 M 1/3  1.0 × 10 0.01 M⊕ −1   13/4 gas r × years. min gas 1 AU      2 M −1 gas −1 r 3  0.5 × 10 years. min M⊕ gas 1 AU 6

(2.3)

(2.4)

The ratio of the two time scales is     5/4  τaero r e −1 M 4/3  2000 . τgrav 0.01 M⊕ 1 AU

(2.5)

When the mass is larger than ∼1025 g, gravitational gas drag is more effective than aerodynamical gas drag (Ward 1993, Tanaka and Ida 1999). So far, N-body simulation including gravitational gas drag has not been done, because aerodynamical gas drag is more effective for leftover planetesimals, and dynamical friction from the leftover planetesimals sufficiently damps the eccentricities of protoplanets unless final accretion stage is considered (e.g., Kokubo and Ida 1996, 1998, 2000). N-body simulation of final accretion stage (e.g., Chambers and Wetherill 1998, Agnor et al. 1999, Chambers 2001) also neglected gravitational gas drag. However, as shown in Eq. (2.4), the gravitational gas drag from even a significantly depleted gas disk is still important. Note that the effect of dynamical friction is formally incorporated into gravitational gas drag, because the dependence on planetary mass is the same. Actually, the damping time scale derived from Chandrasekhar’s dynamical friction formula (Eq. (5.4)) is reduced to τgrav given by Eq. (2.2) by replacing velocity dispersion and surface density of leftover planetesimals with sound velocity and surface density of a gas disk, as shown in Appendix.

46

KOMINAMI AND IDA

Also note that a functional form of τgrav would change when evkep  cs or gas solid , where solid is the smoothed-out surface density of solid component including protoplanets. When evkep  cs , τgrav is longer for larger e (Papaloizou and Larwood 2000, Artymowicz 1993). When gas solid , τgrav may be longer than that given by Eq. (2.4), unless transfered energy from protoplanets to gas is effectively dissipated. However, detailed dependence is yet to be clarified. Hence, we adopt the damping time scale rather than gas as a parameter, to make clear how a damping mechanism affects terrestrial planet formation. The dependence on τgrav of M and r would be important for final outcome of planet accretion (e-dependence is a future issue), so that the parameter τdamp we use in the present paper is defined by  τgrav = τdamp

M 0.2M⊕

−1 

r 1 AU

2 ,

(2.6)

where τdamp is the damping time scale for a body with M = 0.2M⊕ at 1 AU. A mass of 0.2M⊕ is used for mass scaling, according to initial condition of our simulations (see below). If the simple formula (2.4) is applied, τdamp is explicitly given by  τdamp = 2500

gas min gas

−1 years.

(2.7)

3. CALCULATION METHOD

3.1. Initial Condition We used the result of Kokubo and Ida (1998, 2000) to set up initial conditions. In most runs (standard set), 15 protoplanets with mass of 0.2M⊕ are placed with semimajor axis separation of 6–10 Hill radius, where Hill radius rH is defined by  rH =

f GD

v − v gas =− , τgrav

(2.8)

where v and v gas are the velocity of a protoplanet and the gas. We here assume the gas motion is a noninclined circular Keplerian motion. Equation (2.8) corresponds to Chadrasekhar’s dynamical friction formula. (Eq. (5.1)) with the above mentioned replacement (see Appendix), which may support the approximate inclusion of the gravitational gas drag as Eq. (2.8). As mentioned above, τdamp is a parameter. In many cases of our calculations, τdamp is kept constant throughout a run to clarify the effect of the damping. We also have done some runs with τdamp increasing with time, corresponding to the depletion of disk gas. Through orbit averaging following calculation by Adachi et al. (1976), the force f GD results in damping of the eccentricity e and the inclination i as 1 de −1 = −τgrav e dt 1 di = −(2τgrav )−1 . i dt

1/3

 r  0.007

M 0.2M⊕

1/3 r



solid M  0.1 min solid

3/2 

r 1 AU

3/4 

a 10rH

3/2 M⊕ ,

We confirmed the above relations by direct orbital integration of one protoplanet with f GD .

(3.2)

min where solid is a solid surface density in the minimum mass disk model, given by 7(r/1 AU)−3/2 gcm−2 . In order to simplify the problem and save computation time, we set all protoplanets with equal mass (0.2M⊕ ) in the standard set. Smoothed-out surface density of solid material in our model is



solid

a = 15 10rH

−1 

r 1 AU

−2

[gcm−2 ],

(3.3)

which has rather steeper radial gradient than that in the minimum mass disk model. Since protoplanets are spread out quickly by mutual perturbations, the effective surface density is more simmin ilar to solid = 7(r/1 AU)−3/2 gcm−2 . Calculation starts from the point when the orbital crossing starts. The result of Kokubo and Ida (2000) shows that the eccentricities of protoplanets are about 10−3 −10−2 , which means that the protoplanets have almost circular orbits. However, the protoplanets will eventually start orbital crossing by long-term mutual distant perturbations on a time scale (τcross ) depending on their orbital separation, mass (Chambers et al. 1996), initial eccentricities (Yoshinaga et al. 1999), and how much gas there

TABLE I Maximum and Minimum Semimajor Axis of the Initial Condition a(rH )

amin (AU)

amax (AU)

6 7 7.5 8 9 10

0.734 0.697 0.679 0.662 0.629 0.597

1.363 1.435 1.472 1.511 1.592 1.675

(2.9) (2.10)

(3.1)

and random angular distribution. The minimum and maximam values of the semimajor axis for the initial conditions are indicated in Table I. The result of Kokubo and Ida (2000) shows that the mass distribution of the protoplanets made through oligarchic growth has radial dependence,

2.2. Drag Formula We include the effects of the gravitational gas drag in N-body simulation by directly adding drag force ( f GD ) to equations of motion as

2M 3M

47

TERRESTRIAL PLANET FORMATION

is around the protoplanets (Iwasaki et al. 2001, 2002). The state we are investigating is after the onset of orbital crossing. It is rational to start the calculation with relatively high eccentricities (∼10−2
rH /a), supposing the state where the eccentricities has increased and orbital crossing is ready to start. We also did some calculations with smaller initial eccentricities. We found the results do not change except for existence of certain length (τcross ) of orbitally stable stage, because all bodies’ eccentricities increase almost simultaneously, not particular ones. Iwasaki et al. (2001, 2002) show that the eccentricities of the protoplanets secularly increase and orbital crossing starts if the damping time scale (τaero or τgrav ) is longer than 0.1–0.3τcross where τcross is the crossing time scale in the gas-free case (otherwise, orbital crossing is inhibited by the drag). In this case, the crossing time τcross is the same as that in the gas-free case. In the standard set of runs (24 runs), initial mass is 0.2M⊕ and the eccentricity damping timescale τdamp , which corresponds to the amount of gas, is constant throughout a run to see the effects of the drag clearly. We set the initial orbital separation as a. The list of the calculations of the standard set are shown in Table II. Three initial angular distribution types are identified as type a, b, and c. Angular distribution is given randomly. Each distribution is made of a different set of random number. The range of a is 6–10rH , and the range of the damping time scale is 105 –108 years. We also did the calculations of τdamp = 1015 years

and of those with no gas to see what happens in the case with extremely dissipated gas. We denote a run with τdamp = 10α years, a = βrH , and angular distribution k (k = a, b, c) as Runαβk . Each calculation has time span of 107 years because disk gas would not exist on a time scale much longer than 107 years. In addition to the standard set, we also performed three runs starting from 30 protoplanets with 0.1M⊕ (set B; Section 4.1) and five runs with time-dependent τdamp , which corresponds to disk gas depletion (set C; Section 4.3). In set B, initial solid is the same as a corresponding run in standard set. We do not see any significant difference in formed Earth-sized planets. Calculations with time-dependent τdamp in set C are more realistic. Here we assume exponentially increasing τdamp which corresponds to exponential decay of gas through Eq. (2.7) (for details, see Section 4.3). Initial conditions of distribution of protoplanets are the same as those in the standard set except for smaller initial eccentricities and inclinations. Earth-like planets are also formed with suitable choice of decay time scale of gas . The results of the runs in the standard set help us understand the results in set C. 3.2. Orbital Integration We integrate orbits with fourth-order Hermite scheme (Makino and Aarseth 1992) and hierarchial individual timestep

TABLE II List of Simulations of the Standard Set with Initial Conditions, and Mass and Eccentricity of Largest Planet Formed after 107 Years Simulation

τdamp (yr)

a(rH )

M1 (M⊕ )

e1

M2 (M⊕ )

e2

M1 ≥ 0.8M⊕

e1 ≤ 0.04

Run56a Run58a Run59a Run510a Run68a Run69a Run610a Run6.57.5a Run76a Run76b Run77a Run77b Run78a Run78b Run78c Run79a Run79b Run710a Run710b Run86a Run87a Run156a Run1510a RunGF9a

1 × 105 1 × 105 1 × 105 1 × 105 1 × 106 1 × 106 1 × 106 3 × 106 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 1 × 108 1 × 108 1 × 1015 1 × 1015 gas free

6 8 9 10 8 9 10 7.5 6 6 7 7 8 8 8 9 9 10 10 6 7 6 10 9

0.6 0.6 0.4 0.6 0.6 0.4 0.8 0.8 1.2 1.2 1.8 1.2 0.6 0.8 1.4 1.0 1.4 1.4 1.0 1.0 0.8 1.0 1.2 1.0

0.00077 0.00040 0.00015 0.00012 0.00010 0.00034 0.00013 0.00013 0.0047 0.0024 0.0053 0.052 0.087 0.036 0.013 0.014 0.0081 0.012 0.0098 0.13 0.11 0.057 0.049 0.18

0.4 0.4 0.4 0.4 0.6 0.4 0.6 0.8 0.8 0.8 0.8 0.4 0.6 0.6 0.4 0.6 0.6 0.4 0.6 0.8 0.6 0.8 0.6 0.8

0.00048 0.00043 0.00015 0.00011 0.00010 0.00034 0.00014 0.00013 0.0044 0.018 0.0037 0.058 0.087 0.034 0.023 0.016 0.020 0.052 0.025 0.14 0.099 0.076 0.098 0.19

No No No No No No Yes Yes Yes Yes Yes Yes No Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes

Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes No No Yes Yes Yes Yes Yes Yes No No No No No

Note. a,b,c labeled in the simulation number indicates the initial angular distribution type. Angular distribution is given randomly. Each distribution is made of different set of random number. M1 , M2 are the largest and the second largest of the final planets in each run. e1 , e2 are their eccentricities, respectively. If there are more than one second largest planets, the average is taken within the same mass. If the number of the largest mass is two or more, M1 and M2 are equal to that largest mass and e1 and e2 are the average eccentricity of those largest masses.

48

KOMINAMI AND IDA

(Makino 1991). The equation of motion of particle k is  GMj v k − v gas dv k G M =− rk − (r j − r k ) − . 3 3 dt |r k | |r j −r k | τgrav j=k

The physical radius of a protoplanet is determined by its mass and internal density as (3.4)

The first term is the gravity from the Sun. The second term is the mutual gravity between the protoplanets. And the last term is the gravitational drag from disk gas. We include the dependence of the drag force on the mass and the radial location of a ptotoplanet as Eq. (2.6). When protoplanets collide, perfect accretion is assumed. In other words, after the collision, a new body has total mass of the two colliding protoplanets, their coordinates are at the center of mass at the collision, and the velocity is determined so that total momentum is conserved before and after the collision.

 rp =

3 M 4π ρp

1/3 .

(3.5)

The internal density ρp is set to be 3 gcm−3 . 4. RESULT

4.1. Numerical Results The results of Run56a , Run76a , Run156a are shown in Figs. 1a, 1b, and 1c, respectively. Initial eccentricities and inclinations are ∼0.01. Figure 1 shows the time evolution of the orbits.

FIG. 1. Orbital evolution of the protoplanets. Dotted lines are apocenters and pericenters. Solid lines are semimajor axis. Thichness of the solid lines min . In (b), τ 7 is proportional to the mass. In (a), τdamp = 105 years, which is equivalent to a gas disk with ∼0.025 gas damp = 10 years, which corresponds to −4 min 15 gas ∼ 2.5 × 10 gas . In (c), τdamp = 10 years, which means the gas is almost completely depleted.

TERRESTRIAL PLANET FORMATION

FIG. 2. The mass, eccentricity, and semimajor axis of planets surviving at t = 107 years. The orbital elements are averaged over the last 5 × 105 years in (a) and (c), 6 × 105 years in (b). The area of the circles is proportional to the mass. The open circles are the result of our calculations, and the filled ones are the present terrestrial planets in our Solar System; Earth and Venus are at 1.0 AU and 0.72 AU. (a) the result of τdamp = 105 years, (b) τdamp = 107 years, and (c) τdamp = 1015 years.

Solid lines represent semimajor axis a. The dotted lines represent pericenter and apocenter; the amplitude between pericenter and apocenter is 2ea, which is proportional to eccentricity. The thickness of the solid lines represents the square root of planet mass. In Run156a (Fig. 1c), orbital crossing has not been finished even after 107 years, although some planets have grown to the size of Earth or Venus. Figure 2c shows the masses, eccentricities, and semimajor axes of surviving planets at t = 107 years. The orbital elements are averaged over the last 5 × 105 years. The area of the circles is proportional to the mass of a planet. The largest one at a = 0.65 AU and the second largest one at a = 1.05 AU have the mass of 1.2M⊕ and 0.8M⊕ . Although planets as massive as Earth or Venus are made, their eccentricities are too large compared to those of present Earth and Venus (0.017 and 0.0068) or their time-averaged free eccentricities ∼0.03; both of the largest and the second largest planets have e  0.06. The high eccentricities are the remnant of higher eccentricities during orbital crossing (e ∼ 0.1–0.2) pumped up by secular perturbations and close encounters between the protoplanets; also see analytical estimate in Section 4.2. Collisional damping is not enough to diminish e to ∼0.03. The amount of gas is not enough, either. In the other runs with τdamp = 108 and 1015 years, similar results are obtained regardless of the initial radial separation a, as shown in Fig. 3 and Table II. Figure 3 shows plots

49

similar to Fig. 2 for other representative runs. Table II shows data of all the runs. These results are comparable to the results of Chambers and Wetherill (1998) and Agnor et al. (1999). If we continue the integration on a time scale
108 years, the small bodies that are still orbitally crossing would be accreted by the large bodies and only two or three bodies would remain with still large eccentricities. Figure 1a is the result of Run56a . In this case, the gravitational gas drag damps the eccentricities before the protoplanets can grow by coagulations during the orbital crossing. The time span of the orbital crossing is significantly shortened by the damping of the eccentricities. Only a few coagulation events occur during first 6 × 105 years. As e is damped, the planets are isolated and τcross quickly increases (Yoshinaga et al. 1999). Further coagulation between the protoplanets does not occur at least until 107 years. As a result, numerous small planets remain. In other words, the planets do not grow enough. Figure 2a shows that most planets have mass smaller than 0.4M⊕ except for one body with 0.6M⊕ , although the eccentricities are small enough to be <0.01. Other runs with τdamp = 105 years but different initial radial separation, Run58a , Run59a , Run510a , resulted in the same kind of trend (Fig. 3 and Table II). On the other hand, when τdamp = 107 years, planets with about the size of Earth and small eccentricities (
0.03) are formed. As shown in Fig. 1b (result of Run76a ), orbital crossing lasts for about the first 3 × 106 years. As shown in Fig. 2b, the largest planet formed at a = 0.8 AU and the second largest one at

FIG. 3. The same plots as Fig. 2 for other representative runs with various τdamp , a and initial angular distributions. Runs in the top row are with strong gravitational drag (τdamp = 105 years), in the middle row, τdamp = 107 years, and in the bottom row, τdamp = 108 (left and middle) and 1015 (right) years. Present terrestrial planets are indicated as Solar System (figure with filled circles).

50

KOMINAMI AND IDA

a = 1.05 AU have the mass of 1.2M⊕ and 0.8M⊕ . After the isolation of the planets, the eccentricities of the surviving planets are damped gradually, spending the rest of the time of calculation (Fig. 1b). The final eccentricities of the largest and the second largest planets, averaged over the last 6 × 105 years, are 0.0058 and 0.0044, which are smaller than the eccentricities of the present Earth and Venus. There are also small mass and high eccentricity planets in innermost and outermost regions (Fig. 2b). Since the gravitational drag is weak for a small mass planet, the high eccentricities are not damped enough during the calculation time. These small bodies have experienced no collision. The innermost (outermost) one is scattered by another body to the region. Since there is no body in this region, the small scattered one tends to be isolated. Just after the scattering (t = (2 − 3) × 106 years in Fig. 1b), it is still crossing the orbit of the perturbing body, because the eccentricity is pumped up. However, the eccentricity is gradually damped by the gas and the scattered planet becomes isolated. This process is essentially the same as “orbital repulsion” described by Kokubo and Ida (1995). In many runs with τdamp = 107 years, we find similar small bodies that experienced no collision or only one collision in innermost and outermost regions. We also did three calculations starting with smaller mass planets (set B), which are more consistent with Eq. (3.2). These three simulations are started with 30 protoplanets with mass of 0.1M⊕ and different initial orbital separation (a). One result is shown in Figs. 4 and 5. In this simulation, the inner and outer edge, and τdamp is the same as those in Run77a . a is 4.2rH . As shown in Figs. 4 and 5, Earth-sized planets with low eccentricities are formed (the largest mass is 1.1M⊕ with e ∼ 0.008, and the second largest mass is 0.8M⊕ with e ∼ 0.02). In this case, bodies in innermost and outermost regions tend to be smaller

FIG. 4. Orbital evolution of the protoplanets. Here we start with 30 protoplanets with 0.1M⊕ . a = 4.2rH and τdamp = 107 years. The inner edge and the outer edge of the disk is 0.697 AU and 1.435 AU, respectively.

FIG. 5.

The final planets of Fig. 4.

than those formed in the standard set. These small bodies tend to have larger eccentricities because of weaker gravitational drag. These planets could be analogous to Mercury and Mars in our Solar System, although the number of remaining small bodies is rather large (also see discussion below). Note that Eq. (3.2) shows masses consistent with Mercury and Mars, in the inner and outer regions, respectively. In order for the orbital repulsion mechanism to work, the inwardly (outwardly) scattered planet must not meet another large body in the region into which the planet is scattered. In our calculations, bodies are placed in limited regions so that the scattered planets do not meet other bodies. This “edge effect” makes the orbital repulsion effective. Alexander and Agnor (1998) and Agnor et al. (1999) also found a similar trend because their disks are also truncated at 0.5 AU and ∼1.5 AU. For a Mercury-like planet, this edge effect may be reasonable because the mass of a protoplanet resulted by oligarchic growth decreases with decrease in r as shown in Eq. (3.2). To make a Mars-like planet, bodies beyond ∼1.5 AU must be cleared before the orbital crossing stage. If the orbital crossing is triggered by formation of giant planets, which are formed before the orbital crossing stage the asteroid belt can be quickly cleared by the perturbation of the giant planets (Chambers and Wetherill 2001). However, clearing bodies between 1.5 and 2.0 AU may take longer time scale (
108 years; Chambers and Wetherill 2001). Hence, in order to clarify the formation of Mars, it would be necessary to consider formation of terrestrial planets, Jupiter and Saturn, and the evolution of the asteroid belts simultaneously. These trends, caused by different τdamp , are also present in the results from various initial radial separation a, and different initial angular distribution as shown in Fig. 3 and Table II. Figure 6 shows mass and eccentricity of the largest and the second largest bodies averaged over all runs with the same τdamp (the plot of τdamp = 1015 years also includes the run without gas, RunGF9a ). Since the number of runs used for averaging the masses and the eccentricities of runs of τdamp = 108 is only 2,

TERRESTRIAL PLANET FORMATION

51

cretion will not occur without excitation of eccentricities. In reality, τdamp increases with time as the disk gas is depleted (Section 4.3). When τdamp becomes
0.3τcross , eccentricities are pumped up to start orbital crossing (Iwasaki et al. 2002). If the stage of τdamp 107 years lasts long enough (
107 years), planets that resemble the result of Run76a will be formed. If such a stage is short and the gas dissipates quickly, the results would be similar to the result of Run156a . Note that direct conversion of τdamp ∼ 107 years with Eq. (2.7) min is gas ∼ gas /4000. Although Eq. (2.7) may include some uncertainties in the case with gas solid , this conversion indicates that a residual disk with very small fraction of the minimum mass disk is enough to diminish the eccentricities of Earthsized planets. This is consistent with estimate by Agnor and Ward (2002) that analytically studied damping of eccentricities of planets by gravitational drag in the presence of distant perturbation between planets after the orbital crossing stage (they did not consider orbital crossing and collision). FIG. 6. Average mass and eccentricity of the largest and the second largest planets formed at the end of the calculation (1 × 107 years) of the standard set. Solid line with filled circles is for the mass and the dotted line with open circles is for the eccentricity. The average is taken over the runs of each damping time scale (The plot of τdamp = 1015 is the average of the runs with τdamp = 1015 and the run without gas (RunGF9a )). When τdamp = 107 years, the average mass reaches to ∼M⊕ , while the eccentricity is still below 0.03.

and τdamp = 1015 is only 3, the dip of the averaged mass of τdamp = 108 can be said as not relevant. It clearly indicates that calculations with τdamp < 106 years only form small planets, and with τdamp > 108 years, the final eccentricities are too large. When τdamp = 107 years, Earth-like planets with M ∼ M⊕ and e
0.03 are formed. Note that these runs do not reproduce the present terrestrial planets completely in the point that extra small bodies often remain in the middle regions. The average number of remaining planets is five or six in the cases with τdamp = 107 years. The calculations with each of the initial protoplanet mass being smaller (0.1M⊕ ) show similar results. Further accretion will not occur in Run76a after the calculation time (107 years) for following reason. The minimum planet separation is ∼19rH . The time (τcross ) for orbital instability to occur when orbital separation is 19rH is ∼1014 years if the results by Chambers et al. (1996) is extrapolated. This time scale is longer than the age of Solar System. The planets resulted in other Run7βγ also seem to be orbitally stable, although some small bodies could be orbitally excited to collide with large bodies or ejected on a longer time scale. Inclusion of Jupiter and Saturn perturbations might reduce number of final planets to three to five (see discussion in Section 5). The final planet separation in Run56a ranges from 8 to 13rH . This corresponds to τcross ∼106 years (Chambers et al. 1996). However, this time scale is longer than the eccentricity damping time scale, which is 105 years in this case. Eccentricities of the planets will be damped by the surrounding gas rather than pumped up by the planets’ mutual gravitational interaction. Ac-

4.2. Analytical Estimate Why are planets about the size of the Earth and small eccentricities successfully formed for τdamp ∼ 107 years? The results by Chambers and Wetherill (1998) and Agnor et al. (1999) showed that the orbital crossing lasted for 107 –108 years to make Earth-sized planets. From the observation, it is inferred that gas is likely to remain in the disk only for ∼107 years (e.g., Strom et al. 1993, Zuckerman et al. 1995). These might suggest that the eccentricities of the surviving planets cannot be damped. However, in their simulations, planets as large as Earth are formed on a time scale of ∼106 –107 years. The remaining 107 –108 years are required for small bodies with high eccentricities to collide with the Earth-sized planets. If the small amount of gas exists, it damps the high eccentricities of the small bodies to isolate the bodies and/or to accelerate their accretion by the Earth-sized bodies. Then the orbital crossing time would be shorter than the lifetime of the disk gas. Chambers (2001) that included dynamical friction also reported the shortening of orbital crossing period. Although orbital crossing time is shortened, enough coagulation events (10 events in Fig. 1b) occurred to form planets as massive as Earth or Venus. After formation of Earth-sized planets, the decrease of the eccentricities can happen within the disk lifetime. This analysis explains why similar results would be obtained for τdamp ∼ 107 years. Furthermore, the observationally inferred disk lifetime ∼107 years is a time scale for depletion of a major fraction of a disk. The required small fraction of disk gas can persist on a longer time scale. If gas decreases exponentially, the decay of gas by a factor 10−4 –10−3 takes seven to nine e-folding times. Hence, the amount of gas required for τdamp ∼ 107 years may persist for several times 107 years. We briefly comment on the orbital separation and the mass of planets after orbital crossing. The following argument explains why a few Earth-sized planets with orbital separation 0.2–0.5 AU are usually formed for τdamp >106 years. The width of the feeding

52

KOMINAMI AND IDA

zone a of a planet is written as

the estimated eccentricity would be

a ∼ 2ae.

(4.1)

e ∼ 0.25.

(4.5)

Thus, the mass and the separation between the planets are estimated as

The mass of the planet M is M ∼ 2πaa solid ,

(4.2)

where solid is the surface density that is obtained by smoothing out the protoplanets. solid is determined by initial conditions of the protoplanets. We need to estimate the eccentricity e to get a and M. Orbital evolution is devided into three stages in our calculations that succeeded to form Earth-like planets. The excitation and damping mechanism of the eccentricities depend on the stages. In the first stage, orbital crossing does not occur yet. In the second stage, the orbital crossing starts and accretion of protoplanets takes place. The third stage is the isolation of the planets. After the accretion, surrounding gas damps the eccentricities and surviving planets are isolated with each other. In the first stage (also in the third stage), since there is no collision, the only damping mechanism is the gravitational drag. Once the instability is triggered, the second stage begins; the accretion starts. The pumping is caused by relatively close encounters. Since the bodies collide with each other, collisional damping also diminishes the eccentricities. Collisional damping is more effective than gravitational interaction in the second stage. The collisional damping time scale is estimated by mean collision time to be 105 –106 years, while the time scale of gravitational drag in the case we are concerned with is ∼107 years. The damping force is dominated by collisional damping in this stage. Accretion in the second stage determines a and M of the surviving planets. The equilibrium velocity dispersion (equivalently, the eccentricity e) is determined by comparing the collisional cross section and the scattering cross section (e.g., Safronov 1969, Aarseth et al. 1993) as the velocity with Safronov parameter ∼1,  veq ∼

GM rp

1/2 ,

(4.3)

where rp is the physical radius the planet. We recognize that secular perturbations (accumulation of successive noncrossing perturbations) may be more important than close encounters (Chambers and Wetherill 1998); however, we here consider close encounters for simplicity. (Note that secular pertubations here are different from secular resonance discussed in later section.) The eccentricity is the ratio of veq and the Kepler velocity vkep as veq e∼ ∼ vkep



M a M r p

1/2 .

(4.4)

When M ∼ M⊕ and the density of the planets ρp is 3 gcm−3 ,

a ∼ 2ea ∼ 0.5 AU M ∼ 2πaa solid ∼ M⊕ .

(4.6) (4.7)

The numerical values are estimated with density closer to ∼10 gcm−2 , which is comparable to the solid surface density of the minimum mass disk model. Although we start with higher solid (Eq. (3.3)), rapid diffusion reduces solid to ∼10 gcm−2 . If the collision events occur before e is fully pumped, a and M are slightly smaller. In our calculations, we actually found e = 0.1–0.2. We do not find any orbital crossing after first orbital crossing is finished. Right after the second stage, the eccentricities are still high and the protoplanets are not sufficiently isolated, then τcross is relatively small. The gas drag gradually damps the eccentricities if the planets tentatively become isolated by the collision and close encounters are no longer available to pump up the eccentricities. The decrease of the eccentricities results in significant lengthening of τcross (Yoshinaga et al. 1999), which accounts for the inhibitation of further orbital crossing in our simulation. 4.3. The Case of the Dissipating Disk We showed that if the orbital crossing begins when the disk has gas of τdamp ∼ 107 years and if such gas stays in the disk for ∼107 years, Earth-sized planets with low eccentricities can be formed. Direct conversion of τdamp ∼ 107 years to gas min surface density is ∼ a few ×10−4 gas , which is well below −2 solid (∼10 gas ). Although Eq. (2.7) may not be correct enough in the case of gas solid , this means that a residual disk with only small fraction of initial mass is enough to diminish the eccentricities of the Earth-sized planets. Since in the minimum mass disk model, solid ∼ 10−2 gas , it may be also possible to make Earth-like planets by including dynamical friction force due to the leftover planetesimals with total mass as massive as that of protoplanets, which may explain the result in Chambers (2001) that the eccentricities of the finally formed planets are smaller than those in the case without dynamical friction. However, the eccentricities are not small enough to be comparative to the eccentricities of Earth and Venus. This may be because the leftover planetesimal sizes he set are so large that viscous stirring inhibits efficient dynamical friction (e.g., Ida 1990) or gap opening around protoplanets in a planetesimal disk weakens dynamical friction (Ida and Makino 1993, Papaloizou and Larwood 2000). More smaller planetesimals may be needed for sufficient dynamical friction. We will need further investigation on the size and spatial distribution of the leftover planetesimals.

53

TERRESTRIAL PLANET FORMATION

So far, we have fixed τdamp throughout a run. In reality, τdamp increases with time as gas decreases. We have also carried out five simulations with time-dependent τdamp as  τdamp (t) = τdamp (0) exp

t tgas

 ,

(4.8)

where tgas is depleting time scale of the gas. If we use Eq. (2.7), even in the stages of gas solid , Eq. (4.8) corresponds to gas surface density gas depleting exponentially with time as gas (t) = gas (0) exp(−t/tgas ). tgas may correspond to observationally inferred disk depletion time ∼107 years. The result of one of the five runs is shown in Figs. 7 and 8. The initial mass of each protoplanet and the separations are the same as Run78a . The eccentricities and inclinations are given randomly with the magnitude of ∼0.0001, which is small enough to prevent orbital crossing initially. In this calculation, we used tgas = 3 × 106 years and τdamp (0) = 2.5 × 105 years. The latter min min corresponds to gas (0) = 0.01 gas . Because 0.01 gas is still large enough to suppress the orbital crossing as shown below, the calculation being done until the gas amount diminishes to this value would merely show quite stable orbital evolution; e and i min are kept small. If we start the calculation from gas (0) = gas , the time is shifted by 1.4 × 107 years. Figure 7 shows the orbital evolution. Like Fig. 1, the pericenter, apocenter, and the semimajor axis are plotted versus time. During the first 4 × 106 years, the eccentricities are kept low by the gravitational drag. If there is no gas, orbital crossing begins at t  1 × 106 years (τcross  1 × 106 years). The gas is suppressing the

FIG. 7. Orbital evolution of the protoplanets. The same plots as Figs. 1a, 1b, 1c. Except that τdamp is increasing with time as τdamp (t) = τdamp (0) exp (t/tgas ), corresponding to the depletion of disk gas. In this calculation, τdamp (0) = 2.5 × min ), and t 6 105 years ( gas (0) = 0.01 gas gas = 3 × 10 years. The initial number of the protoplanets, their individual mass, and the separation are the same as Run 78a . The initial eccentricities and inclinations are about 10−3 .

FIG. 8.

The final planets of Fig. 7.

eccentricities during the first 4 × 106 years. At t  4 × 106 years, τdamp increases to 1.0 × 106 years ( gas decreases to 2.6 × min ), and the drag becomes weak enough to allow or10−3 gas bital crossing to start. The result is consistent with the condition to start orbital crossing, τdamp a
0.3τcross , obtained by Iwasaki et al. (2002), although they calculated only constant τdamp cases. Orbital crossing lasts for 3.5 × 106 years. The eccentricities are damped after the orbital crossing stage ends until the end of simulation at 1 × 107 years (τdamp  7 × 106 years, which cormin years). The planets’ mass responds to gas = 3.6 × 10−4 gas 7 and eccentricities at t = 1 × 10 years are shown in Fig. 8. Like Fig. 2, filled circles represent the present terrestrial planets, and the open circles are the simulation results. The largest planet has the mass of 0.8M⊕ and the eccentricity of 0.0017 at a = 0.94 AU. The second largest planet’s mass is 0.6M⊕ , the eccentricity is 0.0024 at 0.72 AU. Hence, the planets analogous to Earth and Venus are formed, although Mercury and Mars analogues are not produced in this case. Whether the final planets are stable or not depends on the separations. Two small planets at ∼1.7 AU have orbital separation of 0.37rH and are likely to merge. Among the other orbital separations, the smallest was 20rH , which will take at least ∼1014 years for orbital crossing to begin (Chambers et al. 1996). Since this time scale is much longer than the age of the Solar System, the other larger planets are orbitally stable. Other simulations have different initial τdamp and initial orbital separations. Four of the simulations are carried out with the same tgas min and the same initial gas surface density ( gas (0) = 0.01 gas ) but min . with different a, and one of them is with gas (0) = 0.03 gas In four out of five runs, planets with ∼M⊕ and low eccentricity (
0.03) are formed. If the initial orbital separation is too wide, τcross becomes so long that the amount of gas there when the orbital crossing begins is not enough to suppress the eccentricities of the planets. However, if tgas is longer, this problem may be overcome. We need more runs with different tgas and different functional forms of time dependence of τdamp . These four results that succeeded to form Earth-sized low eccentricity planets have

54

KOMINAMI AND IDA

similar final separation width; bodies with mass of ∼M⊕ tend to be separated with a distance of
20rH . This shows that even if the gas completely decays, the planets are stable. Also, the number of the final planets tend to be more than number of the present terrestrial planets. Those small planets near the Earthsized planets may coagulate with the large planets if there had been giant planets to pump up the eccentricities. The effect of giant planets are the issue of future studies. 5. CONCLUSION AND DISCUSSION

We have investigated the final accretion stage of terrestrial planet formation, through N-body simulation, including the eccentricity damping caused by tidal interaction with a remnant gas disk (“gravitational gas drag”). The final stage would be dominated by giant impacts among the protoplanets formed by runaway and oligarchic growth. N-body simulations by Chambers and Wetherill (1998) and Agnor et al. (1999) showed that the giant impacts result in planets of about the size of Earth, while the eccentricities of the planets tend to be significantly higher (∼0.1) than time-averaged eccentricities of Earth and Venus (∼0.03). Collisional damping is not strong enough. In order to diminish the eccentricities, some kind of drag force that is effective even after collision events are over and protoplanets are isolated from each other should be incorporated. The inclusion of strong drag force may prevent the orbital crossing of the protoplanets in the first place, not allowing the planets to grow sufficiently. However, we showed that the inclusion of “gravitational gas drag,” damping due to tidal interaction with a protoplanetary gas disk (Ward 1989, 1993, Artymovicz 1993) can successfully reproduce Earth-sized planets with eccentricities
0.03 if we consider the drag exerted by a mostly depleted gas disk. Aerodynamical gas drag is too weak in this case. Dynamical friction by leftover planetesimals may play a similar role (Section 4.2; Chambers 2001). In most runs (standard set), we start from 15 protoplanets with equal mass (0.2M⊕ ) and several Hill radius orbital separations and integrate the orbits for 107 years. The gravitational gas drag force is directly included in the equation of motion for orbital integration of the protoplanets. In the standard set, the damping time scale τdamp is fixed throughout the calculation. When τdamp < 106 years, since the drag force is too strong, it damps the eccentricities so quickly that the orbital crossing does not last long enough to allow the planets to grow to the size of Earth. When τdamp > 108 years, orbital crossing lasts throughout the calculation time and Earth-sized planets can be formed. However, the drag is too weak to diminish the eccentricities of the planets. The eccentricities are as large as ∼0.1. When 106 years < τdamp < 108 years, the orbital crossing stage lasts for ∼ a few ×106 years. Some planets grow to the size of Earth during orbital crossing period. The eccentricities are damped for the rest of the calculation time, resulting in the final values
0.03. Usually, the largest planets form around 1 AU. Also some small planets that have relatively high eccentricities

are formed at innermost (∼0.5 AU) and outermost (∼1.5 AU) part of the calculation area, which is partly due to the “edge effect” that initial distribution of protoplanets are limited to the region between ∼0.5 AU and ∼1.5 AU, and also reflects the mass dependence of the drag force. These small bodies have experienced no collision or only one collision. They are scattered by another body to the innermost or outermost regions and are isolated by the eccentricity damping (similar to “orbital repulsion” mechanism; Kokubo and Ida 1995). These planets may be analogous to Mars and Mercury in our Solar System. τdamp ∼ 107 years would correspond to the gas surface density min ), ( gas ) of 10−4 –10−3 times the minimum mass disk model ( gas if we use Eq. (2.7). If remnant gas with only small fraction of surface density in the minimum mass model exists in the disk during the orbital crossing, a planetary system similar to terrestrial planets can be formed. In reality, τdamp increases with time according to depletion of disk gas. We did five calculations with time-dependent τdamp as τdamp (t) = τdamp (0) exp(t/tgas ) with tgas = 3 × 106 years. When the protoplanets are formed through runaway and oligarchic growth on a time scale ∼105 –106 years, it is possible to suppose the gas still exists in the disk with the full amount of the minimum mass model. This amount of gas prevents the eccentricities of the protoplanets from being pumped up. Eventually, the gas is depleted to the amount that the eccentricities of the protoplanets can be elevated by mutual gravitational interaction and/or the secular resonance from the giant gaseous planets. The excitation time scale by the secular resonance is τres ∼ 106 years (Ito and Tanikawa 1999). The result by Iwasaki et al. (2002) predicts that orbital crossing starts when τdamp ∼ 0.3 min(τcross , τres ) (where τcross is the crossing time without the resonance). This condition min would correspond to gas ∼ 10−3 −10−2 gas . As shown in many runs in set C, a few planets with mass of ∼M⊕ are formed. The gas that still exists in the disk with min gas ∼ 10−4 −10−3 gas diminishes the eccentricities of the planets until the gas is depleted completely from the disk. The final eccentricities of the Earth-sized planets may be
0.03. The diminishing of the eccentricities makes τcross lengthen even more (Yoshinaga et al. 1999), making the planetary orbits become more stable. Although we need to carry out more detailed calculations with more realistic initial distribution of the protoplanets as well as with time-dependent damping force according to depletion of disk gas, the model of terrestrial planet formation with gravitational drag by disk gas is promising. There are some cosmochemical evidences that the Earth– Moon impact and Earth’s core formation occurred ∼2.5–6.5 × 107 years after the start of Solar System (Halliday et al. 2000). min . Such Orbital crossing would start when gas ∼ 10−3 −10−2 gas decay requires several e-folding time if we assume exponential decay. Since observationally inferred disk lifetime ∼107 years may correspond to e-folding time, our results suggest formation time scale of the Earth would be several tens to a hundred million years, which may be consistent with the cosmochemical evidences.

55

TERRESTRIAL PLANET FORMATION

In our calculation, the effect of the giant gas planets is not included. The secular resonance can highly increase the eccentricities of the protoplantes. If one of the protoplanets gets caught around the resonant points, the eccentricities would be pumped up, causing others to increase the eccentricities gradually (Ito and Tanikawa 1999). However, the resonances are only effective in the narrow range. As the planets grow, their distribution becomes sparse and the possibility that one of the planets would get trapped in them becomes low. So, the effect of the resonances on the grown planets can be considered as negligible. In fact, the secular resonances exist in the terrestrial planet region now, but they do not affect the orbits of the present terrestrial planets. However, the secular resonances could clean excessive small bodies that are often found in the results of our calculations; we have to carry out calculations including Jupiter and Saturn in future study. The eccentricities of planets are damped by disk-planet interaction. The same interaction also damps semimajor axes (Goldreich and Tremaine 1980, Ward 1986, 1989, 1997. If the gas is considered in the calculation, such an effect should be considered as well as the damping of the eccentricities. However, in the stage of orbital crossing of the protoplanets when min gas ∼10−4 –10−3 gas , the time scale for the migration to occur is longer than the calculation time we are considering in this study. It is acceptable to ignore the migration effect in our calculation. When the protoplanets are just formed, however, the amount of gas present may be enough to make the bodies migrate by significant distance. The distribution of the protoplanets may vary. This issue is left for future study. APPENDIX: DYNAMICAL FRICTION TIME SCALE Consider a planet with mass M and velocity v gravitationally interacting with a swarm of planetesimals with spatial density ρ, velocity dispersion vdis , and bulk velocity v bulk (vbulk vdis ). If v vdis , Chandrasekhar’s dynamical friction formula   v − v bulk 8 π dv , − ln G 2 ρ M 3 dt 3 2 vdis

(5.1)

where ln  expresses the effect of distant perturbations (e.g., Chandrasekhar 1943, Binney and Tremaine 1987). For the system of a planet and planetesimals, ln  can be approximated to ln  ∼ 3 − 5

(5.2)

(Stewart and Ida 2000). Since ρ ∼ pl /(2vdis ) ( pl : the surface density of the planetesimals) and 2 r/M , Eq. (5.1) reads as G = vkep dv v − v bulk , =− dt τDF

(5.3)

where the damping time scale (τDF ) of the dynamical friction is τDF =

3 2π ln 



M M



M solid r 2



vdis vkep

4

−1 kep .

(5.4)

Note that Eq. (5.4) is reduced to Eq. (2.2) except for a numerical factor by relacing vdis and pl with cs and gas . It may suggest that dynamical friction and gravitational gas drag are essentially the same dynamical process.

ACKNOWLEDGMENTS We thank John Chambers and an anonymous refree for useful comments on the first version manuscript. J. K. appreciates the help by Hiroshi Daisaka on simulation code. We also thank Kazunori Iwasaki and Hidekazu Tanaka for sharing their results before publishing, and Kiyoshi Nakazawa for encouragements.

REFERENCES Aarseth, S. J., D. N. C. Lin, and P. L. Palmer 1993. Evolution of planetesimals. II. Numerical simulations. Astrophys. J. 403, 351. Adachi, I., C. Hayashi, and K. Nakazawa 1976. The gas drag effect on the elliptical motion of a solid body in the primordial solar nebula. Prog. Theor. Phys. 56, 1756–1771. Agnor, C. B., and W. R. Ward 2002. Damping of terrestrial planet eccentricities by density wave interactions with a remnant gas disk. Astrophys. J., in press. Agnor, C. B., R. M. Canup, and H. F. Levison 1999. On the character and consequences of large impacts in the late stage of terrestrial planet formation. Icarus 142, 219–237. Alexander, S. G., and C. B. Agnor 1998. N-Body simulations of late stage planetary formation with a simple fragmentation model. Icarus 132, 113–124. Artymowicz, P. 1993. Disk-satellite interaction via density wave and the eccentricity evolution of bodies embedded in disks. Astron. J. 419, 166–180. Barge, P., and R. Pellat 1991. Mass spectrum and velocity dispersion during planetesimal accumullation. I. Accretion. Icarus 93, 270–287. Brouwer, D., and A. J. J. van Woerkom 1950. The secular variations of the orbital elements of the principal planets. Astron. Pap. Am. Ephemeris. Naut. Alm. 13, 81–107. Chambers, J. E. 2001. Making more terrestrial planets. Icarus 152, 205–224. Chambers, J. E., and G. W. Wetherill 1998. Making the terrestrial planets: N-body integrations of planetary enbryos in three dimensions. Icarus 136, 304–327. Chambers, J. E., and G. W. Wetherill 2001. Planets in asteroid belt. Meteoritics & Planetary Science 36, 381–399. Chambers, J. E., G. W. Wetherill, and A. P. Boss 1996. The stability of multi-planet systems. Icarus 119, 261–268. Goldreich, P., and S. Tremaine 1980. Disk-satellite interactions. Astrophys. J. 241, 425–441. Greenberg, R., J. Wacker, W. K. Hartmann, and C. R. Chapman 1978. Planetesimals to planets: Numerical simulations and collisional evolution. Icarus 35, 1–26. Halliday, A. N., D.-C. Lee, and S. B. Jacobsen 2000. Tungsten isotopes, the timing of metal–silicate fractionation, and the origin of the Earth and Moon. In Origin of the Earth and Moon (R. M. Canup and K. Rightor, Eds.), pp. 45–62. Univ. of Arizona Press, Tucson. Hayashi, C. 1981. Structure of the solar nebula, growth and decay of magnetic fields and effects of magnetic and turbulent viscosities on the nebula. Prog. Theor. Phys. Suppl. 70, 35–53. Ida, S., and J. Makino 1992a. N-body simulation of gravitational interaction between planetesimals and a protoplanet I. Velocity distribution of planetesimals. Icarus 96, 107–120. Ida, S., and J. Makino 1992b. N-body simulation of gravitational interaction between planetesimals and a protoplanet II. Dynamical friction. Icarus 98, 28–37.

56

KOMINAMI AND IDA

Inaba, S., H. Tanaka, K. Nakazawa, G. W. Wetherill, and E. Kokubo 2001. High-accuracy statistical simulation of planetary accretion: II. Comparison with N-body simulation. Icarus 149, 235–250. Ito, T., and K. Tanikawa 1999. Stability and instability of the terrestrial protoplanet system and their possible roles in the final stage of planet formation. Icarus 139, 336–349. Iwasaki, K., H. Tanaka, K. Nakazawa, and H. Emori 2001. The gas drag effect on the orbital instability of a protoplanet system. Publ. Astron. Soc. Japan 53, 321–329. Iwasaki, K., H. Emori, K. Nakazawa, and H. Tanaka 2002. Orbital stability of a protoplanet system under the drag force proportional to the random velocity. Publ. Astron. Soc. Japan, submitted. Kokubo, E., and S. Ida 1995. Orbital evolution of protoplanets embedded in a swarm of planetesimals. Icarus 114, 247–257. Kokubo, E., and S. Ida 1996. On runaway growth of planetesimals. Icarus 123, 180–191. Kokubo, E., and S. Ida 1998. Oligarchic growth of protoplanets. Icarus 131, 171–178. Kokubo, E., and S. Ida 2000. Formation of protoplanets from planetesimals in the solar nebula. Icarus 143, 15–27. Makino, J. 1991. A modified Aarseth code for GRAPE and vector processors. Publ. Astron. Soc. Japan 43, 859–876. Makino, J., and S. J. Aarseth 1992. On a hermite integrator with Ahmad-Cohen scheme for gravitational many-body problems. Publ. Astron. Soc. Tpn. 44, 141–151. Nagasawa, M., H. Tanaka, and S. Ida 2000. Orbital evolution of asteroids during depletion of solar nebula. Astron. J. 119, 1480–1497. Papaloizou, J. C. B., and J. D. Larwood 2000. On the orbital evolution and growth of protoplanets embedded in a gaseous disk. MNRAS 315, 823–833. Safronov, V. S. 1969. Evolution of the protoplanetary cloud and Formation of the Earth and planets. Nauka, Moscow. [In Russian; English translation: NASA TTF-677, 1972]. Spaute, D., S. Weidenschilling, D. R. Davis, and F. Marzari 1991. Accretional evolution of a planetesimal swarm: I. A new simulation. Icarus 92, 147–164.

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 timescales 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. Tanaka, H., and S. Ida 1999. Growth of migrating protoplanet. Icarus 139, 350–366. Tanaka, H., T. Takeuchi, and W. R. Ward 2002. Three-dimensional interaction between a planet and an isothermal gaseous disk. I. Corotation and Lindblad resonances and planet migration. Astrophys. J. 565, 1257–1274. Ward, W. R. 1986. Density waves in the solar nebula: Differential Lindblad torque. Icarus 67, 164–180. Ward, W. R. 1989. On the rapid formation of giant planet cores. Astrophys. J. 345, L99–L102. Ward, W. R. 1993. Density wave in the solar nebula: Planetesimal velocities. Icarus 106, 274–287. Ward, W. R. 1997. Survival of planetary systems. Astrophys. J. 482, L211. Weidenschilling, S. J. 1977. The distribution of mass in the planetary system and solar nebula. Astrophys. Space Sci. 51, 153–158. Weidenschilling, S. J., D. Spaute, D. R. Davis, F. Marzari, and K. Ohtsuki 1997. Accretional evolution of a planetesimal swarm. 2. The terrestrial zone. Icarus 128, 429–455. Wetherill, G. W. 1992. An alternative model for the formation of the asteroid belt. Icarus 100, 307–325. Wetherill, G. W., and S. R. Stewart 1989. Accumulation of a swarm of small planetesimals. Icarus 77, 330–357. Wetherill, G. W., and S. R. Stewart 1993. Formation of planetary embryos— Effects of fragmentation, low relative velocity, and independent variation of eccentricity and inclination. Icarus 106, 190. Yoshinaga, K., E. Kokubo, and J. Makino 1999. The stability of protoplanet systems. Icarus 139, 328–335. Zuckerman, B., T. Forveille, and J. H. Kastner 1995. Inhibitation of giant-planet formation by rapid gas depletion around young stars. Nature 373, 494–496.