INSTABILITY-DRIVEN DYNAMICAL EVOLUTION MODEL OF A PRIMORDIALLY FIVE-PLANET OUTER SOLAR SYSTEM

In the last two decades, discoveries of the Kuiper Belt (Jewitt & Luu 1993), as well as planets orbiting stars other than the Sun (Mayor & Queloz 1995), have supplied the centuries-old quest to understand the formation of the solar system with fresh constraints and insights into the physical processes at play. Among a multitude of newly proposed formation scenarios, the "Nice" model (Tsiganis et al. 2005; Gomes et al. 2005; Morbidelli et al. 2005) is particularly notable, as it has attained a considerable amount of success in reproducing the various observed features of the solar system. Within the context of the scenario envisioned by the Nice model, giant planets start their post-nebular evolution in a compact, multi-resonant configuration, and following a brief period of dynamical instability, scatter onto their current orbits (Morbidelli et al. 2007; Batygin & Brown 2010; Walsh et al. 2011).

The first success of the Nice model lies in its ability to quantitatively reproduce the observed orbits of the giant planets, as well as their dynamical architecture (i.e., secular eigenmodes of the system; Tsiganis et al. 2005; Morbidelli et al. 2009). Simultaneously, the brief instability, inherent to the model, provides a natural trigger to the late heavy bombardment (Gomes et al. 2005), as well as a transport mechanism for emplacement of dynamically "hot" Kuiper Belt objects from inside ∼35 AU (Levison et al. 2008). Meanwhile, it has been recently demonstrated that survival of a dynamically "cold" primordial population between Neptune's current 3:2 and 2:1 exterior mean-motion resonances (MMRs) is fully consistent with a Nice model-like evolution of the planets, implying an in situ formation of the cold classical population of the Kuiper Belt (Batygin et al. 2011). Finally, the presence of Jupiter's and Neptune's Trojan asteroids has been attributed to chaotic capture of planetesimals during the instability (Morbidelli et al. 2005; Nesvorný & Vokrouhlický 2009).

All successful realizations of the Nice model to date have been comprised exclusively of the four currently present giant planets. However, there exists no strong evidence that suggests that additional planets were not present in the solar system at the epoch of the dispersion of the nebula. In fact, theoretical arguments, presented by Goldreich et al. (2004), point to a possibility of initially forming as many as five ice giants, three of which get subsequently removed via ejections (see however Levison & Morbidelli 2007). The dynamical sensibility of such a scenario is further strengthened by the fact that a considerable fraction of standard Nice model simulations result in an ejection of an ice giant after an encounter with at least one of the gas giants.

In this Letter, we explore an instability-driven dynamical evolution of a five-planet system (two gas giants + three ice giants) with an eye toward identifying a pathway toward reproduction of a solar-system-like dynamical architecture. In principle, the realm of possibility available to this study is enormous. Consequently, rather than performing a comprehensive parameter search, here we limit ourselves to systems that contain an additional Uranus-like planet, with the aim of presenting a few proof-of-concept numerical experiments. The plan of the Letter is as follows. In Section 2, we describe our numerical setup. In Section 3, we show that in our model, the planetary orbits, their secular eigenmodes, as well as various populations of the Kuiper Belt are approximately reproduced. We conclude and discuss our results in Section 4.

The numerical setup of the simulations performed here was qualitatively similar to those presented by Batygin & Brown (2010) and Batygin et al. (2011). Particularly, the five giant planets were initialized in a compact, multi-resonant initial condition, surrounded by a massive planetesimal disk that extended between its immediate stability boundary and 30 AU.

Two of the three ice giants were taken to have the same mass as Uranus and were initially placed next to Saturn and as the outermost planet, respectively. The middle ice giant was taken to have Neptune's mass. In all simulations, Jupiter and Saturn started out in a 3:2 MMR, in accord with the results of hydrodynamical simulations of convergent migration of the planets in the solar nebula (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Pierens & Nelson 2008). The ice giants were also sequentially assembled into first-order MMRs by applying dissipative forces, designed to mimic the presence of the nebula (Lee & Peale 2002). Following resonant locking, each assembled multi-resonant initial condition was evolved in isolation for 10 Myr, as an immediate test of orbital stability. This procedure yielded a total of 81 stable multi-resonant initial conditions.

The search for adequate dynamical evolutions was performed in two steps. First, we evolved 10 permutations of each initial condition, with planetesimal disks composed of N = 1000 planetesimals. Disk masses were chosen randomly between M^min_disk = 25 M_⊕ and M^max_disk = 100 M_⊕. The density profiles followed a power-law distribution, Σ∝r^k, where the power-law index, k, was chosen randomly between k_min = 1 and k_max = 2. The planetesimals were initialized on near-circular orbits (e ∼ sin i ∼ 10⁻³). To reduce the already substantial computational cost, self-gravity of the planetesimal swarm was neglected.

Subsequently, we eliminated all initial conditions that did not yield any final systems that were comprised of four planets, reducing the number of viable initial conditions to 25. Then, an additional 30 permutations of these initial conditions were integrated with disks composed of N = 3000 planetesimals (but otherwise identical to those described above). Each integration was performed using the mercury6 integration software package (Chambers 1999) and spanned 50 Myr.³ The calculations were performed on Caltech's PANGU super-computer.

After their completion, simulations that were deemed successful were reintegrated with the use of tracer simulations (see Levison et al. 2008; Batygin et al. 2011), to address the dynamical evolution of a locally formed population of KBOs. In particular, each run was supplemented with an additional disk of massless particles that resided in the cold classical region of the Kuiper Belt (i.e., between the final exterior 3:2 and 2:1 MMRs of Neptune).

Out of the 810 integrations that were initially performed, 214 (∼25%) cases featured an ejection of a single ice giant, yielding a system composed of four planets. Perhaps unsurprisingly, in most cases the ejected planet is the ice giant that neighbors Saturn. Of the 750 simulations that were performed following the elimination of initial conditions, 33 evolutions resulted in orbits reminiscent of the solar system. Specifically, we searched for solutions where the Saturn–Jupiter period ratio exceed 2, while the final semi-major axes of the ice giants were within 3 AU of their observed counterparts. No strong requirements were placed on the planetary eccentricities and inclinations.

It is noteworthy that evolving the giant planets onto solar-system-like orbits is insufficient for a simulation to be deemed successful for indeed, there are additional constraints that must be satisfied. The first orbital constraint is the reproduction of the secular architecture of the system. The secular orbital angular momentum exchange (i.e., eccentricity evolution) of a planetary system containing N secondaries can be approximately represented as a superposition of N eigenmodes, each corresponding to a fundamental frequency of the system (see Murray & Dermott 1999; Morbidelli et al. 2009). Physically, the maximum eccentricity that a given planet attains in its secular cycle is equal to the sum of all of its corresponding eigenmode amplitudes. In the context of the traditional Nice model, some difficulty has been noted in correctly reproducing the dynamical character of Jupiter's and Saturn's eccentricity evolution. Particularly, it has been shown that smooth passage of Jupiter and Saturn through the 2:1 MMR has a tendency to underexcite the g₅ eccentricity eigenmode as well as their mutual inclinations (Morbidelli et al. 2009). This difficulty can be overcome sometimes by invoking a close encounter between an ice giant and the gas giants. Consequently, we have checked, using Fourier analysis (see, for example, Morbidelli et al. 2009), the relative strength of the g₅ and g₆ eigenmodes in all simulations whose orbital end state resembled the solar system. We did not restrict the success criteria of our simulations to include the correct reproduction of the mean eccentricities of the planets and in some cases, the mean final eccentricities (namely those of the ice giants) exceeded their observed counterparts by as much as a factor of ∼2. This is, however, likely an artifact of the coarse representation of the planetesimal disk and the resulting dynamical friction, that we employed in our calculations and should not be viewed as a major drawback.

In total, we found 10 cases (corresponding to eight different initial conditions) where their amplitudes are satisfactorily reproduced. Specifically, in these 10 cases, the amplitude of Jupiter's g₅ eigenmode exceeds that of the g₆ eigenmode, while the amplitudes are roughly the same for Saturn. We made no attempt at quantitatively matching the pair's inclination eigenmodes; however, their reproduction does not appear to be problematic (Morbidelli et al. 2009). We also examined the amplitudes of the secular eigenmodes of Uranus and Neptune. Generally it appears that the dynamical architecture of the ice giants is set in an essentially random manner, depending on the particular encounter history. Consequently, we decided to not use ice-giant secular architecture as a distinctive property in our analysis.

Successful formation of the Kuiper Belt is another important constraint of the Nice model. Levison et al. (2008) have shown that the excited populations of the Kuiper Belt are naturally emplaced from inside ∼35 AU during the instability. Given that this aspect of the planetary evolution is not particularly different between the four-planet and the five-planet scenarios, there is little reason to speculate that the dynamical pathway for formation of the resonant, scattered, and hot classical populations of the Kuiper Belt will be inhibited. The same is likely to be true for the chaotic capture of Trojan asteroids. The cold classical population of the Kuiper Belt, however, is a different story.

A series of observational dissimilarities between the cold classical population and the rest of the Kuiper Belt (e.g., uniquely red colors, Trujillo & Brown 2002; Lykawka & Mukai 2005; strongly enhanced wide binary fraction, Stephens & Noll 2006; Parker & Kavelaars 2010) suggest that the cold classicals formed in situ, and maintained dynamical coherence despite Neptune's temporary acquisition of high eccentricity and inclination (a characteristic orbital feature of the cold classical population is inclination that does not exceed ∼5° (Brown 2001)). In a recent study, Batygin et al. (2011) showed that local formation of the cold classicals is fully consistent with an instability-driven evolution of the planets, given favorable conditions during the instability. Particularly, Batygin et al. (2011) required the apsidal precession and nodal recession rates of Neptune to be comparatively fast to prevent secular excitation of the cold classical orbits, in addition to a sufficiently small apohelion distance of Neptune, to avoid orbital excitation due to close encounters.

In a five-planet scenario, the retention of unexcited orbits of the cold population can be jeopardized by the ejecting planet. This possibility served as a premise for recalculation of the dynamical evolution of test particles in the cold classical region, with the aid of tracer simulations. We found that in only three of the ten simulations where the secular eigenmodes were successfully reproduced, a primordial cold classical population of the Kuiper Belt also retained unexcited orbits. Consequently, it appears that the retention of an unexcited cold belt is not the norm of the five-planet scenario. This is not surprising, since dynamical excitation can only be avoided if ejection and close encounter timescales are sufficiently short. From a quantitative point of view, our simulations imply that dynamically cold orbits can only be sustained if the ejecting ice giant spends ∼10⁴ years or less crossing the classical Kuiper Belt region. As a result, we can expect that the incorporation of yet additional ice giants into the primordial solar system (as suggested by Goldreich et al. 2004) will further diminish the chances of reproducing the cold classical Kuiper Belt, since all ejecting planets would have to do so very rapidly.

Two of the three successful simulations discussed above are presented in Figures 1–3. Specifically, Figure 1 shows the orbital evolutions of the runs, while Figures 2 and 3 show the eccentricity and inclination distributions of planetesimals in the Kuiper Belt region, respectively. The red and blue dots, shown in Figures 2 and 3, depict the orbital distribution of emplaced and local planetesimal populations, respectively, at t = 50 Myr. The red and blue triangles depict the planetesimals whose orbits remain stable at t = 500 Myr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The starting multi-resonant initial condition of the simulation presented in panels (A) of the figures is one where Saturn and the first ice giant are locked in a 2:1 MMR, while both pairs of ice giants are locked in 5:4 MMRs. The planetesimal disk in this simulation was comprised of N = 3000 particles and contained a total of 26 M_⊕. In this evolution, the local population of test particles suffers numerous short close encounters with the ejecting ice giant, yielding a more excited and depleted cold classical population, compared to that of run (B). Note also that at t = 50 Myr, the ice-giant eccentricities are considerably greater than that of Uranus and Neptune. These high eccentricities do not get damped away by dynamical friction in the following 500 Myr of dormant evolution. Consequently, in this simulation, the inner edge of the cold belt gets dynamically depleted over the following 500 Myr. It is furthermore noteworthy that another simulation that originated from the same initial condition reproduced the eigenmodes of the system correctly, although the primordial cold Kuiper Belt in this integration was entirely destroyed by close encounters.

The starting multi-resonant initial condition of the simulation presented in panels (B) of the figures is one where Saturn and the first ice giant are locked in a 3:2 MMR, the inner pair of ice giants is locked in a 2:1 MMR, and the outer pair of the ice giants is locked in a 4:3 MMR. In this simulation, the disk consisted of N = 1000 particles and had a cumulative mass of 42 M_⊕. Incidentally, the frequency spectrum of the eccentricity vectors of Jupiter and Saturn produced in this simulation, matches that of the real Jupiter and Saturn exceptionally well, signaling a nearly ideal reproduction of the secular eigenmodes. Particularly, the simulation yields (in the notation of Murray & Dermott 1999) e^(sim)₅₅/e₅₆^(sim) = 2.28 and e^(sim)₆₅/e₆₆^(sim) = 0.51, whereas the solar system is characterized by e₅₅/e₅₆ = 2.81 and e₆₅/e₆₆ = 0.68. The scaling of the eigenvectors are also well reproduced: e^sim₅₅ = 0.0465, e₆₆^sim = 0.067, while the solar system has e₅₅ = 0.0442, e₆₆ = 0.0482 (Morbidelli et al. 2009). In this simulation, the cold Kuiper Belt suffers only a single short encounter with the escaping ice giant, allowing for the orbits (inclinations in particular) to remain dynamically cold.

In this Letter, we have presented a successful realization of the Nice model that starts out with five planets. The numerical experiments presented here explicitly show that such an evolution is plausible since the resulting four-planet systems can closely resemble the solar system. Particularly, in both simulations presented here, the secular architecture of the outer solar system is well reproduced. Furthermore, the demonstrated survival of a local, primordial cold classical Kuiper Belt suggests that all constraints that can be matched with a four-planet model can also be matched with a five-planet model to an equal degree of satisfaction.

It is noteworthy that ejection was not always necessary in our simulations to generate a four-planet system. In a handful of runs (one of which successfully reproduced the secular eigenmodes, but not the cold Kuiper Belt), one of the ice giants ended up merging with Saturn. In principle, such a scenario may help explain Saturn's enhanced metallicity in comparison with Jupiter, although here again the explanation is not unique (see Stevenson 1982 and references therein).

In a traditional realization of the Nice model, the rate of successful reproduction of the secular eigenmodes is rather low, i.e., ∼10% of the integrations for a favorable initial condition (Batygin & Brown 2010). This is in part because an ice-giant/gas-giant encounter often leads to an ejection of the ice giant, leaving behind only three planets. Thus, the need for an ice-giant/gas-giant encounter in the orbital history of the solar system is in itself motivation for a five-planet model.

The statistics of simulations presented in this work suggest that a five-planet model is neither more nor less advantageous. Recall that the probability of ending up with only four planets is 214/810 ∼ 25%. The probability of reproducing the secular eigenmodes of Jupiter and Saturn is 10/750 ∼ 1.5%. Naively, this yields an overall probability of success of only ∼0.4%. However, it is important to keep in mind that the characteristic outcomes are generally dependent on initial conditions⁴ and runs that originated from the initial condition presented in panels (A) correctly reproduced the secular eigenmodes in 2 out of 30 simulations. This statistic is similar to the four-planet model. Finally, the 1/30 probability also of retaining an unexcited cold classical Kuiper Belt puts the five-planet model and the four-planet model on equal footing in terms of success rate (Batygin et al. 2011). That said, it is important to note that this success rate is only characteristic of the particular five-planet model that we have constructed. In other words, it is likely that if one allows the mass of the ejected planet to also be a variable parameter, tuning of the initial state may in principle lead to a more frequent reproduction of the solar system.

The results presented in this work imply that the solar system is one of many possible outcomes of dynamical evolution and can originate from many possible initial conditions. As a result, the possibility of having an extra planet initially present in the system, yet its ejection leaving no observable signature, erases any hope for construction of a deterministic model for solar system evolution. The forward-process-like nature of the Nice model is not surprising, given that the solar system exhibits large-scale chaos, characterized by Lypunov times that are comparable to orbital timescales, during the instability. Moreover, the similarity between the orbital architectures of simulations whose outcomes were deemed unsuccessful in this work and those of extra-solar planetary systems further confirms that planet–planet scattering is likely to be the physical process responsible for shaping the orbital distribution of planets (Jurić & Tremaine 2008). Consequently, we conclude that an instability-driven dynamical history remains a sensible choice as a baseline scenario for solar system's early dynamical evolution.

We thank Alessandro Morbidelli, Hal Levison, and Peter Goldreich for useful conversations. We thank Naveed Near-Ansari for operational help with the PANGU supercomputer. K.B. acknowledges support from NASA NESSF Graduate Fellowship.