THE ASTEROID BELT AS A RELIC FROM A CHAOTIC EARLY SOLAR SYSTEM André Izidoro1,2, Sean N. Raymond1, Arnaud Pierens1, Alessandro Morbidelli3, Othon C. Winter4, and David Nesvornỳ5 1 Laboratoire d’astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France; izidoro.costa@gmail.com 2 Capes Foundation, Ministry of Education of Brazil, Brasília/DF 70040-020, Brazil 3 University of Nice-Sophia Antipolis, CNRS, Observatoire de la Côte d’Azur, Laboratoire Lagrange, BP 4229, F-06304 Nice Cedex 4, France 4 UNESP, Univ. Estadual Paulista—Grupo de Dinâmica Orbital & Planetologia, Guaratinguetá, CEP 12.516-410, São Paulo, Brazil 5 Department of Space Studies, Southwest Research Institute, 1050 Walnut St., Suite 300, Boulder, CO 80302, USA Received 2016 July 12; revised 2016 September 16; accepted 2016 September 19; published 2016 December 6 ABSTRACT The orbital structure of the asteroid belt holds a record of the solar system’s dynamical history. The current belt only contains ∼10−3 Earth masses yet the asteroids’ orbits are dynamically excited, with a large spread in eccentricity and inclination. In the context of models of terrestrial planet formation, the belt may have been excited by Jupiter’s orbital migration. The terrestrial planets can also be reproduced without invoking a migrating Jupiter; however, as it requires a severe mass deficit beyond Earth’s orbit, this model systematically under-excites the asteroid belt. Here we show that the orbits of the asteroids may have been excited to their current state if Jupiter’s and Saturn’s early orbits were chaotic. Stochastic variations in the gas giants’ orbits cause resonances to continually jump across the main belt and excite the asteroids’ orbits on a timescale of tens of millions of years. While hydrodynamical simulations show that the gas giants were likely in mean motion resonance at the end of the gaseous disk phase, small perturbations could have driven them into a chaotic but stable state. The gas giants’ current orbits were achieved later, during an instability in the outer solar system. Although it is well known that the present-day solar system exhibits chaotic behavior, our results suggest that the early solar system may also have been chaotic. Key words: chaos – minor planets, asteroids: general – planets and satellites: gaseous planets – planets and satellites: terrestrial planets 1. INTRODUCTION The distribution of asteroids strongly constrains planet formation models. While the terrestrial planets’ orbits are nearly circular and coplanar, the orbital eccentricities of asteroids are excited, filling parameter space from e=0 to 0.3, and inclination i=0° to 20°. The asteroid belt’s total mass is also only ~ -10 3 Earth masses. There are two basic views on how the inner solar system was built, with different implications for the asteroid belt. In one view, the asteroid belt contained a few Earth masses in solid material but was rapidly depleted and excited by dynamical mechanisms. Gravitational scattering of asteroids by a popula- tion of moon- to Mars-sized planetary embryos originally in the belt can promote significant depletion and excitation of the belt (Wetherill 1978, 1992, 1986; Chambers & Wetherill 1998; Agnor et al. 1999; Petit et al. 1999; Chambers 2001; Petit et al. 2001, 2002, pp. 711–723; O’Brien et al. 2007). One problem with this scenario in the context of terrestrial planet formation is that Mars analogs produced in these simulations are far more massive than the actual one (e.g., Raymond et al. 2004; O’Brien et al. 2006; Raymond et al. 2006; Morishima et al. 2008; Raymond et al. 2009; Izidoro et al. 2013; Lykawka & Ito 2013; Fischer & Ciesla 2014; Izidoro et al. 2014). The Grand Tack scenario (Walsh et al. 2011)—which invokes the inward then outward migration of Jupiter through the asteroid belt region—removes enough mass beyond 1 au to explain why Mars is much smaller than Earth and to sculpt the asteroid belt in a way that will become consistent with its current structure via subsequent dynamical evolution (Roig & Nesvorný 2015; Deienno et al. 2016). In the opposite view, the asteroid belt had low mass even at early times (Izidoro et al. 2015b; Levison et al. 2015; Moriarty & Fischer 2015; Drazkowska et al. 2016), and Jupiter and Saturn did not migrate across the asteroid belt. In this framework, a primordial low-mass asteroid belt should be far less dynamically excited than the observed one (Izidoro et al. 2015b), and what remains to be explained is the belt’s dynamical excitation (for a recent review, see Morbidelli et al. 2015, pp. 493–507). In this paper we propose a novel mechanism for explaining the dynamical excitation of the asteroid belt. The mechanism relies on the chaotic evolution of Jupiter’s and Saturn’s orbits at early times. In Section 2 we present an example for the origin of chaos in Jupiter’s and Saturn’s early orbits. In Section 3 we present our results for the chaotic excitation of the asteroid belt. In Section 4 we more fully address the possible origins of chaos in the giant planets’ orbits, and present alternative scenarios to trigger chaos. In Section 5 we discuss the implication of our results for models of solar system formation. Finally, in Section 6 we briefly summarize our findings. 2. AN EXAMPLE OF CHAOS IN JUPITER’S AND SATURN’S EARLY ORBITS The excitation of the asteroid belt took place after the gaseous protoplanetary disk had dissipated, yet it is during the disk phase that the gas giants’ orbits could have changed most dramatically due to orbital migration. Embedded in the gaseous disk, Jupiter and Saturn systematically migrate into mean motion resonance (MMR), where their orbital periods are related by a ratio of small integers (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Pierens & Nelson 2008; D’Angelo & Marzari 2012; Pierens et al. 2014). The most common are the 3:2 and 2:1 MMRs. The present-day orbits of the giant planets The Astrophysical Journal, 833:40 (18pp), 2016 December 10 doi:10.3847/1538-4357/833/1/40 © 2016. The American Astronomical Society. All rights reserved. 1 mailto:izidoro.costa@gmail.com http://dx.doi.org/10.3847/1538-4357/833/1/40 http://crossmark.crossref.org/dialog/?doi=10.3847/1538-4357/833/1/40&domain=pdf&date_stamp=2016-12-06 http://crossmark.crossref.org/dialog/?doi=10.3847/1538-4357/833/1/40&domain=pdf&date_stamp=2016-12-06 were achieved later, after the gaseous disk was gone, during a dynamical instability (Nesvorný & Morbidelli 2012). In our model, the asteroids were excited between the dissipation of the disk and the instability. In hydrodynamical simulations, Jupiter’s and Saturn’s migration typically leads to deep capture in resonance, with orbits characterized by regular motion. However, very small perturbations may push them into chaos (Sándor & Kley 2006; Batygin & Morbidelli 2013). Perturbations come from (i) dispersal of the gaseous disk and the corresponding loss of its damping (Papaloizou & Larwood 2000; Cresswell & Nelson 2008); (ii) gravitational forcing from the ice giants, both during their inward migration (Izidoro et al. 2015a) and after the dissipation of the gaseous disk; and (iii) gravitational scattering of small remnant planetary embryos and planetesimals in the giant planet region. We performed a suite of numerical experiments to show that seemingly trivial perturbations can trigger chaos in the gas giants’ orbits. Figure 1 shows one simulation in which a regularly evolving configuration of Jupiter and Saturn (the JSREG simulation), locked in 2:1 resonance, became chaotic as a result of the ejection of a Mars-sized embryo from the system. To perform this simulation we used the Mercury integrator (Chambers 1999). The system was composed by the central solar mass star, the fully formed Jupiter and Saturn in the 2:1 MMR, and a Mars-mass planetary embryo. Jupiter was initially at 5.25 au, Saturn at ∼8.33 au, and the planetary embryo at 12 au. The gas giants’ eccentricities were initially 0.025 and their mutual orbital inclination 0°.5. The planetary embryo started with zero orbital eccentricity and inclination. We used an integration time step of 100 days and assumed that the gaseous disk was already fully dissipated. While the planetary embryo in Figure 1 was only 1/4000 the combined mass of the gas giants, it triggered chaos. The gas giants’ orbits were chaotic but dynamically stable for long timescales, consistent with a late dynamical instability that re- arranged their orbits (Tsiganis et al. 2005; Morbidelli et al. 2007; Levison et al. 2011). The perturbed giant planets even remained in 2:1 resonance with modest (but chaotic) eccentricities and inclinations (see for example Figure 2). Depending on the nature and strength of the perturbation, chaos was generated in up to 100% of our simulations starting from a regular, resonant configuration. In fact, stochastic forcing from turbulent density fluctuations during the disk phase may have pushed the planets out of (deep) resonance (Adams et al. 2008; Lecoanet et al. 2009; Pierens et al. 2011), onto nearby orbits where the density of chaos is high (see the Appendix). This simulation (Figure 1) represents a simple proof of concept; a regular orbital configuration of Jupiter and Saturn can easily be converted into a chaotic one. The perturbation required must be strong enough to transition the system to a new dynamical state but not so strong as to make the system dynamically unstable. The magnitude of the perturbation in the simulation from Figure 1 is entirely plausible, as it is likely that some leftovers remained when the protoplanetary disk dissipated. Indeed, the so-called “late veneer” represents geochemical evidence that planetary leftovers remained scattered throughout the inner solar system after the Moon- forming impact on Earth, long after the dissipation of the disk (e.g., Day et al. 2007; Walker 2009; Bottke et al. 2010; Jacobson et al. 2014). In Section 4 we present additional scenarios for generating chaos in Jupiter’s and Saturn’s early orbits. 3. CHAOTIC EXCITATION OF THE ASTEROID BELT We now turn our attention to showing how the asteroids’ orbits may have been excited by chaos in the orbits of Jupiter and Saturn. We used the Symba (Duncan et al. 1998), Swift (Levison & Duncan 1994), Mercury (Chambers 1999), and Rebound (Rein & Spiegel 2014; Rein & Tamayo 2015) integrators to perform our simulations of the belt excitation. Asteroids were modeled as massless test particles. The integration time step in all our simulations was at most 1/20 of the orbital period of the innermost body in the system. We stress that we are aware of an issue with Symba5 which compromises the performance of the integrator in the case where planets have close encounters every orbital period. When this kind of evolution takes place it degrades the symplectic nature of the integrator, causing large errors. This may be the case for example when Jupiter and Saturn evolve in a compact resonant configuration (e.g., 3:2 MMR). To make sure the chaos observed in our simulations has no numeric origin we tested an ample number of integrators over the different simulations presented here. In Rebound we test both the WHFAST and IAS15 integrators and in Mercury we performed tests with the Bulirsch–Stoer and “Hybrid” integrators. Figure 2 shows the dynamical evolution of Jupiter and Saturn in the JSREG and JSCHA simulations, which we use as fiducial cases to illustrate and contrast mechanisms of excitation of the asteroid belt. These two simulations were generated from almost identical initial orbital arrangements; the Figure 1. Onset of chaos in a characteristic dynamical simulation. Jupiter and Saturn started in the JSREG configuration, locked in 2:1 MMR with low eccentricities (starting eccentricities of 0.025 for both planets) and exhibiting regular motion. The evolution of their orbital eccentricities/inclinations are shown in the top/bottom panel. Gravitational scattering of a Mars-mass embryo initially at 12au triggered chaos in the giant planets’ orbits. The planets remained in resonance, as confirmed by libration of the 2:1 MMR critical angle l l v- -2 Sat Jup Jup, on chaotic but long-term stable orbits. lSat, lJup, and vJup are Saturn’s mean longitude, Jupiter’s mean longitude and longitude of pericenter, respectively. The initial conditions of this simulation are provided in the Appendix. 2 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. only difference is that Jupiter and Saturn’s eccentricities were each set to 0.025 in the JSREG simulation (resulting in regular motion) and to 0.03 in the JSCHA simulation (which triggered chaos). In both simulations, Jupiter and Saturn’s semimajor axes are initially 5.4 and about ∼8.57 au, respectively. Their mutual initial orbital inclinations are 0°.5. Their argument of pericenter and longitude of the ascending node are set zero. Jupiter’s mean anomaly is initially zero and Saturn’s mean anomaly is initially 180°. Figure 2 shows the evolution of the gas giant’s period ratio, eccentricities, orbital inclinations, and an angle associated with the 2:1 MMR between the planets. In both simulations, the planets are in 2:1 resonance since the critical angle f l l v= - -22 Sat Jup Jup librates around zero degree, where lJup and lSat are the mean longitudes of Jupiter and Saturn, respectively. The critical angle f l l v= - -21 Sat Jup Sat circulates, wherevSat is Saturn’s longitude of pericenter. Thus, the planets are not in apsidal corotation resonance—in which the planets undergo apsidal libration as well as libration of both resonant arguments (Michtchenko et al. 2008)—as shown by the circulation of the angles v v-Sat Jup (Figure 2, bottom panels). Yet the eccentricity and inclination evolution of Jupiter and Saturn are quite different in the the two simulations. In particular, there are larger variations of eccentricity and inclination over 100Myr in the chaotic case. These variations are linked with the precession of longitudes, which translate to shifting resonances and gravitational perturbations in the belt. We now show how, once the gaseous disk was gone, chaos in Jupiter’s and Saturn’s early orbits may have excited the asteroids’ orbits even if the belt’s primordial mass was very low (comparable to its current mass). The asteroid belt is speckled with resonances, locations where there is an integer match between the characteristic orbital frequencies of asteroids and the giant planets. MMRs are located where an asteroid’s orbital period forms an integer ratio with a planet’s period (Jupiter in this case). In secular resonances (SRs), a quantity related to the precession of an asteroid’s orbit matches one of the giant planets. The most important SRs in the main belt are the n6 and n16 resonances, where an asteroid’s apsidal and nodal precession rate (or frequency), respectively, match that of Saturn (Froeschle & Scholl 1989; Morbidelli & Henrard 1991). When Jupiter’s and Saturn’s orbits evolve in a regular fashion, MMRs and SRs are stationary, so asteroids in certain parts of the belt are excited whereas asteroids in other parts are not (Morbidelli & Henrard 1991). Figure 2. Dynamical evolution of regular Jupiter and Saturn in the JSREG and JSCHA simulations. In both cases the planets are in 2:1 MMR (the critical angle f l l v= - -22 Sat Jup Jup associated with the 2:1 MMR librate around zero degree while the f l l v= - -21 Sat Jup Sat circulates). lSat and vSat are Saturn’s mean longitude and longitude of pericenter. lJup and vJup are Jupiter’s mean longitude and longitude of pericenter. The initial conditions of these two simulations are provided in the Appendix. 3 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. When the giant planets evolve chaotically, their orbital alignments undergo stochastic jumps, i.e., they may precess at many different frequencies and their orbits may even have its direction of apsidal precession temporarily reversed. The location of SRs within the belt undergo corresponding jumps (MMRs are much less sensitive to such variations). When a resonance jumps to the location of a given asteroid its orbit is significantly excited on a relatively short (~104 6– year) timescale. Figure 3 shows how the n6 and n16 SRs pump the eccentricity and orbital inclinations, respectively, of two asteroids. Note that the time intervals shown in both plots of Figure 3 were purposely chosen to clearly illustrate the effect of the respective resonances. The lifetime of a particle in the belt depends on the chaotic evolution of the gas giants and the particle’s initial orbit. We discuss and show in the Appendix that the effects of the chaotic excitation can also eject particles from the system and even empty parts of the belt. In addition to n6 and n16, other resonances also play a role in pumping asteroids’ eccentricities and/or inclinations (see the Appendix). These include resonances resulting from the linear combination of principal secular frequencies (nonlinear secular resonances), MMRs, secondary resonances (linear combination between the main secular and short period frequencies associated with orbital period), and Kozai resonances. Another factor is that, unlike MMRs, which are linked to a given orbital radius, the locations of SRs are a function of an asteroid’s proper inclination and eccentricity (Froeschle & Scholl 1989; Morbidelli & Henrard 1991). We can understand the chaotic excitation of asteroids using a Fourier analysis of Saturn’s longitude of pericenter vSat (Figure 4; top panels). When Jupiter and Saturn undergo regular motion, the power spectrum of vSat is peaked, dominated by characteristic frequencies (in particular the g6 frequency at ∼1/62000 -yr 1, which controls the location of the n6 SR) and their harmonics (Laskar 1990, 1993; Michtchenko & Ferraz-Mello 1995). The peaked nature of the power spectrum indicates that Saturn’s precession frequency—and thus the location of the n6 and other resonances—is fixed. This explains why eccentricities and inclinations of asteroids are only strongly excited at specific locations in the belt (Figure 4). In contrast, the power spectrum of a simulation in which Jupiter’s and Saturn’s orbits evolve chaotically (the JSCHA simulation) shows a broad band of frequencies instead of a few strong peaks. In this case the SR n6 jumps across the entire belt because vSat precesses at many different frequencies due to Saturn’s chaotic interactions with Jupiter. Note from the frequency analysis in Figure 4 that in the JSREG simulation the longitude of pericenter of Saturn precesses (slowly) positively while the longitude of pericenter of Jupiter precesses (quickly) backwards. If planets were in apsidal corotation, that would imply both their longitudes of pericenter would precess in the same direction (negative). Because asteroids suffering secular perturbation always precess positively, n6 could not exist in the belt if the asteroids’ and Saturn’s longitude of pericenter precess in different directions. During the gas giants’ chaotic evolution the n6 SR jumps across the entire asteroid belt. To estimate the radial extent of the jumps of n6 we first compute the precession frequency that an asteroid would have with Jupiter and Saturn in 2:1 resonance. Using linear secular theory (see Murray & Dermott 1999), this frequency is given by ⎡ ⎣ ⎢⎢ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎤ ⎦ ⎥⎥ p a a = +   A n M M a a b M M a a b 1 8 , 1 Ast Jup Ast Jup 2 3 2 1 Jup Sat Ast Sat 2 3 2 1 Sat ( ) ( ) ( ) ( ) ( ) where nAst and aAst are the mean motion and semimajor axis of the asteroid. M , MJup, and MSat are the solar mass, Jupiter’s mass and Saturn’s mass, respectively. aJup and aSat are Jupiter’s and Saturn’s semimajor axes, respectively. b3 2 1( ) is the Laplace coefficient which is computed as a function of aJup and aSat, which are given by a = a a Jup Ast Jup a = a a Sat Ast Sat If we assume that vSat precesses with frequencies between ∼10−5 and~ -10 year4 the location of n6 varies from about 1.4 to 3.5 au (Figure 5). The asteroids’ eccentricities are thus Figure 3. Chaotic excitation of the orbits of two asteroids (modeled as massless test particles) in a dynamical simulation. Left: the eccentricity of an asteroid at 3 au increases rapidly during a short interval while it is locked in the n6 SR. n6 occurs when the angle v v-Sat Ast librates, where vSat and vAst are the longitudes of pericenter of Saturn and the asteroid, respectively. Right: the inclination of an asteroid at 2.7 au is pumped during temporary capture in the n16 SR. n16 occurs when the angle W - WSat Ast librates, where WSat and WAst are the longitudes of ascending node of Saturn and the asteroid, respectively. The asteroids’ semimajor axes remain roughly constant. 4 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. excited across the entire main belt. An analogous process acts to pump the asteroids’ inclinations. In the JSCHA simulation, although n16 is much stronger and wider than in the JSREG simulation, it does not jump across the entire belt. However, other resonances contribute to exciting the asteroids’ inclina- tions across the entire belt (see the Appendix). The timescale for chaotic excitation of the full asteroid belt is a few million years to hundreds of millions of years depending on the evolution of the gas giants (additional examples are given in the Appendix). The surviving asteroids broadly match the observed distribution (Figure 4; right-hand side, middle and bottom panels). In simulations that successfully excited the asteroid belt, two conditions were typically observed. First, Jupiter’s and Saturn’s orbits were chaotic in both their eccentricities and inclinations (Barnes et al. 2015). The hardest aspect of the asteroid belt to reproduce is its broad inclination distribution (Izidoro et al. 2015b). Second, Jupiter’s eccen- tricity was not too high. In simulations in which Jupiter’s eccentricity remained much larger than its current value of ∼0.05 for longer than 100Myr, parts of the belt were emptied. Finally, although we have illustrated this mechanism with Jupiter and Saturn in 2:1 MMR, we observe chaotic excitation in a number of resonant (or near-resonant) orbits, including the 2:1, 3:2, 7:4, and 5:3 (see the Appendix). Figure 4. Fourier analysis and excitation of the asteroid belt by the gas giants in two N-body simulations. In both cases, Jupiter and Saturn are locked in 2:1 MMR, with starting semimajor axes of 5.25 au and ∼8.33 au. In the JSREG simulation (left panels) the gas giants’ initial eccentricities are 0.025 and their orbits exhibit regular motion. In the JSCHA simulation (right) the giants’ initial eccentricities are 0.03 and their orbits are chaotic. Their mutual orbital inclination is initially 0°. 5. Each system was integrated for 136 Myr, and to perform the Fourier analysis we used an output time step of 2 years; the top panels show the power spectra forvSat for each case. The middle and bottom panels show the dynamical excitation of the main belt in the two simulations, using a snapshot at 40 Myr. Some asteroids in the JSCHA simulation have larger eccentricities and inclinations than those observed. This does not mean that our results are inconsistent with the present-day asteroid belt. Highly eccentric and inclined objects are removed from the system during the later evolution of the solar system: during the dynamical instability between the giant planets (Morbidelli et al. 2010) and over the subsequent 3.8 billion years (Minton & Malhotra 2010). These simulations also do not take into account the gravitational influence of the growing terrestrial planets, which may remove a large fraction of dynamically overexcited asteroids. In the JSCHA simulation, asteroids with low orbital inclination are also observed between 1.8 and 2 au. Equally, these objects do not exist in the real belt today. It is highly likely that these objects will also be removed from this region during the accretion of the terrestrial planets and by the effects of SRs when Jupiter and Saturn reach their current orbits (e.g., n16 is at ∼1.9 au today; Froeschle & Scholl 1989). Figure 5. Precession frequency of asteroids at different distances from the star calculated from linear secular theory in the dynamical system including Jupiter and Saturn. The gas giants are initially as in the JSREG simulation. The n6 and n16 locations in the JSREG simulation corresponds to the intersection between the A curve and the g6 and s6 frequencies, respectively. Frequencies corresponding to 10−5 to 10−4 yr−1 are shown to estimate how much n6 jumps across the belt in the JSCHA simulation. 5 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. 4. PATHS TO CHAOS IN JUPITER’S AND SATURN’S EARLY ORBITS We have argued that Jupiter’s and Saturn’s orbits may have evolved chaotically at early times. In this section we further justify this argument. We first map the prevalence of chaotic motion in the phase space in the orbits of Jupiter and Saturn. Next we perform a series of numerical experiments to mimic the orbital migration and resonant capture of Jupiter and Saturn while they were embedded in the gaseous disk. After the disk’s dissipation, a large number of simulations exhibited chaotic behavior. Finally, we show the long-term dynamical stability between the giant planets evolving in a chaotic resonant configuration. 4.1. A Map of Chaos in Jupiter’s and Saturn’s Orbits To get a sense of the presence of chaos across the phase space of Jupiter’s and Saturn’s orbital configuration, we performed about 9000 simulations to build a dynamical map for a wide range of orbital period ratios between Jupiter and Saturn. The MEGNO (Mean Exponential Growth factor of Nearby Orbits; Cincotta & Simo 2000) chaos indicator is a powerful tool used to identify chaos in dynamical systems. Chaotic orbits are characterized by a large MEGNO value (á ñ Y 2) while regular or quasi-periodic orbits are associated with á ñ Y 2 (e.g., Cincotta & Simo 2000). Our initial conditions were set as follows. Jupiter was placed at 5.25 au while Saturn’s semimajor axis was sampled from about 6.25 to 8.6 au (period ratio between 1.3 and 2.1). Their mutual inclination was sampled randomly from 0° to 0°.5. The eccentricity of Jupiter was randomly selected between 0 and 0.01. Saturn’s initial eccentricity ranges from 0 to 0.1. Note that although the initial eccentricity of Jupiter in this set of simulations is smaller than the corresponding initial values in Figures 1 and 2, this quantity does not remain constant over time. Jupiter and Saturn’s (secular) interaction leads to eccentricity oscillations that imply that Jupiter’s eccentricity may reach values comparable to or even larger than 0.025–0.03. Angular orbital elements of both planets were all sampled randomly between 0° and 360°. Simulations were integrated for 50Myr using the REBOUND integrator (Rein & Tamayo 2015) computing the MEGNO value of each of these dynamical states. Figure 6 shows a dynamical map of the behavior of Jupiter’s and Saturn’s orbits at different orbital separations. The simulations were integrated for fifty million years and the results are color-coded by MEGNO value. Black regions are potentially unstable; they show orbital configurations where Jupiter and Saturn undergo close encounters. Orange regions exhibit chaotic motion (with á ñY 4⪆ ) while blue regions encloses regular motion (with á ñ =Y 2). Chaotic regions are generally confined to orbital period ratios and are associated with specific MMRs. There is a broad chaotic region near the 2:1 MMR and narrower regions close to the 5:3 (period ratio of 1.66) and 7:4 (period ratio of 1.75) resonances, with some chaos present just exterior to the 3:2 resonance (period ratio of 1.5). The dynamical map in Figure 6 shows that chaos is common in the phase space available for Jupiter’s and Saturn’s orbits. Yet their actual orbits at early times were not chosen at random. Rather, the gas giants’ orbital configuration was generated by interactions between the growing planets and the gaseous protoplanetary disk, in particular by a combination of orbital migration (Lin & Papaloizou 1986; Ward 1997) and eccen- tricity and inclination damping (e.g., Papaloizou & Larwood 2000; Bitsch et al. 2013). The dynamical evolution of Jupiter and Saturn in the gas disk depends on the disk properties (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Zhang & Zhou 2010; Pierens & Raymond 2011). To study chaos in Jupiter’s and Saturn’s orbits in the context of orbital migration we perform N-body simulations using artificial forces to mimic the effects of the gas and also pure hydrodynamical simula- tions. We present these results next. 4.2. Chaos in Jupiter’s and Saturn’s Orbits in the Context of Orbital Migration In our simulations of Jupiter and Saturn migrating in a gaseous disk, Jupiter was initially at 5.25 au and Saturn was initially placed exterior to the 2:1 MMR with Jupiter (beyond 8.33 au). Following Baruteau et al. (2014, p. 4293), Jupiter was assumed to migrate in a type-II mode with the migration timescale computed as ⎛ ⎝⎜ ⎞ ⎠⎟n = ´t a M m 2 3 min 1, , 2m,Jup Jup 2 disk Jup ( ) where aJup and mJup are Jupiter’s semimajor axis and mass, respectively. n is the gas viscosity and Mdisk is the gas disk mass. We modeled the disk viscosity using the standard “alpha” prescription given by n a= c Hs (Shakura & Sunyaev 1973), where cs is the isothermal sound speed and H is the disk scale height. In our simulations a = 0.002 and the disk aspect ratio is h∼0.07. To account for the damping of eccentricity and inclination on Jupiter’s orbit we assume the following relationship between migration timescale and Figure 6. A dynamical map of chaos in Jupiter’s and Saturn’s orbits as a function of their initial orbital period ratio. The vertical axes show the initial eccentricity of Saturn (Jupiter’s initial eccentricity was randomly chosen between zero and 0.01). The horizontal axes show the initial period ratio between Jupiter and Saturn. Each dynamical state is color-coded showing the MEGNO value after 50 Myr of integration. Black is used to denote orbits where Jupiter and Saturn’s mutual distance gets smaller than 1au. Because we use a symplectic integrator, these orbits would not be solved properly in this case. Bluish colors show regular or quasi-periodic motion while orange-ish colors show chaotic orbits. 6 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. eccentricity/inclination damping (Crida et al. 2008): = =t t t K. 3e i m,Jup ,Jup ,Jup ( ) In our simulations, we generally adopted the typical value of K=10, although we also performed simulations with K=1 and K=100. To mimic the migration of Saturn in the gas disk, for simplicity we use the type-I migration/damping approach (Papaloizou & Larwood 2000; Tanaka et al. 2002; Tanaka & Ward 2004; Cresswell & Nelson 2006, 2008). Since Saturn’s gap is not fully open, this is an acceptable approximation. We stress that our goal here is only to have convergent and smooth migration of Jupiter and Saturn such that we can access the plausibility of chaos origin in this kind of simulation. The initial gas surface density at the location of Saturn is S = - -r900 g cmSat Sat 0.5 2. To implement migration, eccentri- city, and inclination damping on Saturn, we use the following formulas: ⎜ ⎟ ⎜ ⎟⎛ ⎝ ⎞ ⎠ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝ ⎞ ⎠ ⎛ ⎝ ⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟b = + S + - W- t M m M a h r 2 2.7 1.1 1 1 , 4 m g er h er h k,Sat 2 2 1.3 5 1.1 4 1( ) ( ) ( ) ⎛ ⎝ ⎜⎜ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎞ ⎠ ⎟⎟ = - + + t t e h r e h r e h r i h r 0.780 1 0.14 0.06 0.18 , 5 e,Sat wave 2 3 2 ( ) and ⎛ ⎝ ⎜⎜ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎞ ⎠ ⎟⎟ = - + + t t i h r i h r e h r i h r 0.544 1 0.3 0.24 0.14 6 i,Sat wave 2 3 2 ( ) where ⎜ ⎟ ⎜ ⎟⎛ ⎝ ⎞ ⎠ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎛ ⎝ ⎞ ⎠= S W- t M m M a h r 7 g kwave 2 4 1 ( ) and M , a, m, i, and e are the solar mass, planet’s semimajor axis, planet’s mass, orbital inclination, and eccentricity, respectively. r is the planet’s heliocentric distance. Sg and b are the gas disk surface density and gas surface profile index at the planet’s location, respectively. In our simulations, in the case of Saturn, b was calibrated from hydrodynamical simulations. The synthetic accelerations to account for the effects of the gas on the planet were modeled as = -a v t 8 m mig ,pla ( ) = -a v r r r t 2 . 9e e 2 ,pla ( ) ( ) = -a k v t , 10i z i,pla ( ) where k is the unit vector in the z-direction.The label “pla,” which appears in Equations (8)–(10), takes the form “Jup” or “Sat” to refer to accelerations applied to Jupiter and Saturn, respectively. To ensure the robustness of this method, we compared selected simulations with true hydrodynamical simulations using similar disk setups and they agree in terms of the final period ratio, eccentricities, and orbital inclination of Jupiter and Saturn. In simulations in which the ice giants were also included, they were assumed to migrate inward in the type-I mode described above. We performed 1000 simulations that differed only in the initial surface density of the disk. In each case the gas disk was assumed to dissipate exponentially in 1Myr (tgas=100 kyr). After the gas disk dissipated, simulations were integrated for another 50Myr. To compute how chaotic the orbits of the giant planets are in each simulation we again used the MEGNO chaos indicator incorporated in the REBOUND integrator (Rein & Tamayo 2015). Figure 7 shows the results of these simulations. The vertical axis shows the relative initial surface density of the disk and the horizontal axis shows the final orbital period ratio between Jupiter and Saturn. Each final dynamical state of our simulation is represented with a circle whose color shows the MEGNO value of the final dynamical state of Jupiter and Saturn. The experiment presented in Figure 7 produced more regular configurations of Jupiter and Saturn than chaotic ones. In agreement with pure hydrodynamical simulations (Zhang & Zhou 2010; Pierens et al. 2014), Jupiter and Saturn typically park in either 3:2 (in high-mass disks) or 2:1 (low-mass disks) MMR. We did not find chaos in simulations where Jupiter and Saturn ends in 3:2. All our instances of chaos were related to Saturn and Jupiter’s period ratio being close to 2. We observe chaos in about 1%–25% of these simulations depending on the disk parameters and K-value (a parameter that defines the ratio between the migration timescale and eccentricity/inclination damping; see Equation (3)). We expect that chaotic configura- tions for Jupiter and Saturn in 3:2 may be more likely for larger Figure 7. Final period ratio of Jupiter and Saturn migrating in N-body simulations including the effects of a 1D gas disk. The vertical axis shows the relative surface density of the disk. The gas disk lasts 1 Myr. After the gas dissipation, the orbits of the giant planets are numerically integrated for another 50 Myr. The color code shows the MEGNO value of each system at the end of the simulation. In this set of simulations, there is no case where Jupiter and Saturn suffer close encounters (there are no black points in the figure, which is different from Figure 6). 7 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. eccentricities than our migration and damping prescriptions allow (see Figure 6). We performed 15 hydrodynamical simulations using the GENESIS code (de Val-Borro 2006). None of the simulations with Jupiter and Saturn migrating in the disk generated chaos in Jupiter’s and Saturn’s orbits. However, both hydrodynamical simulations and our N-body simulations with synthetic forces are extremely idealized. A number of factors could change the outcome. First, the protoplanetary disks in these simulations are typically laminar and perfectly axisymmetric. Second, the gas disk’s dissipation in numerical simulations is in general poorly modeled with an exponential density decay over the entire disk. Third, at the end of the gas disk phase, planets should evolve in a sea of planetesimals and leftover building blocks of their own process of formation. In these simulations we did not include any external perturbation in the system (but see Section 2). Fourth, the migrating ice giants (and/or their building blocks; Izidoro et al. 2015a) represent a further source of perturbations. In fact, the results shown in Figures 6 and 7 are complementary in nature. The setup of the simulations in Figure 7 favor deep capture in resonance and regular motion between Jupiter and Saturn while that in Figure 6 (because of the random selection of angles) may preferentially put the dynamical state slightly off the libration center and favor mostly chaotic orbits. The existence of many chaotic config- urations in Figure 6 demonstrates the prevalence of chaos and even if the gas giants’ orbits behaved regularly immediately after the dissipation of the gaseous disk, there are a number of processes that could have rendered them chaotic. 4.3. Long-term Dynamical Stability of the Gas Giants in a Chaotic Configuration The present-day orbits of the giant planets are thought to have been achieved by an instability in the giant planets’ orbits that occurred long after the dissipation of the gases (Tsiganis et al. 2005; Morbidelli et al. 2007; Levison et al. 2011). Our model is entirely consistent with a late (∼500 Myr-later) instability in the giant planets’ orbits, whatever the exact configuration of Jupiter and Saturn (see Batygin et al. 2012; Nesvorný & Morbidelli 2012; Pierens et al. 2014). In fact, our model is also consistent with an earlier giant planet instability(Kaib & Chambers 2016), as long as there is a sufficiently long interval during which the belt can be chaotically excited. This interval is roughly longer than 10Myr in most of the simulations we have run, but can be as short as 2 Myr. To illustrate the possibility of a long-term stability between the giant planets evolving in chaotic orbits, we performed N-body simulations in which Jupiter, Saturn, Uranus, and Neptune migrate in a disk. The prescription for migration used in these simulations is analogous to that explained before in this text. Isothermal type-I migration and damping are also applied to the ice giants. Here, we again used those damping timescales defined in Equations (4)–(7) and accelerations given by Equations (8)–(10). Figure 8 shows a long-term stable dynamical evolution of the gas giants in chaotic orbits. At the end of the gas disk phase, Jupiter and Saturn are locked in 2:1 MMR, Saturn–Uranus in 2:1 MMR, and Uranus–Neptune in 3:2 MMR. The system is dynamically stable for over 500Myr. 5. DISCUSSION Chaotic excitation of the asteroid belt represents a novel mechanism for solving a longstanding problem in planetary science. Yet there are a number of questions that arise when considering whether this mechanism is consistent with our current vision of the solar system’s global evolution. In this section we address a number of questions. 5.1. Chaotic Excitation versus the Grand Tack The “small Mars” problem highlights the fact that a mass deficit is needed from ~1 4 au– to explain Mars’ small mass relative to Earth’s(Wetherill 1996; Raymond et al. 2009). In the Grand Tack model, this deficit is generated by Jupiter’s long-distance orbital migration from several au to 1.5 au, then back out to beyond 5 au(Walsh et al. 2011; Jacobson & Morbidelli 2014; Raymond & Morbidelli 2014; Brasser et al. 2016) By explaining the asteroid belt’s orbital structure, our results revive models in which the asteroid belt was initially very low mass (Izidoro et al. 2015b; Levison et al. 2015). While standard disk models typically invoke a smooth mass distribution within disks(for example, the minimum-mass solar nebula model of Figure 8. Example of the genesis of chaos in the giant planets’ orbits and their long-term stability. Top: orbital migration of Jupiter, Saturn, Uranus, and Neptune embedded in the gaseous protoplanetary disk. At the end of the gas disk phase, the planets are locked in a chain of MMRs (2:1, 2:1, 3:2). Bottom: long-term evolution of the Saturn-to-Jupiter orbital period ratio, and the gas giants’ eccentricities and inclinations over the next 500 Myr. This configuration is consistent with a later instability in the orbits of the giant planets (Tsiganis et al. 2005; Levison et al. 2011; Nesvorný & Morbidelli 2012). 8 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. Weidenschilling 1977 and Hayashi 1981 generates a smooth disk from discrete planets), it remains unclear whether the planetary building blocks embedded in these gaseous disks should really follow a smooth distribution. In fact, models often find that planetesimals preferentially grow in special locations within the disk, such as at pressure bumps (e.g., Johansen et al. 2014, p. 547). Some models naturally create confined rings of planetesimals within broad disks (Surville et al. 2016). There also exist mechanisms to systematically drain solids from certain areas of the disk, thereby creating localized depletions and enhancements (e.g., Levison et al. 2015; Moriarty & Fischer 2015; Drazkowska et al. 2016). Our model thus forms the basis of an alternative to the Grand Tack model. Within the context of this model, Mars’ small mass can be explained by a broad mass depletion between Earth’s and Jupiter’s orbits. The asteroid belt, which could not be stirred by the dynamical effects of local embryos, was chaotically excited by Jupiter and Saturn. The next step is to search for ways to distinguish between the Grand Tack and this new chaotic model. Tests may be based on observations of small bodies in the solar system or geochemical measurements. Alternately, the models may be differentiated by more detailed studies of the underlying physical mechan- isms involving planetary orbital migration and pebble accretion. Like the Grand Tack, our model assumes that Jupiter and Saturn migrated during the gas disk phase into a resonant configuration, mostly likely to 3:2 or 2:1 MMR (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Pierens & Nelson 2008; D’Angelo & Marzari 2012; Pierens et al. 2014). However, the scale of radial migration of Jupiter and Saturn in our scenario could be less specific and much smaller than that in the Grand Tack scenario, where Jupiter and Saturn migrated inward then outward. In our case—since it is quite unlikely that Jupiter and Saturn grew already in resonance —only an inward convergent phase of migration between the giant planets is needed to put them into a resonant configuration. The Grand Tack and the chaotic excitation models are not contradictory models. In the Grand Tack model, some level of chaotic excitation could have operated between Jupiter’s and Saturn’s two-phase migration and the Nice model instability. However, since the asteroid belt is already sufficiently dynamically excited after the Grand Tack (Deienno et al. 2016) it would be unnecessary to invoke chaotic excitation. We also note that both the Grand Tack and chaotic excitation mechanisms may operate with Jupiter and Saturn in either the 3:2 or 2:1 MMR (see also Pierens et al. 2014). Of course, the 3:2 resonance has been much more carefully studied for the Grand Tack, and our results suggest that the 2:1 resonance may be favored with regards to chaotic excitation. Yet we caution that Jupiter and Saturn’s configuration during the disk phase may not necessarily differentiate between the two models. 5.2. A Chaotic Young Solar System? It is entirely reasonable to imagine a chaotic young solar system. Several exoplanetary systems with planets on near- resonant orbits are thought to be chaotic, such as 16 Cyg B (Holman et al. 1997), GJ876 (Rivera et al. 2010), and Kepler 36 (Deck et al. 2012). The present-day solar system is well known to be chaotic (Laskar 1989, 2003; Sussman & Wisdom 1992). The orbits of the terrestrial planets undergo chaotic diffusion on a timescale of a few million years (Laskar 1989; Batygin et al. 2015). It is unknown whether the present-day giant planets are in a chaotic configuration; an accurate determination is precluded by uncertainties in our knowledge of the planets’ orbital positions (Michtchenko & Ferraz-Mello 2001; Hayes 2008). It is often assumed that the early solar system was characterized by regular motion of the planets (Brasser et al. 2013). Our work suggests that this may not have been the case, and that the structure of the asteroid belt is a signpost that Jupiter’s and Saturn’s early orbits were in fact chaotic. We showed that a small nudge could trigger chaos in Jupiter’s and Saturn’s early orbits (Figure 1). Could another small nudge or accumulated effects of scattering planetesimals make the system regular again? It is important to note that when the scattering of planetesimals takes place, not only damping of eccentricity but also radial migration should be observed. This process would lead to the divergent migration of Jupiter and Saturn (and ice giants) and not necessarily sinking toward the resonant center and consequently regular orbits. Typically, a late dynamical evolution is envisioned in models of solar system evolution (Gomes et al. 2005; Levison et al. 2011) which implies Jupiter and Saturn hanging out in or very near resonances for hundreds of Myr after gas disk phase. What is required for chaotic excitation to work is simply a sufficient time interval in this phase while the giant planets remain in a chaotic state. This also depends on the chaotic configuration of Jupiter and Saturn and the setup over the disk of planetesimals (Nesvorný & Morbidelli 2012). Obviously, it would be computationally challenging to perform a systematic analysis to identify all dynamical configurations of Jupiter and Saturn that could excite the belt. However, we indeed found that non-resonant or temporary resonant configurations can also excite the belt. We also recognized that it is not clear why previous classical simulations of terrestrial planet formation did not find similar effects in the belt. One possibility is that the chaotic effects on asteroids have been erased or mitigated by the typical presence of large planetary embryos in the belt in classical simulations (Chambers 2001; O’Brien et al. 2006; Raymond et al. 2006; Izidoro et al. 2014, 2015b). Also, these previous simulations have preferentially considered Jupiter and Saturn near their current orbits or in 3:2 MMR (in almost circular orbits; see for example Raymond et al. 2009). In our most successful simulations of belt excitation Jupiter and Saturn have eccentricities of about 0.03–0.05; the latter is consistent with results from hydrodynamical simulations (Pierens et al. 2014). 5.3. The Absence of “Fossilized” Kirkwood Gaps The Kirkwood gaps in the asteroid population are created by MMRs with Jupiter. The most prominent is the gap created by Jupiter’s 3:1 MMR centered at 2.50 au in the present-day belt. The chaotic excitation of the asteroid belt likely took place before a late instability in the giant planets’ orbits (the so-called Nice model; Tsiganis et al. 2005; Morbidelli et al. 2007; Levison et al. 2011). If a late instability shifted Jupiter’s orbit inward by ∼0.2 au (Tsiganis et al. 2005) then there should exist a fossilized gap in the asteroid belt at about 2.6 au, just exterior 9 The Astrophysical Journal, 833:40 (18pp), 2016 December 10 Izidoro et al. to the current Kirkwood gap associated with the 3:1 resonance. No such gap exists (see Section 6.4 of Morbidelli et al. 2010). We caution that the data set used to study fossilized gaps (or lack thereof) are relatively sparse, containing just 335 large asteroids (with absolute magnitude