Excitation and Depletion of the Asteroid Belt in the Early Instability Scenario

The modern asteroid belt's (AB) structure starkly contrasts with that of the terrestrial and giant planet systems in that it contains less than ∼5 × 10⁻⁴ M_⊕ of material on dynamically excited (large eccentricities and inclinations) orbits (DeMeo & Carry 2013; Kuchynka & Folkner 2013). Though hundreds of thousands of constituents have been observed, around half of the main belt's (MB) mass is concentrated in just four asteroids (in order of descending mass: Ceres, Vesta, Pallas, and Hygiea). Furthermore, the composition of the belt is far from homogeneous. The inner MB is primarily composed of silicate-rich, moderate-albedo S-types, whereas the belt's outer regions are dominated by carbonaceous, low-albedo C-types (Chapman et al. 1975; Gradie & Tedesco 1982; Gradie et al. 1989; Bus & Binzel 2002; DeMeo & Carry 2013). However, these two populations overlap substantially, and C-types account for around two thirds of all large (D > 50 km) bright (H < 9.7) asteroids (Campins et al. 2018).

Although the present-day AB is well constrained, the primordial belt is not. In particular, there is a huge disparity in the assumed total mass in the early belt. The primordial AB's total mass as inferred from disk models (Hayashi 1981; Bitsch et al. 2015) was around ∼1–5 M_⊕. However, what counts from a dynamical point of view is the asteroidal mass in planetesimals and planetary embryos. Simulations of planetesimal formation in evolving disks have found a range of outcomes (e.g., Carrera et al. (2017) and Dra̧żkowska & Dullemond 2018), including the possibility of forming a ring of planetesimals around ∼1 au, but none in the AB (Dra̧żkowska et al. 2016). The low-mass AB model (Izidoro et al. 2015; Raymond & Izidoro 2017b) discussed below makes the assumption that very few planetesimals originated in the belt.

Because the primordial mass distribution is not well characterized, studying the long-term evolution of the MB can be challenging (Bottke et al. 2015). Meteorites believed to originate from Vesta indicate that it differentiated and formed a crust just a few megayears (Shukolyukov & Lugmair 2002) after calcium aluminum-rich inclusion formation, which are the oldest samples known in the solar system. Thus, the most massive asteroids grew to their present size when the solar system was in its infancy and have not grown substantially larger since. This is primarily the result of the belt's high degree of dynamical excitation greatly lengthening accretion timescales.

Over gigayear timescales, collisional grinding and fragmentation tend to push the belt's size distribution toward smaller diameters (Bottke et al. 2005b). This seems consistent with the large number of known collisional families, particularly in the inner MB (Bottke et al. 2006a; Walsh et al. 2013). In fact, Dermott et al. (2018) argued that nearly all asteroids in the inner MB are members of collisional families. Using constraints for the formation of Vesta's two enormous craters (the 505 km Rheasilvia crater on the south pole and the underlying 395 km Veneneia crater; Schenk et al. 2012), Bottke et al. (2005a) calculated the probability that these impacts occurred during the complete collisional evolution of the MB. If the two basins were formed in the last ∼2 Gyr, the probability of both events happening given the current MB size distribution would be ∼1% (Bottke et al. 2015). The lack of other such basins on Vesta would therefore imply that the belt has always had a low mass. Thus, from a probabilistic standpoint, it seems unlikely that the AB has lost a substantial amount of mass over the last several gigayears. However, this implication would not hold if Vesta was implanted into the belt from a different region of the solar system (Bottke et al. 2006b; Mastrobuono-Battisti & Perets 2017; Raymond & Izidoro 2017a).

Mean motion resonances (MMRs) and secular resonances are the dominant mass loss mechanisms in the modern AB. MMRs occur when an asteroid's orbital period is in integer ratio with that of another body, in this case Jupiter. Secular resonances are the result of an object's longitude of perihelion ($\dot{\varpi }$) or longitude of ascending node ($\dot{{\rm{\Omega }}}$) precession frequency equaling one of the solar system's dominant eigenfrequencies. Because the ν₅ and ν₆ secular resonances overlie several dominant MMRs with Jupiter (4:1, 3:1, 5:2, 7:3, 2:1), asteroids in these regions are quickly excited onto planet-crossing orbits (Morbidelli et al. 1995). In addition to depleting the MB, these processes also form the Kirkwood gaps in the belt's orbital distribution (Petit et al. 2001; O'Brien et al. 2007; Deienno et al. 2016). In spite of these loss mechanisms, Minton & Malhotra (2010) concluded that the MB has lost only 50% of large asteroids since attaining its current dynamical state. However, the simulations of Minton & Malhotra (2010) only modeled large asteroids as test particles. This means that around two to three orders of magnitude worth of AB depletion must be accounted for before the end of the planet formation epoch (for recent summaries on the evolution of the AB, consult Bottke et al. 2015 and Morbidelli et al. 2015).

Most models for terrestrial planet formation account for depletion at the 99%–99.9% level with some combination of the following mechanisms: primordial depletion, giant planet influence, and embryo excitation. In the classical model of terrestrial planet formation, a population of ∼100 Moon-to-Mars massed planet-forming embryos extend throughout the inner terrestrial disk and MB region (Wetherill 1992; Chambers & Wetherill 1998; Chambers 2001, 2007; Chambers & Wetherill 2001; Raymond et al. 2009). The orbital excitation provided by a population of embryos can result in substantial mass loss (Petit et al. 2001; O'Brien et al. 2007). However, these standard initial conditions fail in that they systematically produce overmassed Mars analogs and Mars-to-Earth massed planets in the AB (Chambers 2001; Raymond et al. 2009; Morbidelli et al. 2012). Furthermore, given the strong radial dependence of embryo growth (e.g., Kokubo & Ida 1996, 1998, 2000), it is possible that large embryos never existed in the young belt. Thus, the most compelling solutions to the so-called "small Mars problem" are those that also deplete the primordial MB and sufficiently mix the radial distribution of S- and C-types. Here, we summarize the various classes of models (see Raymond et al. 2018 for a review of the models' assumptions).

1. Low-mass AB model: Izidoro et al. (2014) showed that the inner solar system could be consistently replicated if built from a steep initial radial mass distribution. These initial conditions, wherein the primordial Mars-forming and AB regions never contained a substantial amount of mass (Izidoro et al. 2015), are largely consistent with modern pebble accretion simulations of embryo and planetesimal formation (Levison et al. 2015; Dra̧żkowska et al. 2016). Raymond & Izidoro (2017a) provided the explanation for the MB's compositional dichotomy by suggesting that planetesimals near the growing giant planets could be destabilized via aerodynamic drag and scattered into the empty primordial belt. Thus, the MB's current population of C-types originated in the outer solar system prior to being implanted throughout the AB. The major weakness of the low-mass AB model lies in the initial conditions and whether such a steep initial mass distribution profile is realistic.

2. The "Grand Tack" hypothesis: By assuming that the terrestrial planets formed out of a narrow annulus of material extending from ∼0.7–1.0 au, Hansen (2009) consistently replicated the terrestrial planets' mass distribution, in particular the large mass ratios between neighboring planets (Mercury/Venus and Mars/Earth). Walsh et al. (2011) provided a dynamical mechanism for these initial conditions by proposing that, during the solar system's gas disk phase, Jupiter and Saturn migrated in and out of the inner solar system (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Pierens et al. 2014). Jupiter's presence serves to shepherd primordial MB objects onto orbits where they are scattered out of the region. Thus, the MB was already substantially depleted during the gas disk phase. Furthermore, Deienno et al. (2016) investigated the effects of the Nice Model (an orbital instability among the giant planets; the leading evolutionary model for the outer solar system Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005) on the post-Grand Tack MB distribution and found that the final orbital structure was largely consistent with the current AB. However, the high inclination orbital parameter space of the post-instability belt in that work was overpopulated. Although the Grand Tack scenario succeeds at explaining the radial mixing of C- and S-types, the outward migration mechanism is highly dependent on the supposed disk structure and gas accretion rates (D'Angelo & Marzari 2012; Pierens et al. 2014; Raymond & Morbidelli 2014).

3. An early instability: The classical Nice Model (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005) was originally timed in conjunction with the Late Heavy Bombardment, an inferred delayed spike in the lunar cratering record around ∼400 Myr after the planets formed. Recent evidence (Boehnke & Harrison 2016; Zellner 2017; Morbidelli et al. 2018; Nesvorný et al. 2018) has called the Late Heavy Bombardment's existence into doubt. Because the instability need not be tied to a specific time, moving its occurrence earlier might prevent the disruption of the fully formed terrestrial planets. Indeed, planet ejections and collisions are common in simulations of a delayed Nice Model instability (Brasser et al. 2009; Agnor & Lin 2012; Brasser et al. 2013; Kaib & Chambers 2016). Clement et al. (2018) showed that timing the instability ∼1–10 Myr after gas disk dispersal substantially limits the mass and formation time of Mars. This provides a natural explanation for the differences in the inferred geological formation times of Earth (Kleine et al. 2009) and Mars (Dauphas & Pourmand 2011). Although the coadded simulations of Clement et al. (2018) matched the broad orbital structure of the MB, and depleted the region at greater than the ∼95% level, they began with unrealistic populations of large embryos and planetesimals (with the smallest simulation particles being more than an order of magnitude larger than the entire current AB mass). In a study similar to this one, Deienno et al. (2018) studied the effect of an early "Jumping Jupiter" style instability on a primordial terrestrial disk constructed of 10,000 massless test particles. The final MB orbital distribution of Deienno et al. (2018) was a good match to the real one, and the authors argued that their inferred net depletion of ∼90% would be consistent with the low-mass AB framework of Izidoro et al. (2015) and Raymond & Izidoro (2017b). However, the Deienno et al. (2018) simulations only considered the Jumping Jupiter style instability (Brasser et al. 2009; Nesvorný et al. 2013). In addition, the integrations used massless test particles and relied on a primordial low-mass AB to match the current low mass of the MB.

Here, we expand on the early instability framework of Clement et al. (2018) with detailed simulations of the scenario's consequences in the AB. Our work differs from that of Deienno et al. (2018) in that we model the MB with fully self-gravitating bodies and include control cases with test particles. Additionally, we assume a primordial AB mass (∼2 M_⊕) consistent with the value derived from disk models. Because each simulated instability is highly chaotic, and therefore inherently unique, our study investigates a range of instabilities. By placing the specific dynamical state of the solar system's giant planets on our spectrum of simulated instabilities, our study seeks to infer the approximate range of MB depletion possible in an early Nice Model instability. Thus, although Deienno et al. (2018) analyzed one particular instability with an outcome very similar to the modern solar system, our work scrutinizes a range of instabilities (none of which is a perfect match to the solar system) and looks for trends.

We begin by selecting eight simulations from Clement et al. (2018) that best replicated both the inner and outer solar system structure in accordance with the success criteria from that work. Using the naming conventions of Clement et al. (2018), these simulations are n1/1 Myr/6, n2/10 Myr/14, n2/10 Myr/44, n2/1 Myr/40, n1/10 Myr/15, n1/1 Myr/46, n1/10 Myr/29, and n2/10 Myr/25 (number of additional primordial ice giants/instability delay time/run number; henceforward referred to as runs 1–8, respectively). We commence our study with each of these systems evolved up to the instability time. In Clement et al. (2018), instability simulations were created by first separately evolving and then combining two independent sets of integrations. The first was a set of giant planet resonant configurations (Lee & Peale 2002) interacting with an external disk of primordial Kuiper Belt objects (Nesvorný & Morbidelli 2012; Deienno et al. 2017) up until the instability time. The second followed the evolution of a 5 M_⊕ terrestrial-forming disk of 100 embryos and 1000 planetesimals in the presence of a static Jupiter and Saturn (locked in a 3:2 MMR). Snapshots of the terrestrial disks at various phases of evolution (1.0 or 10.0 Myr for the systems concerned in this article) were then added to the giant planet simulations and integrated through the Nice Model instability for an additional 200 Myr. It is at this combination time that we select the initial conditions for the simulations of this article.

Accurately modeling the AB with numerical integrators (Duncan et al. 1998; Chambers 1999; Rein & Tamayo 2015) is difficult because of the vast disparity between the sizes of modern asteroids and the belt's extrapolated primordial mass. Simulating the primordial belt using Ceres-sized particles would require close to 10,000 objects, exceeding the capabilities of the conventional integrators used to study terrestrial planet formation. Doing so using particles the size of Hygiea would involve over 100,000 individual primordial objects. Thus, most authors are forced to model the detailed structure of the MB with massless test particles, or approximate its early state with unrealistically massive planet embryos and planetesimals. Improving simulation resolution with test particles can be a useful tool for studying complex orbital dynamics in models where the MB is already heavily depleted in the gas disk phase because of the giant planet's influence (Walsh et al. 2011; Izidoro et al. 2015; Raymond & Izidoro 2017a, 2017b). In this low surface density limit, the AB can be well characterized with test particles. However, the early instability scenario as described in Clement et al. (2018) assumes no prior depletion in the AB at the start of the giant impact phase. Because of this higher AB surface density, characterizing the levels of depletion in a fully self-gravitating AB with somewhat realistic mass resolution is an important step in validating the viability of the early instability scenario.

To study the evolution of the MB in sufficient detail, we select the GENGA (Gravitational Encounters with GPU (Graphics Processing Unit) Acceleration) simulation package (Grimm & Stadel 2014). GENGA is based on the Mercury hybrid integrator (Chambers 1999) and runs on Nvidia GPUs. For each system, we remove all terrestrial-forming objects in the AB region (a > 2.0 au) and replace them with 3000 primordial asteroids, each four times the current mass of Ceres. Thus, the total MB mass is ∼1.8 M_⊕. Our selection of total initial mass is derived from the average pre-instability mass of the a > 2.0 au region in runs 1–8.

To create our 3000 particle ABs, we assume a surface density profile proportional to r^−3/2 (Birnstiel et al. 2012). Eccentricities and inclinations are drawn from near-circular and coplanar Gaussian distributions (σ_e = .02 and σ_i = .2°). The remaining orbital elements are selected randomly from uniform distributions. Next, we integrate our model AB in the presence of a static Jupiter and Saturn in a 3:2 MMR (a_j = 5.6 au, a_s = 7.6 au) for 1.0 Myr using a 50 day time-step; 305 accretion events occur during this initial phase of evolution, with the largest object growing to five times its initial size. These lengthy accretion times are roughly consistent with the results of semianalytical studies of runaway growth and embryo accretion (Kokubo & Ida 1996, 1998). However, we do note that such oversimplified initial conditions might influence our results. Indeed, constructing the primordial belt strictly out of large bodies might overestimate the effects of planetesimal scattering and dynamical spreading in the disk. Furthermore, our simulations are not designed to simulate fragmenting and hit-and-run collisions. Because we are studying MB dynamics in the high surface density limit, not accounting for the dynamical damping effects of such collisions (Chambers 2013) might artificially increase the orbital excitation in our fully evolved systems.

The eccentricity and inclination structure of our self-stirred initial MB is plotted in Figure 1. We perform 18 separate, 200 Myr integrations (Table 1) of this AB embedded within our respective Nice Model instability simulations using a 6 day time-step as follows (note that each individual simulation contains over 4000 objects including the forming inner terrestrial disk, 3000 asteroids, unstable giant planets, and primordial Kuiper Belt):

• (1)  
  Eight systems (runs 1–8) are integrated with all simulation objects treated as fully self-gravitating.
• (2)  
  The same eight systems (now denoted runs 1a–8a) are integrated with the asteroids treated as test particles, and all other objects as fully self-gravitating. Simulations presented in Levison et al. (2011) indicated that disk self-gravity is important in driving the evolution of the young Kuiper Belt. By performing two sets of simulations both with and without asteroid self-gravity, we are able to test whether accounting for self-gravity affects the overall depletion or orbital structure in the MB.
• (3)  
  For systems 1 and 2, we perform an additional batch of integrations (runs 1b and 2b) where we embed four Mercury-massed planet embryos within the MB. This allows us to probe whether embedded embryos are necessary for the success of the early instability scenario (Clement et al. 2018 found that much of the depletion in the scenario was due to embryo excitation).

Zoom In Zoom Out Reset image size

Traditionally, the success criteria considered by dynamical studies of the AB's early evolution fall into three general categories: depletion, radial mixing of C- and S-types, and orbital structure (specifically producing a dynamically excited MB without overpopulating the high inclination parameter space above the ν₆ secular resonance). In the subsequent three sections, we address each of these constraints individually.

3.1. Depletion

Because each instability is inherently unique, it is important to place the solar system outcome on the spectrum of instabilities our systems experience (and corresponding levels of depletion). It is important to note here that none of our final giant planet configurations simultaneously replicate all the important qualities of the modern outer solar system (e.g., success criteria A-D in Nesvorný & Morbidelli 2012 and E in Deienno et al. 2017). Of our 18 simulations, only seven finished with the correct number of giant planets. One simulation finished with five giant planets, and the instability of another (run 4a) was so violent that Saturn was ejected. However, as shown directly in Deienno et al. (2018), the perturbative effects of the ice giants on the MB are relatively weak compared to that of Jupiter and Saturn (Morbidelli et al. 2009). Therefore, systems with significantly different ice giant evolutionary schemes might still be good analogs to the solar system if their Jupiter/Saturn systems are similar. The ratio of Saturn's orbital period to that of Jupiter is ∼2.48, and eight of our systems have period ratios between 1.95–2.65 (note, N_GP = 4 and P_S/P_J < 2.8 are success criteria A and D in Nesvorný & Morbidelli 2012 and Deienno et al. 2017, respectively).

To compare the final orbital structure of our various Jupiter/Saturn systems with the solar system (besides run 4a that lost Saturn), we begin by calculating the normalized angular momentum deficit (AMD, Equation (1)) for each system (Laskar 1997). AMD quantifies the degree to which a system of orbits differs from that of a coplanar, circular system. Thus, a simulation that finishes with an AMD_JS close to that of the modern solar system would possess a Jupiter/Saturn pair with combined dynamical excitation similar to the actual planets. We plot the relationship between AB depletion and AMD_JS in the top panel of Figure 2. There is significant scatter in the data because the instability is a highly chaotic event. However, the overall trend is very clear. Systems that finish with a more dynamically excited Jupiter/Saturn pair also experience greater depletion in the AB. Given these results, the systems that best match the solar system's value for AMD_JS experience one to three orders of magnitude worth of depletion in the AB.

$$\begin{eqnarray}&&\mathrm{AMD}=\displaystyle \frac{{\displaystyle \sum }_{i}{m}_{i}\sqrt{{a}_{i}}[1-\sqrt{(1-{e}_{i}^{2})}\cos {i}_{i}]}{{\displaystyle \sum }_{i}{m}_{i}\sqrt{{a}_{i}}}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

However, many of our systems with values of AMD_JS similar to that of the actual solar system are poor analogs for other reasons. In particular, the secular architecture of the solar system is dominated by the highly excited g₅ mode (success criteria C in Nesvorný & Morbidelli 2012). For a complete discussion of the solar system's secular structure, consult Brasser et al. (2009) and Morbidelli et al. (2009). To compare the secular landscapes of our simulated systems with the solar system, we calculate e₅₅ (the amplitude of the Jupiter's g₅ eigenmode; e_{55, SS} = 0.044) and e₅₆ (the amplitude of Saturn's forcing on Jupiter's eccentricity; e₅₆ =0.016) for each integration via Fourier analysis (Šidlichovský & Nesvorný 1996). The bottom two panels of Figure 2 show the dependence of AB depletion on the values of e₅₅ and e₅₆. Again, in spite of significant data scatter and the small number statistics, both plots show the same general trend of increased secular amplitudes correlating with larger AB depletion fractions. The relationship between e₅₅ and AB depletion is perhaps the clearest of the three panels in Figure 2. Multiple previous studies have concluded that the sufficient excitation of e₅₅ is supremely important in driving the secular evolution of the solar system (Morbidelli et al. 2009; Nesvorný & Morbidelli 2012; Clement et al. 2018). Therefore, it is promising that all but one of the systems in our study that sufficiently excite e₅₅ to within a factor of two of its present value experience at least two orders of magnitude worth of depletion in the AB. This depletion, which is related to orbital excitation, is thus largely driven by e₅₅ transmitting excitation to the MB via large forced vectors during the instability phase. Our results therefore indicate that the early instability scenario should be able to account for AB depletion at the 99%–99.9% level. As an example, run 3 finishes with N_GP = 4, P_S/P_J = 2.3, AMD_JS = 2.7 × AMD_{JS, SS}, e₅₅ = 0.62 × e_{55, SS}, and 0.13% of AB objects remaining. Furthermore, achieving this level of AB depletion and successfully forming the inner solar system are not mutually exclusive traits in our simulations. To illustrate this, Figure 3 plots the final planets from run 2a, compared with the current solar system. The stochastic nature of planet formation makes it impossible to select a perfect solar system analog out of just 18 simulations. In run 2a, the eccentricities of Earth and Venus are too large, the first three outer planets are underexcited, Mercury is overmassed, Neptune's semimajor axis is incorrect, and an extra Mars analog forms near the inner edge of the AB. However, we still plot this result to demonstrate how the run depletes 91.2% of the AB's mass without totally inhibiting the formation of the terrestrial planets.

Zoom In Zoom Out Reset image size

It is also important to point how challenging it is for N-body simulations of the Nice Model instability to simultaneously replicate all facets of the outer solar system (Nesvorný & Morbidelli 2012; Deienno et al. 2017; Gomes et al. 2018). In particular, few simulations are able to sufficiently excite AMD_JS and e₅₅ while keeping P_S/P_J low. Nesvorný & Morbidelli (2012) performed over 5000 instability simulations using a variety of different giant planet configurations and only 5% of the most successful batch simultaneously met constraints for e₅₅ and P_S/P_J (success criterion C and D in that work). Within our sample of 18 instabilities, only one (run 3) successfully excites the orbits of Jupiter and Saturn while maintaining P_S/P_J < 2.8. Though it is encouraging that this run depletes the MB at close to the 99.9% level, our results should still be taken in context with the fact that all our other runs that experience three orders of magnitude worth of MB depletion (Table 2) finish with P_S/P_J > 2.8.

Table 2.  Table of Results for Our 18 Runs

Run	AMDJS	PS/PJ	MAB, f/MAB,fi	ν6 Ratio
1	0.048	2.35	0.23	0.31
2	0.097	1.93	0.34	2.31
3	2.33	2.77	0.0014	inf
4	74.3	21.5	0	N/A
5	1.72	3.94	0.005	inf
6	0.057	2.08	0.030	4.25
7	47.7	17.9	0	N/A
8	0.093	3.94	0.17	0.50
1a	1.03	5.16	0.15	0.60
2a	0.096	2.53	0.088	2.00
3a	5.84	4.12	0.0033	3.00
4a	N/A	N/A	0.036	0.91
5a	0.074	2.41	0.14	1.45
6a	0.78	4.37	0.0064	0.67
7a	21.29	6.19	0	N/A
8a	12.7	8.02	0.022	1.67
1b	0.31	2.70	0.13	1.00
2b	0.0086	2.23	0.18	inf

Note. The columns are as follows: (1) the run number, (2) the AMD of the Jupiter/Saturn system, (3) Jupiter and Saturn period ratio, (4) the run's AB depletion fraction, and (5) the ratio of asteroids above to below the ν₆ ratio. Only asteroids not on planet-crossing orbits are considered.

Download table as:  ASCIITypeset image

Though the differences in depletion trends between our fully self-gravitating and test particle ABs are not statistically significant, the null result is still important. Nearly all dynamical studies of early depletion in the AB only consider test particles (O'Brien et al. 2007; Deienno et al. 2016, 2018; Izidoro et al. 2016). As discussed in Section 2, test particles provide a good approximation in the low surface density limit. However, because we assume no prior depletion in the AB, including the effects of asteroid self-gravity could be important when studying the total early depletion.

3.2. Radial Mixing

It is inherently difficult for N-body studies of terrestrial planet formation to explain the compositional dichotomy between S- and C-types because there is no consensus as to where each type originated from in the disk. However, nucleosynthetic differences between ordinary and carbonaceous chondrite groups (Burkhardt et al. 2011) indicate that they formed at different radial distances (Kruijer et al. 2017). Early studies of planet formation in the inner solar system assumed that the different populations exist today because of a primordial division (commonly referred to as a "snow line") in the material (Wetherill 1992; Chambers & Wetherill 1998). Grimm & McSween (1993) proposed that S-types formed inside ∼2.7 au because of the higher probability of capturing energy from radionuclides (e.g., ²⁶Al). In the same manner, C-types formed exterior to ∼2.7 au where ²⁶Al is extinct. However, whether ice exists at a snow line depends on the drift of small particles (Ciesla & Cuzzi 2006). Studies of disk dynamics seem to indicate that snow lines tend to migrate inward (Lecar et al. 2006; Kennedy & Kenyon 2008; Martin & Livio 2012). Further complicating the problem, this migration can be blocked by the growing massive planets (Lambrechts et al. 2014), and thus the condensation temperature at any given point in the disk may not be correlated with the presence of volatiles (Morbidelli et al. 2016). Alternatively, Raymond & Izidoro (2017a) proposed that the AB could have already been populated with material from the outer disk (∼5–20 au) via scattering events during the giant planet growth phase.

We adopt a simple model to test how the giant planet instability mixes material from different radial origins in the AB. We assume that our simulations begin with the outer AB (a > 2.8 au) highly populated with C-types, and the inner MB dominated by S-types. In Figure 4 we plot the final distributions of asteroids originating from inside and outside 2.8 au for our three simulation sets. Though this figure is biased toward simulations with less depletion, our fully interacting runs provide the best match to the current two distributions. The majority of the total belt material is in the form of C-types, and the interior population of S-types is significantly mixed with the outer C-types. Specifically, 67% of all remaining asteroids in these runs are considered C-types (originate outside 2.8 au), as compared to 60% of all modern asteroids with D > 50 km; 43% of those same asteroids that reside in the inner MB (a < 2.8 au) are C-types (compared with 49% in the solar system). Furthermore, our runs with more violent instabilities (and large depletion fractions) where P_S/P_J finished greater than 2.8 display a similar mixing trend. In these runs, 76% of all surviving asteroids are C-types. Thus, our mixing results are loosely consistent across the range of our instabilities. Though this trend is similar for all three simulation sets, the test particle runs overdeplete both the inner MB and total S-type inventory. Contrarily, the embryo runs scatter a significant number of C-types inward, but fail to scatter the primordial S-types outward.

Zoom In Zoom Out Reset image size

In general, the test-particle-only runs (1a–8a) overdeplete the inner MB of S-types because there is no dynamical friction present to save material with a < 2.5 au from loss when the ν₁₆ and ν₆ secular resonances sweep through the belt. Indeed, our test particle simulations lose an average of 29% more mass in the inner MB (2.0 < a < 2.5 au) in the 100 Kyr after the first ice giant ejection than in our other simulations. This is consistent with other studies of resonant sweeping in the AB (Walsh & Morbidelli 2011). However, it is also possible that this is the result of small number statistics and the limited number of simulations used in our study. When disk self-gravity is included (runs 1–8), the overall distribution is a much better match to that of the actual solar system.

Though the populations of S- and C-types in runs 1b and 2b are not as well mixed as those in runs 1–8, the small number statistics makes it difficult to say whether the difference is significant. Indeed, runs 1b and 2b both experience uncharacteristically weak instabilities when compared with our other runs (see Figure 2). Both only excite e₅₅ to less than half the modern value and finish with an AMD_JS less than a third that of AMD_{JS, SS}. Because the giant planets' resonant perturbations were less strong in these simulations, it makes sense that the final MB populations are less mixed than those of the other simulation sets.

Nevertheless, the radial mixing results of our interaction set are encouraging and in decent agreement with the modern solar system's distribution (though this assumes perfect zoning at the beginning of our runs, and is biased toward runs with less depletion). Furthermore, because the Kirkwood gaps have yet to fully shape the radial structure of the MB at t = 200 Myr, we need not expect the shapes (in particular the intergap peaks) of the distributions to match exactly (Morbidelli et al. 1995; Petit et al. 2001; O'Brien et al. 2007; Deienno et al. 2016, 2018).

3.3. Orbital Structure

Most studies of the MB's detailed structure attempt to match a simulated belt's a/e and a/i distributions with that of the actual MB. The actual MB orbital distribution is well characterized in a/e space by a uniform distribution of excited orbits that inhabit the entire allowed parameter space (that is to say all excited orbits that are not planet crossing). However, the present-day MB's a/i distribution is peculiar in that it is largely devoid of high inclination asteroids above the ν₆ secular resonance. The inability to match the ratio of simulated asteroids in the inner MB above to below this resonance with that of the solar system (often cited as ∼0.08; calculated by taking the 100 biggest asteroids with a < 2.8 au) is a common pitfall of many studies of MB dynamics (Walsh & Morbidelli 2011; Deienno et al. 2016). However, the MB's radial mass distribution profile is far from uniform because half of its mass is concentrated in just four asteroids. From the standpoint of mass distribution, the MB's structure seems to be better characterized by just a handful of massive bodies embedded in a sea of objects several orders of magnitude smaller (many of which are likely the products of collisional grinding; Bottke et al. 2005a, 2015; Dermott et al. 2018). Indeed, the ν₆ ratio calculated for all known asteroids with a < 2.8 au with respect to number of asteroids (0.04) is a significantly smaller than when it is calculated with respect to mass (0.14). The reason for this is that the third largest asteroid, Pallas, has an inclination above the ν₆ resonance (33°).

In the sense that the AB's structure is dominated by a small number of asteroids, the third largest of which is above the ν₆ resonance, our simulations are largely successful. Because our initial asteroids have masses within an order of magnitude or so of the actual larger asteroids in the MB, it is encouraging that most of our simulations deplete the MB down to just a few of such objects. Furthermore, the MB structure is also dominated by multiple large collisional families that formed after the epoch of terrestrial planet formation (Bottke et al. 2015; Dermott et al. 2018). Thus, it is supremely difficult to correlate the current belt's orbital and size distributions with those at t = 200 Myr (Bottke et al. 2005a).

Nevertheless, our simulated belt's eccentricity and inclination distributions are fairly similar to the real distribution. In Figures 5–7, we compare the orbital distributions for all bright (H < 9.7) asteroids larger than 50 km with our results from run 1, and all of our interaction runs where Jupiter and Saturn's period ratio finished below 2.8. We select run 1 in this figure because it finished with the most MB objects of any interaction run with P_S/P_J < 2.8 and N_GP = 4. Thus, we provide an example of a single successful distribution without coadding results. However, this plot should be taken in context with the fact that run 1 finishes with a low AMD_JS,SS, and significantly underexcited e₅₅ and e₅₆ amplitudes. Furthermore, the coadded plot is biased toward similar simulations with lower depletion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In keeping with previous works, our simulations consistently struggle to replicate the modern AB's structure about the ν₆ secular resonance. Walsh & Morbidelli (2011) reported a ratio of ∼5.2 in a smooth migration scenario. Deienno et al. (2016) evaluated a post-Grand Tack version of the late Nice Model and determined a ratio ∼1.2 (0.10 when removing post-Grand Tack asteroids with i > 20°). Studies of an early instability have shown promise in replicating some aspects of the belt's structure, but still fall short when it comes to the ν₆ ratio. Clement et al. (2018) coadded the remaining asteroids from successful simulations (P_S/P_J < 2.8) and found the ratio to be ∼0.73, but the value was much worse when only considering the inner MB asteroids with a < 2.5 au (∼2.2). In a similar study, Deienno et al. (2018) considered a Jumping Jupiter style instability and found the ratio to be ∼1.3 (0.07 when removing post-instability asteroids with i > 20°). In general, our simulations (Table 2) finish with ratios of similar order to those found in previous early instability studies. Our best ratio occurred in run 1 (0.31). Because the sculpting mechanisms in each instability are unique, it is difficult to extrapolate a correlation between depletion and the final ν₆ ratio.

Two of our runs finished with the entire inner MB population's inclinations above the ν₆ resonance. This occurs when the location of the ν₆ resonance drags more smoothly through the inner MB. Instabilities where the migration of the giant planet's orbits and respective resonances is more chaotic are less likely to excite the entire region uniformly, and they finish with better ratios. This is also the case in the nominal Jumping Jupiter style instability (Brasser et al. 2009; Nesvorný et al. 2013; Deienno et al. 2018). In either scenario, the location of the resonance does not spend a significant amount of time in any one location during the planetary instability.

In spite of several moderately successful runs (e.g., runs 1, 8, 1a, 4a, and 6a) that finish with ν₆ ratios less than unity, reproducing the inner MB's a/i structure is still an outstanding problem for the early instability model. However, the MB's orbital distribution has evolved as the result of family-forming events over the past ∼4 Gyr (Bottke et al. 2015). In fact, Dermott et al. (2018) estimated that ∼85% of inner MB asteroids originated from one of just five original collisional families (Flora, Vesta, Nysa, Polana, and Eulalia). Thus, it is possible that the primordial AB contained very few objects. If Pallas was one of this primordial population of a few, "large" asteroids, and its high inclination (above the ν₆ resonance) is primordial, our results may not be particularly problematic. Indeed, planet formation simulations that include collisional fragmentation (Clement et al. 2018) indicate that collisional fragments tend to preferentially populate the inclination parameter space below the ν₆ secular resonance. Thus, future work on the topic must use similar GPU-accelerated integration schemes to study how collisional fragmentation shapes the young AB. If it turns out that the fragmentation process dramatically alters the primordial belt's size distribution, it is possible that the dynamical friction generated by a diverse population of fragments could save material from loss during the instability. If that were the case, the depletion values discussed in Section 3.1 might be overestimated.

3.4. Seeded Embryos

3.5. The Formation of Veneneia

In this article we presented 18 N-body simulations of the effect of an early Nice Model instability (Tsiganis et al. 2005; Clement et al. 2018) on a 3000 particle AB. In 10 of our simulations, the entire AB is modeled as fully self-gravitating. Our simulations show that the early instability scenario can deplete a ∼2 M_⊕ primordial AB (commensurate with the minimum mass solar nebula; Weidenschilling 1977) at the ∼99%–99.9% level. The level of depletion experienced by a particular system is related to the unique evolution of the giant planets within the chaotic orbital instability. Because each instability is unique and it is impossible to know the exact evolution of the giant planets in the young solar system's own instability, our work establishes a spectrum outcomes on which we place the solar system. In particular, depletion increases with higher values of AMD_JS and e₅₅ (the amplitude of Jupiter's g₅ eigenmode). Systems that finish with values of AMD_JS and e₅₅ similar to those of the present solar system tend to deplete the AB by two to three orders of magnitude.

Assuming the outer primordial AB was populated with volatile rich C-type asteroids (either because of a primordial snow line or via scattering during the giant planet growth phase as described in Raymond & Izidoro 2017a), our fully self-gravitating simulations are very successful at reproducing the observed distribution of S- and C-types. Contrarily, radial mixing is not as strong in our simulations that only consider asteroids as test particles. Furthermore, these test-particle-only simulations typically overerode the inner MB region (2.0 < a < 2.5 au). Dynamical friction between asteroids helps to deexcite orbits in the region, and prevent this erosion in our fully self-interacting systems.

In general, our systems adequately replicate the eccentricity and inclination distributions of the modern AB. However, consistent with various previous studies (Walsh & Morbidelli 2011; Deienno et al. 2016; Clement et al. 2018; Deienno et al. 2018), our simulations often fail to replicate the AB's population about the ν₆ secular resonance. However, the MB size distribution has evolved significantly since the epoch of terrestrial planet formation (Bottke et al. 2005a), and its orbital distribution has likely been altered via the formation of collisional families (Dermott et al. 2018). Because the MB's mass profile is dominated by just a few asteroids, one of which (Pallas) is on a highly inclined orbit above the ν₆ resonance, our systems still broadly replicate the distribution of the most massive asteroids.

Our results should be taken in the appropriate context given the sensitivity of the AB's depletion and final structure to the particular dynamics of the giant planet instability. It is impossible to know the exact evolutionary path followed by the giant planets during the Nice Model instability. Furthermore, none of our final giant planet systems is an exact match to the modern solar system. Therefore, our study can only correlate specific dynamical qualities of the present solar system with those of our simulations. Nevertheless, our fully self-gravitating ABs are largely successful at replicating the modern belt for several reasons. These include a number of systems with low ratios of asteroids above to below the ν₆ resonance, accurate radial distributions of S- and C-types, and relatively high levels of MB depletion.

Although our work demonstrates that GPU-accelerated N-body simulations offer an unprecedented opportunity to explore higher-N systems, each of our simulations still requires large amounts of computing resources. Specifically, each simulation in this article required the same number of node hours to complete as 32 conventional (CPU only) integrations of terrestrial planet formation using standard initial conditions (Chambers 2001; Clement et al. 2018). In spite of generating a suite of simulated ABs sculpted by giant planet instabilities, the stochasticity of instabilities and the small number of surviving asteroids prevents us from developing detailed predictions of AB structure as a function of instability characteristics. Doing so will likely require still more extensive sets of instability simulations.

This material is based on research supported by the Chateaubriand Fellowship of the Office for Science and Technology of the Embassy of France in the United States. We thank Rogerio Deienno and Bill Bottke for thoughtful comments and insight. M.S.C. and N.A.K. thank the National Science Foundation for support under award No. AST-1615975. S.N.R. thanks the Agence Nationale pour la Recherche for support via grant No. ANR-13-BS05-0003-002 (grant MOJO) and acknowledges NASA Astrobiology Institutes Virtual Planetary Laboratory Lead Team, funded via the NASA Astrobiology Institute under solicitation No. NNH12ZDA002C and cooperative agreement No. NNA13AA93A. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards No. OCI-0725070 and No. ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications (Bode et al. 2013; Kramer et al. 2015).