We study the effect of a massive planetesimal disk on the dynamical stability of the outer planets in a system representing the early solar system assuming, as has been suggested recently, that these planets were initially locked in a compact and multiresonant configuration as a result of gas-driven migration in a protoplanetary disk. The planetesimal disk is represented by an ensemble of 2000 lunar mass bodies for which the gravitational interaction is calculated *self-consistently* using the Mercury6.5 code. Several initial multiresonant configurations and planetesimal disk models are considered. Under such conditions a strong dynamical instability, manifested as a rapid giant planet migration and planetesimal disk dispersal, develops on a timescale of less than 40 Myr in most cases. Dynamical disk heating due to the gravitational interactions among planetesimals leads to more frequent interactions between the planetesimals and the ice giants, in comparison to models in which planetesimal–planetesimal interactions are neglected. The number of particles used to represent the planetesimal disk has implications for our results, and although our studies represent the first self-consistent calculations of unstable planetesimal-driven migration, our results point toward the need for using more realistic treatments of the planetesimal disk. Finally, in the framework of our model, we discuss the possible implications of our results on the early evolution of the solar system.

## 1.INTRODUCTION

The last phases of the formation of the solar system are generally believed to have involved a period of significant radial migration of the giant planets due to their interaction with a massive, remnant planetesimal disk located outward of the outermost planet (for a review, see Levison et al. 2007 pp 669–683). This planetesimal-driven migration involves the inward migration of Jupiter; the outward migration of Saturn, Uranus, and Neptune; and the dispersal of most of the planetesimals in the disk into orbits at much greater distances from the Sun.

The current scenario has evolved from the original suggestion of Fernandez & Ip (1984), which was then developed by Malhotra (1993), in which the outer planets migrated as they efficiently exchanged angular momentum with the planetesimals in the outer solar system as the planetesimals were dispersed to the Oort Cloud and beyond. Early numerical simulations of this process (Hahn & Malhotra 1999; Gomes et al. 2004) indicated that Jupiter migrated inward typically by a fraction of an AU, while Saturn, Uranus, and Neptune migrated outward, changing their initial semimajor axes by approximately 10%, 50%, and 100%, respectively.

In a series of papers by Gomes et al. (2005), Tsiganis et al. (2005), and Morbidelli et al. (2005), a variant of the migration scenario proposed by Malhotra (1993, 1995) was presented. As part of what is now known as the Nice model for the early evolution of the solar system, an initially very compact configuration of the outer solar system was proposed as the condition at the time the primordial gaseous nebula was dispersed. A distinctive feature of the initial configuration originally proposed in the Nice model was the fact that it was unstable on a timescale of several hundred megayears, with Saturn crossing the 2:1 mean motion resonance (MMR) with Jupiter, leading to a major rearrangement of the outer solar system.

In recent years, the Nice model has undergone major revisions. The current model has incorporated the idea that the initial orbits of the outer planets are in a multiresonant configuration, as expected following an a priori gas-driven migration phase (Masset & Snellgrove 2006, Morbidelli et al. 2007 and Batygin & Brown 2010, among others). Such configurations are known to be stable only if the gravitational interaction among the planets is considered (Batygin & Brown 2010). Second, as argued by Brasser et al. (2007), the stability of the inner solar system during the period of outer planet migration seems to require an abrupt change in Jupiter's orbit (see also Minton & Malhotra 2009). In principle, this change could have resulted from a close encounter between Jupiter and an ice giant that was eventually ejected from the solar system, leading to the so-called "jumping Jupiter" revision of the Nice 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 interaction in the planetesimal disk (an effect not previously considered), Levison et al. (2011) suggested that, starting from a stable multiresonant configuration, self-gravity instabilities slowly develop and can lead to the rapid migration of Neptune into the planetesimal disk on a timescale consistent with the occurrence of the late heavy bombardment (LHB) several hundred megayears after dispersal of the gas disk (see Hartmann et al. 2007).

In this paper we present a series of *N*-body simulations that probably resemble the dynamical evolution the outer solar system experienced after the primordial gaseous nebula was dispersed in order to study in some detail the stability of the system, taking into account the *self-consistent* gravitational interactions among the planetesimals. To our knowledge, with the exception of the work by Levison et al. (2011) and Li et al. (2011), previous studies of planetesimal-driven migration in the outer solar system have considered these planetesimals to interact gravitationally only with the planets but not among themselves. This approximation is only justified on account of the great computational expense of the direct *N*-body simulations required to model this process. In contrast to the study we present, the work of Li et al. (2011) considers the evolution of a system starting from a planetary configuration in which the planets are not in resonance. Such an initial condition leads to stable planetesimal-driven migration such as that studied in Hanh & Malhotra (1999) and Gomes et al. (2004). In agreement with our results, Li et al. (2011) find that gravitational stirring leads to a greater number of planetesimals interacting with Neptune, thus accelerating and extending its migration through the planetesimal disk. In this paper, we present results of the first numerical simulations to study the effect of a self-gravitating disk on the process of unstable planetesimal-driven migration in a fully self-consistent manner.

The paper is organized as follows. In Section 2 we present the details of our model and we describe the methodology used and the cases studied. In Section 3 we present the results of a series of simulations with different model parameters. A review of some of the assumptions made in our study and some of the implications for the scenarios of planetary system evolution are discussed in Section 4. Finally, in Section 5 we summarize our main results and present our conclusions.

## 2.MODEL DESCRIPTION

Models of the process of planetesimal-driven migration generally assume that the process begins in earnest once the gaseous component of the protoplanetary disk is dispersed somehow, at least to the point that dynamical effects of the gas can be neglected. In this study we adopt this scenario in order to make a direct comparison to results from other authors.

We consider the gravitational force of the Sun, the four outer planets with their present day mass, and a massive disk of much smaller bodies representing the remnant planetesimal disk exterior to the planetary region. A novel feature of our simulations, in comparison to previous studies of planetesimal-driven migration, is the inclusion of the gravitational force among the planetesimals themselves in a fully self-consistent manner.

### 2.1.Planetary Configuration

In accordance with the latest versions of the Nice model based on the results of Masset & Snellgrove (2006) and other authors (e.g., Morbidelli et al. 2007, Batygin & Brown 2010 and Zhang & Zhou 2010), we assume that initially the outer planets were locked in MMR with each other. Specifically, in Table 1 we adopt some of the multiresonant configurations described by Batygin & Brown (2010) as a frequent outcome of their study of resonant trapping in gaseous disks. In our notation, a 3J:2S, 4S:3U, 4U:3N configuration is one in which Jupiter and Saturn are in a 3:2 MMR (i.e., Jupiter makes three revolutions around the Sun while Saturn orbits twice), Saturn and Uranus are in a 4:3 MMR and Uranus and Neptune are also in a 4:3 MMR.

**Table 1.**Planet Configuration and Disk Parameters in Simulations

Case | Planet Configuration | M_{disk} | R_{disk} |
---|---|---|---|

J:S , S:U, U:N | M | R | |

A | 3:2, 4:3, 4:3 | 50 | 18 |

B | 3:2, 4:3, 4:3 | 25 | 18 |

C | 3:2, 4:3, 4:3 | 100 | 18 |

D | 3:2, 4:3, 4:3 | 50 | 16 |

E | 3:2, 4:3, 4:3 | 50 | 20 |

F | 3:2, 4:3, 3:2 | 50 | 18 |

G | 3:2, 4:3, 3:2 | 25 | 18 |

H | 3:2, 4:3, 3:2 | 100 | 18 |

I | 3:2, 4:3, 3:2 | 50 | 16 |

J | 3:2, 4:3, 3:2 | 50 | 20 |

K | 3:2, 3:2, 5:4 | 50 | 18 |

L | 3:2, 3:2, 5:4 | 25 | 18 |

M | 3:2, 3:2, 5:4 | 100 | 18 |

N | 3:2, 3:2, 5:4 | 50 | 16 |

O | 3:2, 3:2, 5:4 | 50 | 20 |

Download table as: ASCIITypeset image

The precise value of the main orbital parameters used for these initial conditions are shown in Table 2 (more details can be obtained by contacting the corresponding author). These configurations 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 eccentricity, *e*, of each of the planets considered in case A of our simulations. The behavior of the other multiresonant planetary configurations that we have considered without the planetesimals is very similar. For reference, we include in Table 2 the average value of these parameters corresponding to the 3J:2S, 3S:2U, 4U:3N multiresonant system considered by Levison et al. (2011).

**Table 2.**Initial Values of the Main Orbital Elements for the Planets in the Multiresonant Configurations Used

Planet | a | e | i |
---|---|---|---|

AU | deg | ||

3J:2S, 4S:3U, 4U:3N | |||

J | 5.84724 | 0.00690581 | 0 |

S | 7.83006 | 0.0260594 | 0 |

U | 9.67303 | 0.0163112 | 0 |

N | 11.6361 | 0.0179751 | 0 |

3J:2S, 4S:3U, 3U:2N | |||

J | 5.87084 | 0.0037631 | 0 |

S | 7.99698 | 0.0165349 | 0 |

U | 9.98186 | 0.0168401 | 0 |

N | 13.1645 | 0.0064222 | 0 |

3J:2S, 3S:2U, 5U:4N | |||

J | 5.88018 | 0.00597505 | 0 |

S | 7.88684 | 0.0245716 | 0 |

U | 10.3836 | 0.0305705 | 0 |

N | 12.005 | 0.00827154 | 0 |

3J:2S, 3S:2U, 4U:3N (Levison et al. 2011) | |||

J | 5.42 | 0.0044 | 0.016 |

S | 7.32 | 0.017 | 0.016 |

Ua | 9.61 | 0.053 | 0.044 |

Na | 11.67 | 0.011 | 0.029 |

^{ a}Levison et al. (2011) do not consider Uranus and Neptune specifically, but *generic* 15 *M*_{E} ice giants instead. Indicated values correspond to reported average values for these parameters.

Download table as: ASCIITypeset image

### 2.2.Planetesimal Disk Model

Following Levison et al. (2011), who studied the effect of the planetesimal disk's self-gravity on the dynamical equilibrium of the outer solar system for the first time, we assume that initially the disk extended from a radius located outside of the outermost planet, *R*_{disk}, to an outer radius approximately at the distance of Neptune's current orbital semimajor axis, *R*_{out}. The inner disk radius is taken as a parameter in our model with several values considered as given in Table 1. Note that the Hill radius for Neptune at the corresponding distance in these compact configurations is approximately 0.31 AU, so that all values of *R*_{disk} considered are more than 3.5 Hill radii from the outermost planet; i.e., outside its so-called feeding zone. The outer radius for all cases considered is taken at 31.5 AU, in accordance with the study of Levison et al. (2011).

On account of the computation time required for these simulations, the disk is restricted to *N*_{p} = 2000 particles for most of our calculations. Although our present computational capabilities do not permit considering a much greater number of particles, in the Discussion section we address the important implications of this choice. The mass of each particle is then , where *M*_{disk} is the total disk mass. For the disk mass values adopted in this study (Table 1), the planetesimal mass ranges from ∼0.01 to 0.05 *M*_{E} or between approximately 1 and 4 lunar masses, where *M*_{E} is the mass of the Earth. Although this mass is still a factor of a few greater than the average mass of the dwarf planets, which can probably be considered the largest planetesimals, including a greater number of planetesimals with a self-consistent treatment of gravitational forces, is beyond 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 is also comparable to the mass considered for gravitational perturbers in the study of Levison et al. (2011).

The mass surface density in the planetesimal disk, Σ, is assumed to be initially axisymmetric and to have a power-law dependence on the radial distance from the Sun, , following Levison et al. (2011) and consistent with minimum mass solar nebula models. The value of Σ_{0} is adopted to reproduce the total disk mass in a specific model. The planetesimal disk is assumed to be cold, with eccentricities, *e*, and inclinations, *i* (in radians) chosen randomly from uniform number distributions in the range . The rest of the orbital elements, Ω, *ω*, and *M*, are also chosen from a uniform random distribution within the range .

### 2.3.Numerical Details

We use the `Mercury6.5` code developed by Chambers (1999) to study the dynamics of the four outer planets and the planetesimal disk. It uses a fixed time-step for the symplectic integrator coupled with a Bulirsch–Stoer (B-S) scheme, which is used to follow close encounters among bodies: we set the error tolerance of the B-S integrator to 10^{−12}. The B-S scheme is used whenever the distance between two bodies is less than three Hill radii of the largest body.

In general our simulations were carried out for up to 70 Myr. The relatively short timescale we have considered is related to the main purpose of our study, which is to assess the stability 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 solar system migration simulations (e.g., Hahn & Malhotra 1999). The simulations presented in this paper were carried out on a Linux PC, with OpenSuse 11.4 as the operating system and a `gfortran`-2.4 compiler.

## 3.RESULTS

We have carried out a series of simulations aimed at understanding the conditions under which planetesimal-driven migration, starting from multiresonant planetary configurations, becomes unstable. To this end, several initial conditions for the planets and for the planetesimal disk, were considered (see Table 1). In the following we present our results by considering the effect of varying the disk mass and inner radius for each of the multiresonant initial conditions studied.

For comparison, we have computed cases A, F, and K (described in Table 1) with the same disk properties considered by Levison et al. (2011) for their fiducial disk model: *M*_{E} and *R*_{disk} = 18 AU. However, the planetary configurations we consider (see Tables 1 and 2) differ slightly from the 3J:2S, 3S:2U, 4U:3N multiresonant system that these authors study. Levison et al. (2011) take their initial planetary configuration from Morbidelli et al. (2007) and report an average semimajor axis for Jupiter that is approximately 10% smaller than in our cases. As a consequence, the initial semimajor axis of the outer planet in the configuration considered by Levison et al. (2011), is smaller in all but the 3J:2S, 4S:3U, 4U:3N configuration we consider. This configuration is the most compact we studied and also the one that destabilizes at an earlier time for all cases.

The main effect of planetesimal–planetesimal interactions on the evolution of the system can be seen by comparing Figures 2 and 3 which show the evolution of systems without and with self-gravity in the planetesimal disk, respectively. Both figures show several snapshots of the and phase space domain of a system corresponding to the configuration of case A. Without the effect of self-gravity, Figure 2, the system is not unstable on the timescale considered, 10 Myr. The behavior of such a system is consistent with the results of typical planetesimal-driven migration simulations starting from Nice-model-like initial conditions. If we include the gravitational interaction among planetesimals, the system undergoes a dynamical instability in approximately 5 Myr. Comparing the first and second rows of panels in Figure 3, which illustrate the system's dynamical state initially and after 0.1 Myr, respectively, we see that one of the main effects of the interplanetesimal gravitational interactions is the heating of the disk, i.e., the increase in orbital eccentricities and inclinations also referred to as gravitational stirring. We discuss this process in more detail and its dependence 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 most cases considered) as the number of planetesimals dispersed into orbits capable of having close encounters with the outermost planets increases. The cumulative effect of these interactions, mainly exchanging angular momentum with Neptune and Uranus, leads to a growth of the orbital eccentricity for these planets until the multiresonant state among the giant planets is broken. The instability manifests as a very fast orbital rearrangement of all planets. The outward 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. In several cases studied, the planetary system is strongly unstable, with one or both ice giants ejected from the system, or even catastrophic such as with one of them colliding with the giant planets as indicated by the numerical simulations.

### 3.1.Effect of the Disk Mass

In order to assess the effect of the planetesimal disk mass on the evolution of each of our initial multiresonant configurations, we have carried out simulations with and . In all cases, we consider the same inner and outer disk radii of the fiducial model of Levison et al. (2011). Disk mass is varied by changing the mass of each of the particles representing the planetesimals, keeping the number of planetesimals constant.

The time evolution of the planetary orbital parameters, *a* and *e*, for all three multiresonant configurations we considered are shown in Figures 4 and 5, respectively. In these figures each row corresponds to a different planetary multiresonant initial condition. We find that, as expected, the instability arises earlier for more massive disks in which gravitational stirring increases the amount (and cumulative mass) of planetesimals interacting with the planets.

With the exception of case G, all systems considered are unstable on a timescale shorter than 42 Myr, as a result of the interaction with the self-gravitating planetesimal disk. Case G also becomes unstable but on a slightly greater timescale, less than 70 Myr. It is worth noting that case G corresponds to the combination of: (a) the less massive disk we have considered and (b) the planetary configuration in which Uranus and Neptune are most separated.

### 3.2.Effect of the Disk Inner Radius

We have also studied the effect of the location of the inner radius of the planetesimal disk, *R*_{disk}, on the stability of these systems. To do so, we fix the disk mass at 50 , following Levison et al. (2011), and vary between 16 and 20 AU. The results for all initial planetary configurations as a function of *R*_{disk} are displayed in Figures 6 and 7, which show the evolution of the semimajor axis, *a*, and eccentricity, *e*, respectively, for each of the planets 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–20 AU (a 25% increment), the time for onset of the instability, *t*_{inst}, increases by a factor of about 2–6. This result is reminiscent of the results of Gomes et al. (2005) who find that the initial location of the inner disk radius must be fine tuned to a specific value in order to reproduce the LHB timescale if this phenomenon is to be related to a planetary instability. Once again, as in our simulations with different disk mass, all cases considered with a 50 planetesimal disk, the outer solar system becomes unstable within a timescale of 30 Myr. For comparison, we have carried out simulations of the three cases with different inner radii for the 3J:2S, 4S:3U, 4U:3N configuration and the same initial distribution of planetesimals, but without the effect of planetesimal–planetesimal interactions. In all cases, the system remains stable for a timescale an order of magnitude greater or more than that with self-gravitating planetesimal disks.

### 3.3.Insight into the Instability

An additional feature of the evolution of the unstable systems is the fact that the planets remain locked in resonance until the instability arises. This is evident in Figure 8 which shows the mean motion resonant angles which, for the 3J:2S, 4S:3U, and 4U:3N initial planetary configuration, can be written as:

where is the mean longitude for the *i*th planet, with *M*_{i} being the mean anomaly, the argument of the planet's periapsis, and Ω_{i} the longitude of the ascending node. Also, is the longitude of the orbit's periapsis.

The behavior seen in Figure 8 is qualitatively similar for all systems before the instability develops, namely: (a) the and resonant arguments for Jupiter and Saturn librate about 0° and 180°, respectively, (b) either one of the or resonant arguments for Saturn and Uranus clearly librates while the other argument may circulate, and (c) at least one of the resonant arguments or for Uranus and Neptune librates.

The development of the dynamical instability is similar for all systems studied and involves a gradual increase in the number of planetesimals moving inside the disk inner boundary, *R*_{disk}, and interacting strongly with the planets. This leads to a gradual increase of the orbital eccentricities of Uranus and Neptune to values such that the MMR of these planets is broken and they can have close encounters. This is illustrated in Figure 9, which shows 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 three of the cases we simulated corresponding to disks with and AU. For comparison, we also show the behavior of Neptune's perihelia, Uranus's aphelia, and the eccentricity of their orbits.

As the system approaches instability, the cumulative effect of the interaction of the dispersed planetesimals (as a result of self-gravity) with the ice giants leads to an increase in their eccentricity until their resonant state is broken. At such a point, their semimajor axes generally increase in an abrupt manner (as shown in Figure 9), suggesting that a threshold in the interaction between planetesimals and the ice giants is triggered as the number of these gradually increases. Roughly, for the cases shown in Figure 9, the total mass that has interacted with Neptune at the time of the onset of the instability, *t*_{inst}, is between 2 and 3 Earth masses. This provides a possible explanation to the fact that disks with a mass smaller than the mass of the ice giants are not capable of destabilizing the system on the timescales considered in this study.

Figure 10 shows the evolution of the planets in phase space as the system crosses the instability. The interaction with the planetesimal disk leads to a gradual increase in the eccentricity of the ice giants (related to angular momentum) while keeping the semimajor axis (energy related) more or less constant. Also shown are the boundaries of the resonances of Saturn, Uranus, and Neptune with their interior planet and of Jupiter with the exterior planet. These are estimated by assuming that each planet is the test particle in the context of the restricted three-body problem which, in the so-called pendulum approximations, allows one to estimate the width of a given first-order resonance as described in Holmes et al. (2003). As can be seen in Figure 10, the instability occurs when either Uranus or Neptune evolves toward a dynamical state in which either is not locked in resonance, resulting in a fast disruption of the planetary configuration.

## 4.DISCUSSION

In this section we revisit some of the basic assumptions made in our study and how these affect the interpretation of our results. In addition, we discuss some of the implications of our results on the theory of planetary system evolution.

### 4.1.Self-gravity

The simulations presented above are a first attempt to study, in a self-consistent manner, the role played by the gravitational interaction between planetesimals and planets in a compact and multiresonant initial configuration for the early solar system such as that proposed as part of the Nice model. Our aim is to assess the dynamical stability and evolution of such configurations in the context of a self-consistent model. In particular, we study the conditions for the development of the instability that knocks the planets out of the multiresonant state, as found by Levison et al. (2011), and gives rise to a fast orbital rearrangement of the solar system and the ensuing LHB.

To test the effect of planetesimal–planetesimal interactions, we carried out numerical simulations of some of the initial multiresonant 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* among themselves (as usually done in studies of planetesimal-driven migration). We find that some of these systems are also unstable, but typically on a timescale of hundreds of Myr. Out of the three systems tested in this manner (cases A, F, and K), only one is found to be unstable on a timescale of less than 100 Myr. As we are focusing on the effect of disk self-gravity in this study, we have decided to leave a further exploration of this issue to future contributions.

Our main dynamical result brings into question the presumed stability of multiresonant configurations and suggests that early planetesimal-driven migration and outer disk dispersal are to be expected for a wide range of initial planetary system configurations.

### 4.2.Model Limitations

The important implications of our study suggesting that the initial configuration of the solar system was much more unstable than previously thought, a fact that may affect our understanding of the evolution of the solar system, highlights the importance of self-consistently modeling the dynamics of the disk of planetesimals and planets. The development of such models would provide firmer support for the previous conclusion and validate studies using approximate methods.

The simulations presented in this study are the first effort to study, in a fully self-consistent manner, the process of planetesimal-driven migration of a dynamically unstable initial condition for the early solar system such as that proposed as part of the Nice model. As such, it is important to point out some factors related to our present computational limitations that may affect our conclusions. First and foremost is the mass of the individual particles representing the planetesimals in our study. In our simulations each disk particle has a mass ranging from (0.75–3) g that is times the mass of Eris, the largest of the dwarf planets. Considering that the present mass of the Kuiper Belt is believed to be of the mass assumed to be originally present in the outer planetesimal disk (Gladman et al. 2001) and scaling the present day number of expected dwarf planets (around 10), one can estimate that the original number of dwarf planet-like bodies was in the vicinity of a few times 10^{3}.

It is usually assumed that the size distribution of small bodies in the early outer solar system was roughly similar to that of the present day KB, this means that each of our 2000 particles represents one dwarf planet and a very large number of much smaller bodies similar to the known KBOs and cometary nuclei, which contain most of the mass in the planetesimal disk. It is worth pointing out that there are still significant uncertainties in the present day size distribution of KBOs, particularly at the lower mass end which is still poorly constrained (Bianco et al. 2010). This translates into an even greater uncertainty in the initial size/mass distribution of planetesimals in the early solar system which is still the subject of debate. For example, contrary to the assumption described above, Johansen et al. (2014) propose that the initial mass distribution of planetesimals () scales as ; where *M* is the mass of the planetesimals. In this case, the largest bodies in the distribution make up a large fraction of the total mass in the planetesimal disk.

The issue of how the mass in the planetesimal disk is distributed among the different size bins is an important one in relation to this study as it pertains to the degree to which our present simulations are strictly applicable to the early solar system. It is known that the number of particles, *N*, in an *N*-body simulation affects the dynamical evolution of self-gravitating systems (e.g., Binney & Tremaine 2008). It is likely that a much greater number of particles in our simulations, assuming that the number of small bodies is as predicted by the present day size distribution of KBOs, can have an important effect on the dynamics of the early solar system. The presence of a very large number of small bodies in our simulations would allow for a further distribution of the interaction energy among these by increasing the number of degrees of freedom of the whole system. In this way, the small bodies can behave as a stabilizing sink for the stirring produced by the self-gravity interactions, which leads to an increase in the eccentricity of the planetesimals and the eventual destabilization of the system as these planetesimals interact with the planets. This effect clearly depends on the number of bodies representing the planetesimals as does the timescale for the development of the dynamical instability although the precise manner in which this is affected is not easy to predict.

The effect of the reduced number of particles used to describe the planetesimals in our simulations is sometimes also called numerical disk heating in the context of galactic dynamic *N*-body simulations. In this process the velocity dispersion of the disk particles, , modeled as a diffusion process, grows with time *t* as , where is a diffusion coefficient, *p* a parameter, and a constant (e.g., Lacey 1991; Velázquez 2005). The values of *D*_{p} and *p* are strongly model dependent and they have to be evaluated numerically under specific physical conditions. For example, for the case of galactic disks inside dark matter halos, one obtains with (Velázquez 2005), clearly showing a dependency on the number of particles in the system.

Assuming that the dominant physical effects in the galactic scenarios considered by Velázquez (2005) are similar to the ones in our planetesimal disk, we use a scaling relation for the diffusion coefficient to estimate the timescale for the development of the dynamical instability. Considering that the coefficient *D*_{p} defines the timescale for the growth of *σ*_{e} as a function of the number of particles in the simulation, we proceed in the following manner. The instability essentially arises when a sufficient number (mass) of planetesimals are dispersed onto Neptune crossing orbits (with ) thus leading to their strong interaction with the planet. This happens approximately when the eccentricity dispersion exceeds a critical value, , such that the Neptune crossing fraction of the population of planetesimals is capable of affecting the planetary orbital motion sufficiently enough to break its resonant condition. Assuming that is much greater than the initial eccentricity dispersion, the characteristic timescale for the growth of the eccentricity dispersion is approximately 1/*D*_{p}. Hence, under this assumption, the timescale for the development of the instability scales with the number of particles as .

One must be careful in translating the significance of these results to our specific solar system problem as the results described in Velázquez (2005) involve a series of approximations in the calculation of the total gravitational potential resulting from different components in the galactic problem (e.g., dark matter halo, bulge) which are not present in our case. In addition, it is well known from the study of Nice-model-like initial conditions that the dynamics of unstable systems are strongly chaotic so that the triggering of instability, once some threshold of perturbation for the planetary orbit has been crossed, probably depends sensitively on the details of the interaction of a few particles with the planet. Nevertheless, as we will see below, the results of a small number of simulations we have conducted, varying *N* over almost an order of magnitude, suggest that the scaling discussed above may not be entirely inadequate for our planetesimal disk.

For the specific scenario we are studying, we have analyzed the dependence of the timescale for the development of the dynamical instability, *t*_{inst}, on the number of particles in a simulation. This time approximately corresponds to the time when the eccentricity of planetesimals is excited to the point that they start approaching Neptune (). Figure 11 shows *t*_{inst} as a function of *N* for a series of simulations with the same initial conditions as case C (see Table 1), but with the number of particles ranging from *N* = 500 to *N* = 4000. While a gradual increase in *t*_{inst} as *N* increases is suggested by these results, particularly for , there are strong variations from this tendency as *N* increases up to 4000 that preclude any firm conclusion from being reached. Nevertheless, it is interesting that a simple power law scaling relation , indicated by the dotted line in Figure 11, apparently reproduces the general tendency of this limited number of simulations. This result clearly indicates the need for additional simulations with a greater number of particles. However, these kind of simulations will take an excessive amount of time with our present computational capabilities and must be postponed for future work.

The previous discussion highlights the importance of the effect of the number of particles when modeling the evolution of the planetesimal disk and the outer planets in a self-consistent manner. It also points out the complexity of this dynamical system, particularly in relation to the interplay between the compact planetary configuration (known to be strongly sensitive to minute perturbations) and the interaction among disk particles. Our results indicate, on one hand, the difficulty of assessing simply the likely outcome of the system's evolution resulting from an increase in the number of particles such as the value of *t*_{inst}. On the other hand, there is a need for a thorough and continuing study of this problem, increasing the number of particles as advances in computational technology permit. All this is particularly important to determine the validity of alternative approaches to this problem using approximate treatments such as in the commonly used models where interactions between planetesimal particles are neglected.

### 4.3.Implications for the Early Evolution of the Solar System

In the context of the simulations we have conducted, our results indicate that when gravitational interactions among the planetesimals are considered, several—perhaps most—of the multiresonant planetary configurations can become unstable on a timescale of less than a few tens of Myr. This result may pose a serious difficulty for solar system formation scenarios in which the LHB is explained as resulting from the unstable, fast migration of the outer planets and the consequent perturbation of planetesimal disk as proposed, for example, in the Nice model.

According to our results, the late development of the instability of multiresonant planets with a self-gravitating planetesimal disk on a timescale consistent with the LHB can result from the following.

(a) The disk inner radius being located at a distance much greater than 20 AU when the gaseous component of the solar nebula was dispersed, which marks the beginning of our simulated dynamical evolution. As can be seen from Figure 12, which shows the time for the development of the instability as a function of the inner disk radius of the planetesimal disk, a value of *R*_{disk} much greater than 20 AU is required to have an instability development time of Myr, which is the timescale of the LHB.

This initial location of the inner boundary of the planetesimal disk at such a great distance from the outermost planet poses, in our view, a difficult conundrum. While it is generally believed that as planets form they 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 configuration suggested above, the inner radius of the planetesimal disk is at a distance of more than 30 Hill radii from the outermost planet. It is not evident why planetesimals were cleared or were not formed in the first place only between approximately 12 and 20 AU.

(b) Another possible explanation for obtaining a long *t*_{inst} is that the initial mass of the planetesimal disk was much smaller than is usually considered. Nice model studies have regularly used *M*_{disk} around 35 , and the latest incarnation of the model favors a slightly heavier disk of 50 (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 that required for the LHB.

Using a power-law fit to the resulting timescale for the development of instability as a function of disk mass, determined from our simulations as shown in Figure 13, suggests that is required if the LHB is to be explained in terms of orbital dynamical instability. However, one would expect that very low mass disks may be altogether unable to destabilize the planetary system as the mass and angular momentum contained in the planetesimal disk is smaller than that of any of the outer planets. We can estimate a lower limit for the disk mass required to significantly modify Neptune's orbit by comparing the initial angular momentum of the planet, , with the total angular momentum in the disk, . The requirement that leads to the condition for the planetesimal disk to destabilize the outermost planet.

(c) Finally, as discussed in the previous subsection, the timescale for the development of the dynamical instability depends on the number of particles used in the simulation to represent the planetesimal disk. Using the scaling of *t*_{inst} with *N* as discussed previously, we can write the relation:

To estimate the effect of having a different number of particles in our simulations. If the mass of the planetesimal disk is concentrated in the largest bodies, as suggested by the simulations of Johansen et al. (2014), then the dynamics of the planetesimal disk can probably be well represented by 10^{4}–10^{5} particles. In such a case, according to this relation, the predicted timescale for the development of instability would be less than approximately 10 Myr. By contrast, if the mass in the disk is strongly dominated by the smallest planetesimals in the distribution in such a manner that the number of bodies in small size bins is of the order of or greater, then we expect the disk to become unstable on a timescale of approximately 700 Myr.

A convenient combination of low mass and large inner radius for the planetesimal disk, together with the appropriate size distribution for the planetesimals, can in principle be fine tuned to yield any value desired for *t*_{inst} including one consistent with the LHB. It is important to bear in mind that, as we have argued, the scaling relation we are using corresponds to disk models with single size particle distributions. Estimations done on this basis strictly correspond to systems in which no large planetesimals or dwarf planets have been formed and the disk is sufficiently far away from the orbits of giant planets, which can perturb the dynamics of such bodies and possibly lead to a different scaling relation for the disk stirring as a function of the number of particles.

An important assumption in our model is that the solar system started from a multiresonant, compact configuration as generally proposed in the latest incarnations of the Nice model. Multiresonant configurations are believed to be the natural end product of planet rearrangement due to gas-driven migration in the primeval gaseous solar nebula. The presumed stability of multiresonant configurations is further invoked to explain the long delay for the planetary configuration to undergo a fast dynamical instability and the consequent stirring of the outer planetesimal disk and the LHB (Levison et al. 2011). However, one could envision that the timescale necessary for the formation of all the giant planets and their locking into a multiresonant configuration did not exactly coincide with the timescale for the dispersal of the gaseous component of protoplanetary disks. This is particularly true if gas dispersal is related to photoevaporation of the disk or any other process unrelated to the formation of the planets themselves. In such a case, the early planetary system may be left free of gas so that no further gas-driven migration is possible before the planets reach the proposed compact multiresonant configuration.

The subsequent evolution of these initially compact systems caught outside of full multiresonance 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 similar to 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 far from Jupiter and Saturn, the evolution can be expected to proceed in a relatively smooth manner as studied by Hahn & Malhotra (1999) or Gomes et al. (2004). It must be stressed, however, that in these studies the effect of planetesimal–planetesimal interactions was not considered.

One further possibility consistent with the present study is that the unstable rearrangement of the solar system is not linked to the LHB. In this scenario the solar system underwent an unstable rearrangement shortly after gas was dispersed on a timescale of tens of 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 originally present between Mars and Jupiter (e.g., Brasser & Morbidelli 2011), it is worth noting that such a scenario cannot be entirely ruled out as a particular solution on dynamical grounds, albeit its low probability. This situation is reminiscent of the jumping Jupiter scenario invoked in the current version of the Nice model.

## 5.CONCLUSIONS

We have carried out a series of direct numerical simulations of the the dynamical evolution of the early outer solar system in the context of the Nice-2 model. We study this process, for the first time, in a fully self-consistent manner, i.e., by taking into account the gravitational interaction between all bodies included in the simulations.

Our starting point was a set of initially compact, multiresonant planetary configurations similar to those proposed in the study of Levison et al. (2011). In contrast to previous attempts to model such systems in the framework of an approximate treatment of the *N*-body problem (e.g., Levison et al. 2011), we find that a dynamical instability of the outer planets develops on a short timescale of less than 70 Myr in all the initial configurations we explored. The early development of the instability is directly related to the effect of planetesimal–planetesimal interactions as the resulting dynamical heating of the planetesimal disk accelerates the interaction with the outer ice giants, which leads to the planetary dynamical instability. It should be stressed, however, that this result depends on the number of particles in a simulation, albeit not in an easily predictable manner. In this sense the number of particles in our simulations (at the limit of our computational capabilities) means that the direct application of our results to the early solar system is uncertain. This problem is particularly important if most of the mass in the early outer planetesimal disk is contained in small kilometer-sized bodies, but may not be so critical if the initial size distribution of planetesimals is as predicted by gravitational instability formation scenarios for such bodies (Johansen et al. 2014).

The above result highlights the need for further, more consistent studies of this problem as it has important consequences for our understanding of the dynamics of the early solar system. For example, our results suggest that when the interplanetesimal gravitational interactions are included in a fully self-consistent manner, the dynamical instability of the early solar system invoked in the context of the Nice model to explain the LHB occurs much sooner than the several hundred Myr timescale indicated by the lunar crater record (Tera et al. 1974). Further simulations, particularly those increasing the number of particles representing the planetesimals, are criticial to validate our self-consistent results and those of approximate treatments of the problem of planetesimal-driven migration on which most models to date have been based.

The authors are grateful to Dr. Konstantin Batygin for kindly providing the multiresonant initial conditions used in this study, to Prof. Renu Malhotra for valuable comments on a draft of this paper, and to an anonymous referee for comments that led to an improvement of the paper. M. R. R. acknowledges financial support from DGAPA-UNAM project IN115429, H. A. from CONACyT project 179662, and all the authors from DGAPA-UNAM project IN108914.