Stability of the Outer Planets in Multiresonant Configurations with a Self-gravitating Planetesimal Disk
SSubmitted to ApJ, July, 2013
Stability of the Outer Planets in Multiresonant Configurationswith a Self-gravitating Planetesimal Disk
M. Reyes-Ruiz , H. Aceves , and C. E. Chavez [email protected] ABSTRACT
We study the effect of a massive planetesimal disk on the dynamical stabilityof the outer planets assuming, as has been suggested recently, that these wereinitially locked in a compact and multiresonant configuration as a result of gas-driven migration in a protoplanetary disk. The gravitational interaction among all bodies in our simulations is included self-consistently using the
Mercury6.5 code. Several initial multiresonant configurations and planetesimal disk modelsare considered. Under such conditions a strong dynamical instability, manifestedas a rapid giant planet migration and planetesimal disk dispersal, develops ona timescale of less than 40 Myr in most cases. Dynamical disk heating due tothe gravitational interactions among planetesimals leads to more frequent inter-actions between the planetesimals and the ice giants Uranus and Neptune, incomparison to models in which planetesitmal-planetesimal interactions are ne-glected. On account of the rapid evolution of the multiresonant configurationsobtained with fully self-consistent simulations, our results are inconsistent withthe dynamical instability origin of the Late Heavy Bombardment as currentlyconsidered by the Nice model for the Solar System.
Subject headings: planets and satellites: formation — planets and satellites:dynamical evolution and stability — planet-disk interactions Universidad Nacional Aut´onoma de M´exico, Instituto de Astronom´ıa, Apdo.Postal 106, Ensenada, B.C.22860 M´exico Facultad de Ingenier´ıa Mec´anica y El´ectrica, Universidad Aut´onoma de Nuevo Le´on, Monterrey, NuevoLe´on, 66451, M´exico a r X i v : . [ a s t r o - ph . E P ] J un
1. Introduction
The formation of the Solar System is generally believed to involve, during its last phases,a period of significant radial migration of the giant planets due to their interaction with amassive, remnant planetesimal disk located outward of the outermost planet (for a reviewsee Levison et al. 2007). This planetesimal-driven migration involves the inward migrationof Jupiter, the outward migration of Saturn, Uranus and Neptune, and the dispersal of mostof the planetesimals in the disc into orbits at much greater distances from the Sun.The current scenario has evolved from the original suggestion of Fernandez & Ip (1984),then developed by Malhotra (1993), in which the outer planets migrate as they exchangeangular momentum efficiently with the planetesimals in the outer Solar System while theplanetesimals are dispersed to the Oort cloud and beyond. Early numerical simulations ofthis process (Hahn & Malhotra 1999, Gomes et al. 2004) indicated that Jupiter migratesinward typically by a fraction of an AU, while Saturn, Uranus and Neptune migrate outwardschanging their initial semimajor axis by approximately 10, 50 and 100%, respectively.In a series of papers by Gomes et al. (2005), Tsiganis et al. (2005) and Morbidelliet al. (2005), a variant of the migration scenario proposed by Malhotra (1993, 1995) waspresented. As part of what is now known as the Nice model for the early evolution of the SolarSystem, an initially very compact configuration of the outer Solar System was proposed asthe condition at the time the primordial gaseous nebula was dispersed. A distinctive featureof the initial configuration originally proposed in the Nice model was the fact that it isunstable on a timescale of several hundred megayears, with Saturn crossing the 2:1 meanmotion resonance (MMR) with Jupiter, and leading to a major rearrangement of the outerSolar System.In recent years, the Nice model has undergone major revisions. The model now hasincorporated the idea that in the initial orbits of the outer planets are in a multi-resonantconfiguration, as expected following an apriori gas-driven migration phase (Masset & Snell-grove 2006, Morbidelli et al. 2007 and Batygin & Brown 2010, among others). Such con-figurations are known to be stable if only the gravitational interaction among the planets isconsidered (Batygin & Brown 2010). Secondly, as argued by Brasser et al. (2007), the sta-bility of the inner solar system during the period of outer planet migration seems to requirean abrupt change in Jupiter’s orbit (see also Minton & Malhotra 2009). This change canin principle result from a close encounter of Jupiter with an ice giant that may eventuallybe ejected from the solar system, leading to the so-called “jumping Jupiter” revision of theNice model (Morbidelli et al. 2009b, 2010, Brasser et al. 2009, Walsh & Morbidelli 2011,Agnor & Lin 2012). Finally, by using an approximate treatment for the self-gravitational in-teraction in the planetesimal disk (an effect not previously considered), Levison et al. (2011) 3 –suggested that, starting from a stable multi-resonant configuration, self-gravity instabilitiesslowly develop and can lead to the rapid migration of Neptune into the planetesimal disk, ona timescale consistent with the occurrence of the Late Heavy Bombardment (LHB) severalhundred megayears after gas disk dispersal.In this paper we present a series of N -body simulations that probably resemble thedynamical evolution the outer Solar System had after the primordial gaseous nebula wasdispersed, in order to study in some detail the stability of the system taking into account self-consistently the gravitational interactions among the planetesimals. To our knowledge,with the exception of the work by Levison et al. (2011), previous studies of planetesimaldriven migration have considered these small bodies to interact gravitationally only withthe planets but not among themselves. This approximation is only justified on account ofthe great computational expense of the direct N -body simulations required to model thisprocess. In this paper, we present results of the first numerical simulations to study theeffect of a self-gravitating disk on the process of planetesimal driven migration in a fullyself-consistent manner.The paper is organized as follows. In Section 2 we present the details of our model, wedescribe the methodology used and the cases studied. In § §
4. Finally, in Section 5 we summarize our main results and present ourconclusions.
2. Model Description
Models of the process of planetesimal driven migration generally assume that the pro-cess begins in earnest once the gaseous component of the protoplanetary disk is dispersedsomehow, at least to the point that dynamical effects of the gas can be neglected. In thisstudy we adopt this scenario in order to make a direct comparison to results from other au-thors. We note, however, that preliminary results of an ongoing investigation by our group,suggest that the remnant gaseous component of protoplanetary disks during the so-calledtransition disk phase, can have important consequences on the evolution of the planetesimalpopulation and the outer planets.We consider the gravitational force of the Sun, the 4 outer planets with their presentday mass, and a massive disk of much smaller bodies representing the remnant planetesimaldisk exterior to the planetary region. A novel feature of our simulations, in comparison to 4 –previous studies of planetesimal-driven migration, is the inclusion of the gravitational forceamong the planetesimals themselves in a fully self-consistent manner.
In accordance to the latest versions of the Nice model, based on the results of Masset& Snellgrove (2006) and other authors (e.g. Morbidelli et al. 2007, Nesvorny & Brown 2010and Zhang & Zhou 2010), we assume that initially the outer planets were locked in meanmotion resonance (MMR) with each other. Specifically, we adopt some of the multiresonantconfigurations described in Table 1 of Batygin & Brown (2010) as a frequent outcome oftheir study of resonant trapping in gaseous disks. In our notation, a 3J:2S, 4S:3U, 4U:3Nconfiguration is one in which Jupiter and Saturn are in a 3:2 MMR (i.e. Jupiter makes 3revolutions around the Sun while Saturn orbits twice), Saturn and Uranus are in a 4:3 MMRand Uranus and Neptune are also in a 4:3 MMR.The precise value of the main orbital parameters used for these initial conditions areshown in Table 2 (more details can be obtained by mailing the corresponding author). Theseconfigurations are dynamically stable (without a planetesimal disk) over gigayear timescales.As can be seen in Figure 1, which shows the evolution of the semimajor axis, a , and eccen-tricity, e , of each of the planets considered in case A of our simulations. The behavior of theother multiresonant planetary configurations we have considered without the planetesimalis very similar. For reference, we include in Table 2 the average value of this parameterscorresponding to the 3J:2S, 3S:2U, 4U:3N multiresonant system considered by Levison et al.(2011). Following Levison et al. (2011), who studied the effect of the planetesimal disk self-gravity on the equilibrium of the system for the first time in an approximate manner, weassume that initially the disk extends from a radius located outwards of the outermostplanet, R disk , to an outer radius approximately at the distance of Neptune’s current orbitalsemimajor axis, R out . The inner disk radius is taken as a parameter in our model withseveral values considered as given in Table 1. Note that the Hill radius for Neptune at thecorresponding distance in these compact configurations is approximately 0.31 AU, so that allvalues of R disk considered are more than 3 Hill radii from the outermost planet; i.e. outsideits so-called feeding zone. The outer radius for all cases considered is taken at 31.5 AU, in 5 –concordance with the study of Levison et al. (2011).On account of the computation time required for these simulations, the disk is restrictedto be made of N p = 2000 particles for most of our calculations. Although our presentcomputational capabilities do not permit considering a much greater number of particles,in the Discussion section we address the probable implications of this choice briefly. Themass of each particle is then m p = M disk /N p , where M disk is the total disk mass. For thedisk mass values adopted in this study (Table 1), the planetesimal mass ranges from ∼ M E or between 1 and 4 Lunar masses approximately; where M E is the mass ofthe Earth. Although this mass is still a factor of a few greater than the average mass ofthe dwarf planets, which can probably be considered the largest planetesimals, including agreater number of planetesimals, with a self consistent treatment of gravitational forces, isbeyond the limit of our present computational capabilities.The mass of particles we have considered is comparable to that in the studies of Hahn& Malhotra (1999), which is in the range 0.01 and 0.2 M E , and also of Gomes et al. (2004),who consider planetesimals of approximately 0.02 M E . In both cases however, planetesimal-planetesimal interactions were not considered. The planetesimal mass we have considered isalso comparable to the mass considered for gravitational perturbers in the study of Levisonet al. (2011).The mass surface density in the planetesimal disk, Σ, is assumed to be initially axisym-metric and to have a power law dependence on the radial distance from the Sun, Σ = Σ R − ,following Levison et al. (2011) and consistent with minimum mass solar nebula models. Thevalue of Σ is adopted to reproduce the total disk mass in a specific model. The planetesi-mal disk is assumed to be cold, with eccentricities, e , and inclinations, i (in radians) chosenrandomly from uniform number distributions in the range [0 , − ]. The rest of the orbitalelements, Ω, ω and M , are also chosen from a uniform random distribution within the range[0 , o ]. We use the
Mercury6.5 code, developed by Chambers (1999), to study the dynamicsof the four outer planets and the planetesimal disk. It uses a fixed time-step for the sym-plectic integrator coupled with a Bulirsch-Stoer (B-S) scheme, which is used to follow closeencounters among bodies: we set the error tolerance of the B-S integrator to 10 − . TheB-S scheme is used whenever the distance between two bodies is less than 3 Hill radii of thelargest body. 6 –Our simulations were in general carried out for up to 70 Myr. The relatively shorttimescale we have considered is related to the main purpose of our study, which is to assess thestability of these systems and not to determine the long-term final planetary configuration.The time-step used in our simulations is 186.4 days, as commonly adopted in outer solarsystem migration simulations (e.g. Hahn & Malhotra 1999). The simulations presented inthis paper were carried out in a Linux PC, with OpenSuse 11.4 as operating system and a gfortran -2.4 compiler.
3. Results
We have carried out a series of simulations aimed to understand the conditions underwhich planetesimal-driven migration, starting from multiresonant planetary configurations,becomes unstable. To this end, several initial conditions for the planets and for the plan-etesimal disk, were considered (see Table 1). In the following we present our results byconsidering the effect of varying the disk mass and inner radius for each of the multiresonantinitial conditions studied.For comparison, we have computed cases A, F and K (described in Table 1) with thesame disk properties considered by Levison et al. (2011) as their fiducial disk model; M disk =50 M E and R disk = 18 AU. However, the planetary configurations we consider (see Tables 1and 2), differ slightly from the 3J:2S, 3S:2U, 4U:3N multiresonant system that these authorsstudy. Levison et al. (2011) take their initial planetary configuration from Morbidelli etal. (2007), and report and average semimajor axis for Jupiter approximately 10% smallerthan in our cases. As a consequence, the initial semi-major axis of the outer planet in theconfiguration considered by Levison et al. (2011), is smaller in all but the 3J:2S,4S:3U,4U:3Nconfiguration we consider. This configuration is the most compact we studied and also theone that destabilizes at an earlier time for all cases.The main effect of planetesimal-planetesimal interactions on the evolution of the systemcan be seen comparing Figure 2 and Figure 3 which show the evolution of systems without andwith self-gravity in the planetesimal disk, respectively. Both figures shows several snapshotsof the a − e and a − i phase space domain, of a system corresponding to the configurationof case A. Without the effect of self-gravity, Figure 2, the system is not unstable on thetimescale considered, 10 Myr. The behavior of such system is consistent with the resultsof typical planetesimal-driven migration simulations starting from Nice model-like initialconditions. If we include the gravitational interaction among planetesimals, the systemundergoes a dynamical instability in approximately 5 Myr. Comparing the first and secondrow panels of Figure 3, which illustrate the system’s dynamical state initially and after 0.1 7 –Myr respectively, we see that one of the main effects of the inter-planetesimal gravitationalinteractions is the heating of the disk, i.e. the increase in orbital eccentricities and inclinationsalso referred to as gravitational stirring. We discuss this process in more detail, and itsdependence on the number of particles used in our simulation, in the following section.As can be seen in Figure 3, the disk remains in a quiescent state for a few Myr (for mostcases considered) as the number of planetesimals dispersed into orbits capable of having closeencounters with the outermost planets increases. The cumulative effect of these interactions,exchanging angular momentum with Neptune and Uranus mainly, leads to a growth of theorbital eccentricity for these planets until the multiresonant state among the giant planetsis broken. The instability manifests as a very fast orbital rearrangement of all planets. Theoutward migration of Uranus and Neptune in turn leads to a greater number of planetesimal-planet interactions and even faster migration into the disk until it is mostly dispersed. Inseveral cases studied, the planetary system is strongly unstable with one or both ice giantsejected from the system, or even catastrophic as with one of them colliding with the giantplanets as indicated by the numerical simulations. In order to assess the effect of the planetesimal disk mass on the evolution of each of ourinitial multiresonant configurations, we have carried out simulations with M disk = 25 ,
50 and100 M ⊕ . In all cases, we consider the same inner and outer disk radii of the fiducial modelof Levison et al. (2011). Disk mass is varied by changing the mass of each of the particlesrepresenting the planetesimals keeping the number of planetesimals constant.The time evolution of the planetary orbital parameters, a and e , for all three mul-tiresonant configurations we considered are shown in Figures 4 and 5, respectively. In thesefigures each row corresponds to a different planetary multiresonant initial condition. We findthat, as expected, the instability arises earlier for more massive disks in which gravitationalstirring increases the amount (and cumulative mass) of planetesimals interacting with theplanets.With the exception of case G, all systems considered are unstable on a timescale shorterthan 42 Myr, as a result of the interaction with the self-gravitating planetesimal disk. Case Galso becomes unstable but on a slightly greater timescale, less than 70 Myr. It is worth notingthat case G corresponds to the combination of: (a) the less massive disk we have consideredand (b) the planetary configuration in which Uranus and Neptune are most separated. 8 – We have also studied the effect of the location of the inner radius of the planetesimaldisk, R disk , on the stability of these systems. To do so, we fix the disk mass at 50 M ⊕ ,following Levison et al. (2011), and vary R disk between 16 and 20 AU. The results for allinitial planetary configurations as a function of R disk are displayed in Figures 6 and 7, whichshow the evolution of the semimajor axis, a , and eccentricity, e , respectively, for each of theplanets in our simulations.We find that the time for the onset of the instability depends strongly on the value of R disk . In going from an inner disk radius of 16 AU to 20 AU (a 25% increment), the timefor onset of the instability, t inst , increases by a factor of about 2-6. This result is reminiscentof the results of Gomes et al. (2005) who find that the initial location of the inner diskradius must be fine tuned to a specific value in order to reproduce the LHB timescale, if thisphenomenon is to be related to a planetary instability. Once again, as in our simulationswith different disk mass, all cases considered with a 50 M ⊕ planetesimal disk, the outerSolar System becomes unstable within a time scale of 30 Myr. For comparison, we havecarried out simulations of the three cases with different inner radius, for the 3J:2S, 4S:3U,4U:3N configuration and the same initial distribution of planetesimals, but without theeffect of planetesimal-planetesimal interactions. In all cases, the system remains stable for atimescale an order of magnitude greater, or more, than that with selfgravitating planetesimaldisks. An additional feature of the evolution of the unstable systems, is the fact that theplanets remain locked in resonance until the instability arises. This is evident in Figure 8which shows the mean motion resonant angles which, for the 3J:2S, 4S:3U, and 4U:3N initialplanetary configuration, can be written as: θ = 2 λ J − λ S + (cid:36) S , (1) θ = 2 λ J − λ S + (cid:36) J ,θ = 3 λ S − λ U + (cid:36) U ,θ = 3 λ S − λ U + (cid:36) S ,θ = 3 λ U − λ N + (cid:36) N and θ = 3 λ U − λ N + (cid:36) U , λ i = M i + ω i + Ω i is the mean longitude for the i -th planet, with M i being the meananomaly, ω i the argument of the planet’s periapsis and Ω i the longitude of the ascendingnode. Also, (cid:36) i is the longitude of the orbit’s periapsis.The behavior seen in Figure 8 is qualitatively similar for all systems before the instabilitydevelops, namely: (a) the θ and θ resonant arguments for Jupiter and Saturn librate about0 and 180 o , respectively. (b) Either the θ or θ resonant arguments for Saturn and Uranusclearly librate while the other argument may circulate, and (c) at least one of the resonantarguments θ or θ for Uranus and Neptune librates.The development of the dynamical instability is similar for all systems studied andinvolves a gradual increase in the number of planetesimals moving inside the disk innerboundary, R disk , and interacting strongly with the planets. This leads to a gradual increaseof the orbital eccentricities of Uranus and Neptune to values such that the MMR of theseplanets is broken and they can have close encounters. This is illustrated in Figure 9, whichshows the time evolution of the number of planetesimals in Neptune crossing orbits, i.e.planetesimals with perihelia less than Neptune’s aphelia (at a given time), for 3 of the caseswe simulated corresponding to disks with M disk = 50 M E and R disk = 18 AU. For comparison,we also show the behavior of Neptune’s perihelia, Uranus aphelia and the eccentricity of theirorbits.The cumulative effect of the interaction of the dispersed planetesimals (as a result ofself-gravity) with the ice giants leads, as the system approaches instability, to an increase intheir eccentricity until their resonant state is broken. At such point their semi-major axesgenerally increase in an abrupt manner (as shown in Figure 9), suggests that a thresholdin the interaction between planetesimals and the ice giants is triggered as the number ofthese gradually increases. Roughly, for the cases shown in Figure 9, the total mass that hasinteracted with Neptune at the time of the onset of the instability, t inst , is between 2 and3 Earth masses. This provides a possible explanation to the fact that disks with a masssmaller than the mass of the ice giants are not capable of destabilizing the system on thetimescales considered in this study.Figure 10 shows the evolution of the planets in phase space ( a − e ) as the system crossesthe instability. The interaction with the planetesimal disk leads to a gradual increase in theeccentricity of the ice giants (related to angular momentum) while keeping the semimajoraxis (energy related) more or less constant. Also shown are the boundaries of the resonancesof Saturn, Uranus and Neptune with their interior planet and of Jupiter with the exteriorplanet. These are estimated by assuming that each planet is the test particle in the contextof the restricted 3-body problem which, in the so-called pendulum approximations, allowsone to estimate the width of a given first order resonance as described in Holmes et al. 10 –(2003). As can be seen in Figure 10, the instability occurs when either Uranus or Neptuneevolve toward a dynamical state in which either is not locked in resonance, resulting in afast disruption of the planetary configuration.
4. Discussion
In this section we revisit some of the basic assumptions made in our study and how theseaffect the interpretation of our results. In addition, we discuss some of the implications ofour results on the theory of planetary system evolution.
The simulations presented above are intended to demonstrate that the inclusion of thegravitational interaction between the planetesimals and planets in a compact and multireso-nant initial configuration, as those proposed as part of the Nice model, renders the systemshighly unstable. Implicit in this statement is the idea that, without self-gravity in the plan-etesimal disk, the planets in these systems would remain locked in resonance for a very longtime, as suggested by Levison et al. (2011).To test the above, we carried out numerical simulations of some of the initial multires-onant conditions (specified in Tables 1 and 2) but removing the effect of self gravity, i.e.assuming that the planetesimals interact gravitationally with the planets but not amongthemselves (as usually done in studies of planetesimal-driven migration). We find that someof these systems are also unstable, but on a timescale of several tens to hundreds of Myr.Out of the 3 systems tested in this manner (cases A, F and K), one is found to be unstable ona timescale of less than 100 Myr. As we are focusing on the effect of disk self-gravity in thisstudy, we have decided to leave a further exploration of this issue for future contributions.Our main dynamical result brings into question the presumed stability of multireso-nant configurations and suggests that early planetesimal driven migration, and outer diskdispersal, is to be expected for a wide range of initial planetary system configurations.
The previous results indicate that, when gravitational interactions among the planetes-imals are considered, several, perhaps most, of the multiresonant planetary configurations 11 –can become unstable on a timescale of less than a few tens of Myr. This posses a seriousdifficulty for Solar System formation scenarios in which the LHB is explained as resultingfrom the unstable, fast migration of the outer planets and the consequent perturbation ofplanetesimal disk, as proposed for example in the Nice model.Given the early development of the instability of multiresonant planets with a self-gravitating planetesimal disk, some of the alternative scenarios that could be envisioned forthe LHB are the following.(a) That the disk inner radius was at a distance much greater than 20 AU when thegaseous component of the solar nebula was dispersed, which marks the beginning of oursimulated dynamical evolution. As can be seen from Figure 12, which shows the time forthe development of the instability as a function of the inner disk radius of the planetesimaldisk, a value of R disk much greater than 20 AU is required to have an instability developmenttime of t inst ≈
700 Myr; that is timescale of the LHB.This initial location of the inner boundary of the planetesimal disk, at such a greatdistance from the outermost planet, poses in our view a difficult conundrum. While it isgenerally believed that planets as they form are able to clear material (gas and planetesimals)from an annulus in the vicinity of their location, the so-called feeding zone (Lissauer, 1993)is thought to extend only to a few Hill radii from each planet’s position. In the configurationsuggested above, the inner radius of the planetesimal disk is at a distance of more than 30Hill radii from the outermost planet. It is not evident why planetesimals were cleared, orwere not formed in the first place, only between approximately 12 and 25 AU.(b) Another possible explanation to obtain a long t inst , is that the initial mass of theplanetesimal disk was much smaller than is usually considered. Nice model studies haveregularly used M disk around 35 M ⊕ , and the latest incarnation of the model favors a slightlyheavier disk of 50 M ⊕ (Levison et al. 2011). However, we find that even a disk with 25 M ⊕ can destabilize the system on a timescale more than an order of magnitude shorter than thatrequired for the LHB.Using a power law fit to the resulting timescale for the development of the instabilityas a function of the disk mass, determined from our simulations as shown in Figure 11,suggests that M disk (cid:28) M ⊕ is required if the LHB is to be explained in terms of orbitaldynamical instability. However, one would expect that very low mass disks may be altogetherunable to destabilize the planetary system, as the mass and angular momentum containedin the planetesimal disk is smaller than that of any of the outer planets. We can estimate alower limit for the disk mass required to significantly modify Neptune’s orbit by comparingthe initial angular momentum of the planet, L N ≈ M N V K ( a N ) a N , with the total angular 12 –momentum in the disk, L disk ≈ M disk V K ( R out ) R out . The requirement that L disk > L N , leadsto the condition M disk (cid:38) M E for the planetesimal disk to destabilize the outermost planet.(c) It could also be possible that the Solar System did not start from a multireso-nant, compact configuration as it is generally proposed in the latest incarnations of the Nicemodel. Multiresonant configurations are believed to be the natural end product of planetrearrangement due to gas-driven migration in the primeval gaseous solar nebula. The pre-sumed stability of multiresonant configurations is further invoked to explain the long delayfor the planetary configuration to undergo a fast dynamical instability and the consequentstirring of the outer planetesimal disk and the LHB (Levison et al. 2011). However, as wehave shown in this work, this belief is not entirely justified. In fact, our results suggest thatseveral compact configurations are not stable over a 100 Myr timescale when the effects ofa self-gravitating planetesimal disk, of structural properties similar to those assumed in theframework of the Nice model, are taken into account.One could envision that the timescale necessary for the formation of all the giant plan-ets and their locking into a multiresonant configuration did not exactly coincide with thetimescale for the dispersal of the gaseous component of protoplanetary disks. This is par-ticularly true if gas dispersal is related to photoevaporation of the disk or any other processunrelated to the formation of the planets themselves. In such case, the early planetary sys-tem may be left free of gas so that no further gas-driven migration is possible before theplanets reach the proposed compact, multiresonant configuration.The subsequent evolution of these initially compact systems caught outside of full mul-tiresonance at the time of gas dispersal, is beyond the scope of the present study. However,one would expect that if the configuration is compact enough, the evolution could be similarto that of the original Nice model as presented in Gomes et al. (2005) and Tsiganis et al.(2005). On the other hand, if gas driven migration stops while the ice giants are located farfrom Jupiter and Saturn, the evolution can be expected to proceed in a relatively smoothmanner as studied by Hahn & Malhotra (1999) or Gomes et al. (2004). It must be stressedhowever, that in these studies the effect of planetesimal-planetesimal interactions was notconsidered.(d) Finally, one further possibility consistent with the present study, is that the unstablerearrangement of the solar system is not linked to the LHB. In this scenario the Solar Systemunderwent an unstable rearrangement shortly after gas was dispersed, on a timescale of 10’sof Myr, so that the much later LHB is not related to a planetesimal-driven orbital instability.Although there are strong arguments against alternative scenarios for the origin of the LHB,such as that proposed by Chambers (2007) involving the effect of a small planet originallypresent between Mars and Jupiter (e.g. Brasser & Morbidelli 2011), it is worth noting that 13 –such scenario can not be entirely ruled out as a particular solution on dynamical grounds,albeit of a low probability. This situation, is reminiscent of the jumping Jupiter scenarioinvoked in the current version of the Nice model. The important implications of our study, which may affect our understanding of theevolution of the Solar System, hints toward the importance of careful modeling of the diskof planetesimals in order put on firmer grounds the previous conclusion.Our numerical approach does not allow us, due to computational limitations, to followthe evolution of these planetary systems in a self-consistent way for a LHB-like timescale ofseveral hundred million years. However, a simple extrapolation of our instability results, asa function of disk mass and inner radius, suggests that if the LHB is to be explained in termsof a dynamical instability of the early Solar System, the initial outer planetesimal disk musthave had been much less than 20 Earth masses and/or the disk inner radius must have beenlocated at more than 10 AU from the outermost planet. Either condition has importantimplications for the preceding stages in the evolution of the Solar System. An even moreradical solution to this conundrum is the possibility that the LHB is not directly related toa dynamical instability of the outer planets.Although we consider our calculations to represent the state of the art in the studyof planetesimal-driven migration, on account of their self-consistent treatment, it is worthpointing out some factors related to our present computational limitations that may affectour conclusions. In our study each particle has a mass ranging from (0.75–3) × g thatis ≈ (4 −
17) times the mass of Eris, the largest of the dwarf planets. Considering thatthe present mass of the Kuiper Belt is believed to be ≈ /
500 of the mass assumed to beoriginally present in the outer planetesimal disk (Gladman et al. 2001), and assuming thatthe size distribution of large planetesimals was originally similar to the present one in theKB, one would expect that the original number of dwarf planet-like bodies to be in thevicinity of a few times 10 . This implies that each of the small bodies in our simulationscorresponds to approximately 5 dwarf planet-like bodies originally believed to be present inthe Solar System.The above approximation does not seem too bad, but it is known that the number ofparticles, N , in an N -body simulation affects the dynamical evolution of self-gravitatingsystems (e.g. Binney & Tremaine 1998). It is possible that an order of magnitude increasein the number of small bodies, as predicted by the present day size distribution of KBOs, 14 –can have an important effect on the dynamics of the Solar System. For example, on onehand, the presence of a very large number of small bodies allows for a further distributionof the interaction-energy among these by incrementing the number of degrees of freedom ofthe whole system. In this way, the small bodies may behave as a stabilizing sink, that isprobably dependent on their number, for the development of the instability as is suggestedby the results of this paper. On the other hand, the actual behavior of the disk particles isdependent on the number of particles used to describe it in what is usually called numericaldisk-heating. In this process the velocity dispersion of the disk particles, σ p , modeled as adiffusion process, grows with time t as ∝ ( σ + D p t ) p ; where D p is a diffusion coefficient, p a parameter and σ a constant (e.g. Lacey 1991, Vel´azquez 2005). The values of D p and p are model dependent and they have to be evaluated numerically under specific physicalconditions. However they show a dependency on the number of particles in the system, forexample, D p ∝ N − α , with α ≈ .
5. Conclusions
We have carried out a series of direct numerical simulations of the the dynamical evo-lution of the early outer Solar System in the context of the Nice-2 model. Our aim has beento study this process, for the first time, in a fully self-consistent manner, i.e. by taking intoaccount the gravitational interaction between all bodies included in the simulations.Our starting point was a set of initially compact, multiresonant planetary configura-tions, similar to those proposed in the study of Levison et al. (2011). In difference withprevious attempts to model such systems in the framework of an approximate treatment ofthe N-body problem (e.g. Levison et al. 2011), we find that a dynamical instability of theouter planets develops on a short timescale, of less than 70 Myr in all the initial configura-tions we explored. The early development of the instability is directly related to the effect ofplanetesimal-planetesimal interactions, as the resulting dynamical heating of the planetesi-mal disk accelerates the interaction with the outer ice giants which leads to the planetaryinstability.The above result has important consequences on our understanding of the dynamics ofthe early Solar System. Our results suggest that when the inter-planetesimal gravitationalinteractions are included in a fully self-consistent manner, the dynamical instability of the 15 –early Solar System invoked in the context of the Nice model to explain the LHB, occurs muchsooner than the several hundred Myr timescale indicated by the lunar crater record (Teraet al. 1974). Although further studies are required, particularly increasing the number ofparticles representing the planetesimal disk, these results bring into question the explanationof the LHB in terms of a planetary dynamical instability.The authors are grateful to Dr. Konstantin Batygin for kindly providing the multireso-nant initial conditions used in this study, and to Prof. Renu Malhotra for valuable commentson a draft of this paper. MRR acknowledges financial support from DGAPA-UNAM projectIN115429, HA from CONACyT project 179662, and all the authors from DGAPA-UNAMproject IN108914..
REFERENCES
Binney J. & Tremaine S. 2008, Galactic Dynamics. Princeton U. Press. Princeton, NJBatygin, K. & Brown, M. E. 2010, ApJ, 716, 1323Chambers, J.E. 1999, MNRAS, 304, 793Chambers, J.E. 2007, Icarus, 189, 386Fernandez, J. A. & Ip, W.-H. 1984, Icarus, 58, 109Gladman, et al., 2001, AJ, 122, 1051Gomes, R., Morbidelli, A. & Levison, H. F. 2004, Icarus, 170, 492Gomes, R., Levison, H. F., Tsiganis, K. & Morbidelli, A. 2005, Nature, 435, 466Hahn, J.M. & Malhotra, R. 1999, AJ, 117, 3041Hartmann, W. K., Ryder, G., Dones, L. & Grinspoon, D. 2000, in Origin of the Earthand Moon, eds. R.M. Canup & K. Righter. Tucson: University of Arizona Press.,p.493-512Holmes, E.K., Dermott, S.F., Gustafson, B..S., Grogan, K. 2003, ApJ, 597, 1211Lacey, C.G., 1991, in Dynamics of Disc Galaxies, ed. B. Sundelius. G¨oteborg, Sweeden, p.257 16 –Levison, H.F., Morbidelli, A., Gomes, R. & Backman, D. 2007, in Protostars and PlanetsV, eds. B. Reipurth, D. Jewitt & K. Keil. Tucson: University of Arizona Press., p.669-683.Levison, H. F., Morbidelli, A., Tsiganis, K., Nesvorny, D. & Gomes, R. 2011, AJ, 142, 152Lissauer, J. J. 1993, ARA&A, 31, 129Lissauer, J. J. & de Pater, I. 2013, Fundamental Planetary Science. Cambridge, UK: Cam-bridge University PressMalhotra, R., 1993, Nature, 365, 819Malhotra, R., 1995, AJ, 110, 420Morbidelli, A., Levison, H. F., Tsiganis, K. & Gomes, R. 2005, Nature, 435, 462Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F. & Gomes, R. 2007, AJ, 134, 1790Nesvorny, D. & Morbidelli, A. 2012, AJ, 144, 20Tera, F., Papanastassiou, D.D. & Wasserburg, G.J. , 1974, Earth Planet. Sci. Lett., 22, 1Tsiganis, K., Gomes, R. Morbidelli & A., Levison, H. F. 2005, Nature, 435, 459Vel´azquez H. 2005, RevMAA, 41, 389Zhang, H. & Zhou, J.-L. 2010, ApJ, 714, 532
This preprint was prepared with the AAS L A TEX macros v5.2.
17 –Fig. 1.— Evolution of the orbital elements, a and e , and of the resonant angles given byEqn. (1) for the planets in the 3J:3S,4S:3U,4U:3N multiresonant configuration used in oursimulations, without a planetesimal disk. The orbits are coplanar throughout the evolution. 18 –Fig. 2.— Snapshots of the evolution of the orbital configuration of a system, a versus e (left column) and a versus i (right column), for a disk as in case A described in Table 1 butNOT considering the effect of planetesimal-planetesimal interactions. It illustrates a typicalevolution of systems with Nice-like initial conditions which remain “stable” for 10 Myr.Small gray circles represent the planetesimals, black circles corresponding to the planets arelabeled by their initial. 19 –Fig. 3.— Same as Figure 2 but for a disk with planetesimal-planetesimal interactions corre-sponding to case A described in Table 1. It illustrates a typical unstable evolution of fullyself-gravitating systems which become unstable on a short timescale. 20 –Fig. 4.— Evolution of the planetary orbital semimajor axis for each of the planetesimal diskmasses considered. Columns correspond to different multiresonant planetary configurations.The left column corresponds to cases B, A and C, the middle column to cases G, F and H,and the right-hand column illustrates the for cases L, K and M. 21 –Fig. 5.— Evolution of the planetary orbital eccentricity for each of the planetesimal diskmasses considered. The left column corresponds to cases B, A and C, the middle column tocases G, F and H, and the right-hand column illustrates the for cases L, K and M. 22 –Fig. 6.— Evolution of the planetary orbital semimajor axis for each of the planetesimal diskinner radii considered. Columns correspond to different multiresonant planetary configura-tions. The left column corresponds to cases D, A and E, the middle column to cases I, Fand J, and the right-hand column illustrates the for cases N, K and O. 23 –Fig. 7.— Evolution of the planetary orbital eccentricity for each of the planetesimal diskinner radii considered. The left column corresponds to cases D, A and E, the middle columnto cases I, F and J, and the right-hand column illustrates the for cases N, K and O. 24 –Fig. 8.— Time evolution of the 2 body resonant arguments for each pair of planets asdefined in Eqs. (1) for case A. The top panel corresponds to the 3:2 MMR resonance betweenJupiter and Saturn, the middle panel to the 4S:3U resonance and the bottom panel showsthe arguments of the 4U:3N MMR. In each panel, the black dots denotes the first of the 2arguments shown (e.g. θ in the top panel) and the gray dots show the second argument forthat resonance (e.g. θ in the top panel) . 25 –Fig. 9.— Time evolution of the number of planetesimals scattered onto Neptune crossingorbits for cases A, F and K (black lines) in the top, middle and bottom panel respectively.Also shown is the evolution of the semimajor axis of Neptune in each case, indicating therelation between these 2 processes. 26 –Fig. 10.— Semi-major axis versus eccentricity diagram for each of the planets as the systemapproaches the instability. The solid black lines correspond to the width of the resonancefor each planet with the interior planet. The diamond symbol indicates the last location ofeach planet before the instability breaks out, the square indicates the first location (in a − e space) after the instability onset. 27 –Fig. 11.— Time at which the instability develops, t inst , as a function of initial disk mass.The black line corresponds to the systems starting from the 3:2, 4:3, 4:3 multiresonantconfiguration, the light gray line corresponds to the 3:2, 4:3, 3:2 initial configurations andthe dark gray line to the systems starting from the 3:2, 4:3, 3:2 configuration. 28 –Fig. 12.— Time at which the instability develops, t inst , as a function of initial disk innerradius. The black line corresponds to the systems starting from the 3:2, 4:3, 4:3 multiresonantconfiguration, the light gray line corresponds to the 3:2, 4:3, 3:2 initial configurations andthe dark gray line to the systems starting from the 3:2, 4:3, 3:2 configuration. 29 –Table 1: Planet configuration and disk parameters of each of the cases studiedCase Planet Configuration M disk R disk J:S , S:U, U:N M ⊕ R ⊕ A 3:2, 4:3, 4:3 50 18B 3:2, 4:3, 4:3 25 18C 3:2, 4:3, 4:3 100 18D 3:2, 4:3, 4:3 50 16E 3:2, 4:3, 4:3 50 20F 3:2, 4:3, 3:2 50 18G 3:2, 4:3, 3:2 25 18H 3:2, 4:3, 3:2 100 18I 3:2, 4:3, 3:2 50 16J 3:2, 4:3, 3:2 50 20K 3:2, 3:2, 5:4 50 18L 3:2, 3:2, 5:4 25 18M 3:2, 3:2, 5:4 100 18N 3:2, 3:2, 5:4 50 16O 3:2, 3:2, 5:4 50 20 30 –Table 2: Initial value of the main orbital elements for the planets in the multiresonant, initialconfigurations we consideredPlanet a e i