The Role of Early Giant-planet Instability in Terrestrial Planet Formation

Our N-body simulations start at the time of the gas disk dispersal, t₀. In the standard setup, there are N _pp  = 100 protoplanets and N_pm = 1000 planetesimals, each representing the total mass of 2 M_Earth (i.e., the total mass of solids is M_solid = 4 M_Earth; each protoplanet has the mass M_pp ≃ 1.6 M_Moon and each planetesimal has the mass M_pm ≃ 13 M_Ceres, where M_Earth = 5.97 × 10²⁷ g, M_Moon = 7.35 × 10²⁵ g, and M_Ceres = 9.38 × 10²³ g). The bulk density of protoplanets and planetesimals is set to ρ = 3 g cm⁻³. The initial eccentricities and initial inclinations are chosen to follow the Rayleigh distributions with σ_e = 0.005 and σ_i  = 0.0025. The protoplanets and planetesimals are distributed with the surface density profile Σ(r) = r^−γ , where r is the heliocentric distance, in a disk between r_in and r_out, with γ = 1 in the standard case.

We test three different disk setups. In the first case, both the protoplanet and planetesimal disks have r_in = 0.3 au and r_out = 4 au. This means that ∼1 M_Earth in protoplanets ³ is placed into the asteroid belt region. In the second case, the protoplanets are distributed in a narrow annulus between r_in = 0.7 au and r_out = 1 au, whereas the planetesimal disk follows the same distribution as in the first setup (i.e., r_in = 0.3 au and r_out = 4 au). This would approximate the state of the inner solar system at t₀ in the Grand Tack model (Walsh et al. 2011), except that the planetesimals in the asteroid belt region start on low-eccentricity/-inclination orbits (see Section 3.4 for simulations with excited asteroid orbits). In the third case, protoplanets and planetesimals occupy separate radial intervals: r_in = 0.3 au and r_out = 1.5 au for protoplanets, and r_in = 1.5 au and r_out = 4 au for planetesimals. This case is motivated by the results of Levison et al. (2015), who showed that the protoplanet and planetesimal disks may not overlap at t₀ (also see Walsh & Levison 2019 and Clement et al. 2020a). ⁴ The three cases discussed above are referred to as C98 (after Chambers & Wetherill 1998), W11 (after Walsh et al. 2011), and L15 (after Levison et al. 2015).

We consider six different dynamical models for each protoplanet/planetesimal distribution. In the first model, there are no giant planets whatsoever. This is clearly unrealistic, because the giant planets formed before the terrestrial planets, but it is useful to consider this case as a reference with no external influence on the terrestrial planet formation. In the second model, we start the simulations with the four giant planets (Jupiter to Neptune) on their current orbits. Note that this would be equivalent to assuming that the giant planets evolved to their current orbits before t₀ and there was no planetesimal-driven migration and dynamical instability after t₀. Again, this "static" case is useful to consider as a benchmark for models with planetary migration/instability (Section 2.2). In the third model, we extract the initial configuration of giant planets from the instability Case 1 reported in Nesvorný et al. (2013). In Case 1, five giant planets (Jupiter, Saturn, and three ice giants) start on orbits in a compact resonant chain with the period ratios following the sequence 3:2, 3:2, 2:1, and 3:2 (i.e., Jupiter and Saturn are in 3:2 mean motion resonance). Because we do not consider the outer planetesimal disk, however, the giant planets do not migrate in this model. This would correspond to a case where the terrestrial planet formation was completed prior to the migration/instability of the giant planets. The effects of the giant-planet migration/instability on a fully formed terrestrial planet system were studied in Brasser et al. (2009, 2013), Agnor & Lin (2012), Kaib & Chambers (2016), and Roig et al. (2016). The three models described above are denoted nojov, jovend, and jovini, where "jov" stands for the Jovian planets. The three remaining models are described below.

We take advantage of our previous simulations of the giant-planet migration/instability (e.g., Nesvorný & Morbidelli 2012, hereafter NM12; Deienno et al. 2017) and select three cases. The three models started with fully formed giant planets at t₀. The giant planets were assumed to have orbits in mutual resonances that have been presumably established during the previous stage of convergent migration in a gas disk (Masset & Snellgrove 2001; Pierens & Nelson 2008; Pierens & Raymond 2011). In all three cases, five giant planets began in the resonant chain (3:2, 3:2, 2:1, 3:2), a configuration that was previously shown to work best to generate plausible dynamical histories (NM12; Deienno et al. 2017). The additional ice giant had mass comparable to that of Uranus/Neptune. A massive outer disk (mass M_disk) of planetesimals was placed beyond the outermost ice giant. The dynamical evolution of the giant planets and planetesimals was then tracked using the N-body integrator known as SyMBA (Duncan et al. 1998).

The properties of two of our instability models, hereafter case1 and case2, are illustrated in Figures 1 and 2. The case1 model is identical to Case 1 used in Nesvorný et al. (2013, 2014, 2017), Deienno et al. (2014), and Roig & Nesvorný (2015). This model satisfies many solar system constraints, including the present orbits of the giant planets themselves, the number and orbital distribution of Jupiter Trojans and irregular moons of the giant planets, the dynamical structure of the asteroid and Kuiper belts, etc. (see Nesvorný 2018 for a review). Moreover, Roig et al. (2016) demonstrated, assuming that the terrestrial planets were already in place when the instability happened, that the case1 model could explain the excited orbit of Mercury and the AMD of the terrestrial planet system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our case2 is similar to case1 but the jump of Jupiter (and Saturn) is not as large as in case1 (Figure 2). In case2, the terrestrial planet region is expected to experience a stronger coupling to the giant planets via sweeping secular resonances (e.g., Brasser et al. 2009; Agnor & Lin 2012). Here we follow the general suggestion of Clement et al. (2018) that a broader range of instabilities should be considered, but we do not want to distance ourselves too much from the class of instabilities that work reasonably well for other solar system constraints (Nesvorný 2018). For the third instability model considered here, case3, we refer the reader to Figures 1 and 3 in Deienno et al. (2018), who showed that a more complex dynamical evolution of Jupiter and Saturn during the instability could explain the excited orbits of main belt asteroids (even if asteroids had e ∼ i ∼ 0 at t₀). This contrasts with case1, where the instability itself does not sufficiently excite asteroid orbits (Roig & Nesvorný 2015; Nesvorný et al. 2017). In case1, the orbital excitation of main belt asteroids would have to predate the onset of giant-planet migration/instability.

The final orbital period ratio of Jupiter and Saturn is P_Sat/P_Jup ≃ 2.3–2.5 in the three instability models considered here (e.g., Figures 1 and 2; P_S/P_J = 2.49 in the real solar system). This can be compared to the very broad variation of final P_S/P_J in the models considered in Clement et al. (2018). The final amplitude of Jupiter's eccentric mode is E₅₅ ≃ 0.025–0.042 in our three models (highest in case2, lowest in case1). These values are somewhat lower than E₅₅ ≃ 0.044 of the real Jupiter's orbit.

In all three dynamical models considered here, the early solar system is assumed to have five giant planets (Jupiter, Saturn, and three ice giants). This is because NM12 showed that various constraints, such as the final orbits of outer planets, can most easily be satisfied if the solar system started with five giant planets and one ice giant was ejected during the instability (Nesvorný 2011; Batygin et al. 2012; Deienno et al. 2017). The case with four initial planets requires a massive planetesimal disk to avoid losing a planet, but the massive disk also tends to produce strong dynamical damping and long-range residual migration of Jupiter and Saturn that violate constraints. It is therefore difficult to obtain a plausible evolution history starting with four planets.

A shared property of the selected models is that Jupiter and Saturn undergo a series of planetary encounters with the ejected ice giant. As a result of these encounters, the semimajor axes of Jupiter and Saturn evolve in discrete steps. While the semimajor axis can decrease or increase during one encounter, depending on the encounter geometry, the general trend is such that Jupiter moves inward, i.e., to shorter orbital periods (by scattering the ice giant outward), and Saturn moves outward, i.e., to longer orbital periods (by scattering the ice giant inward). This process leads to the dynamical evolution known as the jumping-Jupiter model, ⁵ where the orbital period ratio of Jupiter and Saturn, P_Sat/P_Jup, changes as shown in Figures 1(b) and 2(b). In the jumping-Jupiter model, the principal coupling between the inner and outer solar system bodies occurs via secular resonances (e.g., g₁ = g₅, g₄ = g₅, g = g₅, g = g₆, s = s₆, where g_j , s_j , s, and g are the fundamental precession frequencies of planetary and asteroid orbits). The secular resonances can cause strong orbital excitation and inhibit the accretional growth at specific radial distances from the Sun (Clement et al. 2018).

NM12 and Deienno et al. (2018) simulations were performed using the symplectic N-body codes SyMBA (Duncan et al. 1998) and Mercury (Chambers 2001), respectively. Here we want to reproduce exactly the orbital evolution of the giant planets in the selected runs and, in addition, include protoplanets and planetesimals in the inner solar system (Section 2.1). This is not easy because the protoplanets and planetesimals need to carry mass, but if they were included as massive bodies in SyMBA or Mercury, they would affect the orbital evolution of giant planets and produce completely different orbital histories. Clement et al. (2018) approached this problem by performing a large suite of N-body simulations where the terrestrial protoplanets/planetesimals, giant planets, and outer disk planetesimals all had mass and gravitationally interacted among themselves. This led to a large collection of results. The cases that at least roughly matched various constraints were then analyzed in more detail. This approach is not possible here where we want to be more selective (Section 2.2). To deal with this issue, we develop a new method that allows us to track the giant-planet orbits taken from the original SyMBA/Mercury simulations and include terrestrial protoplanets/planetesimals that carry mass, interact between themselves, but do not affect the giant planets. This works as follows.

The planetary positions and velocities in the selected instability simulations were recorded with a 1 yr cadence in Nesvorný et al. (2013) and Deienno et al. (2018). Our modified code iSyMBA, where "i" stands for interpolation, then reads the planetary orbits from a file, and interpolates them to any required time subsampling (generally 0.01–0.02 yr, which is the integration time step used here for the terrestrial planets). The interpolation is done in Cartesian coordinates. First, the giant planets are forward-propagated on the ideal Keplerian orbits starting from the positions and velocities recorded by SyMBA at the beginning of each 1 yr interval. Second, the SyMBA position and velocities at the end of each 1 yr interval are propagated backward (again on the ideal Keplerian orbits; planetary perturbations switched off). We then calculate the weighted mean of these two Keplerian trajectories for each planet such that progressively more (less) weight is given to the backward (forward) trajectory as time approaches the end of the 1 yr interval. We verified that this interpolation method produces insignificant errors.

The iSyMBA code accommodates three types of bodies: (1) interpolated giant planets, (2) fully interacting terrestrial protoplanets, and (3) planetesimals that carry mass, affect the terrestrial protoplanets, but do not interact among themselves. The original symplectic structure of SyMBA in Poincaré canonical variables (heliocentric positions, center-of-mass linear momenta) is LD-K-D-K-LD, where LD stands for a linear drift due to the Sun's linear momentum, K for a kick that accounts for planetary perturbations, and D for a drift on the Keplerian orbit. This was modified to GK-LD-TK-D-TK-LD-GK, where GK indicates the kick corresponding to the giant-planet perturbations and TK is the gravitational interaction of the terrestrial planets and planetesimals. The tricky part in iSyMBA is how to define the center of mass in different parts of the algorithm and how to do that consistently for the terrestrial protoplanets/planetesimals and interpolated giant planets. See F. V. Roig et al. (2020, in preparation) for a detailed description of the code. iSyMBA was parallelized with OpenMP directives. Through testing, we established that the code optimally runs, for the setup described in Section 2.1, when ∼8 threads are used (performance saturates for >10 threads). The typical speedup is ∼3–5.

The start of our numerical integrations with iSyMBA (simulated time t = 0) is assumed to coincide with t₀. The simulations are run to t = 0.8–30 Myr (Phase 1; 10 Myr in case1, 30 Myr in case2, and 0.8 Myr in case3). The instability is assumed to happen as in the original simulation, roughly at t_inst = 6 Myr for case1 and t_inst = 20 Myr for case2 (Figures 1 and 2). This is to respect the fact that the instability is triggered during planetary migration, and there is always certain delay between the gas disk dispersal and t_inst. In case3, taken from Deienno et al. (2018), we test t_inst = 0.6 Myr.

The Phase 1 integrations end at t = t₁ ≃ 1.3–1.7 t_inst. Much longer integrations are difficult to achieve with iSyMBA, because the interpolation method has large requirements on the computer memory and disk storage. As the giant planets are orbitally decoupled by the end of Phase 1, however, the continued integrations for t > t₁ (Phase 2) do not need to deal with the stochastic outcomes of encounters between the giant planets; it therefore suffices to use the standard SyMBA code. Phase 2 simulations are run to t₂ = 200 Myr, at which point the terrestrial planet accretion is complete. For the static cases nojov, jovini, and jovend, where no interpolation is needed, we use the standard SyMBA for the whole length of simulations.

Table 1 lists all simulations performed in this work. We have three static models plus three instability models, and three different distributions of protoplanets/planetesimals in the inner solar system. This represents 6 × 3 = 18 cases in total. For each case, we perform 10 simulations where different random seeds are used to generate slightly different initial conditions. This represents 180 base simulations in total. In addition, we set up four additional cases that have the best chance to match constraints and perform 100 additional simulations for each to thoroughly characterize the statistical variability of the results (Sections 3.4 and 3.5). The simulations were done on the NASA Pleiades supercomputer. Each simulation was completed over ∼3 weeks on eight Broadwell cores.

All collisions between protoplanets and planetesimals are assumed to end in inelastic mergers of bodies. Whether the debris generation in giant impacts is important for the terrestrial planet formation is still a matter of debate. On one hand, Mercury's mantle may have been stripped by a hit-and-run collision with a large body (Asphaug & Reufer 2014), and the Moon likely formed from a circumplanetary debris disk generated by a large impact on proto-Earth (see R. M. Canup et al. 2020, in preparation for a review). On the other hand, Deienno et al. (2019) found no change in the terrestrial planet formation results when different energy dissipation schemes were used to deal with collisions. Given that the main goal of our work is to test the effect of giant planets on the terrestrial planet formation (and on the asteroid belt), we prefer to keep things simple and leave a more complex treatment of collisions for future work. The integrations reported here use a 5 day time step, but we verified in several cases that the results are statistically indistinguishable when a 3.5 day time step is used (the orbital period of Mercury is about 88 days, ∼18 and ∼26 times longer than our time steps).

Different measures are used here to to characterize the overall success of simulations. To quantify the radial distribution of planetary mass, Chambers (2001) defined the radial mass concentration, or RMC, as

$$\begin{eqnarray}&&{S}_{{\rm{c}}}=\max \left(\displaystyle \frac{{\displaystyle \sum }_{j}{m}_{j}}{{\displaystyle \sum }_{j}{m}_{j}{\left[{\mathrm{log}}_{10}(a/{a}_{j})\right]}^{2}}\right),\end{eqnarray} \tag{ 1 }$$

where m_j and a_j are the mass and semimajor axis of planet j, the sum goes over all planets in the inner solar system (giant planets not included), and the maximum is taken over a; S_c = 89.9 in the real solar system. The AMD is defined as

$$\begin{eqnarray}&&{S}_{{\rm{d}}}=\displaystyle \frac{{\displaystyle \sum }_{j}{m}_{j}\sqrt{{a}_{j}}(1-\sqrt{1-{e}_{j}^{2}}\cos {i}_{j})}{{\displaystyle \sum }_{j}{m}_{j}\sqrt{{a}_{j}}},\end{eqnarray} \tag{ 2 }$$

where e_j and i_j are planetary eccentricities and inclinations; S_d = 0.0018 in the real solar system. The advantage of these parameters is their objective mathematical formulation and relationship to the dynamical stability of planetary systems (e.g., Petit et al. 2018). Note, however, that they can be computed, and are often reported as such, for simulations that do not produce meaningful results (when the number and/or masses of the final terrestrial planets are clearly wrong). This can bias the statistical interpretation of results because correlations between different measures exist (e.g., Deienno et al. 2019)

For these reasons, we report the S_c and S_d statistics only for the simulations that produced "good" Venus/Earth planets. To this end, we collect all planets with mass M > 1/3 M_Earth between 0.5 and 1.2 au and call them the Venus/Earth analogs. The successful simulations are required to have exactly two Venus/Earth analogs (denoted by indices 1 and 2 below; f_g in Table 1 reports the fraction of these "good" simulations for each model). We also compute the radial separation between good Venus/Earth planets as Δa = a₂ − a₁ (this is a complementary parameter to S_c). In the real solar system, Venus and Earth have Δa = 0.277 au. We also determine, as a complement to S_d, the mean eccentricities/inclinations of Venus/Earth analogs: 〈e, i〉 = (e₁ + i₁ + e₂ + i₂)/4. The real Earth and Venus have 〈e, i〉 = 0.0274.

The values of Δa and 〈e, i〉 measure the radial separation and orbital excitation of the two large planets that form at 0.5–1.2 au. The S_c and S_d parameters are related to that but, unlike Δa and 〈e, i〉, they are influenced, for example, by whether Mars/Mercury form in the simulations, and if so, whether they at least approximately have the right mass (e.g., massive Mars/Mercury on perfect orbits would increase AMD/decrease RMC not because the orbits are too excited/spread but because the planets are too massive). In this work, the planets with a < 0.5 au are called Mercury and the planets with 1.2 < a < 1.8 au are called Mars; they can have any mass M > 0.01 M_Earth with our setup. We also define strict Mercury and strict Mars as a < 0.5 au and 0.025 < M < 0.2 M_Earth for Mercury, and 1.2 < a < 1.8 au and 0.05 < M < 0.2 M_Earth for Mars. We allow for a larger mass range in the case of Mercury to leave space for the possibility of hit-and-run collisions (Chambers 2013; Asphaug & Reufer 2014; Clement et al. 2019a). We pay close attention to the mass of Mars planets obtained in the simulations, because small Mars is thought to be an important diagnostic for the early evolution of the solar system (e.g., Chambers & Wetherill 1998; Chambers 2001; Hansen 2009; Walsh et al. 2011; Walsh & Morbidelli 2011; Jacobson & Morbidelli 2014; Clement et al. 2018).

3.1.1. Model C98

We first briefly comment on C98/nojov to set the stage for the presentation of more successful models. With an extended disk of protoplanets/planetesimals and no giant planets included in simulations, planets with M > 0.1 M_Earth form over the entire radial range 0.3 < r < 4 au (Figure 3). This is clearly unsatisfactory because Mars with a = 1.52 au is the outermost terrestrial planet in the real solar system, and there are no planetary-size bodies in the asteroid belt (2.2 < r < 3.5 au). The masses of the Venus/Earth analogs are somewhat small, indicating that the total mass in solids placed into the inner solar system, M_solid = 4 M_Earth, was not large enough in this model. We increase the total mass of protoplanets/planetesimals in the models described in Section 3.4. Also, the planets that form in C98/nojov at r ∼ 1.5 au are excessively massive, typically ∼2–7 times more massive than the real Mars.

Zoom In Zoom Out Reset image size

The results improve when the giant planets are included in the simulations. With the giant planets on their current orbits (model C98/jovend; Figure 4), ⁶ there is a sharp transition from r < 1.5 au, where massive planets with M > 0.2 M_Earth grow, to r > 1.5 au, where the growth is modest and the final planets have masses 0.02 < M < 0.2 M_Earth. A similar result has been reported in Raymond et al. (2009, their case EJS). The transition is most likely related to the combined effect of protoplanet scattering and secular resonances ν₅, ν₆, and ν₁₆ that excite orbits near a = 2 au. Beyond 2 au, protoplanet scattering, orbital resonances with Jupiter, and the proximity to Jupiter eliminate most protoplanets and planetesimals from the asteroid belt region (f_ast ≃ 0.02; Table 1).

Zoom In Zoom Out Reset image size

In the C98/jovend model, the transition at r = 1.5 au is too sharp with planets near 1.5 au either having M > 0.3 M_Earth or M < 0.05 M_Earth. This leaves Mars in no man's land. In fact, there are only two relatively good Mars planets (small yellow and blue dots in the top panel of Figure 4) in 10 simulations, which is not very satisfactory. The results are better when the giant-planet instability is taken into account (C98/case1 in the bottom panel of Figure 4). In this case, the planetary mass shows a more gradual decrease with the radial distance: from M ≳ 0.5 M_Earth at r = 1 au, to M ∼ 0.1-0.3 M_Earth at r = 1.5 au, to M ∼ 0.02-0.1 M_Earth at r = 2 au. The mass scatter, however, is quite large and only 2 of our 10 C98/case1 simulations produced acceptable matches to the inner solar system (two planets with M > 1/3 M_Earth and 0.5 < a < 1.2 au, Δa < 0.415 au—no more than 1.5 times the real value—and 〈e, i〉 < 0.0548—less than double the real value—and the simulated Mars mass within a factor of 2 to the real mass). The results for C98/case3 are similar (Table 1). We therefore do not see much dependence on the orbital evolution of giant planets during the instability (at least for the class of instabilities investigated here; recall that the case1 and case3 instabilities are quite different, Section 2.2). Interestingly, the C98/case2 model does not produce any strict Mars analogs (planets near 1.5 au are too massive; Table 1). This is probably a consequence of the fact that the giant-planet instability happens relatively late in case2 (t_inst ≃ 20 Myr), and the planets near r = 1.5 au accrete too much mass before this region is affected by the giant planets (Clement et al. 2018).

In summary, the primary effect on the terrestrial planet formation is the presence of the giant planets themselves (if they reach their current orbits early) that creates the transition at 1.5 au. The secondary effect is the evolving orbits of the giant planets during the instability, which makes the transition more gradual and allows for the formation of a small (but not too small) Mars at 1.5 au. For that, however, the giant-planet instability must happen relatively early (before 10 Myr after t₀ or so; Clement et al. 2018). If it would happen late (t > 10 Myr), Mars's growth would not be affected by it.

Interestingly, there are some reasonable Mercury analogs in the C98/case1 model, where 4 out of our 10 simulations produced planets with M  ≲  0.2 M_Earth and a ≃ 0.4 au (bottom panel of Figure 4). Only one such planet was obtained in 10 simulations of the C98/jovend model. This indicates that having small Mercuries in C98/case1 is probably not a consequence of the inner edge of the initial disk (set at r_in = 0.3 au to limit the CPU cost). The orbits of the Mercury analogs are excited to eccentricities e ≃ 0–0.25 and inclinations i ≃ 5°–15°. This is reminiscent of the results of Roig et al. (2016), who showed that the orbit of fully formed Mercury could be excited by sweeping secular resonances (g₁ = g₅ and s₁ = s₇), if giant-planet migration/instability happened relatively late. We therefore believe that the same resonances act here, during the terrestrial planet formation, to excite orbits at r ≃ 0.4 au and stall Mercury's growth.

We find that N_spl = 1.5 and 0.4 planetary-mass bodies (M > 0.01 M_Earth, average per good simulation, Table 1, index "spl" stands for "stranded planets") end up in the asteroid belt in the C98/jovend and C98/case1 models, respectively. These results are therefore more plausible than C98/nojov where N_spl = 4.5. Planetary-mass bodies in the asteroid belt would have a huge influence on the depletion, mixing, and orbital excitation of asteroids (Section 3.2). The fact that N_spl ∼ 1 in the C98/jovend and C98/case1 models is still troublesome, but this can potentially be resolved if (i) the long-term instabilities remove additional bodies from the asteroid belt on gigayear timescales, and/or (ii) the original protoplanet/planetesimal disk in the inner solar system had a steep radial profile (Izidoro et al. 2015; Levison et al. 2015). Here we use γ = 1 (Section 2.1). Simulations with γ = 1.5 are described in Section 3.4. We proceed by discussing the results of the W11 and L15 setups where N_spl ∼ 0 (except if the nojov or jovini models are considered; Table 1).

3.1.2. Models W11 and L15

Figure 5 compares the results for the W11/nojov and W11/jovend models. The protoplanets with M = 0.02 M_Earth start in a narrow annulus in these models (Hansen 2009; Walsh et al. 2011; Walsh & Levison 2016), and the radial mass profile of planets at t = 200 Myr is peaked near the original annulus location. The radial distribution of final planets (0.5 < r < 1.2 au at t = 200 Myr), however, is wider than that of the original annulus (0.7 < r < 1 au at t₀). This happens because there is no active mechanism that would confine protoplanets to the annulus. Instead, as protoplanets grow, they scatter each other, and the general tendency of this relaxation process is to spread bodies over a larger range of radial distances. An undesired consequence of this process is that our simulations typically lead to a relatively large separation between Venus and Earth (Section 3.3.1). This problem could potentially be resolved if the initial protoplanets were more massive (Jacobson & Morbidelli 2014). We test this possibility in Section 3.5.

Zoom In Zoom Out Reset image size

The effect of giant planets on terrestrial planet formation is not as dramatic in the narrow-annulus model (Figure 5) as it is in the C98 model (Figures 3 and 4). There are far fewer planets stranded in the asteroid belt in W11/nojov (N_spl = 1.4) than in C98/nojov (N_spl = 4.5) and including the giant planets in W11 reduces N_spl to zero (except for jovini where N_spl = 0.6). The results for W11/jovend (bottom panel of Figure 5) are relatively good. The Mars planets at 1.2–1.8 au have M = 0.02–0.3 M_Earth, except for one case (large green dot) where the massive planet with M ≃ 0.6 M_Earth formed near 1.42 au.

Model Mercuries with M < 0.2 M_Earth and 0.3 < r < 0.5 au form in both W11/nojov and W11/jovend. In this case, the primary reason for a small planet at r ∼ 0.4 au is the narrow-annulus setup and, consequently, the lack of protoplanets initially available for accretion at r < 0.7 au. Model Mercuries in W11 are protoplanets scattered from r > 0.7 au to r ∼ 0.4 au, where their orbits become stabilized by dynamical friction from planetesimals. The model planets at r ∼ 0.4 au are often a factor of ∼1.5–4 more massive than Mercury, but this may not be an insurmountable problem because Mercury's mass could have been reduced by violent hit-and-run collisions (Chambers 2013; Asphaug & Reufer 2014; Clement et al. 2019a), a process not accounted for here.

The results with and without the instability, W11/case1 in the top panel of Figure 6 and W11/jovend in the bottom panel of Figure 5, are similar in that they show a relatively strong concentration of massive planets near the original annulus location. There are interesting differences as well. Notably, the W11/case1 model produced fewer Mars/Mercury planets than the W11/jovend model. This is probably a consequence of the coupled effect of the narrow annulus and sweeping secular resonances in W11/case1. The results for W11/case3 show the same trend but the statistics of Mars/Mercury planets is better in the W11/case2 model. This probably indicates some mild dependence on the instability case.

Zoom In Zoom Out Reset image size

The L15/case1 model (bottom panel of Figure 6) does not offer any obvious advantage over the W11/case1 model. The final radial mass distribution of planets in L15/case1 is relatively broad, probably because the original protoplanet disk has an edge at r_out = 1.5 au in L15, and there are fewer Mercuries than in W11, probably because there are no planetesimals (i.e., no dynamical friction) to stabilize protoplanets scattered from the original annulus location to r ∼ 0.4 au. We see no important differences between the results obtained for different instability cases, and thus, for the sake of brevity, we do not discuss the case2 and case3 models in detail here (see Section 3.3 for a more detailed comparison).

3.1.3. Comparison of the C98, W11, and L15 Models

The orbits of the terrestrial planets, and their orbital eccentricities and inclinations in particular, are important constraints on the planet/planetesimal formation in the inner solar system. These constraints are best evaluated together with the dynamical structure of asteroid belt, because (i) the main belt asteroids are the only leftover planetesimals that survived in the inner solar system, and (ii) their overall mass and orbital distribution have been influenced both by protoplanet formation processes and the giant-planet instability.

Figure 7 shows the real orbits of four terrestrial planets and all known asteroids with diameters D > 30 km. We only plot the large asteroids because they are less affected by radiation effects and collisions than the small ones; their orbital distribution is thus more diagnostic for processes in the early solar system. Mercury's relatively large orbital eccentricity (mean e ≃ 0.17) and inclination (mean i ≃ 7°) are puzzling aspects of the inner solar system. In the narrow-annulus model, proto-Mercury is scattered from r = 0.7–1 au and stabilized at its current orbit with a ≃ 0.38 au by collisions and/or dynamical friction. In Roig et al. (2016), Mercury can start with e ∼ i ∼ 0, and its orbit becomes excited during the giant-planet instability.

Zoom In Zoom Out Reset image size

The main asteroid belt spans r ≃ 2.2–3.5 au (here we ignore Hildas in the 3:2 resonance with Jupiter, a ≃ 4 au; see Roig & Nesvorný 2015). The orbital eccentricities of asteroids are relatively large—the orbits practically fill the whole region that is dynamically stable. The eccentricities could have been low at t₀ and increased after t₀, for example, during the instability (Deienno et al. 2018; Clement et al. 2019b) and/or by protoplanet scattering (O'Brien et al. 2007; Clement et al. 2019b). The eccentricity distribution could have also been quite broad at t₀, as in the Grand Tack model (Walsh et al. 2011; also see Izidoro et al. 2016), and evolved to a narrower distribution as the high-e orbits were removed over 4.5 Gyr (Roig & Nesvorný 2015; Deienno et al. 2016).

The inclination distribution of asteroids represents a more meaningful constraint. The high-i orbits with i ≃ 20°–30° are generally stable and yet the number of asteroids falls off for i > 20° (Figure 7(b)). This means that a relatively small population of asteroids had i ≃ 20°–30° in the first place, which is an interesting constraint on the Grand Tack model (Deienno et al. 2016) and dynamical evolution of the giant planets after t₀ (Morbidelli et al. 2010; Roig & Nesvorný 2015; Clement et al. 2020b).

Figures 8 and 9 show the orbits of planets and asteroids for the C98 model. Both in the jovend and case1 models, the small planets at r < 0.5 au and r ∼ 1.5 au tend to have more excited orbits than the massive planets at r  ≃  0.5–1.2 au. Obtaining good orbits for Mercury and Mars may thus not be difficult (we would need a larger statistic to be able to evaluate this in detail). Instead, the main problem is that the orbital eccentricities and inclinations of massive planets at 0.5–1.2 au are generally too large (Section 3.3). This problem—traditionally expressed in terms of the total AMD of the terrestrial planets—received much attention since it was first pointed out in Chambers & Wetherill (1998). Increasing the planetesimal population would help because this enhances the effects of dynamical friction on protoplanets (e.g., O'Brien et al. 2006), but Jacobson & Morbidelli (2014) and Levison et al. (2015) argued for a relatively small mass in planetesimals at r ∼ 1 au. Here we split the initial mass equally in planets and planetesimals.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As for the orbital structure of the asteroid belt, we find that the asteroid orbits are strongly but not excessively excited in the C98 model. The planetesimal population at r = 2.2–3.5 au is depleted by factor of ∼50 in C98/jovend and ∼100 in C98/case1 (Table 1), and there are on average 1.5 and 0.4 planets stranded in the asteroid belt at t = 200 Myr in C98/jovend and C98/case1, respectively. In these models, the excitation and depletion of asteroid orbits are mainly caused by the protoplanets that start in the main belt at t₀. Given that the number of protoplanets surviving in the main belt at t = 200 Myr is generally low (some would be removed during the subsequent evolution, a steeper radial profile would limit their number as well; Section 3.4 and O'Brien et al. 2007), the C98 model cannot be ruled out based on this criterion alone. The maximum depletion factor of ∼140 is obtained for C98/case2. Interestingly, at t = 200 Myr, many leftover planetesimals have planet-crossing orbits with a = 1–2 au and orbital inclinations reaching up to i  ≃  40° (Figures 8 and 9). The potential relevance of this impactor population for the Late Heavy Bombardment (LHB) will be addressed in our future work.

The orbital structure of the terrestrial planet system in the W11 model is similar (Figures 10 and 11), with the eccentricities and inclinations of Venus/Earth planets being again somewhat too large on average. In this case, there are no protoplanets in the asteroid belt, and the excitation of asteroid orbits is therefore due to the giant-planet instability itself and, to a lesser degree, due to gravitational scattering from protoplanets whose orbits temporarily overlap with the main belt. The excitation of asteroid inclinations ends up being inadequate in W11/case1 (Figure 10(b)), which is consistent with the results of Roig & Nesvorný (2015), who used the same instability model. In the W11/case1 model, some process would have to excite asteroid orbits before t₀ (e.g., Walsh et al. 2011). Alternatively, for the excitation to happen after t₀, a stronger instability model would have to be adopted (Clement et al. 2018; Deienno et al. 2018). Here, we illustrate this for W11/case3 (Figure 11), where case3 was taken from Deienno et al. (2018). The results obtained for the L15 model are very similar to those obtained for W11, both in what concerns the excitation of the terrestrial planet orbits, and the depletion and excitation of the asteroid belt (Table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There is substantial freedom in how to communicate the results of the terrestrial planet simulations. For example, the S_c and S_d parameters can be reported for a simulation independently of whether the simulations produced one, two, or many terrestrial planets. While this is still useful to demonstrate different trends, it is not ideal to clump the results with different numbers of the terrestrial planets together because various correlations exist. For example, Deienno et al. (2019) pointed out that S_d is correlated with the final number of planets, and the systems with more planets typically have lower S_d and therefore appear to be more plausible in terms of S_d. But, strictly speaking, the results with more than four terrestrial planets are incorrect, and they should not be included in the S_d (and S_c) statistics in the first place.

As Venus and Earth are the chief contributors to S_c and S_d in the real solar system, we find it useful to identify simulations where some sort of reasonable Venus/Earth analogs form and report the results only for those simulations. Here we collect simulated planets with 0.5 < a < 1.2 au and M > 1/3 M_Earth (Section 2.4), call them Venus/Earth, and only consider the runs where the final number of Venus/Earth planets is equal to two. The parameter f_g (second column in Table 1), where the index "g" stands for "good," reports the fraction of simulations in each case that produced exactly two Venus/Earth-analog planets. For example, in the W11 model, f_g = 0.5–0.8, indicating that a satisfactorily large fraction of simulations ended up with two large planets at 0.5–1.2 au.

3.3.1. Radial Mass Concentration

3.3.2. Angular Momentum Deficit

3.3.3. Δa and 〈e, i〉

3.3.4. Mars/Mercury

3.3.5. Asteroid Belt

Figures 12 and 13 summarize the results of simulations described previously. We highlight the simulations that produced reasonable results: M_Mars = 0.054–0.214 M_Earth (within a factor of 2 to the real mass), Δa = 0.139–0.415 au (no less than a half and no more than 1.5 times the real value) and 〈e, i〉 = 0.0137–0.0548 (within a factor of 2 to the real value). For C98, the best results were obtained for case1 where two of 10 simulations simultaneously satisfied these criteria (red triangles in the figures). For W11, the best results were obtained for jovend where three simulations satisfied the criteria (green octagons), and for case1 and case2 (one successful simulation each). For L15, one simulation of jovend looks good in Figures 12 and 13 (blue octagon), but that one ended with very small Earth/Venus analogs (just above our M > 1/3 M_Earth cutoff).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the statistically small fraction of successful simulations, we select two promising models and perform 100 simulations for each. The first model is a slight modification of W11/case1, where planetesimals with a > 1.5 au are placed on orbits with larger eccentricities and larger inclinations. There are two reasons for this. First, the case1 instability, in the absence of other effects, cannot explain the inclination excitation of main belt asteroids (Figure 10), if i ∼ 0 at t₀. Second, for W11 to be a closer proxy for the Grand Tack, the eccentricities and inclinations of planetesimals would have to be substantial at t₀. Specifically, we use the Rayleigh distributions with σ_e = 0.3 and σ_i  = 0.14 for a > 1.5 au (model W11e where "e" stands for "excited"). The second model is a modification of C98/case1. We use a slightly steeper radial profile of protoplanets/planetesimals, γ = 1.5, in an attempt to reduce the number of planets stranded in the asteroid belt (this model is called C98s in the following, where "s" stands for "steeper" profile). We also increase the total initial mass to 5 M_Earth, splitting it equally between protoplanets and planetesimals, because the Venus/Earth analogs were somewhat too small in the original C98 setup (Figure 4).

The results for the W11e/case1 model are documented in Figures 14–16 and Table 1. Overall, the results are similar to those obtained (with smaller statistics) for the W11/case1 model. The radial mass distribution of final planets has the desired profile (Figures 14), but the Mercury analogs are again somewhat too massive. The orbital excitation of Venus/Earth planets in W11e/case1 is slightly larger than in W11/case1 (Table 1), probably because the planetesimals beyond 1.5 au had excited orbits in W11e and could not damp the planetary eccentricities and inclinations as efficiently as in W11 (via a secular coupling). The average Mars mass in W11e/case1 is 0.26 M_Earth, about 2.4 times higher than what would be ideal, but this is slightly biased by the fact that planets at 1.2–1.5 au tend to be more massive (Figure 14). If a stricter selection of Mars analogs were made, the average mass would be lower (e.g., M  ≃  0.2 M_Earth for 1.4 < a < 1.6 au). There are 0.21 strict Mars analogs per good simulation, which is an improvement from W11/case1, where we did not have the needed resolution to identify any strict Mars analogs.

Zoom In Zoom Out Reset image size

Overall, 8% of W11e/case1 simulations produced satisfactory results for the Mars mass, Δa and 〈e, i〉 (as defined by dashed rectangles in Figures 16). This may seem small, but if the probability of satisfying each criterion were 50%, about 13% cases would be expected to satisfy them all (assuming no correlations). From Figure 16(a), we see that the most persistent problem remains matching Δa. All but one satisfactory simulation (within the dashed rectangles in Figures 16) show Δa values that are higher than Δa = 0.277 au of the real planets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The orbital structure of the asteroid belt in W11e/case1 looks great (Figure 15). About 11% of main belt asteroids survive by t = 200 Myr, corresponding to a factor of ∼3 larger dynamical depletion than in the W11/case1 (due to the initially excited orbits of asteroids in W11e). ⁹ Assuming that the asteroid belt was depleted by an additional factor of ∼2 in the 0.2–4.5 Gyr interval (e.g., Minton & Malhotra 2010; Nesvorný et al. 2017), we find f_ast ∼ 0.055 overall. For comparison, Deienno et al. (2016) reported that ∼5% of bodies starting with 1.8 < a < 3.6 au survive in their simulations over 4.5 Gyr.

The current mass of the asteroid belt is ≃5 × 10⁻⁴ M_Earth. Assuming that the total depletion factor over 4.5 Gyr is f_ast = 0.05, we estimate that the initial mass at 2.2–3.5 au and t₀ was 5 × 10⁻⁴/0.05 = 10⁻² M_Earth. This is a factor of ∼100 lower than the mass inferred from the Minimum Mass Solar Nebula (MMSN; Weidenschilling 1977; Hayashi 1981). This means that (1) ∼100 of the population was depleted prior to t₀, which would be consistent with Jupiter's inward-then-outward migration in the Grand Tack model (Walsh et al. 2011), (2) the planetesimal formation at 2.2–3.5 au was relatively inefficient (e.g., Raymond & Izidoro 2017), perhaps due to the early formation of Jupiter (Kruijer et al. 2017), or (3) the MMSN is not a good reference (e.g., observed protoplanetary disks often show radial gaps; Andrews 2020).

Finally, we briefly comment on C98s/case1. This model is not successful for several different reasons. The planetary-mass distribution is similar to that in W11e/case1, but the Mercury analogs are far too massive due to the steeper radial profile of protoplanets (initially extending down to r_in = 0.3 au). About one planetary-mass body is stranded in the asteroid belt at t = 200 Myr in C98s/case1, which is worse than what we found in C98/case1: the slightly steeper radial profile did not help. The orbital excitation of Venus/Earth is too large. In fact, all 41 good simulations in C98s/case1 ended with 〈e, i〉 > 0.03. The average AMD is S_d = 0.0068, about 3.8 times larger that the real value. As for the RMC and radial separation of Venus and Earth, we find average S_c = 24.4 and average Δa = 0.375 au, again significantly larger than the real values. Interestingly, unlike in Figure 16(a), where the Δa values are concentrated near 0.3–0.4 au in W11e/case1, here Δa ranges from 0.2 to 0.6 au, and a small fraction of cases yield the right orbital separation. The smaller variability of Δa (and S_c) in W11e/case1 is most likely related to the initial distribution of protoplanets in the narrow annulus.

Jacobson & Morbidelli (2014) suggested that the RMC can be enhanced in the narrow-annulus model, or equivalently Δa can be reduced, if the late stage of the terrestrial planet accretion started with Mars-size protoplanets at t₀. If the initial protoplanets are more massive, they can grow faster by accretional collisions, and there is presumably not enough time for them to spread radially. To test this effect, we perform 100 additional simulations for both the W11e/case1 and C98s/case1 models, but this time starting with Mars-mass protoplanets at t₀ (20 protoplanets in W11e, 25 protoplanets in C98; these jobs are called W11e/20M and C98s/25M; Table 1). We find that the C98s/25M model does not offer any benefits over the C98 or C98s models. The radial mass distribution of the final planets in C98s/25M is relatively flat, and there are no strict Mars analogs (Table 1). For brevity, we do not discuss C98s/25M here.

We find that the RMC values are generally higher, in agreement with Jacobson & Morbidelli (2014), and the AMD values are generally lower in the W11e/20M model (Table 1). Specifically, we obtain an average S_c = 51.2 in W11e/20M, to be compared with S_c = 46.4 in the standard W11e/case1 and real S_c = 89.9, and S_d = 0.0034 in W11e/20M, to be compared with S_d = 0.0040 in the standard W11e/case1 and real S_d = 0.0018. The differences between the model and real values are still large, but this is influenced by Mercury/Mars. As we discussed in Section 2.4, if Mercury/Mars forms too massive in the simulations, this affects S_c and S_d. This problem is exacerbated in W11e/20M, where the protoplanets are massive at t₀ and Mercury/Mars end up being more massive on average (Table 1). ¹⁰ It may thus seem better to utilize Δa and 〈e, i〉 instead.

The average Δa values are practically the same in the W11e and W11e/20M models (Δa ≃ 0.34 au, Table 1). We therefore do not see much influence of the starting mass of protoplanets on the mean separation of Venus/Earth. Compared to the real Δa  ≃  0.28 au, the average model values are only ∼20% higher (whereas the average model S_c values discussed above are 40% lower that the real S_c = 89.9). The range of Δa values is larger in W11e/20M (≃0.15–0.5; Figure 17(a)) than in W11e (≃0.25–0.4; Figure 16(a)) probably because the accretion of Venus/Earth is more stochastic if protoplanets start more massive. If protoplanets start with roughly the Mars mass, in nearly 20% of simulations of W11e/20M, the strict Mars analogs are scattered from r < 1 au to 1.2–1.8 au and only grow by accreting planetesimals. They end up having masses of 1–1.5 M_Mars (Figure 17). Unlike in Figure 16(a), where the W11e model did not no produce any really close matches to the real terrestrial planets (model Mars was either too small or too big, or model Δa was too large), there are some very good matches in Figure 17(a).

Zoom In Zoom Out Reset image size

The initial mass of the protoplanets affects 〈e, i〉. We obtain 〈e, i〉 = 0.053 in the W11e model and 〈e, i〉 = 0.039 in the W11e/20M model (to be compared with real 〈e, i〉 = 0.027). This is a significant improvement. Moreover, if we only count the simulations where a Mars analog formed in W11e/20M (56% of cases; Figure 17(b)), we find average 〈e, i〉 = 0.032—only 15% higher than the real value. In addition, if we only count the cases with strict Mars analogs (M < 0.2 M_Earth; 17% of cases; Figure 17(b)), then 〈e, i〉 = 0.027, which is identical to the real value. This trend is not as strong in W11e where counting the simulations with a Mars analog gives 〈e, i〉 = 0.046. We therefore endorse the suggestion of Jacobson & Morbidelli (2014) that it is beneficial to initiate the simulations with Mars-mass protoplanets. ¹¹

Finally, we test possible trends with the mass of Venus/Earth analogs. For that, we compute the total mass of the Venus/Earth planets in the good systems and normalize it by the total mass of real Venus and Earth, M_VE = (M₁ + M₂)/M_real with M_real = 1.82 M_Earth. We find that both the W11e and W11e/20M models give average M_VE = 1.07 (i.e., the mass is only 7% higher than the actual value). This is typically caused by a slightly larger mass of the Venus analogs (Figure 14). The radial separation of Venus/Earth planets in W11e/20M correlates with their total mass (Figure 18(a)) but the orbital excitation does not (Figure 18(b)). If we restrict the range to strict Venus/Earth analogs, 0.9 < M_VE < 1.1, and only count systems where a Mars analog formed, we obtain Δa = 0.317 au and 〈e, i〉 = 0.030 in W11e/20M, which is a significant improvement over the values listed in Table 1. We thus see that obtaining the correct planetary masses in the simulations links to better results for Δa au and 〈e, i〉 as well. Such a trend is not expected a priori. The trends for W11e are similar but the overall results are not as good (e.g., 〈e, i〉 = 0.044 for 0.9 < M_VE < 1.1 with a Mars analog).

Zoom In Zoom Out Reset image size

Here we highlight a case from the W11e/20M model (job number 35) that produced a particularly good match to the terrestrial planets (Figures 19 and 20). Venus, Earth, and Mars formed in this simulation with nearly perfect masses and semimajor axis: model 1.017 M_Earth and a = 0.727 au for Venus (0.815 M_Earth and a = 0.723 au), model 1.083 M_Earth and a = 1.034 au for Earth (1.0 M_Earth and a = 1.0 au), model 0.116 M_Earth and a = 1.525 au for Mars (0.107 M_Earth and a = 1.524 au), where the numbers in parentheses list the real values for reference. The separation of Venus and Earth is excellent (Δa = 0.0307 au versus real Δa = 0.0277 au). Here the main difference is that Venus ended with the mass about 20% higher than what would be ideal. The Mercury analog is very massive but that may be related to our assumption of inelastic mergers, as we discussed previously. We therefore do not comment on Mercury below.

Zoom In Zoom Out Reset image size

The model orbits of Venus and Earth are slightly less excited than in the real system, resulting in 〈e, i〉 = 0.0176 (to be compared with real 〈e, i〉 = 0.0274). This is not bad. In fact, it is encouraging to see, for a change, a case with low AMD. The model orbit of Mars has a nearly perfect inclination (mean i = 451 versus real mean i = 441) and somewhat lower eccentricity (mean e = 0.028 versus real mean e = 0.069). We would probably need to perform a larger suite of simulations, possibly with a less massive planetesimal population (Jacobson & Morbidelli 2014), to obtain a better match to Mars's eccentricity. The orbits of main belt asteroids look good (compare Figures 20 and 7). The model orbital inclinations are slightly lower, but this is related to our choice of initial conditions and could easily be remedied by adopting a slightly wider initial distribution.

Zoom In Zoom Out Reset image size

A spectacular feature of the job highlighted here is related to Earth's accretion history and the Moon-forming impact (Figure 21). The initial growth of proto-Earth is fast. After accreting several Mars-class protoplanets, the proto-Earth reaches mass ≃0.52 M_Earth at t ≃ 3 Myr after the start of our simulation. A prolonged stage of planetesimal accretion follows during which the proto-Earth modestly grows to ≃0.57 M_Earth by t ≃ 40 Myr. Shortly after that, at t = 41.3 Myr, an accretional collision between two roughly equal-mass bodies occurs, and the Earth mass shoots up to ≃1.05 M_Earth. ¹² For comparison, Venus grows to ≃85% of its final mass only 3 Myr after the start of the simulation; the remaining ≃15% of mass is supplied to Venus by planetesimals over an extended time interval (∼100 Myr). This supports the suggestion of Jacobson et al. (2017), who argued that Venus has not developed a persistent magnetic field because it did not experience any late energetic impacts that would mechanically stir the core and create a long-lasting dynamo.

Zoom In Zoom Out Reset image size

The accretion growth of Earth shown in Figure 21 would satisfy constraints from the Hf–W and U–Pb isotope systematics, assuming that the degree of metal-silicate equilibration was intermediate between the full and 40% equilibration considered in Kleine & Walker (2017, their Figure 8). The last major accretion event at t = 41.3 Myr corresponds to a very low speed collision (the impact-to-escape speed ratio v_imp/v_esc = 1.01) between two nearly equal-mass protoplanets (the impactor-to-total mass ratio Γ = 0.46). This falls into the preferred regime of Moon-forming impacts investigated in Canup (2012). In particular, the collision of two protoplanets with Γ > 0.4 would lead to the Moon formation from a well-mixed disk and would therefore imply matching oxygen-isotope compositions of Earth and the Moon. The two protoplanets involved in the Moon-forming collision were nearly fully formed and on their pre-impact orbits (a = 0.95 and 1.11 au; radial separation ≃0.16 au) already by t  ≃  3 Myr after t₀. Their orbits intermittently crossed each other during orbital eccentricity oscillations, but it took nearly 40 Myr before the two protoplanets actually collided. A similar case was investigated in DeSouza et al. (2020). This demonstrates that the Moon-forming collision could have happened tens of megayears after the giant-planet instability.