- Email: [email protected]
Numerical simulations of the accretion of Uranus and Neptune
\ PERGAMON
Planetary and Space Science 36 "0888# 480Ð594
Numerical simulations of the accretion of Uranus and Neptune Adrian Bruninia\ \ Julio A[ Fernandezb a
Facultad de Ciencias Astronomics y Geof(sicas\ Universidad Nacional de La Plata "0899# La Plata\ Ar`entina b Departamento de Astronom(a\ Facultad de Ciencias\ I`ua esq[ Mataojo\ 00399 Montevideo\ Uru`uay Received 2 March 0887^ received in revised form 4 November 0887^ accepted 11 November 0887
Abstract New numerical simulations of the accretion process of the outer planets are carried out using an N!body code[ In 10 out of 29 computer runs 1 giant planets form beyond Saturn\ whereas 0 forms in 7 cases and 2 in 0 case[ The time scale of formation is a few 096 years[ This rather short time scale seems to be consistent with the content of non!negligible amounts of hydrogen and helium in Uranus and Neptune\ which could be accreted before the protoplanetary nebula dispersed[ A marked radial drift of the accreting outer planets is found in our simulations in agreement with previous results[ The results also seem to con_rm that this radial migration is an e}ective mechanism to produce orbital con_gurations near mean motion resonances[ Less than 49) of the solid material originally present in the system as planetesimals contributes to the formed protoplanets\ the rest being mostly ejected by Jupiter and Saturn[ The phase of strong scattering of bodies from the UranusÐNeptune region to the inner solar system lasts for a period of ½3×096 years[ Around 09Ð19 Earth masses reach the region of the terrestrial planets\ although the very short time scale for the massive in~ux of material "perhaps shorter than the time scale of Earth formation# might pose some di.culties to the theory that comets were major contributors of water to the Earth[ Þ 0888 Elsevier Science Ltd[ All rights reserved[
0[ Introduction The accretion process of the outer planets was one of the most collective phenomena during the formation of the solar system[ From the region of the terrestrial planets up to its outermost boundary\ the Oort cloud\ the solar system has been strongly in~uenced by this process "Fer! nandez and Gallardo\ 0885#[ Much progress has been done in our understanding of the formation of the outer planets and scattering of the residual mass from the early theoretical research of Safronov "0861#[ Nevertheless\ there are several aspects not yet understood that deserve further study[ The numerical simulations of the macro!accretion pro! cess of Uranus and Neptune carried out by Fernandez and Ip "0873\ 0885#\ have properly accounted for some important features observed in the outer planetary region such as the presence of two massive planets beyond Saturn[ One interesting and unexpected result found by Fernandez and Ip "0873# was that the orbit of Neptune\ and to a lesser extent those of Uranus and Saturn experi! enced an outward drift in the course of their accretion due to the exchange of angular momentum with the inter! Corresponding author[ E!mail] abruniniÝfcaglp[fcaglp[unlp[edu[ar
acting planetesimals[ The angular momentum gained by the orbital expansion of these planets was compensated by a small drift inwards of the massive Jupiter[ The planetary migration may have played an impor! tant role during the late accretion phase of the outer solar system\ being recently claimed as responsible for the presence of Pluto in its eccentric and inclined orbit at 2 ] 1 resonance with Neptune "Malhotra\ 0884#\ and of the global quasi!resonant structure of the outer solar system "Fernandez and Ip\ 0885#[ The numerical models developed by Fernandez and Ip "0873\ 0885# also suggested that the formation of Uranus and Neptune required an initial mass in their accretion zones about 1Ð2 times larger than their combined masses[ The unaccreted solid material was lost to the inner solar system or to the interstellar space[ Therefore\ the accretion of Uranus and Neptune seems to have been very ine.cient during the latest stages\ which may be explained as due to the increasing probability of ejection of interacting planetesimals with protoplanets whose masses exceed a few Earth masses[ The injection of a large number of planetesimals into the inner solar system may have probably played a key role in the formation and evolution of the atmospheres of the terrestrial planets "Oro\ 0850^ Ip and Fernandez\
9921Ð9522:88:, ! see front matter Þ 0888 Elsevier Science Ltd[ All rights reserved[ PII] S 9 9 2 1 Ð 9 5 2 2 " 8 7 # 9 9 0 3 9 Ð 7
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
481
0877#\ and in the collisional evolution of asteroids in the main belt "Brunini and Gil Hutton\ 0886#[ The numerical models previously used to investigate the Ý accretion of the outer solar system were Ý based on Opik|s formalism for two!body interactions "Opik\ 0840#[ These codes are essentially of probabilistic nature\ usually including the dynamics via short range interactions[ The origin of relevant features observed in the outer solar system such as its near!resonant structure must be inves! tigated through more complete numerical models\ includ! ing in a self consistent way Ý long range perturbations[ Another limitation of Opik|s formalism is that it allows bodies to collide if they are on mutual crossing orbits[ This is not a serious limitation whenever relative vel! ocities are high\ so radial excursions may be large[ In this case planetary embryos grow at a moderate rate\ because at large relative velocities\ the gravitational focusing fac! tor is small[ The most important regime of accretion\ the {runaway growth|\ occurs during the _rst stages of the process\ when relative velocities are small enough as to make dynamical friction operative "Wetherill and Stewart\ 0878^ Ida\ 0889^ Ida and Makino\ 0881^ Green! berg et al[\ 0880#[ Dynamical friction will damp the orbital eccentricity and inclination of the largest embryos at the expense of the smaller bodies that get more inclined and eccentric orbits\ thus reaching a kind of equipartition of the energy[ Lower inclinations and eccentricities mean larger collisional cross sections[ In this regime\ collisions are possible between bodies in orbits that do not overlap\ if they can approach to few times the radius of their mutual sphere of action\ that in a non!rotating reference frame is de_ned in the classical way RH a
Mp M=
1:4
0 1
\
"0#
where Mp is the mass of the planet and a its semimajor axis[ According to the most recent simulations of the for! mation of the outer solar system "Fernandez and Ip\ 0885#\ the accretion of Uranus and Neptune would demand time scales of a few 097 years[ These rather short time scales represent an important step if compared with previous results of 098 years or even longer than the age of the solar system "Safronov\ 0861#\ which posed the problem of formation Ý of the outer planets as a puzzle[ Unfortunately the Opik|s treatment will not be able to provide adequate information about time scales of plan! etary formation\ because it is e.cient when the accretion rate is low "in other words\ it is applicable to crossing orbits and not to grazing\ non!intersecting orbits for which the accretion e.ciency is higher#[ It is now a generalized consensus that the mostly gase! ous Jupiter and Saturn had to form before the dispersal of hydrogen and helium from the primordial solar nebula\ which is thought to have occurred much before 09
Myears[ This is supported by recent radio CO obser! vations "Zuckerman et al[\ 0884# and the observation of {naked T Tauri| stars with ages of approximately one million years "Walter et al[\ 0877#\ suggesting that the molecular gas surrounding young solar type stars tends to dissipate rapidly\ over time scales no longer than a few Myears[ The possible presence of a non!negligible amount of hydrogen and helium in Uranus and Neptune\ of the order of 0[5 M$ and 0[0 M$ respectively "Pollack et al[\ 0885#\ suggests that their time scales of formation could have not been much longer than the time scale of dissipation of the gaseous component of the nebula\ perhaps no longer than a few 096 years[ These time scales are however only marginally compatible with the more recent numerical results of the macro accretion of Uranus and Neptune "Fernandez and Ip\ 0885#[ A model includ! ing the runaway growth process\ might substantially reduce the time scales of formation of the outer planets\ bringing them more in line with the previous estimate[ Ý In view of the previously mentioned limitations of the Opik formalism to explain relevant features of the accretion of Uranus and Neptune\ we have decided to investigate this problem by means of a more complete numerical modeling of the accretion process\ based on the numerical integration of the equations of motion that naturally overcomes the limitations mentioned above[ This is now possible due to the recent development of a suitable technique] the simplectic integrators "Wisdom and Holman\ 0880#[ The main details of the numerical algorithm and of the model are described in the next section[ Section 2 is devoted to present the most relevant results of the numerical simulations of the accretion of Uranus and Neptune and Section 3 is left for the con! clusions[
1[ Numerical modeling 1[0[ Algorithm One of the most powerful tools for the numerical inte! gration of the equations of planetary motion is the sim! plectic algorithm recently developed by Wisdom and Holman "0880#[ The details of the implementation of the simplectic integrator are well known and are described there[ However\ as implemented by these authors\ the algorithm is not able to properly account for close plan! etary encounters[ As one of our main interests is to study the interaction between massive planetary embryos and a swarm of bodies of di}erent masses\ we need to intro! duce a way to evaluate close encounters[ One such way was already suggested by Wisdom and Holman "0880#[ When two bodies approach to less than a distance given by their mutual sphere of in~uence\ the motion may be considered as planetocentric\ perturbed by the action of the Sun[ Thus\ the simplectic algorithm may be applied
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
essentially in the same way as the one used to integrate the heliocentric motion[ This is the technique used by Levison and Duncan "0883# in their simplectic code[ Nevertheless\ this procedure has proved to be too expens! ive from the computational point of view when the num! ber of interacting bodies is large "in fact\ it would be completely impracticable on our computational scheme#[ We were thus forced to introduce some simpli_cations in the treatment of the problem[ As in all the previous simulations by Fernandez and Ip "0873\ 0885# the close encounters were computed as two!body scattering vents\ whenever the interacting bod! ies are inside their mutual sphere of action[ This treatment has been proved to be a good starting point for this kind of investigation in these previous studies[ Taking advantage of the symmetry of the two!body encounter\ the planetocentric coordinates and velocities at the end of the encounter may be easily computed by means of the F and G classical functions\ without even solving Kepler|s equation[ At low encounter velocities\ the two!body approxi! mation may not give accurate results for the outcome of any single encounter\ but the average outcome of many encounters have been proved to be in good agreement with detailed orbital integrations "Wetherill and Cox\ 0873#[ Massive three body integrations "Ida\ 0889# have shown that when the relative velocity is less than 9[96Ð9[07 times the escape velocity\ the viscous stirring rate under the solar gravity largely di}ers from the two body scattering formula[ Thus\ in our simulations we have monitored the relative velocity\ as a check of the accuracy and con! _dence of our results[ We have found that during the _rst 094 years\ less than 4) of the encounters are in the low velocity regime\ whereas this value falls to less than 9[4) after few 094 years[ Nevertheless\ we have made one additional run including the solar perturbation during close encounters at the low relative velocity regime[ The results of this simulation\ which are also reported in this paper\ have not shown any noticeable di}erence with the ones carried out with the simpli_ed model[ The prediction of the occurrence of close encounters has been a critical point and has needed a special treat! ment[ Computing at each integration step the mutual distances between N bodies in order to detect all pairs involved in close encounters requires O"N1# operations[ This very poor strategy is in fact not applicable to our computational scheme[ We have used a class of tree code based on recently developed algorithms used to simulate the dynamical evolution of galaxies and star clusters "Hernquist\ 0876#[ Tree structures are e.cient algorithms doing only O"N log N# operations to detect all those pairs involved in close encounters at any time[ The selected criterion to switch from heliocentric to planetocentric motion may introduce some degree of arbitrariness[ In addition\ the use of a relatively long time
482
step of integration may produce the loss of some non {penetrating| encounters[ To overcome these two di.! culties we have implemented the same technique used by Levison and Duncan "0883#] we have de_ned a sphere of 3[4 sphere of action|s radius around each body\ where the motion of the interacting bodies\ if approaching\ were integrated in heliocentric motion but with a shorter time step\ inversely proportional to the relative approaching velocity[ We have adopted this limiting distance based on some numerical experiments "Wisdom\ 0879^ Gladman\ 0882^ Kokubo and Ida\ 0884^ Chambers et al[\ 0885# which have shown that there is a minimum separation between two planets on near circular obits such that they would be stable against close approach[ This algorithm has proved to be very accurate in Gy[ integration of short period comets "Levison and Duncan\ 0883# and it is reminiscent to\ although very much simpler than\ the regularization techniques commonly used to treat close encounters "see e[g[ Stiefel and Shiefele\ 0860 for detailed discussion about regularization#[ Never! theless\ it is important to mention that our integrator "and the one of Levison and Duncan\ 0883# is not exactly simplectic[ It is worth noting\ however\ that the changes in the step!size of integration are sporadic\ thus we do not expect the simplectic behavior of the integrator to be substantially degraded[ Once the protoplanets grow to moderate masses\ the number of planetesimals\ and thus the number of close encounters\ will be small[ In this sense\ we hope to have an accurate representation of the dynamical evolution of the outer protoplanetary system[ For the heliocentric motion\ the time!step of inte! gration was of 9[1 years[ Our own experience has shown that with longer time!steps some close encounters would be lost[ The mutual perturbations between planets\ if present\ were included\ and the perturbing action of the planets over all the other bodies was always considered[ It is worth noting that we have introduced an arti_cial distinction between {particles| and {planets|] {planets| are those bodies massive enough to gravitationally interact with all other bodies of the system "including other planets\# either through long!range perturbations or short!range ones "close encounters#[ {Particles|\ instead\ are perturbed by the planets\ but they do not perturb the planets\ unless during close encounters[ After several collisions\ the mass of some particles may become large enough to be considered planets[ In all the simulations\ we have adopted 0 M$ as the lower limit for mass for a particle to pass to the category of planet[ As mentioned above\ this limit\ although rather arbi! trary\ is based on the following reasons] "0# Physical reasons] 0 M$ is only 09Ð19 times the initial mass of the planetesimals in our simulations[ There! fore\ our limit is consistent with the one used by Beauge et al[ "0883#[ In the outer solar system\ the distant perturbations of a planet of 0 M$ may be
483
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
neglected without substantially modifying the dynamical evolution of the planetesimal swarm[ As it will be clear later\ the presence of small planetary embryos is a common outcome in all our simulations\ specially during the _rst stages of evolution of the system[ However\ during these early stages\ the dynamics is largely dominated by close interactions with the planetesimals\ which are still large in number[ On the other hand\ the planetary embryos grow very fast to larger masses\ therefore the adopted limit of 0 M$ no longer becomes relevant after a very short period of time[ "1# Practical reasons] choosing a smaller mass cut!o} would produce a large number of planetary embryos that would not be possible to handle with our com! putational scheme[ This is because the numerical e}ort to compute mutual interactions between N planets scales as N1[ The introduction of this distinction between planets Ý and planetesimals is one of the di}erences with pure Opik Ý schemes[ In Opik|s algorithms\ the interaction between planets is not taken into account[ This is\ however\ one of the most relevant points for our study[ The computation of close encounters includes many re_nements in the numerical algorithm\ such as the way to decide when the outcome of an encounter is a collision[ This point was solved by considering that a collision occurred when the relative keplerian orbit has a peri! centric distance less then the sum of the geometric radius of the interacting bodies[ In this case\ the collision was considered perfectly inelastic "total accretion of the involved bodies#[ One major di}erence with the previous simulations by Fernandez and Ip "0873\ 0885# is on their adoption of extended atmospheres associated to the plan! etesimals as a way to enhance their e}ective cross sections[ In the present study\ the adopted cross!section is the sum of the real geometrical cross!sections of the interacting bodies[ In few test runs we have considered the presence of extended atmospheres in the same way as in Fernandez and Ip "0873\ 0885#[ The results of these runs do not show any noticeable di}erence[ We conclude this section noting that the only di}er! ences between these simulations and the ones Ý carried out with the stochastic codes based on the Opik algorithm are the predictive nature of the occurrence of close inter! actions and the consideration of planetary perturbations during the heliocentric motion[
integration of the asteroid belt for a Gy time!scale\ it was possible to reproduce the main known features of the asteroidal dynamics "Brunini\ 0887#[ In addition\ the dynamics of mutual interacting bodies was tested by com! paring the long term planetary motion obtained with the simplectic code and with other well known standard integrators such as Bullirsch and Stoer or Radau[ These tests were completely satisfactory[ One of the several procedures that were used to test the behavior of the part of the code that integrates the orbits of planetesimals that su}er close encounters with the planets\ was by studying the circular restricted three body problem\ for which an analytic integral of motion exists] the Jacobi integral[ We have performed a great number of numerical experiments[ The results of one representative case illustrate the accuracy of the method[ Figure 0a displays the relative error in the Jacobi integral for a massless particle that su}ers repeated close encoun! ters with a planet[ In this case\ a Jupiter like planet was placed on a circular orbit at 0 AU from the Sun and the particle was initially placed on an adjacent orbit with low eccentricity and inclination\ in such a way that it undergoes repeated close encounters at low relative velocity[ In this experiment\ the time step was of 0:39 the period of the planet[ Figure 0b shows the distance to the planet in astronomical units[ If the particle falls between the solid and the dotted line\ the time!step is decreased as explained in the preceding section[ If the particle falls
1[1[ Tests of the numerical code Before starting the numerical simulations of planetary accretion\ we have exhaustively tested the numerical code[ The core of the code] the simplectic integrator of the planetary motion\ was extensively tested[ In numerical
Fig[ 0[ The behavior of a test particle on a planet!encountering orbit in the circular restricted three body problem] "a# the relative error in the Jacobi constant^ "b# the relative distance to the planet[
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
below the dotted line the integration is performed in the planetocentric frame as a two body event[ The large observed changes in the Jacobi constant are due to close encounters with the planet[ In hundreds of numerical experiments\ we have not found systematic trends in the errors of the Jacobi constant[ In fact\ the behavior is very similar to the one found by Levison and Duncan "0883#[ The error in the Jacobi constant propagates in a form consistent with a random walk\ and the typical fractional error was 0×09−4Ð4×09−4 after 09 close encounters[ Another analyzed question was if the repeated switches between heliocentric to planetocentric motion could introduce any arti_cial heating in a system of plan! etesimals[ We do not expect such a behavoir because the results of Levison and Duncan "0883# have not shown this[ In fact\ if the solar perturbations during close encounters were included\ any di}erence could be intro! duced by switching Ý the origin of coordinates[ On the other hand\ in Opik|s schemes there are also switches from heliocentric to planetocentric motion[ The only di}erence with our calculation is that in our scheme\ the heliocentric motion is perturbed by the presence of planets[ Nevertheless\ in order to shed light on this ques! tion we have performed a numerical integration of a swarm of 099 Pluto like planetesimals of equal mass distributed between 04 and 24 AU in near circular and near coplanar orbits in the _eld of the four outer planets\ that have been considered with its present masses and orbits[ We have used our code with stepsize of 9[1 years and a full N!body Bullirsch and Stoer integrator with pre! scribed tolerance of 09−09[ During the integration\ which covered a time span of 094 years\ about 399 switches from heliocentric to planetocentric motion were computed[ Figure 1a and b displays the behavior of the r[m[s[ eccen! tricity and inclination of the swarm[ No important di}erences are found between the heating rate in both integrations[ 1[2[ Initial conditions The observational evidence summarized in the intro! duction suggests that the more realistic way to consider the growth of Uranus and Neptune by accumulation of smaller embryos is to assume that Jupiter and Saturn had already reached their present mass\ but that a large number of embryos were still growing at greater distances[ It is currently accepted that when Jupiter and Saturn were fully grown\ embryos of masses ranging from a fraction to several M$ could have been located in the outer solar system "Lissauer\ 0882#[ The macro accretion of this swarm of large bodies probably represented the main contribution to the _nal masses of Uranus and Neptune[ Following this line of reasoning\ in all our numerical simulations we have started with Jupiter and Saturn
484
Fig[ 1[ Comparison between the temporal evolution of the r[m[s[ eccen! tricity "a# and inclination "b# for a swarm of 099 planetesimals in the _eld of the four giant planets\ computed by means of a full N!body integration "open circles# and our simplectic code "triangles#[
already present with their current masses[ The orbital elements of Jupiter were always equal to the present ones\ while for Saturn we took in some cases its current orbital elements\ while in other cases it started somewhat closer to the Sun to make allowance for its outward drift[ In the outer solar system\ a swarm of N bodies of equal mass was distributed in a ~at structure between 01 and 24 AU "in some runs the outer boundary was extended up to 39 AU#[ The adopted radial distribution was dN: daa−g with g 0:1\ which follows some models of the mass distribution in the protoplanetary disk "Wei! denschilling\ 0866#[ We adopted small initial eccentricities and inclinations randomly distributed with maximum values ranging from 9[994 up to 9[0[ The total initial mass of the entire swarm of bodies was de_ned as a factor f of the combined mass of the current Uranus and Neptune[ Because of the mentioned ine.ciency of the accretion process of the outer planets\ we adopted f × 0[ As the maximum number of bodies that can be handled by our computational scheme is limited\ we have per! formed several runs with di}erent initial number N of planetesimals "149\ 499 and 0999#\ in order to check the e}ect of the _nite number of bodies on the main outcomes of the simulations[ As the initial mass in the planetesimal swarm is typically of 1Ð2 times the combined mass of Uranus and Neptune\ the planetesimals have initial masses in the range 09−1Ð09−2 MU\ MU being the present mass of Uranus[ The time span covered by the numerical simulations was of 099 Myears in all the runs[
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
485
Duncan "0883#[ It was _ve time more time consuming than the ones without this e}ect[ We have carefully exam! ined all the outcomes of this simulation[ We have not found any noticeable di}erence with other simulations of the set[ Remarkably\ in all the runs there is at least one massive planet[ From our simulations we found that the accretion process by mutual collisions of the planetesimals usually evolved to the formation of two planets\ though in some cases we obtained one or three large planets[ In those cases where only one massive planet is formed\ its orbit lies beyond 29 AU[ Figure 2 shows the orbital evolution of the proto! Uranus and proto!Neptune for several representative runs[ The orbital migration\ reported in the previous simulations by Fernandez and Ip "0873\ 0885# is also found with our more realistic model\ although we observe a strong dynamical interaction between the protoplanets at early stages of the evolution[ An outward migration is also the case for Saturn[ It has in fact experienced a very strong displacement of several AU[ In most cases\ this strong migration is responsible for the destabilization of the proto!Uranus[ This is a limitation of our model\
2[ Results Table 0 displays the main results of 20 simulations with di}erent initial conditions[ The columns refer to] "0# number of run^ "1# initial number of bodies^ "2# the total mass initially present in terms of the combined masses of Uranus and Neptune " 20[7 M$#^ "3# an estimate of the maximum relative velocity in terms of the circular velocity at each helio! centric distance\ where i9 and e9 are the maximum possible initial inclination "in radians# and eccentricity for the planetesimals in each run^ "4# the maximum number of embryo planets that reach at least one Earth mass during the given run^ "5# the number of planets remaining in the outer solar system at the end of each run^ "6# and "7# the _nal masses of our computed {Uranus| and {Neptune| in terms of Uranus|s mass " 03[43 M$#^ "8# and "09# the time of the last macroaccretional event for proto!Uranus and proto!Neptune in Myears[ The _nal results of the simulations seem to be not strongly dependent on the initial number of plan! etesimals[ Run 20 was made considering the solar per! turbation during close encounters as in Levison and
Table 0 Initial conditions and main results of 20 numerical simulations Run
N
Minitial
ze19¦i19
Nemb
Npla
MU
MN
TU
TN
0 1 2 3 4 5 6 7 8 09 00 01 02 03 04 05 06 07 08 19 10 11 12 13 14 15 16 17 18 29 20
499 499 499 499 499 499 499 499 499 499 499 499 499 0999 0999 0999 0999 0999 0999 0999 0999 0999 0999 0999 0999 149 149 499 0999 0999 499
0[7 0[8 0[7 1[0 1[4 0[6 0[6 1[6 1[9 2[2 1[4 2[9 1[9 1[3 0[4 0[6 0[7 2[9 1[4 2[7 2[9 1[4 2[7 2[9 2[9 1[9 0[4 2[9 1[4 1[4 1[9
9[996 9[996 9[996 9[960 9[960 9[960 9[960 9[903 9[960 9[960 9[960 9[903 9[960 9[949 9[030 9[030 9[903 9[001 9[903 9[917 9[903 9[996 9[903 9[903 9[917 9[903 9[903 9[960 9[960 9[960 9[960
00 09 01 02 6 09 5 7 00 3 3 3 02 12 02 08 7 04 04 01 02 03 07 8 02 00 04 09 08 02 5
2 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 0 1 1 1
0[7 0[1 9[1 9[5 9[3 1[3 * 1[9 9[4 * 9[5 0[4 9[4 9[6 0[1 9[7 9[0 0[4 * 9[5 * 0[7 * 9[0 * 9[5 * * 0[6 9[8 9[6
9[2 0[1 0[0 1[1 0[1 9[3 0[3 0[0 9[6 3[1 2[4 1[0 9[2 0[4 9[4 9[2 0[2 0[6 0[6 0[4 9[4 9[5 0[1 0[9 9[7 9[5 0[6 9[6 0[9 1[9 9[3
56 00 63 68 70 88 * 22 78 * 71 69 21 71 74 57 89 5 * 79 * 39 * 57 * 35 * * 24 14 01
09 08 20 36 35 4 31 12 15 54 50 16 76 59 20 09 58 4 50 69 81 4 50 8 51 00 43 32 36 06 10
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
486
Fig[ 2[ Evolution of the semimajor axis of the outer planets for some of the runs where two massive planets were formed[
imposed by the inclusion of Saturn with its present mass and orbit[ We have also made some runs starting Saturn with a 7 AU without noticeable di}erences "runs 15Ð 29#[ The general rule in our simulation is that Jupiter drifts inwards[ Surprisingly for us\ in a few runs Jupiter also experienced a net outward migration[ This is explained as due to the ejection by Jupiter of a large number of
planetesimals with less angular momentum than they initially had[ The angular momentum lost by these ejected bodies is compensated in these cases by an outward drift of Jupiter[ This e}ect is not observed in the simulations by Fernandez and Ip because they have not considered a collision with the Sun as a possible end!state of a plan! etesimal[ If we followed the evolution of these objects "ignoring the possible collision# they would su}er a close
487
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
encounter with Jupiter at some later time\ being ejected from the system in hyperbolic orbits[ In this situation\ Jupiter|s orbital drift would be always inward\ as was previously found by Fernandez and Ip[ We observe that several protoplanets form during the process of accumulation of the outer planets[ Figure 3 exhibits the radial distribution of the protoplanets in one of the runs that resulted in the formation of two massive planets[ In short time scales\ of the order of few 096 years\ up to 00 protoplanets are competing for a place in the region beyond Saturn[ Many of them accrete to form more massive embryos[ It is worth noting that some large embryos reach semimajor axes in the region of the Kuiper belt "½39Ð59 AU# remaining there for several 095 years[ The dynamical and physical consequences of this phenomenon may be relevant for the evolution of the belt[ A large number of massive protoplanets has\ in some cases\ a destabilizing e}ect over the entire system\ resulting in a strong interaction between proto!Uranus and proto!Neptune[ The large jumps in semimajor axis are associated to encounters with these large proto! planets[ Our results suggest that there is a range of initial masses from which two planets like Uranus and Neptune can be formed\ that is in good agreement with the values previously found by Fernandez and Ip "0873\ 0885#[ There is a minimum amount of initial mass\ in the range of 1Ð2 times the combined masses of Uranus and Neptune\ required to form planets the size of Uranus and Neptune[ For the 11 runs with 1 outer planets\ an average of only ½35) of the solid material initially present was
Fig[ 3[ Evolution of the semimajor axis of the protoplanets formed during one of the runs "open circles#[ We de_ne as {protoplanets| those bodies reaching masses greater than one Earth mass[
incorporated into the planets[ Therefore\ the increasing probability of ejection of interacting planetesimals with the massive protoplanets during the latest stage of accretion is con_rmed by our model[ A large fraction of the bodies scattered by the accreting proto!Uranus and proto!Neptune fall under the gravi! tational control of Saturn and Jupiter[ Once this happens the dynamical evolution of the scattered bodies is char! acterized by a fast random!walk in perihelion distance q in such a way that a fraction of them become Mars crossers before being ejected in hyperbolic or near para! bolic orbits[ Some of them could _nally collide with Jup! iter or Saturn[ We found that 9[0Ð9[4 M$ of solid material from the UranusÐNeptune region was incorporated into Jupiter and into Saturn[ These values are in contrast with the numerical results previously found by Fernandez and Ip "0873#\ where up to 1 M$ of solid material was incor! porated into Jupiter\ although they started their simu! lations with comparable amount of mass[ We have not found any correlation between the initial mass and the total mass incorporated into Jupiter and Saturn[ Figure 4 displays the mean e.ciency\ averaging over all the runs\ of each planet in ejecting planetesimals from the system[ As expected\ Jupiter is the most e.cient ejec! tor\ although there is a non!negligible fraction of bodies ejected by the dynamical in~uence of Saturn\ in the region between ½7 and ½04 AU[ The large fraction of bodies interacting with Saturn "part of the material is ejected and part is handed down to Jupiter# is responsible for its strong orbital migration[ The distribution of perihelion distances of the ejected bodies displayed in Fig[ 4 shows that a large fraction is ejected at q ³ 4 AU\ con_rming our argument about the
Fig[ 4[ Relative distribution of the perihelion distances of those bodies ejected from the system^ mean values for 20 simulations[
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
origin of the outward migration of Jupiter in a few runs[ The distribution of perihelion distances of bodies reach! ing the region of the terrestrial planets shows a steep decrease towards smaller q\ so that only a small fraction of the bodies crossing Mars|s orbit also cross Mercury|s[ As shown in Fig[ 5 for a run with 499 bodies\ the total number of surviving planetesimals decreases expo! nentially with time[ We explained above\ the losses are due to accretion by the growing outer planets "the fast loss regime shown in Fig[ 5 during the _rst Myears#\ or ejection to interstellar space[ The cumulative mass reaching the region of the ter! restrial planets is found to be of the order of 09Ð04 Earth masses\ becoming potential impactors of these four planets[ The evolution of the cumulative mass entering within the zone r ³ 0[4 AU is displayed in Fig[ 6[ The solid line is the average value over our 20 simulations[ The dotted lines are the maximum and minimum values found in the simulations[ The phase of early strong bom! bardment of the terrestrial planets due to scattered plan! etesimals from the UranusÐNeptune region reached its maximum near ½096 years\ lasting\ but at a slower rate\ for ½3×096 years[ This is a much shorter time span than previously found "Fernandez and Ip\ 0885#[ About 19 times the amount of water contained in the oceans could have reached the Earth from this bombardment of icy material "Fernandez and Ip\ 0886#\ but this is a very uncertain result at this point[ The answer depends on whether or not the Earth and the other terrestrial planets were already formed when the in~ux of planetesimals from the outer planetary region reached its peak[ If they were still in the process of accretion most of the accreted volatile material might have been lost in subsequent megaimpacts[ The masses of proto!Uranus and proto!Neptune grow very fast[ Figure 7 displays the growth rate of proto!
Fig[ 5[ Evolution of the number of planetesimals for a run with initially 499 bodies[
488
Fig[ 6[ Cumulative mass\ in M$\ that reached the region of the terrestrial planets[ The solid line is the mean value for 20 simulations[ The dashed lines are upper and lower bounds to the cumulative mass[
Uranus and proto!Neptune for three runs[ We may observe that the time scale of complete planetary for! mation is of few 096 years\ either for Uranus or Neptune[ The slope of the curves exhibits the classical behavior of runaway growth "Ida\ 0889#[ The large radial excursions experienced by the protoplanets\ such as the ones shown in Fig[ 2\ bring them into zones not depleted of plan! etesimals\ thus substantially extending their feeding zones[ Moreover\ the orbits in these undepleted regions were not previously excited by the protoplanets\ and therefore\ the orbital migration seems to be an e.cient way to sustain the regime of runaway growth for long time spans[ Also noticeable in Fig[ 7 is that Neptune generally grows faster than Uranus\ reaching also a big! ger mass[ This fact may be related to the radial excursions that are in general larger for Neptune than for Uranus[ Figure 8 displays\ for one of the runs\ a combination of Figs 2 and 7\ where we illustrate the correlation between the orbital migration and the change in mass[ As previously noted\ the maximum initial values for the orbital eccentricity and inclination varies between 9[994 in some cases and the more realistic 9[0[ We have not noticed important di}erences arising from these di}erent initial conditions\ mainly because even when the initial eccentricity and inclination are very low\ the system evolves very fast to "e\ i# ½ 9[0[ The evolution of the r[m[s[ eccentricity and inclination for one of the runs is shown in Fig[ 09[ In one of the runs we have additionally computed the relative velocity and the distance of closest approach for all the encounters experienced by proto!Neptune[ The average value of the relative velocity in terms of the proto! Neptune orbital velocity is shown in Fig[ 00a[ Each value was computed averaging over 099 encounters[ It is inter! esting to note that the relative velocity increases fast and
599
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
Fig[ 7[ Accretion rate of proto!Uranus and proto!Neptune for some of the runs[
Planetary and Space Science 36 "0888# 480Ð594
Numerical simulations of the accretion of Uranus and Neptune Adrian Bruninia\ \ Julio A[ Fernandezb a
Facultad de Ciencias Astronomics y Geof(sicas\ Universidad Nacional de La Plata "0899# La Plata\ Ar`entina b Departamento de Astronom(a\ Facultad de Ciencias\ I`ua esq[ Mataojo\ 00399 Montevideo\ Uru`uay Received 2 March 0887^ received in revised form 4 November 0887^ accepted 11 November 0887
Abstract New numerical simulations of the accretion process of the outer planets are carried out using an N!body code[ In 10 out of 29 computer runs 1 giant planets form beyond Saturn\ whereas 0 forms in 7 cases and 2 in 0 case[ The time scale of formation is a few 096 years[ This rather short time scale seems to be consistent with the content of non!negligible amounts of hydrogen and helium in Uranus and Neptune\ which could be accreted before the protoplanetary nebula dispersed[ A marked radial drift of the accreting outer planets is found in our simulations in agreement with previous results[ The results also seem to con_rm that this radial migration is an e}ective mechanism to produce orbital con_gurations near mean motion resonances[ Less than 49) of the solid material originally present in the system as planetesimals contributes to the formed protoplanets\ the rest being mostly ejected by Jupiter and Saturn[ The phase of strong scattering of bodies from the UranusÐNeptune region to the inner solar system lasts for a period of ½3×096 years[ Around 09Ð19 Earth masses reach the region of the terrestrial planets\ although the very short time scale for the massive in~ux of material "perhaps shorter than the time scale of Earth formation# might pose some di.culties to the theory that comets were major contributors of water to the Earth[ Þ 0888 Elsevier Science Ltd[ All rights reserved[
0[ Introduction The accretion process of the outer planets was one of the most collective phenomena during the formation of the solar system[ From the region of the terrestrial planets up to its outermost boundary\ the Oort cloud\ the solar system has been strongly in~uenced by this process "Fer! nandez and Gallardo\ 0885#[ Much progress has been done in our understanding of the formation of the outer planets and scattering of the residual mass from the early theoretical research of Safronov "0861#[ Nevertheless\ there are several aspects not yet understood that deserve further study[ The numerical simulations of the macro!accretion pro! cess of Uranus and Neptune carried out by Fernandez and Ip "0873\ 0885#\ have properly accounted for some important features observed in the outer planetary region such as the presence of two massive planets beyond Saturn[ One interesting and unexpected result found by Fernandez and Ip "0873# was that the orbit of Neptune\ and to a lesser extent those of Uranus and Saturn experi! enced an outward drift in the course of their accretion due to the exchange of angular momentum with the inter! Corresponding author[ E!mail] abruniniÝfcaglp[fcaglp[unlp[edu[ar
acting planetesimals[ The angular momentum gained by the orbital expansion of these planets was compensated by a small drift inwards of the massive Jupiter[ The planetary migration may have played an impor! tant role during the late accretion phase of the outer solar system\ being recently claimed as responsible for the presence of Pluto in its eccentric and inclined orbit at 2 ] 1 resonance with Neptune "Malhotra\ 0884#\ and of the global quasi!resonant structure of the outer solar system "Fernandez and Ip\ 0885#[ The numerical models developed by Fernandez and Ip "0873\ 0885# also suggested that the formation of Uranus and Neptune required an initial mass in their accretion zones about 1Ð2 times larger than their combined masses[ The unaccreted solid material was lost to the inner solar system or to the interstellar space[ Therefore\ the accretion of Uranus and Neptune seems to have been very ine.cient during the latest stages\ which may be explained as due to the increasing probability of ejection of interacting planetesimals with protoplanets whose masses exceed a few Earth masses[ The injection of a large number of planetesimals into the inner solar system may have probably played a key role in the formation and evolution of the atmospheres of the terrestrial planets "Oro\ 0850^ Ip and Fernandez\
9921Ð9522:88:, ! see front matter Þ 0888 Elsevier Science Ltd[ All rights reserved[ PII] S 9 9 2 1 Ð 9 5 2 2 " 8 7 # 9 9 0 3 9 Ð 7
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
481
0877#\ and in the collisional evolution of asteroids in the main belt "Brunini and Gil Hutton\ 0886#[ The numerical models previously used to investigate the Ý accretion of the outer solar system were Ý based on Opik|s formalism for two!body interactions "Opik\ 0840#[ These codes are essentially of probabilistic nature\ usually including the dynamics via short range interactions[ The origin of relevant features observed in the outer solar system such as its near!resonant structure must be inves! tigated through more complete numerical models\ includ! ing in a self consistent way Ý long range perturbations[ Another limitation of Opik|s formalism is that it allows bodies to collide if they are on mutual crossing orbits[ This is not a serious limitation whenever relative vel! ocities are high\ so radial excursions may be large[ In this case planetary embryos grow at a moderate rate\ because at large relative velocities\ the gravitational focusing fac! tor is small[ The most important regime of accretion\ the {runaway growth|\ occurs during the _rst stages of the process\ when relative velocities are small enough as to make dynamical friction operative "Wetherill and Stewart\ 0878^ Ida\ 0889^ Ida and Makino\ 0881^ Green! berg et al[\ 0880#[ Dynamical friction will damp the orbital eccentricity and inclination of the largest embryos at the expense of the smaller bodies that get more inclined and eccentric orbits\ thus reaching a kind of equipartition of the energy[ Lower inclinations and eccentricities mean larger collisional cross sections[ In this regime\ collisions are possible between bodies in orbits that do not overlap\ if they can approach to few times the radius of their mutual sphere of action\ that in a non!rotating reference frame is de_ned in the classical way RH a
Mp M=
1:4
0 1
\
"0#
where Mp is the mass of the planet and a its semimajor axis[ According to the most recent simulations of the for! mation of the outer solar system "Fernandez and Ip\ 0885#\ the accretion of Uranus and Neptune would demand time scales of a few 097 years[ These rather short time scales represent an important step if compared with previous results of 098 years or even longer than the age of the solar system "Safronov\ 0861#\ which posed the problem of formation Ý of the outer planets as a puzzle[ Unfortunately the Opik|s treatment will not be able to provide adequate information about time scales of plan! etary formation\ because it is e.cient when the accretion rate is low "in other words\ it is applicable to crossing orbits and not to grazing\ non!intersecting orbits for which the accretion e.ciency is higher#[ It is now a generalized consensus that the mostly gase! ous Jupiter and Saturn had to form before the dispersal of hydrogen and helium from the primordial solar nebula\ which is thought to have occurred much before 09
Myears[ This is supported by recent radio CO obser! vations "Zuckerman et al[\ 0884# and the observation of {naked T Tauri| stars with ages of approximately one million years "Walter et al[\ 0877#\ suggesting that the molecular gas surrounding young solar type stars tends to dissipate rapidly\ over time scales no longer than a few Myears[ The possible presence of a non!negligible amount of hydrogen and helium in Uranus and Neptune\ of the order of 0[5 M$ and 0[0 M$ respectively "Pollack et al[\ 0885#\ suggests that their time scales of formation could have not been much longer than the time scale of dissipation of the gaseous component of the nebula\ perhaps no longer than a few 096 years[ These time scales are however only marginally compatible with the more recent numerical results of the macro accretion of Uranus and Neptune "Fernandez and Ip\ 0885#[ A model includ! ing the runaway growth process\ might substantially reduce the time scales of formation of the outer planets\ bringing them more in line with the previous estimate[ Ý In view of the previously mentioned limitations of the Opik formalism to explain relevant features of the accretion of Uranus and Neptune\ we have decided to investigate this problem by means of a more complete numerical modeling of the accretion process\ based on the numerical integration of the equations of motion that naturally overcomes the limitations mentioned above[ This is now possible due to the recent development of a suitable technique] the simplectic integrators "Wisdom and Holman\ 0880#[ The main details of the numerical algorithm and of the model are described in the next section[ Section 2 is devoted to present the most relevant results of the numerical simulations of the accretion of Uranus and Neptune and Section 3 is left for the con! clusions[
1[ Numerical modeling 1[0[ Algorithm One of the most powerful tools for the numerical inte! gration of the equations of planetary motion is the sim! plectic algorithm recently developed by Wisdom and Holman "0880#[ The details of the implementation of the simplectic integrator are well known and are described there[ However\ as implemented by these authors\ the algorithm is not able to properly account for close plan! etary encounters[ As one of our main interests is to study the interaction between massive planetary embryos and a swarm of bodies of di}erent masses\ we need to intro! duce a way to evaluate close encounters[ One such way was already suggested by Wisdom and Holman "0880#[ When two bodies approach to less than a distance given by their mutual sphere of in~uence\ the motion may be considered as planetocentric\ perturbed by the action of the Sun[ Thus\ the simplectic algorithm may be applied
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
essentially in the same way as the one used to integrate the heliocentric motion[ This is the technique used by Levison and Duncan "0883# in their simplectic code[ Nevertheless\ this procedure has proved to be too expens! ive from the computational point of view when the num! ber of interacting bodies is large "in fact\ it would be completely impracticable on our computational scheme#[ We were thus forced to introduce some simpli_cations in the treatment of the problem[ As in all the previous simulations by Fernandez and Ip "0873\ 0885# the close encounters were computed as two!body scattering vents\ whenever the interacting bod! ies are inside their mutual sphere of action[ This treatment has been proved to be a good starting point for this kind of investigation in these previous studies[ Taking advantage of the symmetry of the two!body encounter\ the planetocentric coordinates and velocities at the end of the encounter may be easily computed by means of the F and G classical functions\ without even solving Kepler|s equation[ At low encounter velocities\ the two!body approxi! mation may not give accurate results for the outcome of any single encounter\ but the average outcome of many encounters have been proved to be in good agreement with detailed orbital integrations "Wetherill and Cox\ 0873#[ Massive three body integrations "Ida\ 0889# have shown that when the relative velocity is less than 9[96Ð9[07 times the escape velocity\ the viscous stirring rate under the solar gravity largely di}ers from the two body scattering formula[ Thus\ in our simulations we have monitored the relative velocity\ as a check of the accuracy and con! _dence of our results[ We have found that during the _rst 094 years\ less than 4) of the encounters are in the low velocity regime\ whereas this value falls to less than 9[4) after few 094 years[ Nevertheless\ we have made one additional run including the solar perturbation during close encounters at the low relative velocity regime[ The results of this simulation\ which are also reported in this paper\ have not shown any noticeable di}erence with the ones carried out with the simpli_ed model[ The prediction of the occurrence of close encounters has been a critical point and has needed a special treat! ment[ Computing at each integration step the mutual distances between N bodies in order to detect all pairs involved in close encounters requires O"N1# operations[ This very poor strategy is in fact not applicable to our computational scheme[ We have used a class of tree code based on recently developed algorithms used to simulate the dynamical evolution of galaxies and star clusters "Hernquist\ 0876#[ Tree structures are e.cient algorithms doing only O"N log N# operations to detect all those pairs involved in close encounters at any time[ The selected criterion to switch from heliocentric to planetocentric motion may introduce some degree of arbitrariness[ In addition\ the use of a relatively long time
482
step of integration may produce the loss of some non {penetrating| encounters[ To overcome these two di.! culties we have implemented the same technique used by Levison and Duncan "0883#] we have de_ned a sphere of 3[4 sphere of action|s radius around each body\ where the motion of the interacting bodies\ if approaching\ were integrated in heliocentric motion but with a shorter time step\ inversely proportional to the relative approaching velocity[ We have adopted this limiting distance based on some numerical experiments "Wisdom\ 0879^ Gladman\ 0882^ Kokubo and Ida\ 0884^ Chambers et al[\ 0885# which have shown that there is a minimum separation between two planets on near circular obits such that they would be stable against close approach[ This algorithm has proved to be very accurate in Gy[ integration of short period comets "Levison and Duncan\ 0883# and it is reminiscent to\ although very much simpler than\ the regularization techniques commonly used to treat close encounters "see e[g[ Stiefel and Shiefele\ 0860 for detailed discussion about regularization#[ Never! theless\ it is important to mention that our integrator "and the one of Levison and Duncan\ 0883# is not exactly simplectic[ It is worth noting\ however\ that the changes in the step!size of integration are sporadic\ thus we do not expect the simplectic behavior of the integrator to be substantially degraded[ Once the protoplanets grow to moderate masses\ the number of planetesimals\ and thus the number of close encounters\ will be small[ In this sense\ we hope to have an accurate representation of the dynamical evolution of the outer protoplanetary system[ For the heliocentric motion\ the time!step of inte! gration was of 9[1 years[ Our own experience has shown that with longer time!steps some close encounters would be lost[ The mutual perturbations between planets\ if present\ were included\ and the perturbing action of the planets over all the other bodies was always considered[ It is worth noting that we have introduced an arti_cial distinction between {particles| and {planets|] {planets| are those bodies massive enough to gravitationally interact with all other bodies of the system "including other planets\# either through long!range perturbations or short!range ones "close encounters#[ {Particles|\ instead\ are perturbed by the planets\ but they do not perturb the planets\ unless during close encounters[ After several collisions\ the mass of some particles may become large enough to be considered planets[ In all the simulations\ we have adopted 0 M$ as the lower limit for mass for a particle to pass to the category of planet[ As mentioned above\ this limit\ although rather arbi! trary\ is based on the following reasons] "0# Physical reasons] 0 M$ is only 09Ð19 times the initial mass of the planetesimals in our simulations[ There! fore\ our limit is consistent with the one used by Beauge et al[ "0883#[ In the outer solar system\ the distant perturbations of a planet of 0 M$ may be
483
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
neglected without substantially modifying the dynamical evolution of the planetesimal swarm[ As it will be clear later\ the presence of small planetary embryos is a common outcome in all our simulations\ specially during the _rst stages of evolution of the system[ However\ during these early stages\ the dynamics is largely dominated by close interactions with the planetesimals\ which are still large in number[ On the other hand\ the planetary embryos grow very fast to larger masses\ therefore the adopted limit of 0 M$ no longer becomes relevant after a very short period of time[ "1# Practical reasons] choosing a smaller mass cut!o} would produce a large number of planetary embryos that would not be possible to handle with our com! putational scheme[ This is because the numerical e}ort to compute mutual interactions between N planets scales as N1[ The introduction of this distinction between planets Ý and planetesimals is one of the di}erences with pure Opik Ý schemes[ In Opik|s algorithms\ the interaction between planets is not taken into account[ This is\ however\ one of the most relevant points for our study[ The computation of close encounters includes many re_nements in the numerical algorithm\ such as the way to decide when the outcome of an encounter is a collision[ This point was solved by considering that a collision occurred when the relative keplerian orbit has a peri! centric distance less then the sum of the geometric radius of the interacting bodies[ In this case\ the collision was considered perfectly inelastic "total accretion of the involved bodies#[ One major di}erence with the previous simulations by Fernandez and Ip "0873\ 0885# is on their adoption of extended atmospheres associated to the plan! etesimals as a way to enhance their e}ective cross sections[ In the present study\ the adopted cross!section is the sum of the real geometrical cross!sections of the interacting bodies[ In few test runs we have considered the presence of extended atmospheres in the same way as in Fernandez and Ip "0873\ 0885#[ The results of these runs do not show any noticeable di}erence[ We conclude this section noting that the only di}er! ences between these simulations and the ones Ý carried out with the stochastic codes based on the Opik algorithm are the predictive nature of the occurrence of close inter! actions and the consideration of planetary perturbations during the heliocentric motion[
integration of the asteroid belt for a Gy time!scale\ it was possible to reproduce the main known features of the asteroidal dynamics "Brunini\ 0887#[ In addition\ the dynamics of mutual interacting bodies was tested by com! paring the long term planetary motion obtained with the simplectic code and with other well known standard integrators such as Bullirsch and Stoer or Radau[ These tests were completely satisfactory[ One of the several procedures that were used to test the behavior of the part of the code that integrates the orbits of planetesimals that su}er close encounters with the planets\ was by studying the circular restricted three body problem\ for which an analytic integral of motion exists] the Jacobi integral[ We have performed a great number of numerical experiments[ The results of one representative case illustrate the accuracy of the method[ Figure 0a displays the relative error in the Jacobi integral for a massless particle that su}ers repeated close encoun! ters with a planet[ In this case\ a Jupiter like planet was placed on a circular orbit at 0 AU from the Sun and the particle was initially placed on an adjacent orbit with low eccentricity and inclination\ in such a way that it undergoes repeated close encounters at low relative velocity[ In this experiment\ the time step was of 0:39 the period of the planet[ Figure 0b shows the distance to the planet in astronomical units[ If the particle falls between the solid and the dotted line\ the time!step is decreased as explained in the preceding section[ If the particle falls
1[1[ Tests of the numerical code Before starting the numerical simulations of planetary accretion\ we have exhaustively tested the numerical code[ The core of the code] the simplectic integrator of the planetary motion\ was extensively tested[ In numerical
Fig[ 0[ The behavior of a test particle on a planet!encountering orbit in the circular restricted three body problem] "a# the relative error in the Jacobi constant^ "b# the relative distance to the planet[
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
below the dotted line the integration is performed in the planetocentric frame as a two body event[ The large observed changes in the Jacobi constant are due to close encounters with the planet[ In hundreds of numerical experiments\ we have not found systematic trends in the errors of the Jacobi constant[ In fact\ the behavior is very similar to the one found by Levison and Duncan "0883#[ The error in the Jacobi constant propagates in a form consistent with a random walk\ and the typical fractional error was 0×09−4Ð4×09−4 after 09 close encounters[ Another analyzed question was if the repeated switches between heliocentric to planetocentric motion could introduce any arti_cial heating in a system of plan! etesimals[ We do not expect such a behavoir because the results of Levison and Duncan "0883# have not shown this[ In fact\ if the solar perturbations during close encounters were included\ any di}erence could be intro! duced by switching Ý the origin of coordinates[ On the other hand\ in Opik|s schemes there are also switches from heliocentric to planetocentric motion[ The only di}erence with our calculation is that in our scheme\ the heliocentric motion is perturbed by the presence of planets[ Nevertheless\ in order to shed light on this ques! tion we have performed a numerical integration of a swarm of 099 Pluto like planetesimals of equal mass distributed between 04 and 24 AU in near circular and near coplanar orbits in the _eld of the four outer planets\ that have been considered with its present masses and orbits[ We have used our code with stepsize of 9[1 years and a full N!body Bullirsch and Stoer integrator with pre! scribed tolerance of 09−09[ During the integration\ which covered a time span of 094 years\ about 399 switches from heliocentric to planetocentric motion were computed[ Figure 1a and b displays the behavior of the r[m[s[ eccen! tricity and inclination of the swarm[ No important di}erences are found between the heating rate in both integrations[ 1[2[ Initial conditions The observational evidence summarized in the intro! duction suggests that the more realistic way to consider the growth of Uranus and Neptune by accumulation of smaller embryos is to assume that Jupiter and Saturn had already reached their present mass\ but that a large number of embryos were still growing at greater distances[ It is currently accepted that when Jupiter and Saturn were fully grown\ embryos of masses ranging from a fraction to several M$ could have been located in the outer solar system "Lissauer\ 0882#[ The macro accretion of this swarm of large bodies probably represented the main contribution to the _nal masses of Uranus and Neptune[ Following this line of reasoning\ in all our numerical simulations we have started with Jupiter and Saturn
484
Fig[ 1[ Comparison between the temporal evolution of the r[m[s[ eccen! tricity "a# and inclination "b# for a swarm of 099 planetesimals in the _eld of the four giant planets\ computed by means of a full N!body integration "open circles# and our simplectic code "triangles#[
already present with their current masses[ The orbital elements of Jupiter were always equal to the present ones\ while for Saturn we took in some cases its current orbital elements\ while in other cases it started somewhat closer to the Sun to make allowance for its outward drift[ In the outer solar system\ a swarm of N bodies of equal mass was distributed in a ~at structure between 01 and 24 AU "in some runs the outer boundary was extended up to 39 AU#[ The adopted radial distribution was dN: daa−g with g 0:1\ which follows some models of the mass distribution in the protoplanetary disk "Wei! denschilling\ 0866#[ We adopted small initial eccentricities and inclinations randomly distributed with maximum values ranging from 9[994 up to 9[0[ The total initial mass of the entire swarm of bodies was de_ned as a factor f of the combined mass of the current Uranus and Neptune[ Because of the mentioned ine.ciency of the accretion process of the outer planets\ we adopted f × 0[ As the maximum number of bodies that can be handled by our computational scheme is limited\ we have per! formed several runs with di}erent initial number N of planetesimals "149\ 499 and 0999#\ in order to check the e}ect of the _nite number of bodies on the main outcomes of the simulations[ As the initial mass in the planetesimal swarm is typically of 1Ð2 times the combined mass of Uranus and Neptune\ the planetesimals have initial masses in the range 09−1Ð09−2 MU\ MU being the present mass of Uranus[ The time span covered by the numerical simulations was of 099 Myears in all the runs[
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
485
Duncan "0883#[ It was _ve time more time consuming than the ones without this e}ect[ We have carefully exam! ined all the outcomes of this simulation[ We have not found any noticeable di}erence with other simulations of the set[ Remarkably\ in all the runs there is at least one massive planet[ From our simulations we found that the accretion process by mutual collisions of the planetesimals usually evolved to the formation of two planets\ though in some cases we obtained one or three large planets[ In those cases where only one massive planet is formed\ its orbit lies beyond 29 AU[ Figure 2 shows the orbital evolution of the proto! Uranus and proto!Neptune for several representative runs[ The orbital migration\ reported in the previous simulations by Fernandez and Ip "0873\ 0885# is also found with our more realistic model\ although we observe a strong dynamical interaction between the protoplanets at early stages of the evolution[ An outward migration is also the case for Saturn[ It has in fact experienced a very strong displacement of several AU[ In most cases\ this strong migration is responsible for the destabilization of the proto!Uranus[ This is a limitation of our model\
2[ Results Table 0 displays the main results of 20 simulations with di}erent initial conditions[ The columns refer to] "0# number of run^ "1# initial number of bodies^ "2# the total mass initially present in terms of the combined masses of Uranus and Neptune " 20[7 M$#^ "3# an estimate of the maximum relative velocity in terms of the circular velocity at each helio! centric distance\ where i9 and e9 are the maximum possible initial inclination "in radians# and eccentricity for the planetesimals in each run^ "4# the maximum number of embryo planets that reach at least one Earth mass during the given run^ "5# the number of planets remaining in the outer solar system at the end of each run^ "6# and "7# the _nal masses of our computed {Uranus| and {Neptune| in terms of Uranus|s mass " 03[43 M$#^ "8# and "09# the time of the last macroaccretional event for proto!Uranus and proto!Neptune in Myears[ The _nal results of the simulations seem to be not strongly dependent on the initial number of plan! etesimals[ Run 20 was made considering the solar per! turbation during close encounters as in Levison and
Table 0 Initial conditions and main results of 20 numerical simulations Run
N
Minitial
ze19¦i19
Nemb
Npla
MU
MN
TU
TN
0 1 2 3 4 5 6 7 8 09 00 01 02 03 04 05 06 07 08 19 10 11 12 13 14 15 16 17 18 29 20
499 499 499 499 499 499 499 499 499 499 499 499 499 0999 0999 0999 0999 0999 0999 0999 0999 0999 0999 0999 0999 149 149 499 0999 0999 499
0[7 0[8 0[7 1[0 1[4 0[6 0[6 1[6 1[9 2[2 1[4 2[9 1[9 1[3 0[4 0[6 0[7 2[9 1[4 2[7 2[9 1[4 2[7 2[9 2[9 1[9 0[4 2[9 1[4 1[4 1[9
9[996 9[996 9[996 9[960 9[960 9[960 9[960 9[903 9[960 9[960 9[960 9[903 9[960 9[949 9[030 9[030 9[903 9[001 9[903 9[917 9[903 9[996 9[903 9[903 9[917 9[903 9[903 9[960 9[960 9[960 9[960
00 09 01 02 6 09 5 7 00 3 3 3 02 12 02 08 7 04 04 01 02 03 07 8 02 00 04 09 08 02 5
2 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 0 1 1 1
0[7 0[1 9[1 9[5 9[3 1[3 * 1[9 9[4 * 9[5 0[4 9[4 9[6 0[1 9[7 9[0 0[4 * 9[5 * 0[7 * 9[0 * 9[5 * * 0[6 9[8 9[6
9[2 0[1 0[0 1[1 0[1 9[3 0[3 0[0 9[6 3[1 2[4 1[0 9[2 0[4 9[4 9[2 0[2 0[6 0[6 0[4 9[4 9[5 0[1 0[9 9[7 9[5 0[6 9[6 0[9 1[9 9[3
56 00 63 68 70 88 * 22 78 * 71 69 21 71 74 57 89 5 * 79 * 39 * 57 * 35 * * 24 14 01
09 08 20 36 35 4 31 12 15 54 50 16 76 59 20 09 58 4 50 69 81 4 50 8 51 00 43 32 36 06 10
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
486
Fig[ 2[ Evolution of the semimajor axis of the outer planets for some of the runs where two massive planets were formed[
imposed by the inclusion of Saturn with its present mass and orbit[ We have also made some runs starting Saturn with a 7 AU without noticeable di}erences "runs 15Ð 29#[ The general rule in our simulation is that Jupiter drifts inwards[ Surprisingly for us\ in a few runs Jupiter also experienced a net outward migration[ This is explained as due to the ejection by Jupiter of a large number of
planetesimals with less angular momentum than they initially had[ The angular momentum lost by these ejected bodies is compensated in these cases by an outward drift of Jupiter[ This e}ect is not observed in the simulations by Fernandez and Ip because they have not considered a collision with the Sun as a possible end!state of a plan! etesimal[ If we followed the evolution of these objects "ignoring the possible collision# they would su}er a close
487
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
encounter with Jupiter at some later time\ being ejected from the system in hyperbolic orbits[ In this situation\ Jupiter|s orbital drift would be always inward\ as was previously found by Fernandez and Ip[ We observe that several protoplanets form during the process of accumulation of the outer planets[ Figure 3 exhibits the radial distribution of the protoplanets in one of the runs that resulted in the formation of two massive planets[ In short time scales\ of the order of few 096 years\ up to 00 protoplanets are competing for a place in the region beyond Saturn[ Many of them accrete to form more massive embryos[ It is worth noting that some large embryos reach semimajor axes in the region of the Kuiper belt "½39Ð59 AU# remaining there for several 095 years[ The dynamical and physical consequences of this phenomenon may be relevant for the evolution of the belt[ A large number of massive protoplanets has\ in some cases\ a destabilizing e}ect over the entire system\ resulting in a strong interaction between proto!Uranus and proto!Neptune[ The large jumps in semimajor axis are associated to encounters with these large proto! planets[ Our results suggest that there is a range of initial masses from which two planets like Uranus and Neptune can be formed\ that is in good agreement with the values previously found by Fernandez and Ip "0873\ 0885#[ There is a minimum amount of initial mass\ in the range of 1Ð2 times the combined masses of Uranus and Neptune\ required to form planets the size of Uranus and Neptune[ For the 11 runs with 1 outer planets\ an average of only ½35) of the solid material initially present was
Fig[ 3[ Evolution of the semimajor axis of the protoplanets formed during one of the runs "open circles#[ We de_ne as {protoplanets| those bodies reaching masses greater than one Earth mass[
incorporated into the planets[ Therefore\ the increasing probability of ejection of interacting planetesimals with the massive protoplanets during the latest stage of accretion is con_rmed by our model[ A large fraction of the bodies scattered by the accreting proto!Uranus and proto!Neptune fall under the gravi! tational control of Saturn and Jupiter[ Once this happens the dynamical evolution of the scattered bodies is char! acterized by a fast random!walk in perihelion distance q in such a way that a fraction of them become Mars crossers before being ejected in hyperbolic or near para! bolic orbits[ Some of them could _nally collide with Jup! iter or Saturn[ We found that 9[0Ð9[4 M$ of solid material from the UranusÐNeptune region was incorporated into Jupiter and into Saturn[ These values are in contrast with the numerical results previously found by Fernandez and Ip "0873#\ where up to 1 M$ of solid material was incor! porated into Jupiter\ although they started their simu! lations with comparable amount of mass[ We have not found any correlation between the initial mass and the total mass incorporated into Jupiter and Saturn[ Figure 4 displays the mean e.ciency\ averaging over all the runs\ of each planet in ejecting planetesimals from the system[ As expected\ Jupiter is the most e.cient ejec! tor\ although there is a non!negligible fraction of bodies ejected by the dynamical in~uence of Saturn\ in the region between ½7 and ½04 AU[ The large fraction of bodies interacting with Saturn "part of the material is ejected and part is handed down to Jupiter# is responsible for its strong orbital migration[ The distribution of perihelion distances of the ejected bodies displayed in Fig[ 4 shows that a large fraction is ejected at q ³ 4 AU\ con_rming our argument about the
Fig[ 4[ Relative distribution of the perihelion distances of those bodies ejected from the system^ mean values for 20 simulations[
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
origin of the outward migration of Jupiter in a few runs[ The distribution of perihelion distances of bodies reach! ing the region of the terrestrial planets shows a steep decrease towards smaller q\ so that only a small fraction of the bodies crossing Mars|s orbit also cross Mercury|s[ As shown in Fig[ 5 for a run with 499 bodies\ the total number of surviving planetesimals decreases expo! nentially with time[ We explained above\ the losses are due to accretion by the growing outer planets "the fast loss regime shown in Fig[ 5 during the _rst Myears#\ or ejection to interstellar space[ The cumulative mass reaching the region of the ter! restrial planets is found to be of the order of 09Ð04 Earth masses\ becoming potential impactors of these four planets[ The evolution of the cumulative mass entering within the zone r ³ 0[4 AU is displayed in Fig[ 6[ The solid line is the average value over our 20 simulations[ The dotted lines are the maximum and minimum values found in the simulations[ The phase of early strong bom! bardment of the terrestrial planets due to scattered plan! etesimals from the UranusÐNeptune region reached its maximum near ½096 years\ lasting\ but at a slower rate\ for ½3×096 years[ This is a much shorter time span than previously found "Fernandez and Ip\ 0885#[ About 19 times the amount of water contained in the oceans could have reached the Earth from this bombardment of icy material "Fernandez and Ip\ 0886#\ but this is a very uncertain result at this point[ The answer depends on whether or not the Earth and the other terrestrial planets were already formed when the in~ux of planetesimals from the outer planetary region reached its peak[ If they were still in the process of accretion most of the accreted volatile material might have been lost in subsequent megaimpacts[ The masses of proto!Uranus and proto!Neptune grow very fast[ Figure 7 displays the growth rate of proto!
Fig[ 5[ Evolution of the number of planetesimals for a run with initially 499 bodies[
488
Fig[ 6[ Cumulative mass\ in M$\ that reached the region of the terrestrial planets[ The solid line is the mean value for 20 simulations[ The dashed lines are upper and lower bounds to the cumulative mass[
Uranus and proto!Neptune for three runs[ We may observe that the time scale of complete planetary for! mation is of few 096 years\ either for Uranus or Neptune[ The slope of the curves exhibits the classical behavior of runaway growth "Ida\ 0889#[ The large radial excursions experienced by the protoplanets\ such as the ones shown in Fig[ 2\ bring them into zones not depleted of plan! etesimals\ thus substantially extending their feeding zones[ Moreover\ the orbits in these undepleted regions were not previously excited by the protoplanets\ and therefore\ the orbital migration seems to be an e.cient way to sustain the regime of runaway growth for long time spans[ Also noticeable in Fig[ 7 is that Neptune generally grows faster than Uranus\ reaching also a big! ger mass[ This fact may be related to the radial excursions that are in general larger for Neptune than for Uranus[ Figure 8 displays\ for one of the runs\ a combination of Figs 2 and 7\ where we illustrate the correlation between the orbital migration and the change in mass[ As previously noted\ the maximum initial values for the orbital eccentricity and inclination varies between 9[994 in some cases and the more realistic 9[0[ We have not noticed important di}erences arising from these di}erent initial conditions\ mainly because even when the initial eccentricity and inclination are very low\ the system evolves very fast to "e\ i# ½ 9[0[ The evolution of the r[m[s[ eccentricity and inclination for one of the runs is shown in Fig[ 09[ In one of the runs we have additionally computed the relative velocity and the distance of closest approach for all the encounters experienced by proto!Neptune[ The average value of the relative velocity in terms of the proto! Neptune orbital velocity is shown in Fig[ 00a[ Each value was computed averaging over 099 encounters[ It is inter! esting to note that the relative velocity increases fast and
599
A[ Brunini\ J[A[ Fernandez : Planetary and Space Science 36 "0888# 480Ð594
Fig[ 7[ Accretion rate of proto!Uranus and proto!Neptune for some of the runs[