# Early Solar System instability triggered by dispersal of the gaseous disk

Beibei Liu, Sean N. Raymond, Seth A. Jacobson

► To cite this version:

Beibei Liu, Sean N. Raymond, Seth A. Jacobson. Early Solar System instability triggered by dispersal of the gaseous disk. Nature, 2022, 604, pp.643-646. 10.1038/s41586-022-04535-1. insu-03678901

HAL Id: insu-03678901

https://insu.hal.science/insu-03678901v1

Submitted on 6 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Early Solar System instability triggered by dispersal of the gaseous disk

https://doi.org/10.1038/s41586-022-04535-1

Beibei Liu<fnref n="1" />,<fnref n="2" />✉, Sean N. Raymond³ & Seth A. Jacobson⁴

Received: 14 December 2020

Accepted: 8 February 2022

Published online: 27 April 2022

<figure><img src="image_2.png" /></figure>

Check for updates

The Solar System's orbital structure is thought to have been sculpted by an episode of dynamical instability among the giant planets<fnref n="1" />⁻⁴. However, the instability trigger and timing have not been clearly established⁵⁻⁹. Hydrodynamical modelling has shown that while the Sun's gaseous protoplanetary disk was present the giant planets migrated into a compact orbital configuration in a chain of resonances<fnref n="2" />,¹⁰. Here we use dynamical simulations to show that the giant planets' instability was probably triggered by the dispersal of the gaseous disk. As the disk evaporated from the inside out, its inner edge swept successively across and dynamically perturbed each planet's orbit in turn. The associated orbital shift caused a dynamical compression of the exterior part of the system, ultimately triggering instability. The final orbits of our simulated systems match those of the Solar System for a viable range of astrophysical parameters. The giant planet instability therefore took place as the gaseous disk dissipated, constrained by astronomical observations to be a few to ten million years after the birth of the Solar System¹¹. Terrestrial planet formation would not complete until after such an early giant planet instability¹²,¹³; the growing terrestrial planets may even have been sculpted by its perturbations, explaining the small mass of Mars relative to Earth¹⁴.

We modelled the dynamical consequences of the dispersal of the Sun's gaseous disk. Stellar photoevaporation dominates the mass loss during this advanced phase, causing the disk to dissipate from the inside out¹⁵,¹⁶. Whereas planets embedded in the disk feel 'two-sided' gravitational torques from both the interior and exterior parts of the disk, planets at the disk's inner edge only interact with the gas exterior to their orbits. As a result of these larger 'one-sided' torques, a planet below the mass threshold for opening a gap will stop migrating inwards at the disk inner edge¹⁷,¹⁸. If the inner edge itself moves outwards owing to disk dispersal, then the planet may subsequently migrate outwards along with it (Methods). This mechanism is termed 'rebound', and was first applied in the context of the magnetospheric cavity on sub-AU scales (AU, astronomical unit) to explain the architecture of close-in super-Earth planets¹⁹,²⁰.

Figure 1 demonstrates an example simulation of a dynamical instability triggered by the disk's dispersal. The expanding edge of the inner disk cavity does not affect all planets equally. Because Jupiter is sufficiently massive enough to open a deep gap around its horseshoe region, the corresponding corotation torque diminishes and the rebound is quenched (Methods). Jupiter then simply enters the cavity as the inner disk edge sweeps by. The one-sided torque is strong enough to expand Saturn's orbit outwards when the disk edge approaches Saturn at $t = 0.6$ Myr (Fig. 1), moving Jupiter and Saturn out of their shared resonance. As Saturn migrates outwards with the expanding cavity, the spacing between the orbits of the outer planets is compressed. The eccentricities of the ice giants increase owing to this dynamical compression. Saturn is left behind and enters the cavity at 9 AU when

$t = 0.65$ Myr. Meanwhile, the innermost ice giant planet becomes so dynamically excited that its orbit crosses Saturn's, and the two planets undergo a close gravitational encounter. This triggers a dynamical instability and the system becomes chaotic: the third ice giant is scattered outwards, whereas the innermost ice giant is eventually ejected into interstellar space at $t = 0.85$ Myr after a series of close encounters with Jupiter. The planets' final orbits are close to those of the present-day Solar System giant planets.

Such a rebound-triggered instability is consistent with the Solar System's orbital architecture. To demonstrate this, we conducted more than 14,000 numerical simulations like the one from Fig. 1, varying three different aspects of the initial conditions (Extended Data Table 1). First, we tested a wide range of plausible starting configurations for the number of ice giants (two, three or four) and their initial orbital resonant states. Second, we used a Monte Carlo method to test the effects of important disk parameters—the onset mass-loss rate $M_{pho}$, the disk dispersal timescale $\tau_d$ and the expansion rate of the inner cavity $v_r$—across the full range of astronomically relevant values. Third, we ran each simulation twice: once including the effect of inside-out disk dissipation (that is, with rebound) and once assuming the disk dissipates smoothly at all radii (that is, without rebound). As a basic check, we used two system-level indicators to test whether our simulated systems are consistent with the global properties of the Solar System: the (normalized) angular momentum deficit (AMD), a measure of the dynamical excitation of the system, and the radial mass concentration statistic (RMC), a measure of the orbital spacing of the system.

<footnote n="1">¹Zhejiang Institute of Modern Physics, Department of Physics & Zhejiang University-Purple Mountain Observatory Joint Research Center for Astronomy, Zhejiang University, Hangzhou, China.</footnote>

<footnote n="2">²Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Lund, Sweden. ³Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, Pessac, France. ⁴Department of Earth and Environmental Sciences, Michigan State University, East Lansing, MI, USA. ✉e-mail: bbliu@zju.edu.cn</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 1 | An early dynamical instability triggered by the dispersal of the Sun's protoplanetary disk.** The initial system consisted of five giant planets: Jupiter, Saturn and three 15 $M_\odot$ ice giants, one of which was ejected into interstellar space during the instability. The curves show the orbital evolution of each body including its semimajor axis (thick), perihelion and aphelion (thin). The black dashed line tracks the edge of the disk's expanding inner cavity. We do not follow the early evolution through the entire gas-rich disk phase, so the onset of disk dispersal is set arbitrarily to be 0.5 Myr after the start of the simulation. The semimajor axes and eccentricities of the present-day giant planets are shown at the right, with vertical lines extending from perihelion to aphelion. The adopted disk parameters are $M_{\text{pho}} = 4 \times 10^{-10} M_\odot \text{ yr}^{-1}$, $\tau_d = 8.6 \times 10^5 \text{ yr}$ and $v_r = 35 \text{ AU Myr}^{-1}$.

When the rebound effect is included in our simulations, the surviving planetary systems fill the AMD-RMC phase space that matches the Solar System (Fig. 2). That space is mostly empty when rebound is not included because dynamical instabilities are much less frequent. More than 90% of systems starting in 3:2 resonances went unstable when rebound was included but only 39% when rebound was ignored. Likewise, 78% of systems with a chain of 2:1 (Jupiter and Saturn) and 3:2 resonances went unstable when rebound was included versus 31% when it was ignored. The rebound-triggered instability occurs across all astronomically relevant disk parameter values (Methods). In these simulations, we had adopted a moderately viscous disk in which Saturn did not open a deep gap. However, the rebound mechanism also generates instability in low-viscosity environments in which Saturn is above the gap-opening mass. In that case, the ice giants' scattering propagates to the gas giants, triggering a system-wide instability at a rate that is only modestly lower than in our fiducial simulations (Methods).

In previous studies$^{1,5,8,9,21,22}$, a primordial planetesimal disk typically contained $20-30 M_\odot$ within 30 AU and played a central role in triggering the instability. In our model, the gas disk is the instability trigger, yet interactions with a putative outer planetesimal disk would further spread out the giant planets' orbits and decrease their eccentricities and inclinations. After the gas disk was fully dissipated, we extended a subset of simulations in a gas-free environment for another 100 Myr including an outer planetesimal disk containing a total of 5, 10 or $20 M_\odot$. In an example with four giant planets (Fig. 3a, b), the rebound-driven instability leaves the system in a configuration that is more compact than the real one. Yet, during the planetesimal disk phase ($t > 10$ Myr) the orbital radius of Uranus and Neptune increased, and the eccentricities of all planets were damped, resulting in a configuration closer to that of the Solar System. An example starting with five giant planets with an outer planetesimal disk of $5 M_\odot$ followed a similar evolutionary path (Fig. 3c, d). In dynamical terms, the rebound-triggered instability increases a giant planet system's level of orbital excitation (and its AMD) and decreases its degree of radial concentration (and RMC), whereas later interactions with the planetesimal disk tend to decrease both the AMD and RMC.

<figure><img src="image_4.png" /></figure>

**Fig. 2 | Metrics for surviving planetary systems of a subsample of our simulations in matching the Solar System.** The simulations on the left included the rebound effect and those on the right did not. Each simulation started with our four present-day giant planets plus one additional ice giant planet. In the top panels, the giant planets were initially placed in a chain of 3:2 orbital resonances. In the bottom panels, Jupiter and Saturn were in a 2:1 resonance and other neighbouring planet pairs were in 3:2 resonances. Each symbol represents the outcome of a given simulation at $t = 10$ Myr. The colour indicates the timing of the instability after the start of gas disk dispersal; pink systems did not undergo an instability (no collision and/or ejection). Diamonds, circles and triangles correspond to systems with five, four and three or fewer surviving planets, respectively. The arrow gives the initial radial mass concentration of the system. The Solar System is marked as a red star for comparison. Comparable figures presenting different subsamples starting from different orbital configurations are included as Extended Data Figs. 1–3. MMR, mean-motion resonance.

The final system architectures provided a better match to the Solar System when planetesimal disks were included (Fig. 4). One challenge for our simulations is adequately exciting Jupiter's eccentricity to its present value of 0.046. This is a systematic problem in simulations of the instability$^{1,3,14,22,23}$. A possible solution is that Jupiter's orbit was already modestly eccentric at the tail end of the gaseous disk phase$^{10,23}$. We do not attempt to explain the Kuiper Belt's architecture in this work, as the triggering mechanism is not the central aspect for establishing these small body populations. The chaos of the instability erases the dynamical memory of the initial triggering—a defining feature of chaos—and the dissipating gas in the dispersal phase only plays a minor role in damping the random velocities of small bodies once they get excited. Thus, results regarding existing models of small body evolution after giant planet instability hold regardless of the triggering mechanism.

A rebound-triggered instability at the time of disk dispersal fills an important gap in Solar System chronology. Observations of the frequency of disks in star clusters of different ages find that the typical disk lifetime is a few to 10 Myr (ref. $^{11}$). The giant planet instability was initially invoked as a delayed event to coincide with the 'late heavy bombardment'$^{21}$. However, recent re-appraisal of the cosmochemical constraints indicates that there was likely no late spike ('terminal
<figure><img src="image_5.png" /></figure>

**Fig. 3 | Dynamical evolution of giant planets in both gas disk dispersal phase and gas-free, planetesimal disk phase.** The initial system consisted of four giant planets in 2:1 resonances (upper) or five giant planets in a combined 2:1 and 3:2 resonances (lower). The left panels show the orbital evolution of each body including its semimajor axis, perihelion and aphelion. The black dashed line tracks the edge of the disk's expanding inner cavity. The onset of

disk dispersal is set arbitrarily to be 0.5 Myr after the start of the simulation. The planetesimal disks of $10M_⊕$ and $5M_⊕$ are implemented after 10 Myr in the above two configurations. The right panels provide the corresponding system's RMC and normalized AMD at $t = 0$ yr, 10 Myr and 100 Myr with a planetesimal disk of $5M_⊕$ (pink), $10M_⊕$ (brown) and $20M_⊕$ (purple), respectively. The Solar System is marked as a red star for comparison.

cataclysm') in the bombardment rate²⁴. Instead, constraints from a binary Jupiter Trojan⁶ and ages of meteoritic inclusions⁷ indicate that the instability took place no later than around 20–100 Myr after the birth of the Solar System. An instability within 10 Myr would have perturbed the final assembly of the terrestrial planets, and an early instability may explain a number of features of the inner Solar System including the large Earth-to-Mars mass ratio and the dynamical excitation of the asteroid belt¹⁴.

Our model provides a generic trigger for dynamical instability linked with the observed timescale for disk dispersal¹¹. Early models relied on fine-tuning the distance between the ice giants and outer planetesimal disk or the degree of self-interaction between planetesimals to match the assumed late timing of the instability⁵,²¹. More recent studies with more self-consistent outer planetesimal disks systematically find shorter instability timescales, but with broad distributions that extend to around 100 Myr and an uncertainty in the triggering mechanism

<figure><img src="image_6.png" /></figure>

**Fig. 4 | Final orbits of the giant planets in simulations that included both the gaseous disk and an outer planetesimal disk.** Only the simulations that finished with four planets are shown. Simulations with and without

planetesimal disks are plotted in circles and triangles, respectively. The Solar System is marked as a star.
itself⁸,⁹. Given its robustness to disk parameters (Extended Data Figs. 4 and 5), an early, rebound-triggered instability is essentially unavoidable for the Solar System. Such a trigger also implies that the planets' final architecture depends only weakly on the mass of the outer planetesimetal disk (Fig. 4), which may have been less massive than previously thought. For instance, the cratering record on outer Solar System moons does not appear to match predictions from models of the collisional evolution of the massive primordial Kuiper Belt needed to trigger the giant planets' instability²⁵,²⁶.

The rebound effect may explain why dynamical instabilities appear to be nearly ubiquitous in exoplanetary systems²⁷. The broad eccentricity distribution of giant exoplanets can be matched if 75–95% of all giant planet systems we see are the survivors of dynamical instabilities²⁸. Although rebound only affects planets with masses below the gap-opening mass (such as ice giants), microlensing studies find that such planets are extremely common²⁹. The instabilities in ice giant systems spread to nearby gas giants at a high probability (Methods). Likewise, the orbital period ratio distribution of close-in small planets is consistent with the vast majority of systems having undergone dynamical instabilities³⁰. Thus, the rebound effect during disk dispersal may be a nearly universal process affecting not just our Solar System but planetary systems across the Galaxy.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04535-1.

1. Tsiganis, K., Gomes, R., Morbidelli, A. & Levison, H. F. Origin of the orbital architecture of the giant planets of the Solar System. *Nature* **435**, 459–461 (2005).

2. Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F. & Gomes, R. Dynamics of the giant planets of the Solar System in the gaseous protoplanetary disk and their relationship to the current orbital architecture. *Astron. J.* **134**, 1790–1798 (2007).

3. Batygin, K., Brown, M. E. & Betts, H. Instability-driven dynamical evolution model of a primordially five-planet outer Solar System. *Astrophys. J. Lett.* **744**, L3 (2012).

4. Nesvorný, D. Dynamical evolution of the early Solar System. *Ann. Rev. Astron. Astrophys.* **56**, 137–174 (2018).

5. Levison, H. F., Morbidelli, A., Tsiganis, K., Nesvorný, D. & Gomes, R. Late orbital instabilities in the outer planets induced by interaction with a self-gravitating planetesimetal disk. *Astron. J.* **142**, 152 (2011).

6. Nesvorný, D., Vokrouhlický, D., Bottke, W. F. & Levison, H. F. Evidence for very early migration of the Solar System planets from the Patroclus–Menoetius binary Jupiter Trojan. *Nat. Astron.* **2**, 878–882 (2018).

7. Mojzsis, S. J., Brasser, R., Kelly, N. M., Abramov, O. & Werner, S. C. Onset of giant planet migration before 4480 million years ago. *Astrophys. J.* **881**, 44 (2019).

8. Quarles, B. & Kaib, N. Instabilities in the early Solar System due to a self-gravitating disk. *Astron. J.* **157**, 67 (2019).

9. de Sousa, R. R. et al. Dynamical evidence for an early giant planet instability. *Icarus* **339**, 113605 (2020).

10. Pierens, A., Raymond, S. N., Nesvorný, D. & Morbidelli, A. Outward migration of Jupiter and Saturn in 3:2 or 2:1 resonance in radiative disks: implications for the Grand Tack and Nice models. *Astrophys. J. Lett.* **795**, L11 (2014).

11. Williams, J. P. & Cieza, L. A. Protoplanetary disks and their evolution. *Ann. Rev. Astron. Astrophys.* **49**, 67–117 (2011).

12. Jacobson, S. A. et al. Highly siderophile elements in Earth's mantle as a clock for the Moon-forming impact. *Nature* **508**, 84–87 (2014).

13. Kleine, T. & Walker, R. J. Tungsten isotopes in planets. *Ann. Rev. Earth Planet. Sci.* **45**, 389–417 (2017).

14. Clement, M. S., Kaib, N. A., Raymond, S. N. & Walsh, K. J. Mars' growth stunted by an early giant planet instability. *Icarus* **311**, 340–356 (2018).

15. Alexander, R., Pascucci, I., Andrews, S., Armitage, P. & Cieza, L. in *Protostars and Planets Vol. VI* (eds. Beuther, H. et al.) 475–496 (Univ. Arizona Press, 2014).

16. Ecolano, B. & Pascucci, I. The dispersal of planet-forming discs: theory confronts observations. *R. Soc. Open Sci.* **4**, 170114 (2017).

17. Masset, F. S., Morbidelli, A., Crida, A. & Ferreira, J. Disk surface density transitions as protoplanet traps. *Astrophys. J.* **642**, 478–487 (2006).

18. Romanova, M. M. et al. 3D simulations of planet trapping at disc-cavity boundaries. *Mon. Not. R. Astron. Soc.* **485**, 2666–2680 (2019).

19. Liu, B., Ormel, C. W. & Lin, D. N. C. Dynamical rearrangement of super-Earths during disk dispersal. I. Outline of the magnetospheric rebound model. *Astron. Astrophys.* **601**, A15 (2017).

20. Liu, B. & Ormel, C. W. Dynamical rearrangement of super-Earths during disk dispersal. II. Assessment of the magnetospheric rebound model for planet formation scenarios. *Astron. Astrophys.* **606**, A66 (2017).

21. Gomes, R., Levison, H. F., Tsiganis, K. & Morbidelli, A. Origin of the cataclysmic Late Heavy Bombardment period of the terrestrial planets. *Nature* **435**, 466–469 (2005).

22. Nesvorný, D. & Morbidelli, A. Statistical study of the early Solar System's instability with four, five, and six giant planets. *Astron. J.* **144**, 117 (2012).

23. Clement, M. et al. Born eccentric: constraints on Jupiter and Saturn's pre-instability orbits. *Icarus* **355**, 114–122 (2021).

24. Zellner, N. E. B. Cataclysm no more: new views on the timing and delivery of lunar impactsors. *Origins Life Evol. Biosphere* **47**, 261–280 (2017).

25. Zahnle, K., Schenk, P., Levison, H. & Dones, L. Cratering rates in the outer Solar System. *Icarus* **163**, 263–289 (2003).

26. Singer, K. N. et al. Impact craters on Pluto and Charon indicate a deficit of small Kuiper belt objects. *Science* **363**, 955–959 (2019).

27. Raymond, S. N., Izidoro, A. & Morbidelli, A. In *Planetary Astrobiology* (eds Meadows, V. et al.), 287–324 (Univ. Arizona Press, 2020).

28. Raymond, S. N., Armitage, P. J. & Gorelick, N. Planet-planet scattering in planetesimetal disks. II. Predictions for outer extrasolar planetary systems. *Astrophys. J.* **711**, 772–795 (2010).

29. Suzuki, D. et al. The exoplanet mass-ratio function from the MOA-II Survey: discovery of a break and likely peak at a Neptune mass. *Astrophys. J.* **833**, 145 (2016).

30. Pu, B. & Wu, Y. Spacing of Kepler planets: sculpting by dynamical instability. *Astrophys. J.* **807**, 44 (2017).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2022
# Methods

## Disk model

The evolution of a protoplanetary disk is driven by internal gas viscous stresses and external stellar ultraviolet/X-ray photoevaporation<sup>15,16</sup>. At early times the influence of stellar photoevaporation is negligible and the disk evolves viscously in a quasi-steady state such that $\dot{M} = 3\pi v\Sigma$, where $\dot{M}$ is the gas disk accretion rate, $\Sigma$ is the gas surface density, $v = \alpha c_s H$ is the disk viscosity<sup>31</sup>, $\alpha$ is the viscous efficiency parameter and $c_s$ and $H$ are the gas disk sound speed and scale height, respectively. At late times, after the disk accretion rate $\dot{M}$ drops below the stellar photoevaporation mass-loss rate $M_{pho}$, disk dissipation is dominated by the photoevaporative wind originating from stellar high-energy radiation<sup>32,33</sup>.

When the thermal energy of the disk gas is greater than its gravitational bounding energy beyond a threshold disk radius $R_{in}$, the disk gas escapes as a wind. Gas inside $R_{in}$ is rapidly accreted on the local viscous timescale. Hence, an inner disk cavity is opened and the main outer disk is optically thin to the direct radiation from the central star<sup>32</sup>. As such, the photoevaporating mass-loss rate is given by<sup>15</sup>

$$ \dot{M}_{pho} \approx 10^{-9} \left( \frac{\Phi}{10^{41} \text{ s}^{-1}} \right)^{1/2} \left( \frac{R_{in}}{5 \text{ AU}} \right)^{1/2} M_{\odot} \text{yr}^{-1}, \quad (1) $$

where the ionizing flux of the young central star $\Phi$ is order of $10^{41}-10^{42} \text{ s}^{-1}$ and the initial size of the inner cavity $R_{in}$ is 1–5 AU. An initially smaller size cavity does not qualitatively change our results, and the determining factor is $M_{pho}$ when the inner disk edge sweeps at Jupiter’s location. The mass-loss rate depends on the primary stellar incident spectrum (extreme ultraviolet, far ultraviolet or X-ray), grain species, abundance and disk chemistry. The fiducial onset mass-loss rate $\dot{M}_{pho}$ of $10^{-9} M_{\odot} \text{yr}^{-1}$ in equation (1) is taken from extreme-ultraviolet-driven stellar radiation models<sup>32</sup>. If the photoevaporation is instead dominated by X-ray radiation, $\dot{M}_{pho}$ can be higher<sup>33</sup>. We treat $\dot{M}_{pho}$ as a free parameter and vary it from $3 \times 10^{-9}$ to $3 \times 10^{-10} M_{\odot} \text{yr}^{-1}$ in the parameter study.

In this photoevaporation-driven disk dispersal phase, the gas surface density can be written as

$$ \Sigma(t) = \Sigma_{pho} \exp\left[-\frac{t}{\tau_d}\right] = \frac{\dot{M}_{pho}}{3\pi v} \exp\left[-\frac{t}{\tau_d}\right], \quad (2) $$

where $\Sigma_{pho}$ is the gas surface density when the stellar photoevaporation becomes dominant and $\tau_d$ is the gaseous disk dispersal timescale. Noticeably, two different timescales are related to gas disk dissipation. First, the disk lifetime is inferred to be 1–10 Myr with a median value of 3 Myr (ref. <sup>34</sup>). Second, as only a small fraction of young stars have evolved disks caught in the transition between disk-bearing and disk-less states, disk clearing is expected to be very rapid, typically an order of magnitude shorter than the disk lifetime<sup>35,36</sup>. Based on these observational constraints, it is reasonable to assume that $\tau_d$ is in a range between 0.1 and 1 Myr.

We assume an optically thin, flaring disk in this dispersing phase<sup>37</sup>, and the corresponding temperature and aspect ratio are

$$ \begin{aligned} T &= \left( \frac{L_*}{16\pi\sigma_{SB}r^2} \right)^{1/4} = T_{\text{in}} \left( \frac{r}{R_{\text{in}}} \right)^{-1/2} \quad \text{and} \\ h &= \frac{H}{r} = \sqrt{\frac{k_B T}{\mu m_p G M_*}} = h_{\text{in}} \left( \frac{r}{R_{\text{in}}} \right)^{1/4}, \end{aligned} \quad (3) $$

where $M$, and $L$ are the stellar mass and luminosity, $r$ is the disk radial distance, $G$ is the gravitational constant, $\sigma_{SB}$ is the Stefan–Boltzmann constant, $k_B$ is the Boltzmann constant, $\mu$ is the gas mean molecule weight and $m_p$ is the proton mass. A previous study<sup>37</sup> calculated that $h = 0.033$ at 1 AU by adopting $L_* = 1L_\odot$. However, during the pre-main-sequence period (stellar age less than a few megayears), the Sun underwent

gravitational contraction and was more luminous (around $(1-10)L_\odot$, along the Hayashi track) than its main-sequence stage. Hence, we choose a higher luminosity for a solar-mass pre-main-sequence star of $L_* = 4L_\odot$ and assume it remains constant during the relatively rapid disk dispersal phase. In this circumstance, the disk aspect ratio is 0.04 at 1 AU. We adopt $R_{in} = 5$ AU and, therefore, $h_{in} = 0.06$ and $T_{in} = 180$ K, where the subscript indicates the quantity evaluated at $R_{in}$ when the stellar photoevaporation dominates. Based on the quasi-steady-state assumption, we rewrite the gas surface density as

$$ \Sigma(t) = \frac{\dot{M}_{pho}}{3\pi v} \exp\left[-\frac{t}{\tau_d}\right] = \Sigma_{\text{in}} \left(\frac{r}{R_{\text{in}}}\right)^{-1} \exp\left[-\frac{t}{\tau_d}\right], \quad (4) $$

where $\Sigma_{\text{in}} = 4 \text{ g cm}^{-2}$ for the fiducial disk parameters of $\dot{M}_{pho} = 10^{-9} M_{\odot} \text{yr}^{-1}$ and $\alpha = 0.005$.

The cavity spreads from inside out as disk gas disperses. We assume that the cavity expands with a constant speed $v_r$ for simplicity. The cavity expansion rate should be principally determined by the stellar ionizing flux, and the size and mass of the disk and, hence, is not a fully independent variable. Nonetheless, when the dispersal time and disk mass are fixed, $v_r$ directly reflects the size of the disk. To a first-order approximation, for a disk with an outer edge $R_{out}$ of 30 AU and dispersal timescale of 0.5 Myr, $v_r \approx R_{out}/\tau_d = 60$ AU Myr<sup>-1</sup>. Owing to the uncertainties in $R_{out}$ and $\tau_d$, $v_r$ could plausibly span two orders of magnitude from 10 to $10^3$ AU Myr<sup>-1</sup>. We neglect the effect of disk self-gravity in this study, which might induce secular resonances with the gas giant planets to sweep through the inner Solar System and cause the orbital excitation of asteroids and terrestrial planets; understanding how the rebound instability impacts the inner Solar System is a focus for future work.

In brief, the final photoevaporation-driven disk dispersal can be described by three key model parameters: the onset mass-loss rate when photoevaporation dominates $\dot{M}_{pho}$, the disk dispersal timescale $\tau_d$ and the inner disk cavity expansion speed $v_r$.

## Planet–gas disk interaction

Here we summarize the formulas of planet migration torques used in our study based on ref. <sup>19</sup>. When the planet is far from the inner edge of the disk, it feels the torque from the disk gas on both sides. The two-sided torque $\Gamma_{2s}$ is adopted from equation (49) of ref. <sup>38</sup>:

$$ \frac{\Gamma_{2s}}{m_p(r_p \Omega_p)^2} = -2.3 q_d \frac{q_p}{h^2}, \quad (5) $$

where $\Omega_k$ is the Keplerian angular velocity, $r_p$ is the distance from the planet to the central star and $q_d \equiv \Sigma_{\text{g}} r_p^2 / M$ and $q_p \equiv m_p / M$ are the mass ratio between the local gas disk and star, and the mass ratio between the planet and star, respectively. When the planet is at the disk edge, only the one-sided torque $\Gamma_{1s}$ exists:

$$ \frac{\Gamma_{1s}}{m_p(r_p \Omega_p)^2} = C_{hs} q_d \left(\frac{q_p}{h^3}\right)^{1/2} + C_l q_d \frac{q_p}{h^3}, \quad (6) $$

where the first and second terms on the right-hand side of equation (6) are the corotation and Lindblad torque components, respectively, with $C_{hs} = 2.46$ and $C_l = -0.65$. In the above equations, the corotation torques are expressed for planets on circular orbits. The saturation of the corotation torque due to non-zero eccentricity is accounted for by adopting equation (13) of ref. <sup>19</sup>. To derive the corotation torque in equation (6), we assume that, at the disk inner edge, the gas removal time $t_{\text{removal}}$ is faster than the gas libration time in the planet horseshoe region $t_{\text{lib}}$. We show that this is a justified treatment as follows. First, at the inner disk edge the gas removal time should be no longer than the viscous diffusion time $t_{\text{vis}}$; otherwise the gas would be accumulated. Thus, we have $t_{\text{removal}} \lesssim t_{\text{vis}} \approx x_{hs}^2/v$, where $x_{hs} \approx r_p/\sqrt{q_p/h}$ is the half-width of the planet horseshoe region. The gas libration time for a planet can be written as
$t_{\text{lib}} \approx 8\pi r/(3\Omega_k x_{\text{hs}})$. Then $t_{\text{removal}} < t_{\text{lib}}$ is required for the one-sided corotation torque in equation (6). The above condition is satisfied as long as $t_{\text{vis}} < t_{\text{lib}}$, and one can obtain that $q \lesssim (8\pi/3\alpha)^{2/3}h^{7/3} \approx 2.4 \times 10^{-4}$ by adopting $\alpha = 0.005$ and $h = 0.07$ (at the location of the 3:2 mean-motion resonance with Jupiter). In other words, to fulfil the condition of large amplitude one-sided corotation torque, the planet needs to be no more massive than Saturn in our disk model.

For a general situation, the total torque can be expressed by interpolating the torques between the two regimes:

$$ \Gamma = f f_{1s} + (1-f) \Gamma_{2s}, \quad (7) $$

where $f = \exp[-(r_p - r_{\text{in}})/x_{\text{hs}}]$ is a smooth fitting function and $f=1$ and $f=0$ refer to the planet at the disk edge and far away from the disk edge, respectively. The disk torque is added into the equation of the planet motion in a cylindrical coordinate³⁹:

$$ \frac{dv_{\theta}}{dt} = \frac{\Gamma}{m_{p}r_{p}}, \quad (8) $$

$$ \frac{dv_{r}}{dt} = - \frac{v_{r}}{t_{\text{ecc}}}, \quad (9) $$

$$ \frac{dv_z}{dt} = - \frac{v_z}{t_{\text{inc}}}. \quad (10) $$

Both orbital migration and eccentricity/inclination damping are included (see equations (22) and (23) of ref. ³⁹), where $t_{\text{ecc}} = t_{\text{inc}}$ is assumed. Although the above equations are derived analytically, we note that the termination of planet migration at the inner disk edge due to one-sided torques is also obtained in hydrodynamic simulations¹⁷,¹⁸.

For planet outward migration with the retreating disk, the one-sided corotation torque needs to be larger than the one-sided Lindblad torque. Thus, from equation (6) the planet needs to satisfy the condition $q_p < (C_{hs}/C_L)^2 h^3$. The above torque formulae are derived in the linear regime. However, when the planet is massive to clear the local surrounding gas in its horseshoe region⁴⁰, the disk feedback is non-trivial and equations (5) and (6) are no longer applicable. Because the positive corotation torque diminishes owing to gap formation, the planet beyond the gap-opening mass fails to undergo outward migration. The gap opening requires the planet's Hill sphere $R_H \equiv (m_p/3M_\odot)^{1/3}r$ to be larger than the disk scale height $H$. Therefore, $q_p < 3h^3$ is needed for torques in the linear regime. In addition, dedicated hydrodynamic simulations indicate that the planet needs to fulfil this criterion as well⁴¹:

$$ \frac{3 H}{4 R_H} + \frac{50}{q_p R_e} > 1, \quad (11) $$

where the Reynolds number $R_e = r^2 \Omega_k / \nu$. We call $q_{p,c}$ the maximum non-gap opening mass ratio obtained from equation (11). To summarize, the planet-to-star mass ratio needs to satisfy the condition

$$ q_p < q_{\text{gap}} = \min \left[ (C_{hs}/C_L)^2 h^3, 3h^3, q_{p,c} \right] \quad (12) $$

for planet outward migration at the inner disk edge.

For our fiducial disk model, Saturn, Uranus and Neptune are in the linear type I torque regime. Jupiter, being more massive, is in the type II gap-opening regime, and rebound fails to operate for Jupiter. As the timescale of type II migration is much longer than that of type I migration, and our study merely focuses on the time associated with rapid disk dispersal, we neglect the migration of Jupiter for simplicity. We assume that the giant planets have reached their present-day masses before the onset of final disk dispersal, so the above assessment of gap opening is based on their full masses. Nonetheless, our model holds as

long as proto-Jupiter has completed its main gas accretion and became more massive than Saturn before rebound operates.

We note that the above migration condition requires a disk with moderately high viscosity and aspect ratio. The simulations are performed with this fiducial disk setup unless otherwise stated. A too low $\alpha$ or $h$ will also cause Saturn to open a deep gap. Then, the picture changes as both Jupiter and Saturn undergo slow type II migration. We also conduct a subset of simulations in a disk with low viscosity and aspect ratio. The influence of these two parameters is investigated in Methods section 'Low-viscosity disks'.

## Planet-planetesimal disk interaction

The outer planetesimal disk exchanges angular momentum with the giant planets, resulting in the expansion of the planet's orbits with damped eccentricities and inclinations⁴². Such a planet-planetesimal disk interaction is often considered to be a trigger for the late giant planet instability, which played a crucial role in shaping the final architecture of the Solar System¹, typically taking place a few hundreds of megayears after the formation of the Solar System²¹. However, the orbits of the fully formed inner terrestrial planets are probably destabilized by such a late instability⁴³,⁴⁴, motivating the consideration of an earlier instability¹⁴.

During the gas-rich disk phase, a planetesimal with radius $R_{\text{plt}}$ experiences aerodynamic gas drag and its orbital decay timescale can be expressed as

$$ t_{\text{drag}} \approx 1.2 \times 10^8 \text{ yr} \left( \frac{R_{\text{plt}}}{10 \text{ km}} \right) \left( \frac{\Sigma_0}{10^3 \text{ g cm}^{-2}} \right)^{-1} \left( \frac{h_0}{0.04} \right)^{-1} \left( \frac{\rho_0}{1.5 \text{ g cm}^{-3}} \right) \left( \frac{a}{30 \text{ AU}} \right)^{9/4}, \quad (13) $$

where $\Sigma_0$ and $h_0$ are the gas surface density at 1 AU and $\rho_0$, is the internal density of the planetesimal. We find that as $t_{\text{drag}}$ is much longer than the gas disk lifetime, planetesimals in the proto-Kuiper Belt with radii larger than 10 km experience negligible radial drift during the gas disk phase.

## Numerical methods

We perform numerical simulations using a modified version of the publicly available N-body code HERMIT4⁴⁵ to study the evolution of multiplanet systems during gas disk dispersal. The code includes the planet-gas disk interaction by implementing the previously mentioned torque recipes¹⁹,³⁹. In addition, we run extended simulations to study the effect of a planetesimal disk on giant planet orbital evolution in a gas-free environment. These simulations are conducted separately using the open-source N-body code MERCURY with a hybrid symplectic and Bulirsch-Stoer integrator⁴⁶.

## Gas disk study

We investigate the evolution of planetary systems during the final gas disk dispersal phase in a statistical manner. The initial disk and planet conditions are listed in Extended Data Table 1. We consider three different planet configurations: all planets in nearly 2:1 resonances, 3:2 resonances, and a combination of 2:1 and 3:2 resonances. Different initial numbers of planets are also explored: $n = 4, 5$ and $6$. For each planetary configuration, we perform 1,000 simulations by Monte Carlo sampling the disk properties ($M_{pho}$, $\tau_d$ and $v_r$). Importantly, we consider both simulations with and without rebound to evaluate the efficacy of this mechanism. Planets feel the classic type I torques as in equation (5) when rebound is absent, whereas they feel torques including the one-sided components as in equations (6) and (7) when rebound is present.

We start the initial planet period ratios 5% higher than the exact resonant states. To further set up the initial conditions, we integrate the planets for 0.5 Myr with only the migration of the outermost planet
turned on, and the gas surface density unchanged. This ensures that the planets go into the desired resonant chains. The eccentricities and inclinations of the planets are assumed to follow Rayleigh distributions, where the scale parameters $e_0 = 2i_0 = 10^{-3}$. The other orbital phase angles are randomly selected between 0° and 360°. After the initial 0.5 Myr integration, we turn on the migration for all planets, and the disk starts to deplete according to the sampled disk properties. All parameter study simulations are terminated at 10 Myr when gas disks are fully dissipated.

Our goal is to demonstrate our new instability trigger, so we use two broad dynamical indicators to show that our simulated systems are indeed consistent with the global properties of the Solar System without attempting to match each detailed constraint. The first is the normalized angular momentum deficit$^{47}$: AMD = $\sum m_j \sqrt{\bar{a}_j (1 - \cos(i_j) \sqrt{1 - e_j^2}) / \sum m_j \sqrt{\bar{a}_j}}$, where $m_j$, $a_j$, $e_j$ and $i_j$ are the mass, semimajor axis, eccentricity and inclination of each giant planet $j$. The AMD increases with increasing orbital eccentricities and inclinations, because eccentric or inclined orbits have a lower $\bar{h}$-projected angular momentum than circular, coplanar ones at the same semimajor axes. The second indicator is the radial mass concentration$^{48}$: RMC = $\max(\sum m_j / \sum m_j [\log_{10}(a/a_j)]^2)$, which measures the degree of radial mass concentration in a given system, with higher RMC corresponding to a more tightly packed system.

Compared to Fig. 2, we show the effect of rebound on different initial potential giant planet architectures in Extended Data Figs. 1–3. The model parameter setups and numerical outcomes can be found in Extended Data Table 1. In addition, the AMDs of systems with initially four, five and six planets as functions of the disk parameters are presented in Extended Data Figs. 4 and 5. For instance, the upper panel of Extended Data Fig. 4 provides the outcome of simulations in which there were initially four giant planets in a chain of 2:1 resonances (run_A4R in Extended Data Table 1). For this setup, the Solar System analogues (black dots) are likely to form when $\dot{M}_{pho}$ is lower than $10^{-9} M_{\odot} \text{yr}^{-1}$.

## Planetesimal disk study

In addition to the above simulations only considering a gas disk, we also run simulations to account for an outer planetesimal disk. Such a disk is expected to continuously exchange angular momentum with the giant planets on a much longer timescale, motivating our inclusion of a planetesimal disk only after the gas disk is entirely depleted. We perform new sets of extended gas-free simulations for another 100 Myr using the MERCURY code, in which the initial orbital information of giant planets is adopted from the previous Hermite simulations at $t = 10$ Myr. For the purpose of illustration, we only perform a limited number of such extended simulations rather than extensive explorations in a multiparameter space.

We adopt the test particle approach to reduce the computational cost. The disk contains 1,000 test particles and each particle represents a swarm of real planetesimals at similar positions and velocities. The particles feel the gravitational forces from the planets, but the interactions between them are neglected. We assume that the planetesimal disk extends from 20 to 30 AU, with a surface density profile of $\Sigma_{\text{pl}} \propto r^{-1}$. The total disk mass is varied from 5, 10 and $20 M_{\odot}$. These test particles are initialized on nearly circular and coplanar orbits, and their eccentricities and inclinations follow Rayleigh distributions with $e_0 = 2i_0 = 10^{-3}$. We also test for numerical convergence using a total number of 500 test particles. The results are consistent between the explored resolutions.

## Secular mode analysis

Here we check the secular dynamics and eccentricity mode of the resulting planetary systems. As Jupiter and Saturn are the dominant mass contributors among Solar System planets, the secular dynamics of the system can be approximately traced by solving the Lagrange–Laplace equations for the Jupiter–Saturn pair$^{49}$. The eccentricities and longitudes of perihelia of Jupiter and Saturn can be written as:

$$ 
\begin{aligned}
e_j \cos \bar{\omega}_j &= M_{55} \cos(g_5 + \beta_5) + M_{56} \cos(g_6 + \beta_6), \\
e_j \sin \bar{\omega}_j &= M_{55} \sin(g_5 + \beta_5) + M_{56} \sin(g_6 + \beta_6), \\
e_s \cos \bar{\omega}_s &= M_{65} \cos(g_5 + \beta_5) + M_{66} \cos(g_6 + \beta_6), \\
e_s \sin \bar{\omega}_s &= M_{65} \sin(g_5 + \beta_5) + M_{66} \sin(g_6 + \beta_6),
\end{aligned} 
\quad (14) $$

where $g_j$ are the eigenfrequencies for the precession of perihelia, $\beta_j$ are the phase angles and $M_{jk}$ are the coefficients of the corresponding eccentricity amplitudes. The subscript $i$ refers to the eight Solar System planets, and 5 and 6 represent Jupiter and Saturn. We note that the perturbations from the other planets are relatively small compared to these from Jupiter and Saturn. Thus, these four amplitudes and two eigenfrequencies can be treated as a good approximation for the secular evolution of Jupiter and Saturn.

A proper excitation of Jupiter’s eccentricity mode, particularly $M_{55}$ to 0.044, is the most difficult property of the gas giant planets’ secular architecture to match; previous numerical simulations such as those of the Nice model or other alternatives generally yield lower values$^{22,23,49}$. Extended Data Figure 6 shows that the amplitude of the $M_{55}$ mode for the simulated planetary systems lies between 0.005 and 0.05. Simulations that only consider evolution in the depleting gas disk generate higher amplitudes compared with those that additionally account for the planetesimal disk. This is simply because the dynamical friction from the outer planetesimal disk continuously damps the eccentricities of the inner giant planets.

## Low-viscosity disks

To address whether the rebound mechanism can trigger a dynamical instability in the low-viscosity disks in which Saturn carves a deep gap, we conducted a new set of simulations focusing on the interactions with the gas disk, that is, no planetesimal disk included. We adopted a layered accretion disk model$^{50}$, where $\alpha$ represents the averaged global disk angular momentum transport efficiency and $\alpha_t$ corresponds to the local turbulent viscosity strength at the disk midplane. Therefore, $\alpha$ sets the gas surface density whereas the gap opening and planet migration are governed by $\alpha_t$. The two values are set to be equal in the fiducial disk model ($\alpha = \alpha_t = 0.005$). Here we also assume $\alpha = 0.005$ but vary $\alpha_t$ to test its influence on the system’s instability rate. We focus on the initial 3:2 resonant configurations of five and four planets (the same as run_B5R and run_B4R), and we specified the disk parameter distributions such that $\dot{M}_{pho}$, $\tau_d$ and $v_r$ are log-uniformed selected from $[10^{-9.5}, 10^{-9}] M_{\odot} \text{yr}^{-1}$, $[10^{5.5}, 10^6]$ yr and $[40, 120]$ AU Myr$^{-1}$, respectively. The adopted midplane $\alpha_t$ and the disk aspect ratio at the onset of the inner disk edge $h_{in}$ are illustrated in Extended Data Fig. 7.

In a low-viscosity disk environment, the outward-sweeping disk edge only directly affects the orbits of the ice giant planets. Both Jupiter and Saturn are in the gap-opening regime, and their slow type II migration is neglected during the rapid gas disk dispersal phase. Although Jupiter and Saturn cannot move out of resonance directly as the disk edge sweeps by, the subsequent instabilities triggered among the ice giants have a chance to propagate to the two larger gas giant planets (see the example in Extended Data Fig. 8).

As shown in Extended Data Table 2, orbital restructuring and system-wide instability are still common outcomes. A high fraction of systems experience instabilities during the first 10 Myr, most of which take place within 1–3 Myr simultaneously with the rapid disk dispersal (Extended Data Fig. 9a). In most circumstances, dynamical instabilities start among the ice giant planets and rapidly propagate to destabilize Jupiter and Saturn’s orbits (Extended Data Fig. 9b). The instability propagation rate is higher than 90% in general. This is true for both initial configurations with four and five planets, which means that the propagation of the instability from the ice giants to Jupiter and Saturn does not rely on the presence of an additional ice giant. Through the above exploration, we conclude that the rebound mechanism is a viable instability trigger even in low-viscosity disks, as long as at least one planet is below the gap-opening threshold.
Data availability

The data that support the plots within this paper and other find-
ings of this study are available at https://github.com/bbliu-astro/
solarsystem-rebound.git.

Code availability

The source code and simulation output for the model used in this study are available on reasonable request from the corresponding authors. The original version of the HERMIT4 N-body code is available on Sverre Aarseth's homepage https://people.ast.cam.ac.uk/sverre/web/pages/ nbody.htm.

31. Shakura, N. I. & Sunyaev, R. A. Black holes in binary systems. Observational appearance. *Astron. Astrophys.* **500**, 33–51 (1973).

32. Alexander, R. D., Clarke, C. J. & Pringle, J. E. Photoevaporation of protoplanetary discs – I. Hydrodynamic models. *Mon. Not. R. Astron. Soc.* **369**, 216–228 (2006).

33. Owen, J. E., Ercolano, B. & Clarke, C. J. Protoplanetary disc evolution and dispersal: the implications of X-ray photoevaporation. *Mon. Not. R. Astron. Soc.* **412**, 13–25 (2011).

34. Haisch, J., Karl, E., Lada, E. A. & Lada, C. J. Disk frequencies and lifetimes in young clusters. *Astrophys. J. Lett.* **553**, L153–L156 (2001).

35. Luhman, K. L., Espaillat, C., Hartmann, L. & Calvet, N. The disk population of the Taurus star-forming region. Astrophys. J. Suppl. **186**, 111–174 (2010).

36. Koepferl, C. M. et al. Disc clearing of young stellar objects: evidence for fast inside-out dispersal. *Mon. Not. R. Astron. Soc.* **428**, 3327–3354 (2013).

37. Hayashi, C. Structure of the solar nebula, growth and decay of magnetic fields and effects of magnetic and turbulent viscosities on the nebula. *Prog. Theor. Phys. Suppl.* **70**, 35–53 (1981).

38. Paardekooper, S. J., Baruteau, C., Crida, A. & Kley, W. A torque formula for non-isothermal type I planetary migration – I. Unsaturated horseshoe drag. *Mon. Not. R. Astron. Soc.* **401**, 1950–1964 (2010).

39. Liu, B., Zhang, X., Lin, D. N. C. & Aarseth, S. J. Migration and growth of protoplanetary embryos. II. Emergence of proto-gas-giant cores versus super Earth progenitors. *Astrophys. J.* **798**, 62 (2015).

40. Lin, D. N. C. & Papaloizou, J. On the tidal interaction between protoplanets and the protoplanetary disk. III. Orbital migration of protoplanets. *Astrophys. J.* **309**, 846 (1986).

41. Crida, A., Morbidelli, A. & Masset, F. On the width and shape of gaps in protoplanetary disks. *Icarus* **181**, 587–604 (2006).

42. Fernandez, J. A. & Ip, W.-H. Some dynamical aspects of the accretion of Uranus and Neptune: The exchange of orbital angular momentum with planetesimals. *Icarus* **58**, 109–120 (1984).

43. Agnor, C. B. & Lin, N. C. On the migration of Jupiter and Saturn: constraints from linear models of secular resonant coupling with the terrestrial planets. *Astrophys. J.* **745**, 143 (2012).

44. Kaib, N. A. & Chambers, J. E. The fragility of the terrestrial planets during a giant-planet instability. *Mon. Not. R. Astron. Soc.* **455**, 3561–3569 (2016).

45. Aarseth, S. J. *Gravitational N-Body Simulations* (Cambridge Univ. Press, 2003).

46. Chambers, J. A hybrid symplectic integrator that permits close encounters between massive bodies. *Mon. Not. R. Astron. Soc.* **304**, 793–799 (1999).

47. Laskar, J. Large scale chaos and the spacing of the inner planets. *Astron. Astrophys.* **317**, L75–L78 (1997).

48. Chambers, J. E. Making more terrestrial planets. *Icarus* **152**, 205–224 (2001).

49. Morbidelli, A., Brasser, R., Tsiganis, K., Gomes, R. & Levison, H. F. Constructing the secular architecture of the solar system. I. The giant planets. *Astron. Astrophys.* **507**, 1041–1052 (2009).

50. Liu, B., Lambrechts, M., Johansen, A. & Liu, F. Super-Earth masses sculpted by pebble isolation around stars of different masses. *Astron. Astrophys.* **631**, A7 (2019).

**Acknowledgements** B.L. was supported by a start-up grant of the Bairen programme from Zhejiang University, the National Natural Science Foundation of China (grant nos. 12173035 and 12111530175), the Swedish Walter Gyllenberg Foundation and the European Research Council (ERC Consolidator Grant no. 724687- PLANETESYS). S.N.R. is grateful to the CNRS's PNP program.

**Author contributions** S.N.R. and S.A.J. proposed this idea and initiated the collaboration. B.L. examined the feasibility and conducted numerical simulations. S.N.R. drafted the manuscript. All authors contributed to analysing and discussing the numerical results, and editing and revising the manuscript.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s41586-022-04535-1.

**Correspondence and requests for materials** should be addressed to Beibei Liu.

**Peer review information** Nature thanks Kleomenis Tsiganis and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_7.png" /></figure>

**Extended Data Fig. 1** | Metrics for surviving planetary systems of a subsample of our simulations in matching the Solar System, comparable to Figure 2. The simulations on the left included the rebound effect and those on the right did not. Each simulation started with our four present-day giant planets. In the top panels, the giant planets were initially placed in a chain of 3:2 orbital resonances. In the bottom panels, Jupiter and Saturn were initially in a 2:1 resonance and each other neighboring planet pair was in a 3:2 resonance.

Each symbol represents the outcome of a given simulation at *t* = 10 Myr. The color indicates the timing of the instability after the start of gas disk dispersal; pink systems did not undergo an instability (no collision and/or ejection). Circles and triangles correspond to systems with four and three or fewer surviving planets, respectively. The arrow gives the initial radial mass concentration of the system. The Solar System is marked as a red star for comparison.
<figure><img src="image_8.png" /></figure>

**Extended Data Fig. 2** | Metrics for surviving planetary systems of a subsample of our simulations in matching the Solar System, comparable to Figure 2. The simulations on the left included the rebound effect and those on the right did not. Each simulation started with the giant planets in a chain of 2:1 orbital resonances. In the top panels, we initially included only our four giant planets, but in the bottom panels, we added an additional ice giant at the start of the simulation. Each symbol represents the outcome of a given

simulation at $t = 10$ Myr. The color indicates the timing of the instability after the start of gas disk dispersal; pink systems did not undergo an instability (no collision and/or ejection). Diamonds, circles, and triangles correspond to systems with five, four, and three or fewer surviving planets, respectively. The arrow gives the initial radial mass concentration of the system. The Solar System is marked as a red star for comparison.
<figure><img src="image_9.png" /><figcaption>**Extended Data Fig. 3** | Metrics for surviving planetary systems of a subsample of our simulations in matching the Solar System, comparable to Figure 2. Both panels, left and right, included the rebound effect. Each simulation started with our four present-day giant planets plus two additional ice giant planets. In the left panel, the giant planets are in a chain of 3:2 orbital resonances. In the right panel, Jupiter and Saturn were initially in a 2:1 resonance and each other neighboring planet pair was in a 3:2 resonance.</figcaption></figure>

Each symbol represents the outcome of a given simulation at *t* = 10 Myr. The color indicates the timing of the instability after the start of gas disk dispersal; pink systems did not undergo an instability (no collision and/or ejection). Pentagons, diamonds, circles, and triangles correspond to systems with six, five, four, and three or fewer surviving planets, respectively. The arrow gives the initial radial mass concentration of the system. The Solar System is marked as a red star for comparison.
<figure><img src="image_10.png" /><figcaption>**Extended Data Fig. 4** | Outcomes of rebound simulations with initially four planets as a function of disk parameters: onset mass-loss rate, the disk dispersal timescale, and the rate of expansion of the inner cavity. Each</figcaption></figure>

simulation started with our four present-day giant planets in a chain of 2:1 orbital resonances (top panels), a chain of 3:2 orbital resonances (middle panels), or a combination of a 2:1 orbital resonance and an ensuing chain of

3:2 orbital resonances (bottom panels). The color bar corresponds to the system's angular momentum deficit (AMD). The circles with a grey edge color refer to the systems whose planets all survive in the end, while the black dots represent the Solar System analogs, defined as systems with four surviving planets in the correct order and their AMDs and RMCs are within a factor of three compared to those of our Solar System.
<figure><img src="image_11.png" /><figcaption>**Extended Data Fig. 5** | Outcomes of rebound simulations with initially five and six planets as a function of disk parameters: onset mass-loss rate, the disk dispersal timescale, and the rate of expansion of the inner cavity, comparable to Extended Data Fig. 4. Each simulation started with our four present-day giant planets plus one additional ice giant planet in a chain of</figcaption></figure>

2:1 resonances (1st row), a chain of 3:2 resonances (2nd row), or a combination of 2:1 and 3:2 resonances (3rd row), or started with our four present-day giant planets plus two additional ice giants in a chain of 3:2 resonances (4th row), or a combination of 2:1 and 3:2 resonances (5th row).
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 6** | Jupiter's eccentricity mode *M*<sub>55</sub> as a function of the period ratio of Saturn to Jupiter obtained in simulations with and without a planetesimal disk. The simulations without and with planetesimal disks are

plotted in triangles and circles, respectively, and the Solar System is marked as a star. Only systems that finish with four planets are shown here.
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 7 | Gap opening mass as a function of disk aspect ratio and midplane viscous $\alpha_t$.** The background color refers to the gap opening mass criterion from Equation (12), and the grey lines indicate the masses of four Solar System giant planets and 8 M⊕. The color symbols represent the disk setups we have explored in Methods section 'Low-viscosity disks', where the red symbols refer to the circumstances where only Jupiter opens a deep gap

(P0 is the same as the fiducial run in the main text), the magenta symbols correspond to the circumstances where both Jupiter and Saturn are in the gap opening regime, and the orange symbol indicates the circumstance that Jupiter, Saturn, Uranus, and Neptune open gaps while the additional ice giant planet with the lowest mass is in the non-gap opening regime. The values of $\alpha_t$ and disk aspect ratio parameters can be found in Extended Data Table 2.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 8 | An early dynamical instability triggered by the dispersal of the Sun's protoplanetary disk, assuming that the disk has a low viscosity.** The initial system consisted of five giant planets: Jupiter, Saturn, and three ice giants. The curves show the orbital evolution of each body including its semimajor axis (thick), perihelion and aphelion (thin). The black dashed line tracks the edge of the disk's expanding inner cavity. We do not follow the early evolution through the entire gas-rich disk phase, so the onset of disk dispersal

is set arbitrarily to be 0.5 Myr after the start of the simulation. The semimajor axes and eccentricities of the present-day giant planets are shown at the right, with vertical lines extending from perihelion to aphelion. The disk model is adopted from run_B5R_P4 where midplane turbulent strength ($\alpha_t = 10^{-4}$) is 50 times lower compared to the example shown in Fig. 1. The other disk parameters are: $M_{pho} = 1.1 \times 10^{-11} M_\odot$ yr$^{-1}$, $\tau_d = 5.0 \times 10^5$ yr, and $v_r = 42$ AU Myr$^{-1}$.
<figure><img src="image_15.png" /><figcaption>**Extended Data Fig. 9 | Cumulative distributions of delay times across three different suites of simulations.** On the left, the cumulative distribution of the time to the first instability regardless of which planets are involved. On the right, a cumulative distribution of the time delay between when the ice giant</figcaption></figure>

planets undergo orbital instability (typically occurs first), and when the gas giant planets undergo orbital instability. The black, blue and orange curves represent the simulations of *run_B5R_P0*, *run_B5R_P2* and *run_B5R_P4* in Extended Data Table 2.
<table><caption>Extended Data Table 1 | Initial conditions and statistical outcomes of the gas disk parameter study</caption><thead><tr><th>run ID</th><th>period ratio</th><th>position order</th><th>with rebound</th><th>p<sub>Nf=6</sub></th><th>p<sub>Nf=5</sub></th><th>p<sub>Nf=4</sub></th><th>p<sub>Nf&lt;=3</sub></th><th>p<sub>SSA</sub></th></tr></thead><tbody><tr><td>run A4R</td><td>all 2:1</td><td>J, S, U, N</td><td>Yes</td><td>0%</td><td>0%</td><td>73.6%</td><td>26.4%</td><td>7.7%</td></tr><tr><td>run A4N</td><td>all 2:1</td><td>J, S, U, N</td><td>No</td><td>0%</td><td>0%</td><td>100.0%</td><td>0.0%</td><td>89.7%</td></tr><tr><td>run B4R</td><td>all 3:2</td><td>J, S, U, N</td><td>Yes</td><td>0%</td><td>0%</td><td>23.1%</td><td>76.9%</td><td>0.9%</td></tr><tr><td>run B4N</td><td>all 3:2</td><td>J, S, U, N</td><td>No</td><td>0%</td><td>0%</td><td>84.5%</td><td>15.5%</td><td>0.0%</td></tr><tr><td>run C4R</td><td>2:1,3:2,3:2,3:2</td><td>J, S, U, N</td><td>Yes</td><td>0%</td><td>0%</td><td>45.1%</td><td>54.9%</td><td>6.4%</td></tr><tr><td>run C4N</td><td>2:1,3:2,3:2,3:2</td><td>J, S, U, N</td><td>No</td><td>0%</td><td>0%</td><td>77.9%</td><td>22.1%</td><td>58.8%</td></tr><tr><td>run A5R</td><td>all 2:1</td><td>J, S, U, U, N</td><td>Yes</td><td>0%</td><td>69.3%</td><td>26.9%</td><td>3.8%</td><td>1.7%</td></tr><tr><td>run A5N</td><td>all 2:1</td><td>J, S, U, U, N</td><td>No</td><td>0%</td><td>100.0%</td><td>0.0%</td><td>0.0%</td><td>0.0%</td></tr><tr><td>run B5R</td><td>all 3:2</td><td>J, S, U, U, N</td><td>Yes</td><td>0%</td><td>7.1%</td><td>29.4%</td><td>63.5%</td><td>1.9%</td></tr><tr><td>run B5N</td><td>all 3:2</td><td>J, S, U, U, N</td><td>No</td><td>0%</td><td>7.1%</td><td>29.4%</td><td>63.5%</td><td>1.9%</td></tr><tr><td>run C5R</td><td>2:1,3:2,3:2,3:2</td><td>J, S, U, U, N</td><td>Yes</td><td>0%</td><td>22.3%</td><td>40.6%</td><td>36.1%</td><td>2.5%</td></tr><tr><td>run C5N</td><td>2:1,3:2,3:2,3:2</td><td>J, S, U, U, N</td><td>No</td><td>0%</td><td>69.1%</td><td>10.7%</td><td>20.2%</td><td>3.3%</td></tr><tr><td>run B6R</td><td>all 3:2</td><td>J, S, U, U, N, N</td><td>Yes</td><td>3.1%</td><td>19.1%</td><td>39.7%</td><td>38.1%</td><td>1.4%</td></tr><tr><td>run C6R</td><td>2:1,3:2,3:2,3:2</td><td>J, S, U, U, N, N</td><td>Yes</td><td>15.4%</td><td>30.1%</td><td>37.9%</td><td>16.6%</td><td>0.3%</td></tr></tbody></table>

The first to fourth columns correspond to the name of the run, planet period ratio, position order (from inner to outer) and the option of simulations including rebound or not, where J, S, U, N are short for the planet with the mass of Jupiter, Saturn, Uranus, and Neptune, respectively. The disk parameters $M_{pho}$, $\tau_{d}$, and $v$ are log-uniformed selected from [10$^{-9.5}$, 10$^{-8.5}$] M$_{\odot}$ yr$^{-1}$, [10$^{5}$, 10$^{6}$] yr, and [20, 200] AU/Myr, respectively. The fifth to eighth columns show the probability of systems with a final number of planets $N_f=6, 5, 4$, and $\le 3$, respectively. The ninth column represents the probability of forming Solar System analogs, defined that systems survive with four planets in the right position order and their AMDs and RMCs are within a factor of three compared to the Solar System. The orbits of the current giant planets are adopted from Table 1 of Ref. $^{22}$.
<table><caption>Extended Data Table 2 | Instability statistics including low-viscosity disks</caption><thead><tr><th>run ID</th><th>[α<sub>t</sub>, h<sub>in</sub>]</th><th>position order</th><th>p<sub>ice</sub></th><th>p<sub>JS</sub></th></tr></thead><tbody><tr><td>run B5R P0</td><td>[5×10<sup>-3</sup>, 0.06]</td><td>J, S, U, U, N</td><td>93.8%</td><td>83.5%</td></tr><tr><td>run B5R P1</td><td>[10<sup>-3</sup>, 0.057]</td><td>J, S, U, U, N</td><td>93.3%</td><td>84.6%</td></tr><tr><td>run B5R P2</td><td>[5×10<sup>-4</sup>, 0.03]</td><td>J, S, U, U, N</td><td>23.3%</td><td>21.9%</td></tr><tr><td>run B5R P3</td><td>[2×10<sup>-4</sup>, 0.05]</td><td>J, S, U, U, N</td><td>61.6%</td><td>56.3%</td></tr><tr><td>run B5R P4</td><td>[10<sup>-4</sup>, 0.04]</td><td>J, S, U, U, N</td><td>42.6%</td><td>40.0%</td></tr><tr><td>run B5R P5</td><td>[10<sup>-4</sup>, 0.024]</td><td>J, S, U, P, N</td><td>95.4%</td><td>69.4%</td></tr><tr><td>run B4R P0</td><td>[5×10<sup>-3</sup>, 0.06]</td><td>J, S, U, N</td><td>77.8%</td><td>68.2%</td></tr><tr><td>run B4R P1</td><td>[10<sup>-3</sup>, 0.057]</td><td>J, S, U, N</td><td>75.2%</td><td>64.3%</td></tr><tr><td>run B4R P2</td><td>[5×10<sup>-4</sup>, 0.03]</td><td>J, S, U, N</td><td>23.3%</td><td>22.1%</td></tr><tr><td>run B4R P4</td><td>[10<sup>-4</sup>, 0.04]</td><td>J, S, U, N</td><td>10.0%</td><td>9.2%</td></tr></tbody></table>

The first and second columns provide the name of the run and the adopted α<sub>t</sub> and h<sub>in</sub>. The third column lists the planet position order, where in run B5R P5 the mass of one ice giant is chosen to be 8 M<sub>⊕</sub>. The instability probability is given by p<sub>ice/JS</sub>, where the subscript ice or JS represents that the dynamical instability occurs for ice giant planets or the instability spreads to the Jupiter and Saturn pair.