An early dynamical instability among the Solar System’s giant planets triggered by the gas disk’s dispersal

Zhejiang Institute of Modern Physics, Department of Physics & Zhejiang University-Purple Mountain Observatory Joint Research Center for Astronomy, Zhejiang University, 38 Zheda Road, Hangzhou 310027, China Department of Astronomy and Theoretical Physics, Lund Observatory, Box 43, SE–22100, Sweden Laboratoire dAstrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy SaintHilaire, 33615 Pessac, France Department of Earth and Environmental Sciences, Michigan State University, East Lansing, MI 48824, USA *To whom correspondence should be addressed; Email: bbliu@zju.edu.cn

giant planet instability 13,14 , thereby avoiding the possibility of de-stabilizing the terrestrial planets 15 and naturally accounting for the small mass of Mars relative to Earth and the mass depletion of the main asteroid belt 16 .
The evolution of gaseous protoplanetary disks is jointly controlled by viscous-driven accretion and wind-driven mass-loss 17 . In the early phase, the disk predominantly undergoes viscous evolution, as gas accretes at a considerable rate onto the central star for a few Myr. Meanwhile, the disk gas is photo-evaporated at a low rate by UV and X-ray photons from the central star 18 .
Once the disk's surface density drops sufficiently that the disk accretion rate 9 M is lower than the photoevaporating mass-loss rate 9 M pho , the disk begins to dissipate from the inside-out 19 . The disk is truncated at a characteristic radius R in ; the gas interior to R in rapidly viscously drains onto the star, and an inner hole is formed. At this point the main outer disk is exposed to direct stellar radiation. Gas-free inner cavities have been directly detected in several disks that appear to be in the act of dispersing 20 . As photo-evaporation continues, the inner cavity expands outward until the disk is gone. While disks may also be evaporated from the outside-in by high-energy photons in dense stellar environments 21 , inside-out dispersal appears to be an ubiquitous process in disk clearing 18 .
This dispersal can be characterized by three physical parameters: the onset mass-loss rate when photoevaporation takes over 9 M pho , the final gas disk dispersal timescale τ d , and the inner cavity expansion rate v r .
A dispersing protoplanetary disk has a direct effect on the orbital evolution of embedded planets. Angular momentum transfer between the planets' orbits and the gas leads to orbital migration, the shrinking or expansion of the planets' orbital radii 22 . Meanwhile, drag forces act to damp their orbital eccentricities and inclinations 23 . As the disk starts to disperse and an inner cavity is formed, planet migration at the inner disk edge is different than for planets fully embedded within the disk 24,25 . Embedded planets undergo "two-sided" gravitational torques from both the interior and exterior parts of the disk. In contrast, planets at the disk's inner edge only interact with gas on the outside of their orbits. As a result of larger-amplitude Lindblad and corotation torques, a planet below the mass threshold for opening a gap stops migrating inward at the disk inner edge 25,26 . If the disk inner edge is itself moving radially outward due to disk dispersal, then the planet may subsequently migrate outward along with the expanding inner edge (see Supple-mentary Information). This mechanism is termed "rebound", and was first applied in the context of the magnetospheric cavity on sub-AU scales to explain the architecture of close-in super-Earth planets 24,27 . In contrast, a planet massive enough to open a gap in the disk migrates interior to the cavity and its orbits is not affected as the cavity expands.
Here we show that the inside-out dispersal of the Sun's gaseous disk likely triggered a dynamical instability in the giant planets' orbits. Such an instability has been invoked to explain a number of the Solar System's features such as the giant planets' orbits 1,3 and those of Trans-Neptunian objects 28 , Jupiter's co-orbital Trojan asteroids 29 , and the giant planets' irregular satellites 30 . Previous dynamical studies have almost universally relied on interactions between the giant planets and an outer planetesimal disk to trigger the instability on timescales of tens to hundreds of Myr 5,7,9,31 . In the framework of the disk dispersal-driven instability, the primordial outer planetesimal disk -whose existence is inferred from the present-day population of Trans-Neptunian objects 32 -did not trigger the instability but still contributed to shaping the system's orbital architecture. This simulation starts from a system with five giant planets, including Jupiter and Saturn and three ice giant planets, one of which was removed during the instability 33 . The planets' initial orbits reflect an earlier phase of growth and migration, such that each pair of neighboring planets is in mean motion resonance: Jupiter and Saturn are in the 2:1 resonance and all other planet pairs are in 3:2 resonances. The system was integrated taking into account the mutual gravitational forces between the planets and the effect of the dispersing gaseous disk.
The expanding edge of the inner disk cavity does not affect all planets equally. Because Jupiter opens a deep gap around its horseshoe region and the corresponding corotation torque diminishes, the rebound is quenched (see Supplementary Information) and Jupiter simply enters the cavity as the inner disk edge sweeps by. When the disk edge approaches Saturn at t"0.6 Myr, the one-sided torque is strong enough to expand its orbit outward. Jupiter and Saturn exit their resonant state, as Saturn migrates outward with the expanding cavity, pushing into the orbits of the outer planets. The eccentricities of the ice giants become excited by this dynamical compression.
Saturn is left behind and enters the cavity at "9 AU at t"0.7 Myr. Meanwhile, the innermost (blue in Figure 1) ice giant planet becomes excited to the point 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 outward, while the innermost ice giant is eventually ejected towards interstellar space at t«0.8 Myr after a series of close encounters with Jupiter. After the disk's inner edge has swept past all of the planets, their final orbits are close to those of the present-day Solar System giant planets ( Figure 1).
It is worth noting that our simulations neglected the outer disk of planetesimals. Such a disk played a central role in triggering the instability and shaping the planets' final orbits in previous work 1,5,7,9,31,33 , which invoked a mass of "20 Earth masses within 30 AU. Within the context of a rebound-triggered instability, the giant planets' orbits are basically decoupled from the mass in an outer planetesimal disk. Gravitational interactions with a remnant planetesimal disk would modestly expand the orbital radii of Uranus and Neptune and decrease their eccentricities on a 10-100 Myr timescale.
A rebound-triggered instability is consistent with the Solar System's orbital architecture. To demonstrate this we carried out 1.4ˆ10 4 dynamical simulations like the one from Figure 1. We varied three different aspects of these simulations (see Table 1). First, we tested a wide range of plausible starting configurations for the giant planets' orbits, including varying the number of ice giants (2, 3, or 4) and their initial resonant states. Second, we used a Monte Carlo method to test the effects of the three physical parameters -the onset mass-loss rate 9 M pho , the disk dispersal timescale τ d , and the rate of expansion of the inner cavity v r -across the full range of astronomically-allowed values. Third, we ran each simulation twice; once including the effect of inside-out disk dissipation (with rebound) and once assuming the disk to dissipate smoothly at all radii (without rebound).
It is an ongoing challenge to use numerical simulations of the instability to concurrently match all of the Solar System's dynamical constraints 1,3,4,16,33,34 . Our goal is to demonstrate our new instability trigger, so we used 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 35 : where m j , a j , e j and i j are the mass, semimajor axis, eccentricity, and inclination of each giant planet. The AMD increases with increasing orbital eccentricities and inclinations, because eccentric or inclined orbits have a lower angular momentum than circular, coplanar ones at the same semimajor axes. The second indicator is the radial mass concentration 36 : RMC " max p ř m j { ř m j rlog 10 pa{a j qs 2 q, which measures the degree of radial mass concentration in a given system, with higher RMC corresponding to a more tightly packed system. , roughly 40% were unstable when rebound was ignored but more than 90% were unstable when it was included. In systems with a chain of 2:1 (Jupiter and Saturn) and 3:2 resonances (bottom panels), about 30% were unstable when the rebound effect was ignored 34 whereas roughly 80% went unstable when rebound was included. In addition, when rebound is taken into account, the distribution of final systems is more shifted from the initial low AMD and high RMC states towards final high AMD and low RMC states. A significant fraction of systems end up with an architecture similar to the Solar System (red star in Figure 2) with an instability timescale shorter than 10 Myr. This is not observed when rebound is not considered. The rebound-triggered instability occurs across all astrophysically-relevant disk parameter values. Figure 3 shows the distribution of outcomes of our simulations across the three parameters ( 9 M pho , τ d , and v r ) for one set of starting conditions, with five giant planets in a chain of 2:1 and 3:2 resonances (run C5R in Table 1). For this setup the systems with five surviving planets (black edgecolor) are likely to form when the onset mass-loss rate 9 M pho is low. We note that for the above  and three ice giants, each of which is 15 M ' and one of which (blue) was ejected into interstellar space during the instability. The curves show the orbital evolution of each body including its perihelion, semimajor axis, and aphelion (thin, thick, and thin curves, respectively). 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 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: 9 M pho "4ˆ10´1 0 M d yr´1, τ d "8.6ˆ10 5 yr, and v r "35 AU Myr´1.    Table 1). The color bar corresponds to the system's angular momentum deficit (AMD). The circles with black edgecolor refer to the systems whose planets all survive in the end, while the circles with magenta edgecolor represent the Solar System analogs, defined as systems with four surviving planets in the right order and their AMDs and RMCs are within a factor of three of our Solar System's values. Comparable figures starting from different orbital configurations are included in the Supplementary Information.

Methods
Disk model. The evolution of protoplanetary disks is thought to be driven by a combination of internal stress from the gas' viscosity and external stellar UV/X-ray photoevaporation 17,18,45 . At early times the influence of stellar photoevaporation is negligible and the disk evolves viscously in a quasi-steady state such that 9 M "3πνΣ, where 9 M is the gas disk accretion rate, Σ is the gas surface density, ν"αc s H is the gas turbulent viscosity 46 , α 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 9 M drops below the stellar photoevaporation mass-loss rate 9 M pho , disk dissipation is dominated by the photoevaporative wind originated from stellar high-energy radiation 19,47 .
When the thermal energy of disk gas is greater than its gravitational bounding energy beyond a threshold disk radius R in , the disk gas escapes as a wind. Gas interior to R in is rapidly accreted through a short, 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 19 . Under such a circumstance, the photoevaporating mass-loss rate is given by 18 9 where Φ represents the ionizing flux from the central star and R in is the initial size of inner cavity.
We note that photoevaporation is a complicated process coupled with radiative transfer and disk hydrodynamics. The mass-loss rate depends on the primary incident spectrum of the central star (EUV, FUV, or X-ray), grain abundance, and disk chemistry. The fiducial value of the onset massloss rate 9 M pho of 10´9 M d yr´1 in Eq. (1) is taken from models of EUV-driven stellar radiation 19 . If the photoevaporation is instead dominated by X-ray radiation, the mass-loss rate can be even higher 47 . We vary 9 M pho from 3ˆ10´9 to 3ˆ10´1 0 M d yr´1 in the parameter study.
In this photoevaporation-driven disk dispersal phase, the gas surface density can be written where Σ pho is the gas surface density when the stellar photoevaportation becomes dominant and τ d is the disk dispersal timescale. Noticeably, two different timescales are related to the gas disk dissipation. First, the disk lifetime is inferred to be 1´10 Myr with a median value of "3 Myr 48, 49 .
Second, since a small fraction of young stars have evolved disks being a transition between diskbearing and diskless, the disk clearing is expect to be very rapid, typically an order of magnitude shorter than the disk lifetime [50][51][52] . Based on these observational constraints, it is reasonable to assume that the disk dispersal timescale τ d is in a range between 0.1 and 1 Myr.
In this paper we assume a optically thin stellar irradiation dominant flaring disk 53 in which the disk temperature and aspect ratio can be expressed as where M ‹ and L ‹ are the stellar mass and luminosity, r is the disk radial distance, G is the gravitational constant, σ SB is the Stefan-Boltzmann constant, k B is the Boltzmann constant, µ is the gas mean molecule weight and m p is the proton mass. Ref 53 calculated that h " 0.033 at 1 AU by adopting L ‹ "1 L d . However, during the pre-main-sequence (age of ă a few Myr), the Sun underwent gravitational contraction and was far more luminous ("1´10 L d , along the Hayashi track) than its main-sequence stage. Hence, we choose a higher luminosity for a solar-mass pre-mainsequence star of L ‹ "4 L d and assume it remains constant during 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 notation X in indicates the quantity X evaluated at the inner disk edge R in when the stellar photoevaporation dominates. Based on the quasi-steady state assumption, we rewrite the gas surface density as where Σ in "4 g cm´2 for the fiducial disk parameters of 9 M pho "10´9 M d yr´1 and α"0.005.
As disk gas disperses from inside-out the inner cavity spreads outward. For simplicity, we adopt a constant v r to represent the cavity expansion speed. This quantity is not well constrained from a observational perspective. The cavity expansion rate should be principally determined by stellar ionizing flux, size and mass of the disk, and hence, it is not a full-independent variable. Nonetheless, when the dispersal time and disk mass are fixed, the expansion rate directly reflects the size of the disk. To a first order approximation, for a disk with an outer edge v r "R out {τ d "60 AU Myr´1. Due to the uncertainties in R out and τ d , the cavity expansion speed could plausibly span two orders of magnitude from 10 to 10 3 AU Myr´1.
In brief, the final photoevaportation-driven disk dispersal can be featured by three key model parameters: the disk dispersal timescale τ d , the onset mass-loss rate when photoevaporation dominates 9 M pho , and the inner disk cavity expansion speed v r .
Planet migration. We summarize the formulas of planet migration torques based on ref. 24 . When the planet is far from the inner edge of the disk, it feels the torque from disk gas on its both sides.
The two-sided torque is adopted from Eq. (49) of ref. 54 , which reads where Ω K is the Keplerian angular velocity, and q d "Σ g r 2 p {M ‹ and q p "m p {M ‹ are the mass ratio between the local gas disk and star and mass ratio between the planet and star, respectively. When the planet is at the disk edge, only the one-sided torque exists, which is given by where the first and second terms on the right side of Eq. (6) are the corotation and Lindblad torque components, and C hs "2.46, C L "´0.65 are the corresponding coefficients. In order to derive the corotation torque in Eq. (6), we assume that at the disk inner edge, the gas removal time t removal is faster than the gas liberation time in the planet horseshoe region t lib . We demonstrate 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 vis ; otherwise the gas would be accumulated. Thus, we have t removal Àt vis »x 2 hs {ν, where x hs "r p a q p {h is the half-width of the planet horseshoe region. The gas liberation time for a planet can be written as t lib "8πr{p3Ω K x hs q. Then t removal ăt lib is required for one-sided corotation torque in Eq. (6). The above condition is satisfied as long as t vis ăt lib , and one can obtain that qÀ8π{3α 2{3 h 7{3 »2.4ˆ10´4 by adopting α"0.005 and h " 0.07 (at the location of the 2:1 MMR with Jupiter). In other words, the planet needs to be no more massive than Saturn to fulfill the above condition. Therefore, the large amplitude one-sided corotation torque condition holds for our considered planet and disk circumstances.
For a general situation, the total torque can be expressed by interpolating the torques in the above two regimes where f " exp r´pr p´rin q{x hs s is a smooth fitting function, r p is the planet disk location, and f "1 (f "0) refers to the planet at the disk edge (far away from the disk edge). Although the above equations are derived analytically, we note that the termination of planet migration at the inner disk edge due to one-sided torque induced from a steep density transition is also obtained from hydrodynamic simulations 25,26 .
In order to satisfy the planet outward migration at the inner disk edge, the one-sided corotation torque needs to be larger than that the one-sided Lindblad torque. Thus, from Eq. (6) the planet needs to satisfy q p ăpC hs {C L q 2 h 3 . On the other hand, the above formulas are derived by assuming that the planet torques are in the linear-regime. When the planet is massive enough to progressively clear the local surrounding gas in its horseshoe region, an annular gap can be formed. In this case the disk feedback is non-trivial, and Eqs. (5) and (6) are not applicable anymore. Because the positive corotation torque diminishes due to gap formation, the planet beyond the gap-opening mass fails to undergo outward migration. The gap-opening first requires the planet's Hill sphere R H "pm p {3M ‹ q 1{3 r larger than the disk scale height H. Therefore, q p ă3h 3 is needed for torques in the linear regime. Dedicated hydrodynamic simulations further indicated that the planet mass needs to fulfill 55 3 4 where R e is the Reynolds number r 2 Ω K {ν. We name q p,c as the maximum non-gap opening mass ratio obtained from Eq. (8). To summarize, the planet-to-star mass ratio needs to satisfy q p ă minrpC hs {C L q 2 h 3 , 3h 3 , q p,c s for planet outward migration at the disk edge.
For our adopted disk model, Saturn, Uranus and Neptune are in the linear type I torque regime. More massive Jupiter is in the type II gap-opening regime, and rebound fails to operate for Jupiter. Since the timescale of type II migration is much longer than that of the type I and our study merely focuses on the time span of rapid disk dispersal, we neglect the migration of Jupiter for simplicity.
Numerical methods. In order to study the evolution of multi-planet systems during gas disk dispersal, we carried out numerical simulations using a modified version of the publicly available N-body code HERMIT4 56 . The code is designed for accurate calculations of both gravitational interactions among multi-bodies and their tidal interactions with the protoplanetary disk gas 24,57 . The gravitational forces among planets are integrated with a time-symmetric scheme and the Burdet-Heggie regularization, while the planet-disk interaction is taken into account by implementing the torque recipes mentioned previously. The code enables fast and reliable integration of the orbital evolution of multiple planets associated with close-encounters and collisions for millions of orbits, well suited for this study.
Parameter study. We conduct a high number of numerical simulations to study the evolution of planetary systems during the final gas disk dispersal phase. The initial disk and planet conditions are listed in Table 1 we have considered both simulations with and without rebound to evaluate the efficacy of this mechanism. Planets feel the classical type I torques 54 when rebound is absent, whereas one-sided torque component in Eq. (6) is accounted for when rebound is present. For instance, Figure 2 demonstrates the results of run B5R/B5N/C5R/C5N in Table 1.
The initial eccentricities and inclinations of the planets are assumed to follow Rayleigh distributions, where e 0 "2i 0 "10´3 is the scale parameter. The other orbital phase angles are randomly distributed between 0˝and 360˝. We start the initial planet period ratios 5% percent higher than the exact resonance values. For the first 0.5 Myr, we only turned on the migration of the outer most planet, and remain the cavity and gas surface density unchanged. This ensures the planets get into the desired resonant chains. After that we turn on the mentioned migration prescription and let the disk starts to deplete. Since we only focus on the early instabilities, simulations are run until t"10 Myr. It is noteworthy that the assumed values of initial eccentricities barely matter the subsequent evolution since the orbits would quickly get regulated by resonance trapping.
Similar to Figure 2, Figure 4 illustrates the RMC and AMD distributions among surviving planetary systems for initial four-planet systems starting with 3:2 resonances and a chain of 2:1 and 3:2 resonances. Figure 5 is for initial four and five planetary systems starting with 2:1 resonances, Table 1: Initial conditions for a parameter study. The second, third and fourth columns correspond to the initial number, period ratio and the order of the system (from inner to outer), where J, S, U, N are short for the planet with the mass of Jupiter, Saturn, Uranus and Neptune, respectively.
The disk parameters 9 M pho , τ d and v r are adopted from log-uniform distributions.
run ID number period mass the AMDs of systems with initial four, five and six planets as a function of disk parameters are presented in Figures 7, 8 and 9, respectively. The model parameter setups and numerical outcomes can be found in Table 1 and Table 2.

Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.           The color bar corresponds to the system's angular momentum deficit AMD. The circles with black edgecolor refer to the systems whose planets all survive in the end, while the circles with magenta edgecolor represent the Solar System analogs, defined as systems with four surviving planets with right orders and their AMDs and RMCs are within a factor of three compared to those of our Solar System.

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