Born eccentric: constraints on Jupiter and Saturn's pre-instability orbits
Mattthew S. Clement, Sean N. Raymond, Nathan A. Kaib, Rogerio Deienno, John E. Chambers, Andre Izidoro
AAccepted for publication in Icarus
Preprint typeset using L A TEX style AASTeX6 v. 1.0
BORN ECCENTRIC: CONSTRAINTS ON JUPITER AND SATURN’S PRE-INSTABILITY ORBITS
Matthew S. Clement , Sean N. Raymond , Nathan A. Kaib , Rogerio Deienno , John E. Chambers & Andr´eIzidoro Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road, NW, Washington, DC 20015, USA Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, all´e Geoffroy Saint-Hilaire, 33615 Pessac, France HL Dodge Department of Physics Astronomy, University of Oklahoma, Norman, OK 73019, USA Southwest Research Institute, 1050 Walnut St. Suite 300, Boulder, CO 80302, USA Department of Earth, Environmental and Planetary Sciences, MS 126, Rice University, Houston, TX 77005, USA
ABSTRACTAn episode of dynamical instability is thought to have sculpted the orbital structure of the outer solarsystem. When modeling this instability, a key constraint comes from Jupiter’s fifth eccentric mode(quantified by its amplitude M ), which is an important driver of the solar system’s secular evolution.Starting from commonly-assumed near-circular orbits, the present-day giant planets’ architecture liesat the limit of numerically generated systems, and M is rarely excited to its true value. Herewe perform a dynamical analysis of a large batch of artificially triggered instabilities, and test avariety of configurations for the giant planets’ primordial orbits. In addition to more standard setups,and motivated by the results of modern hydrodynamical simulations of the giant planets’ evolutionwithin the primordial gaseous disk, we consider the possibility that Jupiter and Saturn emergedfrom the nebular gas locked in 2:1 resonance with non-zero eccentricities. We show that, in sucha scenario, the modern Jupiter-Saturn system represents a typical simulation outcome, and M iscommonly matched. Furthermore, we show that Uranus and Neptune’s final orbits are determined bya combination of the mass in the primordial Kuiper belt and that of an ejected ice giant. Keywords: giant planets, solar system, planet formation INTRODUCTIONThe realization that the outer planets’ orbits divergeover time (Fernandez & Ip 1984) as a consequenceof small exchanges of angular momentum with nearbyplanetesimals has revolutionized our understanding ofthe solar system’s global dynamical history over the pastseveral decades. Early work invoked migration to ex-plain the orbital and resonant structures of the Kuiperbelt (Malhotra 1993, 1995). A consequence of this mi-gration is an epoch of dynamical instability that rapidlytransforms the outer solar system from its primordiallyconcentrated architecture into the radially diffuse sys-tem of orbits that exists today (Thommes et al. 1999).These ideas were eventually built upon to form a com-prehensive evolutionary framework for the solar system.In this paradigm of a migration-shaped solar system,the Kuiper belt’s modern structure is tied to the giantplanets’ early dynamics (Hahn & Malhotra 2005; Levi-son et al. 2008; Batygin et al. 2011; Nesvorn´y 2015a,b;Kaib & Sheppard 2016) which, in turn, are dependenton the properties of the primordial gas disk (Masset &Snellgrove 2001; Morbidelli et al. 2007; Pierens & Nelson2008; Zhang & Zhou 2010) and first solid bodies (Hahn & Malhotra 1999; Gomes et al. 2004). These ideas led tothe development of the so-called
Nice Model ; a robustlyformulated hypothesis for the solar system’s dynamicalevolution (Tsiganis et al. 2005; Morbidelli et al. 2005;Gomes et al. 2005). In the current framework of theNice Model, the giant planets emerge from the primor-dial gaseous nebula in a compact, resonant configuration(Morbidelli et al. 2007). After some period of time, thecollection of synchronized orbits is cataclysmically de-stroyed when one or more of the planets are perturbedout of resonance (Levison et al. 2011; Nesvorn´y & Mor-bidelli 2012; Deienno et al. 2017; Quarles & Kaib 2019).This departure from resonance triggers an epoch of dy-namical instability that expediently reshapes the outersolar system into its modern form (for a recent review,see Nesvorn´y 2018).While researchers continue to debate specific aspectsof the Nice Model (e.g.: its initial conditions, timingand strength; addressed in greater detail in section 2),its wide acceptance within the field is clearly a resultof the scenario’s consistent ability to reproduce manypeculiar aspects of the solar system in numerical sim-ulations. Among these are the orbital structure of the a r X i v : . [ a s t r o - ph . E P ] S e p Kuiper (Levison et al. 2008; Nesvorn´y 2015a,b) and as-teroid belts (Roig & Nesvorn´y 2015; Deienno et al. 2016,2018; Clement et al. 2019b), the capture and evolutionof trojan asteroids in the outer solar system (Morbidelliet al. 2005; Nesvorn´y et al. 2013, 2018), and certainproperties of the giant planets’ moons (Barr & Canup2010; Deienno et al. 2014; Nesvorn´y et al. 2014a,b). Inparticular, a sequence of events where the giant planetsacquire their modern orbits through a series of plane-tary encounters, rather than via smooth migration, isthe only known model capable of explaining the captureof irregular moons in the outer solar system (Nesvorn´yet al. 2007), and preventing sweeping secular resonancesfrom destroying the asteroid belt (Walsh & Morbidelli2011).The highly chaotic nature of the Nice Model insta-bility makes it difficult to investigate with N-body sim-ulations. While the exact initial configuration of theouter solar system prior to the instability is unknown,the giant planets’ initial orbits can be somewhat con-strained by hydrodynamical simulations that study theirmigration in the gas disk phase (e.g.: Pierens & Ray-mond 2011; Pierens et al. 2014), and by comparing largesuites of instability outcomes with important aspects ofthe modern planets’ architecture (e.g.: Nesvorn´y 2011).However, even the most successful sets of initial condi-tions tested in statistical studies of the instability re-produce the outer solar system in broad strokes just afew percent of the time (Nesvorn´y 2011; Nesvorn´y &Morbidelli 2012; Batygin et al. 2012; Kaib & Cham-bers 2016; Deienno et al. 2017; Clement et al. 2018).The stochastic nature of the giant planets’ early evolu-tion thus presents a significant complication for authorsattempting to break degeneracies between the differ-ent possible primordial giant planet configurations andglobal disk properties.In this paper we present a robust dynamical analysisof the solar system’s instability similar to previous worksby Nesvorn´y & Morbidelli (2012) and Deienno et al.(2017). In particular, we summarize the problems withthe common assumption that Jupiter and Saturn werecaptured in to a 3:2 mean motion resonance (MMR),and systematically test whether the 2:1 is a viable al-ternative as proposed by Pierens et al. (2014). We alsoconsider the possibility that the giant planets alreadypossessed moderate eccentricities ( ∼ MOTIVATION2.1.
Background
The 3:2 MMR
The original Nice Model simulations presented in Tsi-ganis et al. (2005) configured the primordial outer solarsystem in a somewhat ad hoc manner. Specifically, thegas giants’ initial orbits were assigned such that Jupiterand Saturn would migrate through their mutual 2:1 res-onance and trigger the instability. However, such ini-tial conditions are at odds with results from investiga-tions of giant planet growth and evolution in the gasdisk phase (particularly foundational were the works ofMasset & Snellgrove 2001; Morbidelli et al. 2007; Mor-bidelli & Crida 2007). Hydrodynamical simulations ofthe giant planets migrating in gaseous disks find that thegiant planets are most likely to emerge from the primor-dial nebula in a mutual resonant configuration (Zhang& Zhou 2010; Pierens & Raymond 2011; D’Angelo &Marzari 2012; Izidoro et al. 2015). Moreover, theseinitial conditions are largely consistent with observedresonant chains of giant exoplanets (perhaps most fa-mously by GJ 876: Rivera et al. 2010) and gaps in proto-planetary disks presumably induced by growing planets(e.g.: Bae et al. 2019). Thus, authors investigating theNice Model must first determine which particular chainof resonances the giant planets were born in.Early hydrodynamical simulations indicated ratherconvincingly that capture within the 3:2 MMR is theonly possibility for a Jupiter-Saturn-like mass configura-tion (Morbidelli et al. 2007; Pierens & Nelson 2008). Asa result, an overwhelming number of dynamical investi-gations dedicated to the study of the Nice Model overthe past decade or so have almost exclusively consid-ered Jupiter and Saturn in a primordial 3:2 MMR (e.g.:Nesvorn´y 2011; Nesvorn´y et al. 2013; Roig & Nesvorn´y2015; Roig et al. 2016; Kaib & Chambers 2016; Dei-enno et al. 2016, 2018; Clement et al. 2018, 2019b). Ingeneral, instabilities originating from the 3:2 Jupiter-Saturn resonance are inherently more violent that thosethat begin with the planets in a 2:1 MMR (Nesvorn´y& Morbidelli 2012). This is an obvious consequence ofthe fact that placing the solar system’s two most mas- Note that it may not be a strict requirement that the giantplanets be in resonance prior to the instability (as was found inthe pebble accretion model of Levison et al. 2015). Indeed, sucha scenario has been shown to be more successful at generatingthe giant planets’ modern obliquities (Brasser & Lee 2015) thanthe standard resonant version of the Nice Model (Vokrouhlick´y& Nesvorn´y 2015). We speculate further about this alternativescenario in section 5. sive planets in closer proximity to one another leads tostronger gravitational encounters within the instability.As Uranus or Neptune are often ejected in dynamicalsimulations of such a violent instability (Morbidelli et al.2009a), versions of the departure-from-resonance NiceModel adaptation (in contrast to the resonance-crossingscenario originally envisioned in Tsiganis et al. 2005)tend to invoke either a fifth or sixth primordial giantplanet. In successful models, the additional ice giantsare ejected from the system during the instability, thusleading to a higher fraction of simulations finishing withfour giant planets (Nesvorn´y 2011).As each individual numerically generated instabilityfrom a given resonant chain is unique, the planets’ par-ticular evolution within the instability (in addition tothe initial conditions) is also important to constrain.The Kuiper belt’s final structure (Nesvorn´y 2015a,b;Nesvorn´y & Vokrouhlick´y 2016) is somewhat sensitiveto the instability’s timing (mainly Neptune’s migrationhistory. However Volk & Malhotra 2019, showed thatNeptune’s eccentricity is important as well). In con-trast, the asteroid belt (Roig & Nesvorn´y 2015; Izidoroet al. 2016; Deienno et al. 2016, 2018) and terrestrialplanets’ (Brasser et al. 2009; Agnor & Lin 2012; Kaib &Chambers 2016) survival is particularly tied to the vio-lent and chaotic evolution within the instability. Pow-erful secular resonances chaotically traversing regions ofthe inner solar system during the instability can lead toplanet ejections, collisions, and levels of dynamical exci-tation inconsistent with that of the modern inner solarsystem. To avoid this problem, Brasser et al. (2009)proposed that Jupiter and Saturn’s semi-major axesmust strictly evolve in a step-wise manner. While thisso-called “jumping Jupiter” type of evolution can ade-quately reproduce the terrestrial planets’ orbits (Roiget al. 2016) and the asteroid belt’s dynamical structure(Deienno et al. 2016, 2018), the inner asteroid belt’shigh-inclination population can be over-populated as aresult of the violent sweeping of the ν secular resonance(Morbidelli et al. 2010; Walsh & Morbidelli 2011; Minton& Malhotra 2011) if the jump is not ideal. However,many of these high-inclination asteroids can be removedwhen ν reverses direction as Jupiter and Saturn ap-proach the 5:2 MMR (Clement et al. 2020), therebymaking somewhat less-regular instabilities with less-ideal jumps potentially viable. Moreover, such finely-tuned instabilities might also not be required to savethe terrestrial planets from loss and over-excitation (e.g.:Roig et al. 2016) if the event occurs while they are stillforming (Lykawka & Ito 2013; Clement et al. 2018; Ray-mond et al. 2018). In such a scenario, Jupiter’s en-hanced eccentricity (Raymond et al. 2009) and chaoti-cally sweeping secular resonances (Clement et al. 2019b)excite and remove material from the asteroid belt- and Mars-forming regions, while leaving Earth and Venus’growth largely undisturbed (Clement et al. 2019a). Assimulations in Clement et al. (2018) found satisfactoryterrestrial planet outcomes in systems that experienceda variety of jumps, the range of possible viable instabil-ity evolutions might be broadened if the event occurs insitu with the process of terrestrial planet formation (wediscuss the timing of the Nice Model in appendix A).While such a violent scenario (without an ideal jump)is capable of depleting the asteroid belt’s total mass by3-4 orders of magnitude (Clement et al. 2019b), it is stillunclear whether such a powerful depletion event and asuccessful inner solar system are mutually-exclusive re-sults.In summary, for the reasons detailed above, a “jump-ing Jupiter” type instability originating from the 3:2Jupiter-Saturn resonance remains the current consensusversion of the Nice Model. However, adequately excit-ing Jupiter’s eccentricity without exceeding Jupiter andSaturn’s modern orbital spacing is extremely challenging(see further discussion about the primordial 3:2 Jupiter-Saturn resonance’s problem in section 2.2). Therefore,we argue that a detailed investigation of alternative pri-mordial outer solar system configurations, specificallythe 2:1 Jupiter-Saturn resonance, is warranted.2.1.2. The 2:1 MMR
As a result of hydrodynamical studies that largelyconsidered isothermal disks of fixed viscosity (e.g.: Mor-bidelli et al. 2007; Pierens & Nelson 2008) favoring the3:2 Jupiter-Saturn resonance, subsequent study of theprimordial 2:1 MMR’s dynamics and potential insta-bility evolution was largely regarded as a purely aca-demic endeavor (Nesvorn´y & Morbidelli 2012). In gen-eral, instabilities beginning from the 2:1 are less violentthan those born out of the 3:2 (Nesvorn´y & Morbidelli2012; Deienno et al. 2017), and often yield evolution-ary histories for the planets similar to those recorded insimulations of the original Nice Model (Tsiganis et al.2005). Thus, the 2:1 is highly successful at generat-ing a final system of planets with low eccentricities andinclinations, and less successful at exciting Jupiter’s ec-centricity (though the primordial 3:2 MMR struggles inthis manner as well) and matching the modern Jupiter-Saturn period ratio. As discussed in the previous sec-tion, exciting Jupiter’s eccentricity requires a relativelystrong encounter with the ejected ice giant (Morbidelliet al. 2009a). When the planets emerge from the gas diskin the 2:1 MMR’s broader radial spacing, the resultinginstability encounters can often be weaker. While thiscan lead to systematically smaller jumps, it is still rel-atively easy for a system to exceed the modern valueof P S /P J = 2.49 when Jupiter and Saturn migratesmoothly after the instability (Nesvorn´y & Morbidelli2012). As the amount of residual migration is a functionof the remnant planetesimal mass, and the magnitudeof Saturn’s jump is related to the mass of the ejectedplanet, it may be possible to improve the likelihood thatsimulations finish with P S /P J near the modern value byadjusting these two parameters (we explore this in sec-tions 4.6 and 4.7).Hydrodynamical simulations probing different diskthermal and viscosity profiles in Pierens et al. (2014)found that the solar system’s gas giants’ capture in the2:1 MMR is a real possibility in relatively low-massdisks. Additionally, the authors found that outward mi-gration occurs in the 2:1 MMR when the disk viscos-ity parameter is low. Of great interest to our presentmanuscript, Pierens et al. (2014) found several caseswhere Jupiter and Saturn attain relatively high eccen-tricities ( ∼ ∼
40% of systems displaying the appropriateevolution). However, Deienno et al. (2017) also foundthat wide resonant chains beginning with the gas giantsand first ice giant in a 2:1,3:2 chain are highly stable,and typically do not decompose into an orbital instabil- ity. Therefore, it is important to point out that tighterice giant resonances (e.g. 4:3 or 5:4) are a possible pre-condition of the 2:1 Jupiter-Saturn configuration.2.2.
Exciting M An instability simulation where the final giant planets’mean eccentricities closely resemble those of the actualplanets might still be a poor solar system analog if itssecular architecture is incorrect. Indeed, mean eccen-tricities are insufficient to fully characterize the secularstructure of a planetary system because planets’ eccen-tricities oscillate on timescales of order P J M (cid:12) /M J (asJupiter is the most significant perturber in the N-bodyproblem of the solar system, e.g.: Poincare 1892; Laskar1990; Morbidelli et al. 2009b). For a system possessing N massive planetary bodies there will be 2 N fundamen-tal frequencies describing the system’s secular evolution.These frequencies, denoted g i for the eccentricity vectorprecessions and s i for the nodal precession rates, andtheir magnitudes in each planet’s orbit ( M ij ; i, j = 1-8) define the secular evolution of the solar system. AsJupiter and Saturn are the most significant eccentricperturbing bodies in the solar system (see further dis-cussion in appendix B, it is reasonable to seek out anevolutionary model that consistently generates the cor-rect orbital orientation and secular structure for the twoplanets. For the reasons outlined in the introduction,the most successful hypothesis to date is the Nice Model(Tsiganis et al. 2005). Moreover, the arguments laid outin Morbidelli et al. (2009a) convincingly explain why theinstability scenario is the only viable option to explainthe solar system’s unique secular architecture: partic-ularly the large magnitude of M and the fact thatJupiter’s eccentricity lies in the | M | > | M | regime.Morbidelli et al. (2009a) also investigated the possibilityof Jupiter’s eccentricity being excited solely via plane-tary migration, MMR crossings, or three planet dynam-ics (e.g.: the migration of an eccentric Uranus or closeencounters with Uranus). The authors concluded thatno alternative could provide strong enough excitationin as self-consistent of a manner as the Nice Model (webriefly explore the possibility of alternative excitationmechanisms in appendix C).To illustrate the types of final Jupiter-Saturn config-urations that emerge from instabilities beginning withthe gas giants locked in a 3:2 MMR, we analyze twobatches of instability simulations from Clement et al.(2018) and Clement et al. (2019a). These initial con-ditions (table 1) are based on the most successful five(5 GP control ) and six giant planet (6 GP control ) 3:2 con-figurations from Nesvorn´y & Morbidelli (2012). In figure1, we plot the 298 simulations (of 1,800 total) from bothpapers that finish with P S /P J < D ,see section 2.3). While higher values of M are gener- Name N Pln M disk δ r r out a nep Resonance Chain M ice ( M ⊕ ) (au) (au) (au) ( M ⊕ )5 GP control GP control Table 1 . Table of giant planet initial resonant configurations reproduced from Clement et al. (2018) and Clement et al. (2019a).The columns are: (1) the name of the simulation set, (2) the number of giant planets, (3) the mass of the planetesimal diskexterior to the giant planets, (4) the distance between the outermost ice giant and the planetesimal disks inner edge, (5) thesemi-major axis of the outermost ice giant (commonly referred to as Neptune, however not necessarily the planet which completesthe simulation at Neptune’s present orbit), (6) the resonant configuration of the giant planets starting with the Jupiter/Saturnresonance, and (7) the masses of the ice giants from inside to outside. ated regularly in simulations where Saturn is scatteredonto a distant orbit, it is clear from figure 1 that the so-lar system result lies at the extreme of the distributionof possible outcomes in M - P S /P J space for P S /P J < M , M , M and M we find that it is extremely difficult to adequatelyexcite M compared to the other three modes. In thismanner, studies that take the standard success crite-ria of P S /P J < M > M and systemswith P S /P J > GP control set are perhaps the best explanation forthe solar system’s dynamical state (Nesvorn´y & Mor-bidelli 2012; Batygin et al. 2012; Roig & Nesvorn´y 2015;Brasser & Lee 2015; Vokrouhlick´y & Nesvorn´y 2015;Roig et al. 2016; Deienno et al. 2016, 2017). An ex-ample of such a successful, “Jumping Jupiter” (Brasseret al. 2009) type of evolution from Clement et al. (2019a)is plotted in figure 2. In this case, Jupiter’s eccentric-ity is excited by a strong encounter with the ejected icegiant during the instability. We also direct the readerto the three instabilities scrutinized by Nesvorn´y et al.(2013), the one discussed in Batygin et al. (2012), andthe one in Deienno et al. (2018) for other examples ofideal evolutions from the primordial 3:2 Jupiter-Saturnresonance. However it is important to note that thesesuccessful systems lie at the extreme of simulation gen-erated outcomes in M - P S /P J space, and thus at theextreme of the spectrum of systems satisfying the tradi-tional success criteria ( P S /P J < M > -3 -2 -1 | M | P S /P J -3 -2 -1 | M | − | M | | M | > | M | regime -3 -2 -1 M Figure 1 . 386 instability simulations from Clement et al.(2018) and Clement et al. (2019a). The top panel plots thefinal value of M against the Jupiter-Saturn period ratio.The bottom panel depicts | M | − | M | vs P S /P J < for 184simulations that finished in the | M | > | M | regime. Theshaded grey area delimits the region of 2.3 < P S /P J < < M < ing the viability of the 2:1 Jupiter-Saturn resonance andheightened eccentricities: Pierens et al. 2014).2.3. Updated success criteria
Over the past decade, most numerical studies (e.g.:Kaib & Chambers 2016; Deienno et al. 2017; Clementet al. 2018; Deienno et al. 2018; Quarles & Kaib 2019) of
Criteria Nesvorn´y & Morbidelli (2012) Deienno et al. (2017) This Work A N GP = 4 N GP = 4 N GP = 4 B | ∆¯ a/a ss | < .
2, ¯ e < .
11, ¯ i < . ◦ | ∆¯ a/a ss | < .
2, ¯ e < .
11, ¯ i < . ◦ | ∆¯ a/a ss | < . C M > M > | ∆ M ij /M ij,ss | < i, j =5,6), M > M D P S /P J : < < P S /P J : < < P S /P J < E N/A 27 < a N <
29 au @ t inst N/A
Table 2 . Summary of success criteria A - D for our present manuscript (right column), compared with previous work by Nesvorn´y& Morbidelli (2012) (left column) and Deienno et al. (2017) (center column). The Subscripts GP , SS , J , S , and N denotethe giant planets, solar system modern values, Jupiter, Saturn, and Neptune, respectively. Note that, unlike the other successcriteria, criterion B is only evaluated when A is satisfied. q , Q ( a u ) Example successful evolution from primordial 3:2
JupiterSaturn -2 -1 Time (Myr) P S / P J Figure 2 . Example instability evolution from Clement et al.(2019a). Using the nomenclature of that work, this sim-ulation was from set 5
GP/ Myr (indicating it took the5 GP control initial conditions and delayed the instability 1Myr after the start of terrestrial planet formation). Thesimulation finished with P S /P J =2.41 and M =.039. It isimportant to note that this represents a very rare instabilityoutcome of these initial conditions.
More typical final sys-tems eject too many planets and possess over-excited Saturnanalogs on orbits that are too distant. The top panel plotsthe perihelion and aphelion of each planet over the length ofthe simulation. The bottom panel shows the Jupiter-Saturnperiod ratio. The horizontal dashed lines in the upper panelindicate the locations of the giant planets’ modern semi-major axes. The shaded region in the lower panel delimitsthe range of 2.3 < P S /P J < the giant planet instability have invoked the success cri-teria established by Nesvorn´y & Morbidelli (2012) (col-umn one of table 2). Criterion A requires that a systemfinish with 4 giant planets, criterion B stipulates thatplanets’ orbits (in terms of a , e and i ) be close to the realones, and criteria C and D mandate M > P S /P J evolve from < < E for the migration of Neptune, based off thework of Nesvorn´y (2015a,b). As our present study ismost interested in finding the initial giant planet config-urations that best generate the modern Jupiter-Saturnsystem, we leave the analysis of Neptune’s migration his-tory to future work. Moreover, factors such as Neptune’seccentricity evolution, the instability’s strength, and thecold Kuiper Belt’s formation location are also importantfor shaping the trans-Neptunian region (Gomes et al.2018; Volk & Malhotra 2019) Thus, we focus on theoriginal four success criteria of Nesvorn´y & Morbidelli(2012), with several minor modifications described be-low:2.3.1. The outer solar system’s secular architecture
Figure 1 illustrates the difficulty of matching P S /P J and M in a single numerical simulation. Thus, cri-terion C satisfying simulations tend to have final M values that are between 0.022 and the modern value of0.044 (table B2). This, coupled with criterion B al-lowing final average eccentricities up to 0.11, leads to“successful” systems that often inhabit the M > M regime. In this somewhat-typical scenario (40% of our5 GP control and 6 GP control simulations that satisfy crite-ria the C and D of Nesvorn´y & Morbidelli (2012) simul-taneously), Saturn’s eccentricity is over-excited relativeto Jupiter’s; often yielding a semi-major axis jump thatis too large ( P S /P J > M > M ( M = 0.037 and M = 0.023), and thegas giant’s final semi-major axes finish beyond of the5:2 MMR ( P S /P J = 2.58). However, because Saturnand the Uranus analog each have final eccentricities of (cid:39) | M | is still less than | M | in the successful simulation plotted in figure 2,Saturn’s final average eccentricity is still much higher( ∼ C to require that M be greater than M , and that each of the four mag-nitudes in the Jupiter-Saturn secular system (equationB4) be within 50% of their modern values. As Saturn’seccentricity is relatively easy to excite, the additional re-quirements for M and M do not overly constrain oursimulations. Moreover, as Jupiter and Saturn’s eccen-tricities are now fully evaluated by criterion C , and theice giant’s eccentricities and inclinations are somewhateasily damped (e.g.: figure 3) during the post-instabilityresidual migration phase (Nesvorn´y & Morbidelli 2012),we relax criterion B to only scrutinize the giant planet’sfinal semi-major axes.2.3.2. The Jupiter-Saturn period ratio
Recent work proposed that Jupiter and Saturn’s pre-cise migration towards the 5:2 MMR sculpted the innerasteroid belt’s inclination distribution (see full discus-sion in Clement et al. 2020). Furthermore, the gas gi-ants’ ultimate approach to the modern value of P S /P J =2.49 is fossilized in the asteroid belt’s orbital precessiondistribution as an absence objects with g > g > g − (cid:48)(cid:48) /yr . In fact, the only significant clustering of ob-jects with precession rates in this range in the mod-ern belt are members of the high-inclination collisionalfamily (31) Euphrosyne near ∼ D to P S /P J < D satisfyingsystems to those that finish with the planets interior tothe 5:2 MMR. q , Q ( a u ) Example of | M | > | M | and P S /P J > . from primordial 3:2 JupiterSaturn -2 -1 Time (Myr) P S / P J Figure 3 . Example instability evolution from Clement et al.(2019a). Using the nomenclature of that work, this sim-ulation was from set 5
GP/ Myr (indicating it took the5 GP control initial conditions and delayed the instability 1Myr after the start of terrestrial planet formation). The sim-ulation finished with P S /P J = 2.58 and M =.023. The toppanel plots the perihelion and aphelion of each planet overthe length of the simulation. The bottom panel shows theJupiter-Saturn period ratio. The horizontal dashed lines inthe upper panel indicate the locations of the giant planets’modern semi-major axes. The shaded region in the lowerpanel delimits the range of 2.3 < P S /P J < Definition of simulation success
It is challenging to study a chaotic event like theNice Model instability because it is impossible to knowwhether the actual outer solar system is a typical out-come of all possible evolutionary paths that could havebeen followed from the unknown authentic initial con-ditions, or an outlier. Section 4 of Nesvorn´y & Mor-bidelli (2012) provides a good synopsis of this somewhatphilosophical conundrum. Following their example, weanalyze our simulations’ rates of satisfying our successcriteria, A - D , both individually and collectively. Thisis important because our success criteria jointly analyzenearly a dozen different parameters that may or maynot be correlated with one another (table 2). It seemsreasonable to expect that any individual set of success-ful initial conditions might only meet all four criteriasimultaneously a few times, if at all, in ∼ (cid:39) A providedthat the N GP = 4 systems are not systematically differ-ent than the others. Conversely, if the N GP = 3 systemsare the only ones within the batch that satisfy criterion B , C and D , the simulation set would be a failure interms of criterion A . In both scenarios, the total successrate for satisfying all 4 criteria simultaneously might bezero. However, the former case would only possess anull population of totally successful simulations as a re-sult of an over-multiplication of constraints. Thus, ifa large number of systems satisfy 2 or 3 of the criteriawithout preference to a particular combination, and onlynarrowly fail the others (for example P S /P J = 2.51 or M = 0.024), the set of simulations could still representa viable evolutionary path for the solar system. Becauseof this, we focus our analysis on both the success criteriathemselves and how mutually exclusive they are. NUMERICAL SIMULATIONS3.1.
Generating eccentric resonant chains
We perform nearly 6,000 instability simulations us-ing the
M ercury ∼ a ) and eccentricity damping ( ˙ e ) terms (see fullderivation in Lee & Peale 2002). To initially place theplanets in resonance, we establish a form of ˙ a = ka and˙ e = ke/
100 (Batygin & Brown 2010), where k is set toachieve a migration timescale of τ mig (cid:39) M ercury a J < . e on Jupiter or Saturn, or reversing its sign.We find that both planets’ orbits are excited most effi-ciently in 2:1 MMR gas giant configurations by pumpingSaturn’s eccentricity and simultaneously reducing thedegree of eccentricity damping on Jupiter. 3:2 MMR re-quire the opposite prescription. Furthermore, excitingSaturn’s eccentricity while maintaining e J low (closer to (cid:39) e J,o and e S,o . We then ver-ify that the planets are indeed locked in a mutual res-onant chain by checking for libration about a series ofresonant angles using the method described in Clementet al. (2018). For the precise values used to constructeach individual resonant chain, we refer the reader tothe online supplementary data files.Of note, though beyond the scope of this paper,we find that certain configurations of primordial giantplanet eccentricities are more difficult to generate thanothers. In particular, some architectures are far less sen-sitive to the specific migration and eccentricity pump-ing sequence, and require substantially less fine tuningto produce. In general, we find higher- e J /lower- e S con-figurations to be less sensitive to changes in initial con-ditions within the 2:1 MMR, while the opposite com-bination (lower- e J /higher- e S ) is more easily producedwith Jupiter and Saturn in a 3:2 MMR. Indeed Jupiter,being more massive, is more resilient to over-excitationduring the eccentricity-pumping process. As such, whenJupiter and Saturn are placed in the more isolated 2:1MMR, a wider range of values for ˙ e J can be used thancan for ˙ e S . For example, when generating 5 planet, 2:1configurations we find that just a ∼
1% change in ˙ e S canbe the difference between not exciting Saturn at all, ade-quate excitation, and pumping the planet’s eccentricityto the point of ejection. On the other hand, we findthat ˙ e J can be altered by several orders of magnitudeand still yield the desired result. In the case of the 3:2MMR, the strong eccentric forcing of Jupiter on Sat-urn due to the planets’ close proximity makes it easy toover-excite Saturn while pumping Jupiter’s eccentricity.3.2. Triggering the instability
Our goal is to study the relationship between the pri-mordial eccentricities of the gas giants and their post-instability secular architecture. However, the orbits ofour eccentricity-pumped systems of resonant giant plan-ets tend to damp to near-zero eccentricity if they areallowed to interact with the external planetesimal diskfor a significant period of time before fully destabilizing(this occurred within about ∼
10 Myr in our tests. Notethat our model does not include disk self-stirring: Levi-son et al. 2011). In order to prevent this from happen-ing in our simulations, we opt for an artificial instabilitytrigger in order to originate each batch of instabilitiesfrom nearly the same values of e J,o and e S,o . As a con-sequence of this selection, we are unable to constrainour final systems by Neptune’s pre-instability migra-tion history (Nesvorn´y 2015a,b; Deienno et al. 2017). IfNeptune’s primordial migration was indeed significant,our instability simulations will yield a systematic under-estimate of the ice giants’ final semi major axes. How-ever, since we are most interested in finding sets of ini-tial conditions that best construct the modern Jupiter-Saturn system (section 2.3), we relegate a detailed studyof the ice giants’ migration to future work. Additionally,the growing body of work arguing for an early instability(appendix A: Morbidelli et al. 2018; Clement et al. 2018;Deienno et al. 2018; Nesvorn´y et al. 2018) is potentiallyconsistent with our simulated scenario.We utilize the same instability trigger as Nesvorn´y(2011) and Nesvorn´y & Morbidelli (2012). Once eachsystem of giant planets is placed in resonance, we addan external disk of 1,000 equal-mass planetesimals witha total mass of 20.0 M ⊕ (loosely based off Nesvorn´y & Morbidelli 2012) It should be noted that these diskparticles are unrealistically massive (Levison et al. 2011;Quarles & Kaib 2019), and thus our selection of initialmasses is a compromise necessary to limit the computa-tional cost of our study. In future work, we intend to in-vestigate the evolution of these systems using more real-istic primordial planetesimal disks composed of (cid:38) r − (Williams & Cieza2011). Eccentricities and inclinations are chosen fromRayleigh distributions ( σ e = 0.01, σ i = 1 ◦ ; note thatthese might be unrealistically small if, for example, theywere excited earlier by the ice giant’s accretion: Ribeiroet al. 2020). The entire system is then integrated for100.0 Kyr with the M ercury ∼ ◦ to trigger theinstability (Nesvorn´y 2011). If an instability still doesnot ensue after 100 Myr, the simulation is discarded(e.g.: Nesvorn´y et al. 2018). Each system is evolved foran additional 20 Myr after the onset of the instability(as determined by either an ice giant ejection or a stepchange in the Jupiter-Saturn period ratio) in order tocapture some of the remnant planetesimal disk’s inter-actions with the post-instability outer solar system.It is important to recognize that our instability trig-ger is somewhat ad hoc, particularly in light of the factthat the pre-instability migration of the ice giants hasbeen shown to be important for explaining the cap-ture of D-type trojan asteroids by Jupiter (Nesvorn´yet al. 2013) and the Kuiper Belt’s modern orbital struc-ture (Nesvorn´y 2015a,b). Moreover, Uranus and Nep-tune’s survival probabilities (and the corresponding suc-cess rates for criterion A and B ) are boosted when theplanets undergo significant pre-instability migration. Itshould also be noted that 20 Myr is likely insufficient tofully capture the giant planets’ residual migration phase(Nesvorn´y & Morbidelli 2012, for example, integrate for100 Myr after the instability). Because we are most in-terested in studying the Jupiter-Saturn secular system Note that our fiducial disk mass choice of 20.0 M ⊕ mightdrive excessive residual migration in some of our 2:1 chains, thuscausing Saturn to migrate past the 5:2 MMR after the instability(e.g.: Nesvorn´y & Morbidelli 2012; Brasser & Lee 2015). For thisreason, we experiment with different total disk masses in section4.6. e J,o = e S,o =0.05 set: table 3) to a total integration time of 100 Myr.We find that the average change in P S /P J and M over this additional 80 Myr of integration is ∼ ∼ P S /P J < M > ∼ ν reso-nance in the asteroid belt; Clement et al. 2020). Wefurther discuss how this affects our results, in terms ofour criteria for simulation success, in section 2.3.Finally, we remove all Kuiper belt objects from eachfully evolved system and integrate only the remaining gi-ant planets for an additional 100 Myr to ensure that thefinal system is stable. We utilize the outputs of this fi-nal simulation phase to calculate the secular frequenciesand amplitudes via frequency-modulated Fourier Trans-form (ˇSidlichovsk´y & Nesvorn´y 1996, note that this stepis only completed for systems that retain both Jupiterand Saturn). Thus, the results presented in section 4are exclusively from systems that destabilize within 100Myr, and do not eject Saturn. For most of the sets ofinitial conditions we test, (cid:38)
85% of configurations desta-bilize within 100 Myr, and (cid:46)
1% eject Saturn.3.3.
Parameter space tested
Our various batches of numerical simulations are sum-marized in table 3. As discussed by previous authors(e.g.: Batygin & Brown 2010; Nesvorn´y 2011; Nesvorn´y& Morbidelli 2012; Deienno et al. 2017; Clement et al.2018), the parameter space of possible resonant chains,planetesimal disk properties, initial number of giantplanets and primordial eccentricities is extensive. For-tuitously, much of this phase space for low primordialeccentricities (Masset & Snellgrove 2001) has alreadybeen probed by the comprehensive study of Nesvorn´y& Morbidelli (2012). This permits us to simplify ourpresent manuscript by not repeating an analysis of pre-viously investigated parameter space. Instead, we per-form three sets of control simulations utilizing the most successful configurations from Nesvorn´y & Morbidelli(2012), in tandem with the identical disk conditions andcalculation mechanisms described in section 3. Thus,we are able to present a self-consistent comparison ofour present results with those of Nesvorn´y & Mor-bidelli (2012), and the 5 GP control / GP control sets fromClement et al. (2018, 2019a) discussed in section 2.2.The parameter space we seek to study is as follows: (1)a range of moderate primordial eccentricities for the gasgiants ( ∼ We find that, in general, exciting the eccentricitiesof Jupiter and Saturn with the planets in a primordial3:2 MMR (note that this case is not necessarily physi-cally motivated) leads to instability evolutions with ex-cessive jumps, and exceedingly high final values of e S .While primordial eccentricity pumping within the 3:2does boost success rates for criterion C (i.e.: the insta-bility more efficiently excites M , see section 4), thecorresponding success rates for the other 3 criteria aresystematically lower than in the zero-eccentricity controlcase. This is not surprising, as low-eccentricity instabili-ties beginning from the 3:2 are intrinsically more violentthan those from the 2:1 (section 2.1.1). In our scenario,pumping Jupiter’s primordial eccentricity only leads tosystematically stronger encounters during the chaos ofthe instability. In many cases, this leads to a strong scat-tering event between Jupiter and Saturn that essentiallydestroys the entire outer solar system. We were alsounable to develop a procedure for generating resonantchains with e J,o = e S,o within the 3:2 MMR’s tighterorbital spacing. Indeed, the pre-instability eccentricityof Saturn for e J,o ≈ ∼ e J,o to ∼ C and D (re-gardless of final N GP ) is near zero (thus, the distributionin figure 1 moves to the right as a greater fraction of sim-1ulations experience excessively large jumps). Therefore,we restrict our current study to three control cases (table3, two of the most successful five planet configurations,and one six planet chain from Nesvorn´y & Morbidelli2012), and two additional, eccentricity-pumped chainswith e J,o (cid:39) e S,o (cid:39)
Broadly speaking, increasing the gas giants’ primor-dial eccentricities leads to more violent evolutions forsystems with Jupiter and Saturn beginning in a 2:1MMR. Thus, these systems now behave similarly tothe circular, 3:2 cases (section 2.1.1) in that they rou-tinely eject one or more ice giants. To illustrate this,we present one set of instabilities where four giant plan-ets inhabit a 2:1-4:3-4:3 resonance chain ( e J,o ≈ e S,o ≈ . Chains beginning with 2:1,3:2
On the opposite end of the spectrum, as we discussin section 2.1.2, we were unable to consistently generateorbital instabilities in resonant chains where Jupiter andSaturn inhabit a 2:1 MMR, and Saturn and the first icegiant are in a 3:2. Even when the inner ice giant’s meananomaly is shifted, the planets in these widely-spacedconfigurations typically continue to migrate smoothly,and often fall back into resonance, as they interact withthe external planetesimal disk. Therefore, all the reso-nant chains we study that investigate the 2:1 Jupiter-Saturn resonance (table 3) place the innermost ice giantin a 4:3 MMR with Saturn.In summary, the majority of our simulations concen-trate on the primordial 2:1 Jupiter-Saturn resonancewith moderate eccentricities ( e J,o ≈ e S,o ≈ M IG = 16.0 M ⊕ for fiveplanet configurations and M IG = 8.0 M ⊕ for six planetcases. This choice is based off the study of Nesvorn´y &Morbidelli (2012), who find six planet instabilities withmore massive additional planets to be too violent. Wealso investigate the effect of varying M IG with addi-tional simulations in section 4.7. As discussed in section2.3, a major limitation of our work is that it lacks adetailed analysis of the potential ice giant evolutionarypaths. Specifically, we do not test tighter resonances(e.g.: 5:4, 6:5, etc.) or initial conditions derived frommodels of ice giant formation (Ribeiro et al. 2020) andobliquity evolution (Izidoro et al. 2015). Moreover, itis noticeably difficult to ascertain whether our initialconditions are realistic and physically motivated in theabsence of robust hydrodynamical studies that follow the growth and evolution of the entire outer solar sys-tem in the gas disk phase. For instance, it is unclearwhether it is possible for the ice giants to migrate pasta mutual 3:2 MMR and become trapped in any of thetighter first order resonance (specifically the 4:3) oftentested in N-body studies (e.g.: Batygin & Brown 2010;Batygin et al. 2011; Nesvorn´y 2011; Batygin et al. 2012;Nesvorn´y & Morbidelli 2012; Deienno et al. 2017; Gomeset al. 2018; Quarles & Kaib 2019). Indeed, the popula-tion of detected exoplanets locked in the 4:3 resonancehas been argued to be in excess of theoretical predictions(Rein et al. 2012; Matsumura et al. 2017; Brasser et al.2018). To guarantee the capture of two planets in a par-ticular resonance via convergent migration in the circu-larly restricted three body approximation, the smallerobject’s eccentricity must be less than a critical value(see derivations in Henrard & Lemaitre 1983; Petro-vich et al. 2013). For the innermost ice giants’ capturewith Saturn, this maximum eccentricity is ∼ ∼ RESULTSTable 4 gives the percentage of systems in each ofour various simulation batches (table 3) that satisfy oursuccess criteria, A - D (table 2). Our analysis is struc-tured as follows: the subsequent four sections (4.1-4.4)focus on the various resonant chains we test, section4.5 discusses the dependency of our results on e J,o and e S,o , and sections 4.6-4.7 present an additional suite ofsimulations (based off our most successful sets of initialconditions) where we vary the initial masses of the in-nermost ice giants and the planetesimal disk. Becausethe goal of our work is to find sets of initial conditionswhere the solar system is not an outlier in M - P S /P J N pln Resonant Chain e J,o e S,o e IG,o a N,o N sim Table 3 . Summary of resonant chains tested. The columns are as follows: (1) the initial number of giant planets, (2) theresonant chain, beginning with the Jupiter-Saturn resonance, (3-4) the initial eccentricities of Jupiter and Saturn, (5) the initialeccentricities of the ice giants in order of increasing semi-major axis, (6) the initial semi-major axis of the outermost ice giant,and (7) the total number of simulations that produce an instability within 100 Myr and retain both Jupiter and Saturn. In allcases, the planetesimal disk mass is set to 20.0 M ⊕ , the radial offset between Neptune and the disk is to 1.5 au, and Jupiter’sinitial semi-major axis is 5.6 au. All four and five giant planet configurations have ice giant masses of 16.0 M ⊕ . Ice giants insimulations beginning with six giant planets have masses of 8.0, 8.0, 16.0 and 16.0 M ⊕ , with increasing semi-major axis. Notethat, for control configurations (based off the most successful five and six planet systems from Nesvorn´y & Morbidelli 2012, anddenoted here by the “*” symbol) with e J,o and e S,o labeled as 0.0, the actual eccentricities are non-zero, but small. space, we provide a plot of this parameter space for eachof our tested resonant chains in figures 4, 5, 6, 9, 10, 12,13, and 14.4.1.
Control runs: the 3:2 Jupiter-Saturn resonanceand the circular 2:1 case
Nesvorn´y & Morbidelli (2012) favored several differ-ent aspects of both the five planet, 3:2,3:2,3:2,3:2 and3:2,3:2,2:1,3:2 chains’ evolution. For consistency, wepresent a set of simulations for each chain. In gen-eral, the broader spacing of the 3:2,3:2,2:1,3:2 leads to higher success rates for criterion B , and the two setsof initial conditions perform similarly when scrutinizedagainst our other constraints. However, we find thatthe 3:2,3:2,3:2,3:2 case yields the best results in termsof final distributions in M - P S /P J space (figure 4). Asour work seeks to find sets of initial configurations thatbest yield realistic magnitudes of M and low jumps,we focus our analysis on this chain, and the six planet,3:2,4:3,3:2,3:2,3:2 configuration.We find that the rates of success for our control sim-ulations (no eccentricity pumping) are similar to those3 N pln Resonant Chain e J,o e S,o
A B C D ALL
Table 4 . Summary of results for our different simulations (table 3) measured against our various success criteria (table 2).The columns are as follows: (1) the initial number of giant planets, (2) the resonant chain, beginning with the Jupiter-Saturnresonance, (3-4) the initial eccentricities of Jupiter and Saturn, (5) the percentage of systems satisfying criterion A ( N GP =4),(6) criterion B (the planets’ final semi-major axes within 20% of the real ones), (7) criterion C ( | ∆ M ij /M ij,ss | < i, j =5,6), M > M ), (8) criterion D ( P S /P J < reported by Nesvorn´y & Morbidelli (2012) and Clementet al. (2018) (our 5 GP control and 6 GP control sets pre-sented in section 2.2). Compared to these previous stud-ies, our control simulations possess slightly lower successrates for criterion B . This is a consequence of the icegiants’ residual migration being subdued in our simula-tions due to shorter integration times and a lower plan-etesimal disk mass for the five planet case (20 M ⊕ asopposed to ∼
35; we revisit the role of the disk’s mass insection 4.6).All of our simulations beginning from the 3:2 Jupiter-Saturn resonance systematically struggle to satisfy cri-terion C . As discussed in section 2.3, establishing cri-terion C such that any simulation with M > | M | > | M | regime, satisfying the constraint. In-deed, 75% of our five planet, 3:2,3:2,3:2,3:2 chains with-out primordial eccentricity pumping (79% when we in-clude primordial excitation) excite M to greater than0.022 (note that the majority of these are in system’sthat fail criterion A , the corresponding rate reportedby Nesvorn´y & Morbidelli 2012, is lower because theauthors only claim success for criterion C if A and B are met as well). However, only 8% of the runs in oursame control five planet simulation batch satisfy our up-dated version of criterion C that scrutinizes all foureccentric amplitudes of the Jupiter-Saturn system. Ina similar manner, our simulations originating from sixplanet, zero-eccentricity, 3:2,4:3,3:2,3:2,3:2 chains finish4with M > C at a rate of 5%. The mostdifficult amplitude for our 3:2 Jupiter-Saturn resonancesimulations to match is M , with only 12% of systemsfinishing within the appropriate range. Conversely, thenumber of simulations possessing proper values of M , M and M are 32% (note that this is not ∼
70% asdiscussed above by virtue of our new constraint impos-ing a maximum limit on M as well as a minimum),54% and 27%, respectively. Furthermore, 67% of thesystems that finish with adequate values of M are inthe | M | > | M | regime. Thus, while it is importantto avoid over-constraining instability simulations whenstudying the Nice Model statistically, our results areindicative of a strong anti-correlation between the ad-equate excitation of M and the broad replication ofthe complete Jupiter-Saturn system for the primordial3:2 resonance. While high M magnitudes are commonoutcomes in our simulation sets studying these initialconditions, they occur preferentially in systems that ex-perience large jumps (i.e.: those that fail criterion D )and possess overexcited values of M .4.1.2. Our simulation sets investigating the effects of mildprimordial eccentricity pumping in the 3:2 Jupiter-Saturn resonance yield success rates for criteria A - D that are consistently similar or worse to those generatedin the low- e case. In general, exciting the eccentricitiesof Jupiter and Saturn in the 3:2 MMR scenario tends toinject more violence into an already brutal event. In ourbatch studying five planet configurations, this manifestsas low success rates for criteria A and B . Indeed, only1% of these simulations satisfy both constraints simul-taneously. As discussed in section 3, stable 3:2 resonantchains with heightened values of e J,o and e S,o are morechallenging to generate because, in the tighter config-uration, Jupiter’s dynamic excitation readily bleeds tothe other giant planets via stochastic diffusion. In spiteof the fact that we maintain artificial damping forces onthe ice giants throughout the process of generating theseresonant chains, the inner two ice giants typically pos-sess eccentricities around ∼ N GP = 4 contain Uranus and Nep-tune analogs with orbits that are systematically moredistant and excited than the real ones. This effect is less severe in our six planet configurations because theinner two ice giants are less massive (8.0 M ⊕ as opposedto 16.0 M ⊕ ). In this scenario, the more massive outer icegiants that typically go on to become the Uranus andNeptune analogs in criterion A satisfying systems aremore dynamically detached from the excited gas giants,and possess lower eccentricities (e ∼ B .Our simulations indicate that the primordial 3:2Jupiter-Saturn resonance still represents a viable post-formation evolutionary pathway for the solar system(Batygin & Brown 2010; Nesvorn´y 2011; Nesvorn´y &Morbidelli 2012; Batygin et al. 2012; Deienno et al. 2017;Clement et al. 2018). Indeed, our tested six planetconfigurations (with and without primordial excitation)each produce one system that simultaneously satisfiesall 4 of our success criteria. This is also the case forour five planet, 3:2,3:2,2:1,3:2 chain, and our five planetset that includes primordial eccentricity pumping. Aspointed out by previous authors, the challenge with the3:2 version of the Nice Model (section 2.1.1) is thatit is extremely violent. Thus, only ∼
20% of simula-tions experience appropriately “weak” instabilities, andyield correspondingly realistic Jupiter-Saturn period ra-tio jumps (criterion D , this manifests as fewer pointsplotted in figure 4 compared to the corresponding plotsfor our 2:1 simulations). Our results indicate that thereis very little difference in the performance of this subsetof weaker instabilities when e J,o and e S,o are artificiallyelevated. The main difference between the two sets isthat the total size of the population of less-violent out-comes decreases with primordial eccentricity pumping.This is evidenced by a 25% success rate for criterion D in our low- e six planet batch, compared to just 5%in the moderate- e case. Therefore, we find that pri-mordial eccentricity pumping in the 3:2 Jupiter-Saturnresonance does not bring the solar system result closerto the center of the distribution of possible outcomesin M - P S /P J space (figure 4). Instead, it pushes theactual Jupiter-Saturn period ratio farther towards theextreme of possible results.4.1.3. Figure 5 plots the results of our 2:1 control runs in M - P S /P J space. In contrast to the 3:2 cases depictedin figure 4, no simulation in our circular 2:1 batch fin-ishes with M excited to at least the solar system valuewithout exceeding P S /P J = 2.5. This is consistent withthe results of Nesvorn´y & Morbidelli (2012) and Dei-5 -3 -2 -1 Five planets; 3:2,3:2,3:2,3:2 e J , o = . Six planets; 3:2,4:3,3:2,3:2,3:22.0 2.1 2.2 2.3 2.4 2.5 2.6 2.710 -3 -2 -1 e J , o = . | M | P S /P J Figure 4 . Similar to figure 1. M - P S /P J space for our simulations beginning from the 3:2 Jupiter Saturn resonance (ta-ble 3). The left two plots depict five giant planet chains (3:2,3:2,3:2,3:2) and the right plots show runs with six planets(3:2,4:3,3:2,3:2,3:2). The plots on top are control runs without primordial eccentricity pumping. The bottom two panels haveminor eccentricity pumping ( e J = 0.025, e S = 0.05). The red stars indicate modern solar system values. The shaded grey areadelimits the region of 2.3 < P S /P J < D ) and 0.022 < M < C ). While one simulation providesin the five planet, e J = 0.025, e S = 0.05 batch (bottom right panel) provides an excellent match to the Jupiter-Saturn system,this is not a typical outcome of these initial conditions. Indeed, 86% of these simulations finish with P S /P J > enno et al. (2017). In both studies, the authors con-cluded that it is extremely difficult to adequately ex-cite Jupiter’s eccentricity out of the primordial 2:1 res-onance. Indeed, while the individual rates of success foreach of our four success criteria are reasonable for ourcircular 2:1 instabilities (table 4), none are successful atsimultaneously satisfying C and D . Specifically, thesesystems systematically struggle to adequately excite theeccentricities of both Jupiter and Saturn when the insta-bility yields an appropriately small jump ( P S /P J < M to greater than half its modern value, and only 10% pos-sess M magnitudes in excess of 0.016 (half the modernmagnitude). As none of our control, 2:1 simulations aresuccessful at simultaneously satisfying criteria C and D ,we conclude that some degree of primordial excitationis likely an important prerequisite to the viability of any2:1 instability scenario.4.2. In general we find that the primordial 2:1 Jupiter-Saturn resonance with heightened eccentricities is also aviable evolutionary pathway for the outer solar system.In many ways, our 2:1 batches of simulations outperform the 3:2 sets discussed in the previous section. However,we caution the reader that our results should be takenas motivation for follow-on study of the 2:1 resonance,and not as reason to abandon the 3:2. As previouslydiscussed, the 2:1 version of the Nice Model’s (section2.1.2) effects on the solar system’s global dynamics arenot as well-studied as are those of the 3:2 (e.g.: Nesvorn´yet al. 2013; Nesvorn´y 2015a,b; Roig & Nesvorn´y 2015;Roig et al. 2016). In particular, the asteroid belt ismost sensitive to the instability’s particular dynamics(Deienno et al. 2016, 2018) and the precise motions ofthe dominant secular resonances’ locations in belt (Mor-bidelli et al. 2010; Izidoro et al. 2016; Clement et al.2020). Thus, our simulations of the 2:1 instability arelimited in that they only analyze the various sets of ini-tial conditions’ success at replicating the modern orbitalconfiguration of the four giant planets (success criteria A - D ). Future work to fully validate the scenario mustfocus on the evolution of the orbital distributions in theasteroid (e.g.: Roig & Nesvorn´y 2015) and Kuiper (e.g.:Nesvorn´y & Vokrouhlick´y 2016) belts, and the obliquityevolution of the giant planets (Vokrouhlick´y & Nesvorn´y2015; Brasser & Lee 2015).We also find it difficult to finely control the insta-6 P S /P J -3 -2 -1 | M | Circular 2:1 case
Figure 5 . M - P S /P J space for our simulations beginning from the circular 2:1 Jupiter Saturn resonance (table 3). The redstars indicate modern solar system values. The shaded grey area delimits the region of 2.3 < P S /P J < D ) and0.022 < M < C ). bility’s timing, and therefore minimize the amount ofdamping in e J,o and e S,o that occurs prior to the insta-bility. Even with our artificial instability trigger (section3), systems often take a few Myr (the median instabilitytime for our simulation batches varies between ∼ t inst and the finalproperties of our simulated systems that are statisticallysignificant, we remind the reader that the gas giants in asubset (albeit, only a small set) of the systems analyzedin these sections have damped to near-zero eccentricityby the time the instability ensues.Figure 6 plots the results of our 2:1,4:3,3:2,3:2 instabil-ities in M - P S /P J space (we refer to this as our ”loose,”“broad” or “wide” configuration in the subsequent text).We find these wider resonant chains to generally be moresuccessful than our compact (or “tight”) configurationsof five planets (section 4.3). As our preliminary workindicated that even broader chains beginning with 2:1(e.g.: 2:1,3:2,3:2,3:2, see discussion in section 3.3.3) sel-dom degenerate in to an instability (Deienno et al. 2017),the success of the 2:1,4:3,3:2,3:2 chain over tighter con-figurations places a fairly strict constraint on the rangepotentially viable five planet configurations. Interest-ingly, these systems possess rates of success for criteria C ( ∼ D ( ∼ M - P S /P J outcomes for our 2:1 instabilities (red star in figure 6).Thus, while high | M | values occur preferentially in sys-tems that experience a large jump out of the primordial3:2 resonance, properly excited Jupiter analogs occurwith similar frequencies in systems across the full spec-trum of P S /P J outcomes in simulations that begin fromthe 2:1 resonance. This is largely a consequence of M already being excited prior to the instability. Additionally, our looser, 2:1,4:3,3:2,3:2 resonantchains yield systematically higher success rates for cri-teria A (30-40%) and B (20-25%) than our 3:2 controlcases, and our tighter 2:1, five planet chains. In a sense,these initial conditions place the planets in a more dy-namically isolated configuration that already somewhatresembles that of the modern solar system. Therefore,these systems typically undergo instabilities that areweaker and more regularly behaved than those experi-enced in more compact configurations (3:2,3:2,3:2,3:2, In general, the coefficients M ij ( i, j =5,6) of our initial giantplanet configurations are partitioned such that the magnitudes of M and M are approximately twice those of M and M priorto the instability. However, the precise relative values also dependon the whether or not e J,o > e
S,o . -3 -2 -1 e S,o = e J , o = . e S,o = -3 -2 -1 e J , o = . Five planets; 2:1,4:3,3:2,3:2 | M | P S /P J Figure 6 . M - P S /P J space for our simulations with five planets beginning from a 2:1,4:3,3:2,3:2 resonant chain. The uppertwo panels plot simulations begin with e J = 0.025, and the bottom panels show systems where Jupiter’s initial eccentricity is0.05. Similarly, the left two panels plot simulations begin with e S = 0.025, and the panels on the right show systems whereSaturn’s initial eccentricity is 0.05. The shaded grey area delimits the region of 2.3 < P S /P J < D ) and 0.022 < M < C ). e J,o and e S,o ( ∼ M - P S /P J outcomes. Moreover, both cases with e J,o = 0.05are more successful at satisfying criterion C and fullyreplicating the complete Jupiter-Saturn secular systemthan those with e J,o =0.025.Figures 7 and 8 focus on our most successful fiveplanet batch (2:1,4:3,3:2,3:2 with e J,o = e S,o = 0.05).While the solar system values of M ij ( i, j =5,6) roughlyfall within the 1- σ range of outcomes depicted in fig-ure 7, it is clear that M = 0.044 is somewhat closerto the extreme of the distribution for low- M . Con-versely, Saturn’s eccentric modes are reproduced quitefrequently in our simulations. However, our numericalintegrations do not fully capture the planets’ residualmigration phase. Thus, we expect the distribution ofJupiter’s eccentric modes in the top panel of figure 7 tomove slightly down towards the solar system outcome asthe system continues to evolve (by ∼ e J are a consequence of the aforemen-tioned challenges with properly exciting M in weakerinstabilities that eject only one planet, the ice giant’sfinal orbital locations are most sensitive to interactionswith the remnant planetesimal disk (Nesvorn´y & Mor-bidelli 2012). There are two main factors that are re-sponsible for our simulated ice giants’ inability to attaintheir modern semi-major axes. First, our simulationsonly model 20 Myr of the residual migration phase (acompromise necessary to limit the computational costof our work). Second, the amount of residual migrationis a function of the total remnant mass in planetesi-mals; which in turn is related to the initial mass placedin the planetesimal disk. For consistency, our simula-tions strictly consider M KB = 20 M ⊕ disks. However,Nesvorn´y & Morbidelli (2012) found that massive disks( ∼ M ⊕ ) are more successful in more-compact fiveplanet configurations. This is an obvious consequence ofcertain wider configurations with more planets requir-ing less residual migration for the planets to reach theirmodern orbits (of note, Nesvorn´y & Morbidelli 2012,also prefer lighter disks in the wider 3:2,3:2,2:1,3:2). Fur-8 -3 -2 -1 | M | -3 -2 -1 | M | Jupiter -3 -2 -1 | M | -3 -2 -1 | M | Saturn
Figure 7 . Eccentric magnitudes, M ij ( i, j =5,6), for theJupiter-Saturn systems that satisfied criterion D in our batchof simulations beginning with five planets in a 2:1,4:3,3:2,3:2resonant chain and e J,o = e S,o =0.05. The top panel showsthe magnitudes M and M in Jupiter’s eccentricity, thebottom panel depicts the magnitudes M and M in Sat-urn’s. The error bars indicate one standard deviation. Thered stars correspond to solar system values. thermore, the relationship between M KB and simulationsuccess is slightly more complicated because the resid-ual migration phase also damps out the eccentric modesof the Jupiter-Saturn system. Thus, the selection ofinitial M KB can be considered a sort of balancing act.Too much residual migration can over-damp the gas gi-ants and lead to failure of criterion C , while too littlecan result in the ice giants stopping short of their mod-ern semi-major axes and lead to failure of criterion B (Gomes et al. 2004). We explore the consequences ofthe particular choice of M KB further with additionalsimulations in section 4.6.4.3. The results of our tighter configurations of 2:1,five planet instabilities are plotted in figures 9(2:1,4:3,4:3,3:2) and 10 (2:1,4:3,4:3,4:3). While thesesets still produce successful systems (including several E cce n t r i c it y Semi Major Axis (au) I n c li n a ti on ( ◦ ) Figure 8 . Final orbits for systems that satisfied criteria A and B in our batch of simulations beginning with five planetsin a 2:1,4:3,3:2,3:2 resonant chain and e J,o = e S,o =0.05. Thetop panel depicts a/e space and the bottom panel plots a/i space. The respective planets and their simulated analogsare color coded as follows: Jupiter in red, Saturn in gold,Uranus in light blue, Neptune in dark blue. The error barsindicate one standard deviation. The stars correspond tosolar system values. simulations that simultaneously satisfy all four successcriteria; table 4), the rates of success for every suc-cess criterion except C are systematically lower thanfor our looser, 2:1,4:3,3:2,3:2 configuration. In partic-ular, a greater number of these more compact systemstend to devolve into extremely violent instabilities. Thenet result of this is more planet ejections (lower successrates for criterion A ), and larger scattering events (lowersuccess rates for criteria B and D ). As with the primor-dial 3:2 resonance (section 4.1), the ice giants in our2:1,4:3,4:3,3:2 and 2:1,4:3,4:3,4:3 configurations attainhigher eccentricities before the instability, and tend toexperience stronger mutual encounters within the chaosof the instability. To illustrate this, we analyzed theclose encounter histories for Saturn in our tightest andloosest five planet chains with e J,o = e S,o = 0.05. Onaverage, we find that encounters less than 3 Hill radii be-9 -3 -2 -1 e S,o = e J , o = . e S,o = -3 -2 -1 e J , o = . Five planets; 2:1,4:3,4:3,3:2 | M | P S /P J Figure 9 . M - P S /P J space for our simulations with five planets beginning from a 2:1,4:3,4:3,3:2 resonant chain. The uppertwo panels plot simulations begin with e J = 0.025, and the bottom panels show systems where Jupiter’s initial eccentricity is0.05. Similarly, the left two panels plot simulations begin with e S = 0.025, and the panels on the right show systems whereSaturn’s initial eccentricity is 0.05. The shaded grey area delimits the region of 2.3 < P S /P J < D ) and 0.022 < M < C ). -3 -2 -1 e S,o = e J , o = . e S,o = -3 -2 -1 e J , o = . Five planets; 2:1,4:3,4:3,4:3 | M | P S /P J Figure 10 . M - P S /P J space for our simulations with five planets beginning from a 2:1,4:3,4:3,4:3 resonant chain. The uppertwo panels plot simulations begin with e J = 0.025, and the bottom panels show systems where Jupiter’s initial eccentricity is0.05. Similarly, the left two panels plot simulations begin with e S = 0.025, and the panels on the right show systems whereSaturn’s initial eccentricity is 0.05. The shaded grey area delimits the region of 2.3 < P S /P J < D ) and 0.022 < M < C ). Recall that the upper right quadrant is blank as we were unable to generate a stable chain with e S > e J for this configuration . q , Q ( a u ) Example five planet evolution from primordial 2:1
JupiterSaturn -2 -1 Time (Myr) P S / P J Figure 11 . Example instability evolution beginning with fiveplanets in a 2:1,4:3,3:2,3:2 resonant chain (figure 6). The sim-ulation finished with P S /P J = 2.41, M =.034, M =.018, M =.027, M =.052 (all four success criteria are satis-fied). The top panel plots the perihelion and aphelion ofeach planet over the length of the simulation. The bottompanel shows the Jupiter-Saturn period ratio. The horizontaldashed lines in the upper panel indicate the locations of thegiant planets’ modern semi-major axes. The shaded regionin the lower panel delimits the range of 2.3 < P S /P J < tween Saturn and the ice giants are 11% more frequentand 3% closer in the tight batch than the correspondingsimulations beginning from a looser configuration. Anexample evolution for a system that satisfies all four suc-cess criteria from our 2:1,4:3,3:2,3:2, e J,o = e S,o =0.05set is plotted in figure 11. It is clear that, even in themost successful system, Uranus and Neptune are over-excited in the instability.Another way to inspect the differences between ourrespective chains is by analyzing the rates at which agiven batch satisfies a specific subset of our four suc-cess criteria simultaneously. For example, our various2:1,4:3,4:3,4:3 simulation batches have success rates ofonly 1-4% for criterion B . As criterion B can only besatisfied when criterion A is already met, we can easilycompare the ratio of criterion B satisfying simulations toall those that finish with N GP = 4 ( A ) between our var-ious different resonant chains. The two more compactconfigurations already possess lower success rates for cri-terion A (7-19% vs. ∼ ∼ B , compared to (cid:38)
60% ofthose that originated in the looser, 2:1,4:3,3:2,3:2 chain. On closer inspection, we find that the majority of thesecriterion A satisfying systems experience uncharacteris-tically weak instabilities that leave the giant planets in afinal orbital configuration that is too compact. An addi-tional several runs satisfy criterion A , but fail B as a re-sult of a series of scattering events between the ice giantsthat drive the Neptune analog’s semi-major axis into thedistant Kuiper belt. The tendency of these compact sys-tems to only finish with 4 planets in weak instabilitiesis also evidenced by the various percentages of systemsthat meet all but one of our constraints. For instance,only 5% of the 2:1,4:3,4:3,3:2, e J,o = e S,o = 0.05 batchsuccessfully meet both A and B . Of this subset of 10simulations, nine systems experience a small jump andsatisfy D , while only one is strong enough to properly re-produce M ( C ). In summary, while our more compactfive planet configurations do yield successful evolution-ary schemes, they also suffer multiple systematic issuesthat do not affect our wider five planet configurationsas severely. 4.4. Perhaps the most striking difference between our fiveand six planet cases is the six planet sets’ higher successrates for criteria A and D . These two results are not mu-tually inclusive, in that N GP = 3 systems satisfy crite-rion D at roughly equal rates as N GP = 4 systems. Oursix planet instabilities are typically weaker, and tendto behave more consistently than the counterpart fiveplanet runs (note, however, that our six planet configu-rations are rather artificial in terms of our inclusion oftwo additional, M IG = 8.0 M ⊕ planets). While five andsix planet cases are roughly equally successful at excit-ing M ( ∼ M > C orNesvorn´y & Morbidelli 2012), six planet instabilities areless likely to over-excite M . Thus, our six planet runspossess success rates for our updated criterion C that aresystematically the same or better than those of our fiveplanet simulations. This result, and the prevalence of N GP = 4 systems are partially a consequence of the factthat the scattering ice giants have lower masses than inour five planet instabilities. As with our six planet, 3:2control runs (section 4.1), the outer two, more massive,ice giants typically survive the instability and becomeUranus and Neptune analogs. In fact, this is the casein every one of our six planet simulations that satisfyall four of our constraints. As the two outermost plan-ets begin the simulation more dynamically isolated fromthe gas giants, and with lower eccentricities (e (cid:46) M - P S /P J space.1 -3 -2 -1 e S,o = e J , o = . e S,o = -3 -2 -1 e J , o = . Six planets; 2:1,4:3,3:2,3:2,3:2 | M | P S /P J Figure 12 . M - P S /P J space for our simulations with six planets beginning from a 2:1,4:3,3:2,3:2,3:2 resonant chain. The uppertwo panels plot simulations begin with e J = 0.025, and the bottom panels show systems where Jupiter’s initial eccentricity is0.05. Similarly, the left two panels plot simulations begin with e S = 0.025, and the panels on the right show systems whereSaturn’s initial eccentricity is 0.05. The shaded grey area delimits the region of 2.3 < P S /P J < D ) and 0.022 < M < C ). Recall that the upper right quadrant is blank as we were unable to generate a stable chain with e S > e J for this configuration . -3 -2 -1 e S,o = e J , o = . e S,o = -3 -2 -1 e J , o = . Six planets; 2:1,4:3,4:3,3:2,3:2 | M | P S /P J Figure 13 . M - P S /P J space for our simulations with six planets beginning from a 2:1,4:3,4:3,3:2,3:2 resonant chain. The uppertwo panels plot simulations begin with e J = 0.025, and the bottom panels show systems where Jupiter’s initial eccentricity is0.05. Similarly, the left two panels plot simulations begin with e S = 0.025, and the panels on the right show systems whereSaturn’s initial eccentricity is 0.05. The shaded grey area delimits the region of 2.3 < P S /P J < D ) and 0.022 < M < C ). -3 -2 -1 e S,o = e J , o = . e S,o = -3 -2 -1 e J , o = . Six planets; 2:1,4:3,4:3,4:3,3:2 | M | P S /P J Figure 14 . M - P S /P J space for our simulations with six planets beginning from a 2:1,4:3,4:3,4:3,3:2 resonant chain. The uppertwo panels plot simulations begin with e J = 0.025, and the bottom panels show systems where Jupiter’s initial eccentricity is0.05. Similarly, the left two panels plot simulations begin with e S = 0.025, and the panels on the right show systems whereSaturn’s initial eccentricity is 0.05. The shaded grey area delimits the region of 2.3 < P S /P J < D ) and 0.022 < M < C ). A - D . For instance, our looserconfiguration (2:1,4:3,3:2,3:2,3:2) is more likely to expe-rience a small jump (36-43% of systems satisfying cri-terion D ), but less successful at meeting criterion B .Utilizing a more compact initial configuration booststhe probability of success in terms of the planets finalsemi-major axes (20-30% of systems satisfying B ), butreduces the total sample of systems experiencing smalljumps. Thus, the choice of initial configuration is a com-promise between more efficiently limiting the gas giants’jump (with a looser chain, and weaker instability), andplacing the eventual Uranus and Neptune analogs at theright initial semi-major axes to produce a successful finalorbital configuration (a tighter chain). A more compactinitial orientation of the planets also tends to be moresuccessful at adequately exciting the secular eigenmodesof the Jupiter-Saturn system, while looser chains aremore likely to under-excite Jupiter or Saturn’s eccen-tricity (specifically, (cid:38)
75% of these systems under-excite M or M ).The relationship between e J,o , e S,o and simulation suc-cess is also more complicated in our six planet batchesthan those that study five planets. Specifically, higherinitial values of e J,o ( (cid:39) C , as they are more likely to leadto adequately excited values of M . Conversely, slightlylower initial eccentricities for Saturn ( (cid:39) D ).However, it is important to note that this combinationof primordial eccentricities is only produced from a spe-cific combination of disk parameters in hydrodynamicmodels (Pierens et al. 2014, discussed further in section4.5). Figures 15 and 16 plot the final distributions ofsecular eigenmodes and 4 planet system orbits for this“ideal” combination of e J,o = 0.05, e S,o = 0.025 for our2:1,4:3,4:3,3:2,3:2 resonant chain. The overall distribu-tions of the magnitudes, M ij ( i, j = 5.6), is similar tothat of our preferred five planet orientation (figure 7),and a significant improvement from the 3:2 version ofthe Nice Model (e.g.: figure 4). The final planet orbitsin this simulation set also provide a fairly good match tothe actual solar system (figure 16), thus implying thata 20 M ⊕ planetesimal disk provides our six planet con-figuration with sufficiently strong encounters to driveUranus and Neptune towards their modern orbits in theresidual migration phase.Overall, our six planet configurations are similar to thefive planet cases in terms of their ability to reproduce im-portant aspects of the Jupiter-Saturn (e.g. figures 7 and15). However, their tendency to experience small jumps -3 -2 -1 | M | -3 -2 -1 | M | Jupiter -3 -2 -1 | M | -3 -2 -1 | M | Saturn
Figure 15 . Eccentric magnitudes, M ij ( i, j =5,6), forthe Jupiter-Saturn systems that satisfied criterion D inour batch of simulations beginning with six planets in a2:1,4:3,4:3,3:2,3:2 resonant chain and e J,o = e S,o =0.05. Thetop panel shows the magnitudes M and M in Jupiter’seccentricity, the bottom panel depicts the magnitudes M and M in Saturn’s. The error bars indicate one standarddeviation. The red stars correspond to solar system values. and finish with the correct number of planets makes thesix planet, 2:1 instability compelling. Indeed, five of 171simulations in our 2:1,4:3,4:3,4:3,3:2, e J,o = e S,o = 0.05satisfy all four of our success criteria. An example ofsuch a successful simulation is plotted in figure 17.4.5.
Jupiter versus Saturn’s primordial eccentricity
Our results (in comparison to previous study of the2:1: Nesvorn´y & Morbidelli 2012; Deienno et al. 2017)clearly illustrate that primordial eccentricity excitationis necessary for the 2:1 Jupiter-Saturn resonance to beviable (Pierens et al. 2014). It is worthwhile to point outthat our simulations do not indicate a significant depen-dency of successful outcomes on the particular choice of e J,o or e S,o . This is not extremely surprising, given thestochasticity of an event like the Nice Model, it is reason-able to expect small differences in initial conditions tobe largely erased during the violent event. However, it is4 E cce n t r i c it y Semi Major Axis (au) I n c li n a ti on ( ◦ ) Figure 16 . Final orbits for systems that satisfied criteria A and B in our batch of simulations beginning with six planetsin a 2:1,4:3,4:3,3:2,3:2 resonant chain and e J,o = e S,o =0.05.The top panel depicts a/e space and the bottom panelplots a/i space. The respective planets and their simulatedanalogs are color coded as follows: Jupiter in red, Saturn ingold, Uranus in light blue, Neptune in dark blue. The errorbars indicate one standard deviation. The stars correspondto solar system values. readily apparent from table 4 that higher initial valuesof e J,o tend to manifest as improved rates of success forcriterion C (the most difficult constraint to satisfy) aftersystems are evolved through the Nice Model instability.In comparison, Deienno et al. (2017) also studied the2:1,4:3,3:2,3:2 configuration without including primor-dial excitation and reported a 0% success rate for excit-ing M to greater than half the modern value. Thus,it is clear that some degree of primordial excitation forJupiter is critical for achieving successful instability out-comes from the primordial 2:1 Jupiter-Saturn resonance.However, we find that higher eccentricities for Saturnlead to systematically larger jumps, and correspondinglylow success rates for criterion D . Therefore, our simu-lations indicate that a moderate value of e J,o ( (cid:38) e S,o ( (cid:46) q , Q ( a u ) Example six planet evolution from primordial 2:1
JupiterSaturn -2 -1 Time (Myr) P S / P J Figure 17 . Example instability evolution beginning from the2:1,4:3,4:3,4:3,3:2 resonant chain (figure 14). The simula-tion finished with P S /P J = 2.46, M =.029, M =.022, M =.022, M =.068 (all four success criteria are satis-fied). The top panel plots the perihelion and aphelion ofeach planet over the length of the simulation. The bottompanel shows the Jupiter-Saturn period ratio. The horizontaldashed lines in the upper panel indicate the locations of thegiant planets’ modern semi-major axes. The shaded regionin the lower panel delimits the range of 2.3 < P S /P J < notable caveat of this result is that such a combinationof eccentricities is only produced from a specific com-bination of disk parameters in hydrodynamical models.Pierens et al. (2014) found that, in most cases, Saturnattains a higher eccentricity than Jupiter when the plan-ets are captured in the mutual 2:1 MMR. However, theauthors reported an outcome with e J,o > e
S,o for their α = 0.01, f = 0.3 disk (figure 4 in that paper). While thecases we test that consider e J,o < e
S,o do produce suc-cessful simulations (table 4), the outcomes of our moresuccessful sets of initial conditions ( e J,o (cid:38) e S,o ) shouldbe taken in the appropriate context given that hydro-dynamical models typically yield the opposite combina-tion.On the higher side of the range of possible primordialeccentricities for Jupiter, we were unable to generatestable chains with e J,o (cid:38) ∼ M )without exceeding P S /P J = 2.5. The most significantsystematic issues with our simulations are (1) inaccuratefinal ice giant semi-major axes, (2) a preponderance ofstrong instabilities that tend to eject too many planets,and (3) total integration times that are inadequate tofully capture the residual migration phase. As (1) islikely related to our initial selection of planetesimal diskmass and (2) is probably a result of the massive inner-most ice giants beginning the simulation on nearly over-lapping orbits, we conclude our study by varying M KB and M IG in an additional suite of simulations. The fol-lowing two sections discuss the results of this follow-onset of integrations.4.6. Varying the planetesimal disk’s mass
To study the dependence of our results on M KB , weperform four additional batches of 200 simulations basedoff our most successful five (2:1,4:3,3:2,3:2, e J,o = e S,o =0.05, see figures 7 and 8) and six (2:1,4:3,4:3,3:2,3:2, e J,o = 0.05, e S,o = 0.025, see figures 15 and 16) planetconfigurations. We conduct 200 simulations for each setof initial conditions where we utilize M KB = 40.0 M ⊕ ,and 200 runs that study 10.0 M ⊕ disks. These simu-lations are performed in the exact same manner (i.e.:simulation time, time-step, etc.) as our initial simula-tions described in section 3. Table 5 summarizes theresults of these additional simulations, and the percent-ages of systems that satisfy each of our success criteria, A - D As we speculated throughout the previous sections,heavier disk masses significantly boost the probabilitythat five planet systems retain four planets (58% versus31%) with the proper semi-major axes (meet criterion B ; 31% versus 20%). Figure 18 depicts how the finalorbits for these criteria A and B satisfying systems pro-vide a much better match to the actual solar systemthan our corresponding M KB = 20.0 M ⊕ systems (fig-ure 8). In contrast, a lighter, 10.0 M ⊕ primordial Kuiperbelt produced systematically poor solar system analogs.Specifically, none of these systems satisfied criterion C or D .While the success rates for each individual criterion are improved in our simulations considering a heavierdisk (table 5), no system satisfies all four criteria simul-taneously. A major reason for this is that interactionswith the remnant planetesimal disk cause Saturn to mi-grate as well. When an instability is mild enough to re-tain four giant planets, it also leaves behind enough massin the Kuiper belt to drive Saturn past P S /P J = 2.5.15% of all runs in our original, M KB = 20.0 M ⊕ batchsatisfied criterion A and D simultaneously ( N GP = 4and P S /P J < M KB = 40.0 M ⊕ batch finish with Jupiterand Saturn inside of their mutual 5:2 resonance, suchsuccessful outcomes occur preferentially in systems thatexperience violent instabilities, lose too many planets,eject the majority of the primordial Kuiper belt objects,and leave behind a planetesimal disk with a total massthat is insufficient to drive any appreciable residual mi-gration for Saturn. Therefore, the vast majority of sys-tems that satisfy both A and B also fail D . Indeed, onlyone of our simulations considering a heavier planetesimaldisk meets these three constraints simultaneously. In-terestingly, the same phenomenon is responsible for theincreased success of our six planet configurations thatutilize M KB = 10.0 M ⊕ . While these systems achievelower rates of success for criteria A - C than the M KB =40.0 M ⊕ set, the lighter disk tends to curtail Saturn’sresidual migration. This manifests as more systems sat-isfying criterion D , and a larger percentage of simula-tions (3%) satisfying all four success criteria simultane-ously. Though the individual rates of success are worsefor most of our criteria in our M KB = 10.0 M ⊕ sim-ulations, the mutual exclusivities between the variousconstraints are lessened. Specifically, only three of ournominal runs ( M KB = 20.0 M ⊕ ) satisfy the subsets of( A , C , D ) or ( A , B , D ) simultaneously. However, twiceas many of our light disk systems are successful in thismanner; thus leading to a greater number of systemssatisfying all constraints simultaneously.In summary, within our tested parameter space, weare unable to identify an “ideal” value of M KB for ourpreferred five planet configuration. Specifically, the se-lection of initial disk mass represents a compromise be-tween boosting the probability of matching Uranus andNeptune’s modern semi-major axes with a heavier pri-mordial belt, and adequately curtailing Saturn’s residualmigration with a lower initial value of M KB . In contrast,our additional simulations indicate that a lighter primor-dial Kuiper belt ( (cid:39) M ⊕ ) might be advantageous forthe six planet case.4.7. Varying the innermost ice giants’ masses N pln Resonant Chain e J,o e S,o M KB ( M ⊕ ) A B C D ALL
Table 5 . Summary of results for our simulations that vary the initial planetesimal disk mass. The columns are as follows: (1) theinitial number of giant planets, (2) the resonant chain, beginning with the Jupiter-Saturn resonance, (3-4) the initial eccentricitiesof Jupiter and Saturn, (5) the mass of the external planetesimal disk, (6) the percentage of systems satisfying criterion A ( N GP =4), (7) criterion B (the planets’ final semi-major axes within 20% of the real ones), (8) criterion C ( | ∆ M ij /M ij,ss | < i, j =5,6), M > M ), (9) criterion D ( P S /P J < E cce n t r i c it y Semi Major Axis (au) I n c li n a ti on ( ◦ ) E cce n t r i c it y Semi Major Axis (au) I n c li n a ti on ( ◦ ) Figure 18 . Left Panel:
Reproduction of figure 8 for comparison.
Right Panel:
Final orbits for systems that satisfied criteria A and B in our batch of simulations beginning with five planets in a 2:1,4:3,3:2,3:2 resonant chain, e J,o = e S,o =0.05 and a diskwith M KB = 40.0 M ⊕ (the same initial conditions as in the left panel except for M KB being doubled). The top panel depicts a/e space and the bottom panel plots a/i space. The respective planets and their simulated analogs are color coded as follows:Jupiter in red, Saturn in gold, Uranus in light blue, Neptune in dark blue. The error bars indicate one standard deviation. Thestars correspond to solar system values. While multiple previous authors have investigated theeffects of varying M KB (e.g.: Nesvorn´y & Morbidelli2012; Batygin et al. 2012), similar attempts to fine-tune instability statistics by altering the ejected planet’smass ( M IG ) are conspicuously absent from the litera-ture. Nesvorn´y (2011) mentioned experimenting withhigher masses, but eventually concluded that the re- sulting stronger scattering events systematically lead toexcessively violent instabilities. Subsequent work byNesvorn´y & Morbidelli (2012) speculated that it mightbe possible to calibrate instability results by adjusting M IG , and favored a lower mass for the ejected planet( M IG =8.0 M ⊕ ) in 2:1 cases because of the improvedrates of success for criterion D . Indeed, the selection of7 N pln Resonant Chain e J,o e S,o M IG ( M ⊕ ) A B C D ALL
Table 6 . Summary of results for runs that vary the inner ice giant masses. The columns are as follows: (1) the initial number ofgiant planets, (2) the resonant chain, beginning with the Jupiter-Saturn resonance, (3-4) the initial eccentricities of Jupiter andSaturn, (5) the mass of the innermost ice giant(s), (6) the percentage of systems satisfying criterion A ( N GP =4), (7) criterion B (the planets’ final semi-major axes within 20% of the real ones), (8) criterion C ( | ∆ M ij /M ij,ss | < i, j =5,6), M > M ),(9) criterion D ( P S /P J < a Neptune-mass additional planet throughout the ma-jority of the contemporary literature is largely based offthe finding that an encounter with a M IG = 15.0 M ⊕ planet is required to consistently excite Jupiter’s pri-mordially circular orbit (Morbidelli et al. 2009a). Asour scenario considers non-zero initial eccentricities forthe gas giants, alternative values of M IG might prove tobe advantageous. To test this hypothesis, we performan additional array of 800 simulations that vary the in-nermost ice giants’ masses. In the same manner as ourstudy of different M KB values, we take our most suc-cessful five and six planet configurations as a startingpoint (table 6). For five planet architectures, where ournominal simulations evaluate M IG = 16.0 M ⊕ , we per-form 200 simulations that test a mass of 8.0 M ⊕ , andan additional 200 that utilize M IG = 24.0 M ⊕ . In thesix planet case (nominal mass of 8.0 M ⊕ ), we performsimulations testing M IG = 6.0 and 12.0 M ⊕ .The results of this additional batch of integrations aresummarized in table 6. Our new simulations clearly in-dicate that lighter-mass ejected ice giants improve in-stability statistics, while higher values of M IG lead toreduced success rates for criterion A , B and D . Broadlyspeaking, a more massive planet generates a strongerscattering event and typical final M values that arecorrespondingly higher. Indeed, for the five planet case,68% of our systems that consider a more massive addi-tional planet finish with M > M IG = 16.0 and 8.0 M ⊕ , respectively. Oursix planet configurations yield similar results. Specifi-cally, 73%, 47% and 35% of these systems attain M > M IG = 12.0, 8.0 and 6.0 M ⊕ , respectively.While the instabilities considering more massive planetsare often increasingly violent, and therefore have greaterpotential to adequately excite the secular modes of theJupiter-Saturn system, this advantage is outweighed bythese systems’ systematically lower success rates for cri-terion A , B and D (though several systems still manageto satisfy all four criteria simultaneously; table 6). As pointed out in the previous discussion, perhaps themost challenging set of constraints for our primordiallyexcited resonant chains to meet simultaneously are A and D (retain four planets and not exceed P S /P J =2.5). While only a few percent of the majority of oursimulation sets are successful in this regard, 52% of ourfive planet configurations that utilize M IG = 8.0 M ⊕ meet both metrics simultaneously. Thus, by greatlyincreasing the likelihood of Jupiter and Saturn main-taining semi-major axes interior to their mutual 5:2 res-onance, our systems testing the lowest values of M IG yield markedly improved rates of satisfying criteria A - D concurrently. Therefore, future study of the 2:1Jupiter-Saturn resonance with primordial eccentricitypumping should investigate ejected ice giants with lowermasses than those considered in the bulk of our presentmanuscript. DISCUSSIONThe main conclusion of our study is not to disprovethe viability of the primordial 3:2 Jupiter-Saturn reso-nance, or the work of previous studies. Indeed, instabil-ities’ originating from the 3:2 still yield reasonable suc-cess rates (even when scrutinized against our updated,stricter constraints: table 2). However, we have shownthat it is systematically challenging to limit Saturn’ssemi-major axis jump, and post-instability eccentricitywhen the two gas giants are born in closer proximity toone another. This problem is obviously lessoned whenJupiter and Saturn emerge from the nebular disk phasein the wider, 2:1 MMR. If the two planets also acquire el-evated eccentricities during the gas disk phase by virtueof having carved larger gaps in the disk’s radial profile(Pierens et al. 2014), systematic challenges with excitingthe giant planet’s secular modes in 2:1 Jupiter-Saturninstabilities (e.g.: Nesvorn´y & Morbidelli 2012; Deiennoet al. 2017) might also be alleviated.It is worthwhile to acknowledge that, in spite of allefforts made, certain aspects of giant planets’ modernarchitecture remain low-probability outcomes in our fa-vored scenarios. While our work uncovers a potential8scheme for replicating the Jupiter-Saturn system’s secu-lar architecture with a moderately improved likelihood,such successful systems often fail to retain two ice gi-ants and adequately replicate their modern orbits. Forthis reason, many of the Jupiter analogs in the N GP =4 systems plotted in figures 8, 16 and 18 tend to pos-sess M values closer to our lower limit of 0.022 thanthe modern value of 0.044. While it is indeed impor-tant to avoid trading one low-probability outcome (e.g.:matching M ij for i, j = 5,6, and P S /P J in the same sys-tem with the primordial 3:2 Jupiter-Saturn resonance)for another, we argue that a successful scenario for pro-ducing the gas giants’ orbital architecture is perhapsmore compelling for a number of reasons. Most signif-icantly, in simulations where M is matched, our 2:1cases tend to yield significantly improved orbits for Sat-urn (namely lower values e S and P S /P J ) than our best3:2 sets. Conversely, evolutions beginning from the pri-mordial 3:2 Jupiter-Saturn resonance often possess over-excited Saturn analogs when Jupiter’s eccentricity is ad-equately excited. This is evidenced by low success ratesfor criterion C for our 3:2 simulations, and the distribu-tions of eccentric magnitudes for our 2:1 batches plottedin figures 7 and 15. Additionally, we have shown thatsuccessful replication of the ice giant system can be fine-tuned by varying M IG and M KB . The ice giant’s finalsemi-major axes (a problem in some of our sets of initialconditions) are sensitive to the range of migration thatoccurs prior to the instability. As our study utilizes anartificial instability trigger to minimize computing time,we do not evaluate the full spectrum of possible primor-dial migration paths. Furthermore, tighter resonancesbetween the innermost ice giant and Saturn (e.g. 5:4 or6:5) might provide additional fidelity in favorably mod-ifying instability statistics. Finally, the specific numberof primordial ice giants is also a free parameter. Wefind that our most successful outcomes occur in config-urations that maximize the dynamical spacing betweenthe eventual Uranus and Neptune analogs and the chaosthat ensues when Jupiter and Saturn interact with theadditional ice giants. Therefore, it is possible that a con-figuration that includes a larger number of less massiveadditional planets in a compact resonant chain locatedbetween the gas giants and modern ice giants might alsobe successful.As discussed in section 4.5, our work is limited to thestudy of chains with initial gas giant eccentricities (cid:46) P S /P J (cid:39) g = g and g = g secular resonances would be encountered. Thismay not be particularly consequential if the instabilityoccurred early (appendix A: Izidoro et al. 2016; Mor-bidelli et al. 2018; Nesvorn´y et al. 2018; Clement et al.2018; Deienno et al. 2018). However, follow-on stud-ies must throughly investigate the consequences of ourproposed scenario on the inner solar system.Finally, it is clear that our proposed scenario repre-sents somewhat of a paradigm shift; particularly giventhat the original motivation for the Nice Model’s de-velopment was to provide a mechanism for exciting thegiant planets’ eccentricities. If the giant planets wereindeed “born eccentric,” it would be logical to ques-tion whether an epoch of instability is necessary to ex-plain the solar system’s dynamical state. However, westrongly assert that our work should not be interpretedin such a manner. While we have shown that plane-tary encounters within the instability might not be thesole source of the giant planets’ modern eccentricities,the prevalence of irregular moons in the outer solar sys-tem and the asteroid belt’s dynamical structure (amongother aspects: Nesvorn´y 2018) strongly conflict with al-ternative giant planet migration scenarios (i.e. smoothmigration: Morbidelli et al. 2009b; Walsh & Morbidelli2011). CONCLUSIONSWe presented a statistical analysis of more than 6,000simulations of the Nice Model instability. Our work in-vestigated the possibility that Jupiter and Saturn in-habited a 2:1 MMR with inflated eccentricities (Pierenset al. 2014) around the time of nebular gas dispersal.Our investigation of the primordial 2:1 Jupiter-Saturnresonance was partially motivated by a detailed analy-sis of issues with the 3:2 version of the Nice Model typ-ically invoked in the literature (e.g.: Batygin & Brown2010; Nesvorn´y 2011; Batygin et al. 2012; Nesvorn´y &Morbidelli 2012; Deienno et al. 2017; Clement et al.2018) that we argued might be systematic problems.Specifically, there is a strong anti-correlation betweenthe adequate excitation of Jupiter’s M eccentric modeand maintaining the gas giant’s period ratio less than2.5 when the planets’ begin in the more compact, 3:29MMR. We used an updated series of success metrics (ta-ble 2) that fully evaluate all the important eigenmodesof the Jupiter-Saturn secular system to show that the2:1 version of the Nice Model improves the probabil-ity of proper replication of these qualities of the solarsystem. However, a caveat of our findings is that ourfavored combinations of eccentricities for Jupiter andSaturn ( e J,o > e
S,o ) are slightly at odds with the resultsof most hydrodynamical studies (Pierens et al. 2014, al-though we find reasonable success rates for the oppositecombination). While our simulations also indicated aminor anti-correlation between acceptable outcomes forthe gas and ice giants, we used an additional suite of sim-ulations to show that certain instability statistics relatedto Uranus and Neptune can be fine-tuned by varying themasses of the primordial Kuiper belt and innermost icegiants. Though our work shows that the primordial 2:1Jupiter-Saturn resonance is a viable evolutionary pathfor the solar system, future work is still required to fullyvalidate our presumed initial conditions, and robustlyanalyze the consequences of such a scenario on the solarsystem’s fragile populations of small bodies. In particu-lar, follow-on investigation of the giant planets’ instabil-ity evolution with Jupiter and Saturn in a primordial2:1 MMR with enhanced eccentricities must considerlonger integration times ( (cid:38)
100 Myr), higher resolutiondisks ( (cid:38)
Agnor, C. B., & Lin, D. N. C. 2012, ApJ, 745, 143Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJL, 884, L41Barr, A. C., & Canup, R. M. 2010, Nature Geoscience, 3, 164Batygin, K., & Brown, M. E. 2010, ApJ, 716, 1323Batygin, K., Brown, M. E., & Betts, H. 2012, ApJL, 744, L3Batygin, K., Brown, M. E., & Fraser, W. C. 2011, ApJ, 738, 13Beaug´e, C., & Nesvorn´y, D. 2012, ApJ, 751, 119Boehnke, P., & Harrison, T. M. 2016, Proceedings of theNational Academy of Science, 113, 10802Bottke, W. F., Vokrouhlick´y, D., Broz, M., Nesvorn´y, D., &Morbidelli, A. 2001, Science, 294, 1693Brasser, R., & Lee, M. H. 2015, AJ, 150, 157Brasser, R., Matsumura, S., Muto, T., & Ida, S. 2018, ApJL,864, L8Brasser, R., Morbidelli, A., Gomes, R., Tsiganis, K., & Levison,H. F. 2009, A&A, 507, 1053Brasser, R., Werner, S., & Mojzsis, S. 2020, Icarus, 338, 113514Chambers, J. E. 1999, MNRAS, 304, 793—. 2007, Icarus, 189, 386Clement, M. S., Kaib, N. A., Raymond, S. N., Chambers, J. E.,& Walsh, K. J. 2019a, Icarus, 321, 778Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J.2018, Icarus, 311, 340Clement, M. S., Morbidelli, A., Raymond, S. N., & Kaib, N. A.2020, MNRAS, 492, L56 Clement, M. S., Raymond, S. N., & Kaib, N. A. 2019b, AJ, 157,38D’Angelo, G., & Marzari, F. 2012, ApJ, 757, 50Deienno, R., Gomes, R. S., Walsh, K. J., Morbidelli, A., &Nesvorn´y, D. 2016, Icarus, 272, 114Deienno, R., Izidoro, A., Morbidelli, A., et al. 2018, ApJ, 864, 50Deienno, R., Morbidelli, A., Gomes, R. S., & Nesvorn´y, D. 2017,AJ, 153, 153Deienno, R., Nesvorn´y, D., Vokrouhlick´y, D., & Yokoyama, T.2014, AJ, 148, 25Delbo, M., Avdellidou, C., & Morbidelli, A. 2019, A&A, 624, A69Delbo’, M., Walsh, K., Bolin, B., Avdellidou, C., & Morbidelli,A. 2017, Science, 357, 1026Elkins-Tanton, L. T., Burgess, S., & Yin, Q.-Z. 2011, Earth andPlanetary Science Letters, 304, 326Fernandez, J. A., & Ip, W.-H. 1984, Icarus, 58, 109Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005,Nature, 435, 466Gomes, R., Nesvorn´y, D., Morbidelli, A., Deienno, R., &Nogueira, E. 2018, Icarus, 306, 319Gomes, R. S., Morbidelli, A., & Levison, H. F. 2004, Icarus, 170,492Grange, M. L., Nemchin, A. A., & Pidgeon, R. T. 2013, Journalof Geophysical Research (Planets), 118, 2180Hahn, J. M., & Malhotra, R. 1999, AJ, 117, 3041—. 2005, AJ, 130, 2392 Henrard, J., & Lemaitre, A. 1983, Celestial Mechanics, 30, 197Izidoro, A., Morbidelli, A., Raymond, S. N., Hersant, F., &Pierens, A. 2015, A&A, 582, A99Izidoro, A., Raymond, S. N., Pierens, A., et al. 2016, ApJ, 833,40Kaib, N. A., & Chambers, J. E. 2016, MNRAS, 455, 3561Kaib, N. A., & Sheppard, S. S. 2016, AJ, 152, 133Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ,861, 140Laskar, J. 1990, Icarus, 88, 266Laskar, J., Fienga, A., Gastineau, M., & Manche, H. 2011, A&A,532, A89Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature,524, 322Levison, H. F., Morbidelli, A., Tsiganis, K., Nesvorn´y, D., &Gomes, R. 2011, AJ, 142, 152Levison, H. F., Morbidelli, A., Van Laerhoven, C., Gomes, R., &Tsiganis, K. 2008, Icarus, 196, 258Liu, D., Jolliff, B. L., Zeigler, R. A., et al. 2012, Earth andPlanetary Science Letters, 319, 277Lykawka, P. S., & Ito, T. 2013, ApJ, 773, 65Malhotra, R. 1993, Nature, 365, 819—. 1995, AJ, 110, 420Masset, F., & Snellgrove, M. 2001, MNRAS, 320, L55Matsumura, S., Brasser, R., & Ida, S. 2017, A&A, 607, A67Matsumura, S., Thommes, E. W., Chatterjee, S., & Rasio, F. A.2010, ApJ, 714, 194McNally, C. P., Nelson, R. P., Paardekooper, S.-J., &Ben´ıtez-Llambay, P. 2019, MNRAS, 484, 728Mercer, C. M., Young, K. E., Weirich, J. R., et al. 2015, ScienceAdvances, 1, e1400050Merle, R. E., Nemchin, A. A., Grange, M. L., Whitehouse, M. J.,& Pidgeon, R. T. 2014, Meteoritics and Planetary Science, 49,2241Michtchenko, T. A., & Malhotra, R. 2004, Icarus, 168, 237Minton, D. A., & Malhotra, R. 2011, ApJ, 732, 53Mojzsis, S. J., Brasser, R., Kelly, N. M., Abramov, O., &Werner, S. C. 2019, ApJ, 881, 44Morbidelli, A. 2002, Modern celestial mechanics : aspects of solarsystem dynamicsMorbidelli, A., Brasser, R., Gomes, R., Levison, H. F., &Tsiganis, K. 2010, AJ, 140, 1391Morbidelli, A., Brasser, R., Tsiganis, K., Gomes, R., & Levison,H. F. 2009a, A&A, 507, 1041—. 2009b, A&A, 507, 1041Morbidelli, A., & Crida, A. 2007, Icarus, 191, 158Morbidelli, A., Levison, H. F., Tsiganis, K., & Gomes, R. 2005,Nature, 435, 462Morbidelli, A., Nesvorny, D., Laurenz, V., et al. 2018, Icarus,305, 262Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F., & Gomes,R. 2007, AJ, 134, 1790Murray, C. D., & Dermott, S. F. 1999, Solar system dynamicsNesvorn´y, D. 2011, ApJL, 742, L22—. 2015a, AJ, 150, 73—. 2015b, AJ, 150, 68 —. 2018, ARA&A, 56, 137Nesvorn´y, D., & Morbidelli, A. 2012, AJ, 144, 117Nesvorn´y, D., & Vokrouhlick´y, D. 2016, ApJ, 825, 94Nesvorn´y, D., Vokrouhlick´y, D., Bottke, W. F., & Levison, H. F.2018, Nature Astronomy, 2, 878Nesvorn´y, D., Vokrouhlick´y, D., & Deienno, R. 2014a, ApJ, 784,22Nesvorn´y, D., Vokrouhlick´y, D., Deienno, R., & Walsh, K. J.2014b, AJ, 148, 52Nesvorn´y, D., Vokrouhlick´y, D., & Morbidelli, A. 2007, AJ, 133,1962—. 2013, ApJ, 768, 45Nobili, A. M., Milani, A., & Carpino, M. 1989, A&A, 210, 313Norman, M. D., Duncan, R. A., & Huard, J. J. 2006, GeoCoA,70, 6032Novakovi´c, B., Cellino, A., & Kneˇzevi´c, Z. 2011, Icarus, 216, 69Petrovich, C., Malhotra, R., & Tremaine, S. 2013, ApJ, 770, 24Pierens, A., & Nelson, R. P. 2008, A&A, 482, 333Pierens, A., & Raymond, S. N. 2011, A&A, 533, A131Pierens, A., Raymond, S. N., Nesvorny, D., & Morbidelli, A.2014, ApJL, 795, L11Poincare, H. 1892, Les methodes nouvelles de la mecaniquecelesteQuarles, B., & Kaib, N. 2019, AJ, 157, 67Raymond, S. N., Izidoro, A., & Morbidelli, A. 2018, arXive-prints, arXiv:1812.01033Raymond, S. N., O’Brien, D. P., Morbidelli, A., & Kaib, N. A.2009, Icarus, 203, 644Rein, H., Brown, G., & Tamayo, D. 2019, MNRAS, 490, 5122Rein, H., Payne, M. J., Veras, D., & Ford, E. B. 2012, MNRAS,426, 187Ribeiro, R. d. S., Morbidelli, A., Raymond, S. N., et al. 2020,Icarus, 339, 113605Rivera, E. J., Laughlin, G., Butler, R. P., et al. 2010, ApJ, 719,890Roig, F., & Nesvorn´y, D. 2015, AJ, 150, 186Roig, F., Nesvorn´y, D., & DeSouza, S. R. 2016, ApJL, 820, L30Spudis, P. D., Wilhelms, D. E., & Robinson, M. S. 2011, Journalof Geophysical Research (Planets), 116, E00H03Stoer, J., Bartels, R., Gautschi, W., Bulirsch, R., & Witzgall, C.2002, Introduction to Numerical Analysis, Texts in AppliedMathematics (Springer New York)Tera, F., Papanastassiou, D. A., & Wasserburg, G. J. 1974,Earth and Planetary Science Letters, 22, 1Thommes, E. W., Duncan, M. J., & Levison, H. F. 1999, Nature,402, 635Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005,Nature, 435, 459ˇSidlichovsk´y, M., & Nesvorn´y, D. 1996, Celestial Mechanics andDynamical Astronomy, 65, 137Vokrouhlick´y, D., & Nesvorn´y, D. 2015, ApJ, 806, 143Volk, K., & Malhotra, R. 2019, AJ, 158, 64Walsh, K. J., & Morbidelli, A. 2011, A&A, 526, A126Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67Zellner, N. E. B. 2017, Origins of Life and Evolution of theBiosphere, 47, 261Zhang, H., & Zhou, J.-L. 2010, ApJ, 719, 671
APPENDIX A. THE TIMING OF THE INSTABILITYAs we neglect pre-instability migration, our work investigates an instability scenario that occurs rather quickly afternebular gas dissipation. Therefore, it is worthwhile to comment on the various proposed timings for the instability.In its original formulation (Gomes et al. 2005), the Nice Model instability was proposed to have occurred coincident1with the Late Heavy Bombardment (LHB; Tera et al. 1974). Thus, a late instability implies a rather specific timingfor the event; (cid:39)
650 Myr after gas disk dispersal in order to provide a trigger for the cratering spike in the innersolar system as part of a LHB scenario. However, it is worth noting that the timescale of giant planet migration issignificantly shorter than 650 Myr for many sets of initial disk conditions. Specifically, maintaining a system of giantplanets stable in resonance for hundreds of Myr requires larger radial spacings between Neptune and the primordialKuiper belt (Levison et al. 2011; Nesvorn´y & Morbidelli 2012; Deienno et al. 2017). Furthermore, average instabilitydelays in simulations that account for Kuiper belt self-gravity (Quarles & Kaib 2019) are typically less than ∼
100 Myrregardless of initial offset from Neptune. In a similar manner, Ribeiro et al. (2020) found late instabilities to be almostnon-existent in giant planet formation models that successfully generate the large obliquities of Uranus and Neptunevia embryo-embryo impacts (Izidoro et al. 2015).While late instabilities do occur in numerical simulations, albeit less frequently than early destabilizations, authors inrecent years have begun invoking earlier versions of the Nice Model for a variety of reasons. Perhaps most significantly,the very existence of a cratering spike has been called into question as a result of updated isotopic dating methodsused in the analysis of Lunar samples. The refined basin ages (many utilizing Ar/ Ar dataing; e.g.: Boehnke &Harrison 2016) reported in contemporary studies seem to imply a smoother decline in bombardment, rather than aterminal cataclysm (Zellner 2017). These include, for example, ∼ ∼ ∼ ∼ (cid:38) B. SECULAR EVOLUTION OF THE SOLAR SYSTEMThe secular dynamics of a system can be studied when the planets are sufficiently far from MMR by expanding themutual interaction components to order one in mass (Michtchenko & Malhotra 2004): H = H + (cid:15) H (B1)Here, H is the integrable approximation, H is the secular normal form, and (cid:15) represents the order of the perturbation.The standard procedure involves writing the secular Hamiltonian in terms of the modified Delaunay variables, assumingthe planets’ eccentricities and inclinations are small, and expanding in Taylor series. For complete derivations, we directthe reader to chapter 7 of Murray & Dermott (1999) or chapter 7 of Morbidelli (2002). The so-called Lagrange-Laplace2 g g g g Table B1 . Secular precession frequencies ( (cid:48)(cid:48) /yr ) as calculated in Laskar et al. (2011). j \ i Table B2 . Amplitudes, M ij (equation B2), for the dominant eigenmodes g - g (rows) in the orbit in each of the giant planets(columns). All values are reproduced from the work of (Nobili et al. 1989). solution is as follows: e i cos (cid:36) i = (cid:80) j M ij cos ( g j t + β j ) e i sin (cid:36) i = (cid:80) j M ij sin ( g j t + β j )sin i i cos Ω i = (cid:80) j M ij cos ( s j t + β j )sin i i sin Ω i = (cid:80) j M ij sin ( s j t + β j ) (B2)The fundamental frequencies g j and s j ( i = 1-8, with s necessarily set to zero by convention) represent the dominanteigenfrequencies for the precession of perihelia and longitudes of nodes, respectively. The solar system’s frequencies,as well as their respective amplitudes in each planets’ orbit, are typically calculated via Fourier analysis of the outputsof high-accuracy numerical simulations (e.g.: Nobili et al. 1989; Laskar 1990; ˇSidlichovsk´y & Nesvorn´y 1996; Laskaret al. 2011; Rein et al. 2019). Tables B1 and B2 summarize the current values of the dominant eigenfrequencies, g j ,and amplitudes, M ij , in the outer solar system.The precise numerical values of the secular precession frequencies in the outer solar system ( g - g ) are largely afunction of the orbital spacing between the respective planets (Murray & Dermott 1999; Nesvorn´y & Morbidelli 2012).However, the various planets’ eccentric magnitudes, M ij , were acquired via dynamical processes (Brasser et al. 2009;Morbidelli et al. 2009a). Therefore, the dominant secular eigenfrequencies, along with their respective amplitudeswithin the relevant planets’ eccentricities and inclinations must be broadly reproduced in any successful evolutionarymodel . Given the hierarchical nature of the solar system’s mass distribution, one can simplify the equations in B2by studying only the Jupiter-Saturn secular problem for a simple approximation of the solar system’s global seculardynamics (Morbidelli et al. 2009a). Thus, the eccentricity evolution of the two gas giants can roughly be describedby: e J cos (cid:36) J = M cos ( g t + β ) − M cos ( g t + β ) e J sin (cid:36) S = M sin ( g t + β ) − M sin ( g t + β ) e S cos (cid:36) J = M cos ( g t + β ) + M cos ( g t + β ) e S sin (cid:36) S = M sin ( g t + β ) + M sin ( g t + β ) (B3)In this Jupiter-Saturn approximation, the time evolution of each planet’s eccentricity is given by: e J ( t ) = (cid:113) M + M − M M cos (cid:0) ( g − g ) t − β (cid:1) e S ( t ) = (cid:113) M + M + 2 M M cos (cid:0) ( g − g ) t − β (cid:1) (B4)where β = β − β . Thus, noting that M is negative, Jupiter’s eccentricity oscillates between an approximateminimum of M + M (cid:39) M - M (cid:39) M - M = 0.015 and M + M = 0.081, and the characteristic period of these oscillations is ∝ | g − g | − (cid:39) g ; see table B2),this figure illustrates how the secular structure of the solar system is reasonably approximated by the 2 frequenciesand 4 amplitudes of the Jupiter-Saturn system. This is important because Jupiter is the most significant eccentricperturber for all the planets except Mars and Neptune (as quantified by the magnitude of g in each planets’ orbit).3 e J M + M M - M Time (Myr) e S M - M M + M Figure B1 . Eccentricity variations in the orbits of Jupiter and Saturn from a 10 Myr integration of the Sun and 8 planetswith the
Mercury g and g eigenfrequency in each planets orbit: M , M , M , and M . In addition to driving variations in the planets’ eccentricities, the dominant g and g frequencies provide a good firstorder approximation of Jupiter and Saturn’s eccentricity vector precessions (figure B2). C. ALTERNATIVE AVENUES FOR ADEQUATELY EXCITING M C.1.
The effect of approaching the 5:2 MMR
We also explored the possibility that Jupiter’s eccentricity can be further excited after the instability occurs. Indeed,the planets orbits evolve chaotically over Gyr timescales (though this departure from quasi-periodic behavior is mostlyconfined to the orbits of the terrestrial worlds, e.g.: Laskar 1990). It is known that the giant planets continue tomigrate smoothly over some radial range after the instability (Deienno et al. 2017) as they continue to interact withand clear out debris in the vicinity of their orbits. Through this process, the eccentricities and inclinations of Uranusand Neptune can damp rather substantially via secular friction (Nesvorn´y & Morbidelli 2012). While this mechanismtends to damp the eccentric modes of the Jupiter/Saturn system as well, the planets’ mutual secular interactions areknown to enhance in close proximity to the 5:2 MMR (Clement et al. 2020). Figure C3 quantifies the degree to which M can be further excited within this phase of approach. In general, as Jupiter and Saturn’s mutual interactionsweaken with increased radial separation, their eccentric forcing on one another lessens correspondingly. Thus, themagnitude of g in Jupiter’s eccentricity lowers, and so does g in Saturn’s. This results in M rising slightly relativeto M , and M increasing minimally with respect to M . However, the net total change in M is only a few percentover the range of mutual separations depicted in figure C3. We conclude that M likely cannot evolve appreciablyafter the instability via residual migration towards the 5:2 MMR (the subsequent section addresses the possibility ofa resonant crossing, see also Morbidelli et al. 2009b).4 J
300 400 500 600 700 800
Time (Kyr) S Figure B2 . Precession of Jupiter and Saturn’s orbital longitude of perihelia ( (cid:36) ) in a numerical integration of all 8 planets. Thevertical lines are spaced in equal increments of g (solid grey) and g ’s (dashed grey) periods. P Sat /P Jup E cce n t r i c it y M M M M Figure C3 . Evolution of the M ij (i,j=5,6) amplitudes as Jupiter and Saturn period ratio is adjusted. This figure is generatedby performing 3,200 integrations (see Clement et al. 2020) of the modern solar system with the Mercury P S /P J varies from 2.39-2.49, andall other orbital elements are left unchanged. Each system is integrated for 10 Myr, and the secular amplitudes and frequenciesare calculated via Fourier analysis of the simulation time outputs (ˇSidlichovsk´y & Nesvorn´y 1996). The figure is limited tothe range of 2.39 < P S /P J < A late dynamical event
To search for potential avenues for the post-instability excitation of M , we integrate 371 system formed in Clementet al. (2018) for an additional 1 Gyr. These systems, born out of the 3:2 mutual Jupiter-Saturn resonance (table 1),finished with a range of final values for P S /P J and M (e.g.: figure 1). It is important to note that many of thesesystems differ from perfect solar system analogs for a variety of reasons . However, we intentionally study a wide rangeof systems that are similar to the solar system in order to conduct a broad exploration of possible additional excitationmechanisms for M . After the complete inner and outer solar systems are integrated for 1 Gyr with the M ercury ∼
165 Myr) that results in the ejection of one of the ice giants (i.e.: a sixgiant planet system that only lost one planet during the first instability ejects an another ice giant; thus becoming afour giant planet system). When this is the case, Jupiter and Saturn experience a corresponding “jump” in semi-majoraxis, leading to a large absolute change in M . We also observe systems that lose an additional terrestrial planetduring this phase of evolution. Typically, this occurs when the eccentricity an additional ∼ Mars-mass planet thatformed between Mars and the asteroid belt (e.g.: Chambers 2007) is excited to the point where it collides with oneof the other planets or is scattered into the Sun. We find that, when a terrestrial world is lost by either ejection orcollision (often with an Earth or Venus analog), M does not change appreciably. In contrast, the average change in M is 81% when an additional ice giant is ejected late.Large absolute decreases (by at least 50%) in M are about twice as likely as increases (20 versus 12 total systems,respectively) in our systems that do not experience a collision or planetary ejection. Thus, our simulations indicatethat it is substantially more difficult to further excite M after the instability than it is to damp Jupiter’s eccentricity.Moreover, the majority of the larger relative changes in M (up or down in magnitude) are in systems with lowerinitial M magnitudes in the M > M regime. Indeed, the average change in M for any system beginning with M > M > M over 1 Gyr, without ejecting a giant planet. On closer inspection, both of these systemsbegin with Jupiter and Saturn just outside of a MMR (5:2 and 3:1). At some point in the simulation, the gas giantsfall into resonance with each other, leading to chaotic and irregular behavior of the secular eigenmodes.We find similar results for the other amplitudes of the Jupiter-Saturn secular system (equation B4). The magnitudes M and M only change by an average of 18% and 9%, respectively, in systems without a late ice giant ejection.Moreover, large magnitude decreases are around three times more likely than increases for these two modes. However,after carefully analyzing each simulation that experiences a large change in one of the eccentric magnitudes, M ij , wefind that either a late dynamical instability or a giant planet MMR crossing occurred. Thus, we cannot discountthe possibility that Jupiter and Saturn’s eccentricities were under-excited during the instability, and subsequentlyfurther excited up to their modern values by a late dynamical event. However, such a scenario would imply that thefragile dynamical structures in the solar system (e.g.: the terrestrial system, asteroid belt and binary trojans: Kaib &Chambers 2016; Delbo’ et al. 2017; Nesvorn´y et al. 2018) survive an additional violent episode ( >
200 Myr after gasdisk dissipation in our simulations). For instance, the consequences of Jupiter and Saturn having temporarily enteredthe 5:2 MMR are unclear, as is whether or not they could have been perturbed back out of the resonance. MoreoverMojzsis et al. (2019) found that a dynamical instability occurring after ∼∼