© The Authors 2026

1 Introduction

Systems of multiple planets with sizes between Earth and Neptune orbiting with short orbital periods have been discovered in the hundreds by transit surveys and are now known to be one of the most common types of planetary systems in the galaxy (Fressin et al. 2013; Petigura et al. 2013). These Kepler-like systems are increasingly understood to be a cohesive population with a set of correlated defining statistics, describing the physical properties of the planets and their orbital architectures (Weiss et al. 2018; Van Eylen & Albrecht 2015; He et al. 2020). Parameters such as planet multiplicity, orbital spacing, and typical planet sizes are well-understood for planets with radii between 1 and 4 R_⊕ and periods <100 d. Formation models that aim to explain the origin of these systems must therefore focus on reproducing these statistics, rather than matching individual systems.

Briefly, the typical Kepler-like system can be summarized as follows: the planets are either super-Earths with radii between 1 and 1.6 R_⊕, or sub-Neptunes above 2 R_⊕, the latter of which have a large volatile envelope. At the time of formation these super-Earths and sub-Neptunes are a single population, with the radius difference emerging later as a consequence of core mass and insolation-dependent atmosphere loss (Wu 2019; Rogers & Owen 2021; Owen & Schlichting 2024). In terms of orbital architecture, many systems have multiple planets that are closely spaced in period but not near mean-motion resonance, while other systems have just a single transiting planet (Winn & Fabrycky 2015). A small number of systems have many planets locked in pairwise mean motion resonance; these are the “resonant chains”. Producing this diversity of orbital architectures for small planet systems is a major challenge of planet formation theory.

A scenario that incorporates near-universal dynamical instabilities of primordial resonant chains is thus far the most successful model to reproduce the observed planetary census (Pu & Wu 2015; Izidoro et al. 2017). In this model, resonant chains indeed form due to migration within the protoplanetary disk (Cresswell & Nelson 2006, 2008) but are catastrophically disrupted after the dissipation of the disk. During the instability, planetary orbits intersect and planets collide with each other, the star, or are ejected from the system entirely. Eventually, the instability ends and leaves behind a more widely spaced and dynamically excited system. The period-ratio distribution, which was strongly peaked near MMR during the disk phase, becomes smooth and nearly uniform, a close match to the data (Izidoro et al. 2017, 2021). The model also explains the observed eccentricity distributions, including the anomalously eccentric single-transiting systems (Van Eylen et al. 2019), patterns of intra-system uniformity in planet mass and period ratio (Goldberg & Batygin 2022; Lammers et al. 2023), and lower densities and smaller mutual inclinations of near-resonant planets (Leleu et al. 2024).

A major limitation of formation models incorporating N-body integrations is that they are restricted by computational constraints to the early stages of system evolution. Many transiting planet systems have Gyr or even 10 Gyr ages, corresponding to 10¹⁰−10¹² orbits of the innermost planet, a daunting challenge for even the fastest single-threaded integrators. Although substantial evolution takes place within the time covered by these models, much dynamical evolution can occur later. It is therefore not clear whether formation models that extend only to ∼ 100 Myr can be compared directly to data (Emsenhuber et al. 2021a).

A new, extremely powerful dataset is becoming available that promises to bridge this gap. The next generation of transit surveys, in particular TESS (Ricker et al. 2014), have discovered 10−20 multiplanet systems around stars younger than 1 Gyr, and including several younger than 100 Myr (Hu et al. 2025). While still a small minority of the known exoplanet population, these young and adolescent systems are a unique window into the formation environment of small planets, otherwise unobservable with current methods. In a population study, Dai et al. (2024) show that in the youngest systems (<100 Myr), near-resonant configurations are ubiquitous, appearing in >80% of observed multiplanet systems. This fraction drops with a timescale of ∼ 100 Myr to about 20% for mature systems. Similar results have been corroborated by studies that derive stellar ages from galactic kinematics. Hamer & Schlaufman (2024) demonstrate that systems with planet pairs near resonance, particular second-order resonance, have smaller galactic velocity dispersions indicative of lower ages.

On its face, this result appears to be direct confirmation of the “breaking-the-chains” hypothesis, because that idea predicts that systems emerging from their natal disks should be in resonant configurations that later break. On closer look however, many details are unresolved. For one, it has frequently been suggested that an additional mechanism outside of the N-body dynamics is necessary to break enough chains to match the observed distribution. Reshuffling of the system during disk dispersal has been proposed as one candidate process (Liu et al. 2017; Hansen et al. 2024; Pan et al. 2025), but this likely occurs too early to explain the abundance of near-resonance for systems with ages 20−100 Myr. In fact, it is not obvious what astrophysical process can occur over these timescales.

In this work, we demonstrate that resonant chains can unravel on timescales that match the observed decline in near-resonant frequency purely through N-body dynamics. In our model, planets grow out in a narrow ring (Batygin & Morbidelli 2023; Ogihara et al. 2024) and then through Type I migration form resonant chains that resemble observed resonant systems. However, in most of them, chaotic diffusion gradually drives the system away from a fixed point until the resonant protection is lost and the planets undergo a dynamical instability. Importantly, the rate of this diffusion is frequently very slow, with a long tail, reproducing the observed frequency of near-resonance. Previous work has demonstrated that systems that begin resonant but later went unstable can readily match the transiting planet population (Izidoro et al. 2017, 2021; Goldberg & Batygin 2022). Our primary goal in this work is to begin to understand how, why, and when this instability happens.

Section 2 describes our simulation setup and computational methods, and the results of these simulations in relation to the observed exoplanet population are presented in Section 3. We explore the details of the dynamics in Section 4 and describe the limitations of our results. Section 5 summarizes our findings.

2 Methods

We performed a suite of N-body simulations accounting for the growth of planets in the gaseous disk phase, planet-disk interactions, removal of the disk, and post-disk dynamical evolution of the system.

All simulations used the N-body code rebound with the trace integrator (Lu et al. 2024). trace is a hybrid timereversible integrator that switches from the classical WisdomHolman Hamiltonian splitting (Rein & Tamayo 2015) to direct Bulirsch-Stoer (BS) integration during close encounters. In comparison to the widely used MERCURY code (Chambers 1999) (and the similar rebound implementation mercurius (Rein et al. 2019a)), trace is uniquely well-suited to this problem for several reasons. First, the switching distance is dynamically computed at each timestep rather than at the beginning of the integration, as it is in mercurius. While technically breaking symplecticity, this choice guarantees that the switching criteria remains robust even as planets grow or migrate (Hernandez & Dehnen 2023; Lu et al. 2024). Secondly, trace properly resolves close pericenter approaches for highly eccentric orbits, which can arise during dynamical instabilities. Finally, during close encounters, trace uses a much faster BS integrator than the IAS15 integrator used by mercurius, resulting in a major speed-up for very violent systems dominated by close encounters without a meaningful loss of accuracy. These benefits allowed us to perform significantly longer integrations and uncover slower instabilities.

The initial conditions of the simulations are summarized in Table 1. Following Batygin & Morbidelli (2023), who argue that planetary growth should commence from rings in which solids have concentrated, we initialized simulations at time t=0 by placing embryos of mass 0.2 M_⊕ with semi-major axes drawn from a Gaussian distribution with mean 1 au and standard deviation σ around a central star of mass 1 M_⊙ and radius 1 R_⊙. Their eccentricities and inclinations were drawn from a Rayleigh distribution with mode 10⁻³ and 0.5 × 10⁻³, respectively, and their longitudes of pericenter, longitudes of ascending node, and mean longitudes were drawn from a uniform distribution from [0, 2π). All particles were set to be fully active. The timestep was set at the beginning of the integration to be 1/30× the period of the innermost embryo, and periodically adjusted if it decreased below 1/12× the period of a body due to inward migration. All bodies other than the star were set to have a physical radius corresponding to a density of 3g cm⁻³ and, if they approached closer than the sum of their radii, underwent perfect mergers conserving total mass and linear momentum.

We included the effects of Type I migration and damping on the embryos. The gas disk had a surface density of $$\Sigma (r,t)={\Sigma }_{0}exp(-t/{\tau }_{\mathrm{diss}}){\left(\frac{r}{r_0}\right)}^{-p}$$(1) where the power law disk slope is $$p=\frac{1}{2}({\alpha }_{\mathrm{o}\mathrm{u}\mathrm{t}}+{\alpha }_{\mathrm{i}\mathrm{n}})+({\alpha }_{\mathrm{o}\mathrm{u}\mathrm{t}}-{\alpha }_{\mathrm{i}\mathrm{n}})\tanh \left(\frac{r-r_0}{w}\right)),$$(2) Σ₀=1600 g cm⁻², τ_diss=1 Myr, α_out=3/2, r₀=1 au, w= 0.15 au, and α_in differs by simulation (Table 1). A constant aspect ratio was set to h/r=0.05. This parametrization represents a disk with a pressure maximum at 1 au, and a transition from a 3/2 power-law profile in the outer disk to a typically flatter profile in the inner disk. Such a profile is frequently found in simulations of disk evolution with magnetic winds that efficiently carve the inner disk (Suzuki et al. 2016; Kunitomo et al. 2020). The disk was dissipated exponentially on a timescale of 1 Myr, and then at t=5 Myr the disk was rapidly removed on a timescale of 10⁴ yr to represent disk photoevaporation. At t=5.1 Myr the disk and its corresponding forces were removed completely. Each system was then integrated to t=500 Myr under purely Newtonian gravity. The evolution of a typical simulation is shown in Fig. 1.

Embryos felt type I migration forces corresponding to Lindblad and corotation torques, computed using the equations of Paardekooper et al. (2010, 2011). These formulae incorporate a dependence on the local power law index of surface density and temperature and account for saturation of the corotation torque, for which we use a Shakura-Sunyaev turbulent α=10⁻⁴. We also included eccentricity and inclination damping forces according to the prescription of Cresswell & Nelson (2006, 2008), which account for reduction in damping at high inclination and eccentricity. For disks with α_in<α_out, a migration trap at 1 au appears for a specific range of planet masses. The trap can hold planets in the outer system until they grow large enough to saturate the corotation torque, and quickly migrate inwards. While not derived from a specific physical model, this simulation setup is intended to represent a disk with a broad gas pressure bump and dust trap that triggers planetesimal and embryo formation (Flock et al. 2017).

The disk is also believed to be truncated at short radii when magnetic forces from the star are more important than accretion ram pressure. These sharp decreases in surface density also act as planet traps where inward migration is stopped and strong outward torques are possible (Masset et al. 2006). While some previous work prescribed a falloff in surface density at the inner edge and let corotation torques provide outward migration (e.g., Izidoro et al. 2017), Ataiee & Kley (2021) showed that the standard torque prescription based on local power laws fails in these regimes. In fact, they found the trap to be much stronger than predicted, so that one planet at the inner edge is capable of halting the migration a resonant chain of several planets. We use a slightly modified version of their recommended prescription, in which we multiply the local migration rate τ_a=ȧ/a by a coefficient $$f_{\text{trans}}=1+\mathcal{A}(tanh\mathcal{R}-1)-sgn\left({\tau }_a\right)\mathcal{B}·(\mathcal{R}-C)exp(-{\mathcal{R}}^2)$$(3) where ℛ=(r−r_in)/𝒲, 𝒜=0.45, ℬ=20, 𝒞=0.226, 𝒲= 0.012 au, and r_in=0.113 au. This leads to a migration trap very close to 0.1 au, which is in very close agreement with the inner edge of Kepler systems being clustered near 10 d (Mulders et al. 2018). Figure 2 shows a migration map at t=0 of this disk with the three values of α_in used in this work.

3 Results

The vast majority of our simulation runs proceed in a similar fashion. Briefly, the embryos merge with each other near their starting location until they are massive enough to efficiently migrate inwards. Pairs frequently migrate convergently, either because the outer embryo is more massive than the inner or because the inner one is trapped at the inner disk edge. Embryos consistently capture into chains of mean motion resonances, most commonly first-order ones, although higher-order and coorbital resonances also appear. While dynamical instabilities are common in the disk phase, damping and migration rapidly recapture the systems into resonance and in nearly all of our simulation runs, the systems exit the disk phase in a mean-motion resonant configuration.

Following the disk dissipation, the system evolution takes one of two distinct paths. In the quiescent pathway, the resonant chain remains effectively unchanged over the entirety of the simulation, with purely periodic changes in orbital elements from resonant and secular evolution, reaching t=500 Myr with librating resonant angles. On the other hand, many systems experience a dynamical instability sometime between t=5.1 Myr and 500 Myr where the resonances become broken and the system is rapidly chaotic. In most cases, planetary orbits cross and planets collide very soon after the resonant angles begin circulating. The remaining systems are non-resonant, more widely spaced, more eccentric and inclined, and the planets are fewer and more massive. The quiescent pathway is interpreted as producing the rarely observed resonant chains such as TRAPPIST-1 and TOI-178 (Gillon et al. 2017; Leleu et al. 2021a), whereas the unstable track is responsible for the bulk of the known population of small planets in compact systems, which are non-resonant. The ability of a single process to reproduce both populations is the essence of the “breaking-the-chains” hypothesis originally proposed by Izidoro et al. (2017).

3.1 Planet accretion, migration, and resonance capture

The initial phase of planet accretion proceeds very quickly. Embryos quickly excite each other to crossing orbits and grow to several Earth masses. Typically, only 10−20 embryos remain by t=100 kyr. In models with a pressure bump, embryos that reach the minimum mass to saturate the corotation torque migrate inwards quickly, sometimes resonantly entraining smaller planets that remain interior to them. Migration occurs earlier in straight disks, with planet masses set by a balance of local growth rates and the timescale to exit the ring (Batygin & Morbidelli 2023). Migration always terminates at the inner edge located at 0.1 au, where the powerful planet trap prevents the chain from falling further towards the star.

The relative rates of migration for adjacent planets depends on their masses and the local disk profile, but convergent migration is frequently achieved because the inner planet is trapped by the inner edge directly or indirectly. Under these conditions, adjacent planet pairs capture efficiently into mean motion resonances, gradually building a large resonant chain anchored at the inner edge. Dynamical instabilities are repeatedly triggered during the disk phase either because migration dynamically overpacks the system or because of resonant overstability (Goldreich & Schlichting 2014). However, high gas densities in the inner disk allow the broken chain to quickly reform. By the end of the disk phase at t=5.1 Myr, 90−100% of systems are clearly resonant chains. The most common resonances are 5:4, 4:3, and 3:2 for the straight simulations and 4:3, 3:2, and 2:1 for 20flat and 40steep. Perhaps the biggest difference between the simulations is that those with a pressure bump frequently contain adjacent pairs of planets that are non-resonant: approximately 10% of pairs in 20flat and 40steep are further than 5% from a first- or second-order mean motion resonance, compared to <1% in the straight simulations. Additionally, the coupling between mass and migration speed leads to a more pronounced mass ordering in the straight simulations so that the innermost planets are almost always the most massive. In contrast, in simulations with a migration trap, all planets interior to 0.5 au tend to have a similar mass.

3.2 Population statistics

Figure 3 summarizes the statistics of the 20flat synthetic population at the end of the 500 Myr integrations, divided by whether the system experienced an instability after the end of the disk phase. Similar figures for the other simulation suites are presented in Appendix A. In these figures, all planets within 100 d have been included.

Some features are common to all simulation sets. Stable systems have more planets of lower mass, compared to unstable systems, and the period ratios of adjacent planets are dominated by resonant ratios. Most of these systems are extremely flat because their inclinations were damped by the protoplanetary disk and first-order resonances do not involve inclination, although a few have inclinations of ∼ 1^°. Eccentricities have been moderately excited by resonances to ∼ 0.05, a quantity which is set by the planet masses and balance of resonant excitation during migration vs. eccentricity damping in the disk (Terquem & Papaloizou 2019).

Unstable systems have qualitatively different features. In the closely packed small planet regime, systems dynamically relax by way of collisions and mergers, thus reducing the number of planets but increasing the mass of the remaining ones. Prior resonances are systematically broken by the powerful planet-planet interactions during orbit crossing. The only source of energy dissipation is through inelastic mergers, which are insufficient to reform mean motion resonances. Final period ratios are determined by dynamical stability and dependent on planet masses and multiplicity (Pu & Wu 2015). High mass runs thus have more widely spaced systems, from a median period ratio of 1.5 in 10straight to 2.49 in 40straight. Instabilities also leave eccentricities and inclinations in an excited state correlated with planet mass, although inelastic collisions sometimes leave systems anomalously flat or circular (Esteves et al. 2020).

The fraction of systems that go unstable before the end of the 500 Myr integrations differs dramatically between simulation suites, from a low of 32.5% in 10straight to a high of 92.5% in 40steep. Two patterns are clear from Table 2. First, for the same disk profile, simulations with higher M_tot have a higher instability fraction. This is expected because over the regime of interest, the final resonance that planets capture into during migration depends much more strongly on the local gas surface density than on the planet masses (Kajtazi et al. 2023). High mass planets therefore capture into the same resonances as low mass planets, but are much more closely packed dynamically (Petit et al. 2020; Goldberg et al. 2022). The second clear trend is that, for the same M_tot, systems that formed in disks with a migration trap have a higher instability fraction. For example, 57.5% of systems in 20flat go unstable, compared to 35% in 20straight. Likewise, 40steep had a 92.5% unstable fraction compared to 80% in 40straight. We explore the possible origins of this trend in Section 4.1.

3.3 Timing of the instability

The time required for a resonant system to go unstable, and thus transform into a non-resonant system, is of particular importance thanks to recent results of the population characteristics of young systems and an increasing number of clear trends in planetary system architectures with age. Our goal in particular is to determine whether purely N-body dynamics can reproduce the rapid decline in the fraction of systems observed near mean motion resonance on a timescale of approximately 100 Myr (Dai et al. 2024). Figure 4 shows the cumulative distributions for the time of instability of our systems, separated by simulation suite. The “instability” is defined as the first time that a planet with period <100 d at the end of the disk phase had its mass or semi-major axis change by more than 5%, and thus does not include instabilities that occur purely in the outer system that do not propagate to the inner planets.

We found that the instability time distribution is wellmodeled by a Weibull distribution (Weibull 1951), the cumulative function of which is given by $$\mathcal{W}\left(t_{\mathrm{inst}}\right)=1-exp(-{\left(\frac{t_{\mathrm{inst}}-t_0}{\mathcal{T}}\right)}^k).$$(4) Here, t₀ represents the earliest possible time for an instability, defined to be the disk dissipation time t₀=5.1 Myr, 𝒯 is a scale parameter, and k is a shape parameter. While not derived from a specific physical model, the Weibull distribution is widely used in failure rate statistics to describe the time of failure of a complex system made of multiple components¹ and its parameters have a straightforward interpretation. The scale parameter 𝒯 is defined simply as the time at which a 1/e fraction of systems have failed. Furthermore, in this context, one speaks of the hazard function, h(t)=−𝒲^′(t) /(1−𝒲(t)) which describes the instantaneous rate of failure for systems still functioning at time t. For the Weibull distribution, the hazard function takes the form $$h\left(t_{\mathrm{inst}}\right)=\frac{k}{\mathcal{T}}{\left(\frac{t_{\mathrm{inst}}-t_0}{\mathcal{T}}\right)}^{k-1}$$(5) which is decreasing in time for k<1 and increasing for k>1 (for k=1 the Weibull distribution simplifies to the exponential distribution). Essentially, k<1 represents a sort of “infant mortality” wherein many systems fail early, but older ones are increasingly robust and have a reduced likelihood of failure.

We fit the simulated instability time distribution of each simulation suite to a Weibull distribution by maximizing the total likelihood. Systems that remained stable over the duration of the integration were assigned likelihood 1−𝒲(500 Myr). Uncertainties in the distribution parameters and sample draws from the fit were computed by treating the estimated Hessian matrix returned by the BFGS optimizer implemented in scipy as a covariance matrix (Nocedal & Wright 2006).

The Weibull parameters best fits and uncertainties are given in Table 2. Draws from the fit are also shown in Fig. 4 along with the simulation-derived empirical distribution function. The instability timescale 𝒯 spans several orders of magnitude across simulation suites. There is a clear dependence on total mass that is steeper from M_tot=20 M_⊕ to 40 M_⊕ than from 10 M_⊕ to 20 M_⊕, suggesting that these resonant chains span a stability limit driven by the overlap of resonances (Quillen 2011; Petit et al. 2020). The simulations with a pressure bump also have 𝒯 that is a factor of ∼ 5 times smaller than the corresponding ones without.

In contrast, the shape parameter k is exceptionally consistent across our simulations. Indeed, four of the five suites are within one sigma of k=0.37 and 20straight is 1.5 sigma away from that value. The exact interpretation of this uniformity is not clear, but it suggests that an underlying mechanism sets the instability time across all of our simulations and differs only in its strength.

3.4 Simulated observations

In order to directly compare our results to the observed transiting planet sample, we applied a synthetic observation algorithm to determine which planetary systems would be observed by a transiting survey. To do this, for each simulated system, we simulated observations from 1000 random viewing angles and recorded which planets would transit the star using the CORBITS package (Brakensiek & Ragozzine 2016), constructing a large synthetic observation catalog.

Our integrations only continued to 500 Myr due to computation time constraints, but typical Kepler systems are much older. We take a fiducial age of ∼ 5 Gyr. As resonant chains continued going unstable until the end of our simulations, we expect that even fewer will be resonant at late times. We extrapolated the instability fraction to 5 Gyr by sampling the fits to the Weibull distribution and computing the cumulative value at 5 Gyr, taking the median and 1 σ ranges for each simulation suite. Then, we created a new observed catalog by resampling the previous one according to the new unstable fraction. For example, 57.5% of the 20flat simulations went unstable by 500 Myr, but extrapolating the Weibull fit leads to $87.6_{-6.7}^{+7.2}\%$ unstable systems at 5 Gyr. We generate a new observed catalog where 87.6% of systems are drawn from unstable systems and the remainder from the stable systems, while still sampling random viewing angles. Previous work performed similar resampling (Izidoro et al. 2017, 2021), but only in an ad hoc way to match the observed exoplanet demographics rather than deriving it from the instability statistics.

The synthetic observations are particularly important because the stable and unstable systems are observed in different ways. Stable systems are very flat, leading to typically an all-or-nothing discovery outcome depending on their sky-plane inclination. We therefore expect observed resonant systems to be mostly representative of their true configuration, at least in the inner system. In contrast, unstable systems have moderate inclination, so it is uncommon to observe the entire system in transit. Frequently, only a single planet is observed, or a planet is missed in between two observed ones, biasing both the observed multiplicity and period ratio distribution.

Figures 5 9 show the inferred multiplicity and period ratio statistics that would be observed by a transiting survey, split into the stable and unstable systems and the resampled catalog in purple. We also plot the multiplicity and period ratio distributions from candidate and confirmed planets detected by the Kepler mission, taken from the Exoplanet Archive². To compare directly with simulations, only planets with periods <100 d are included and systems with a giant planet (defined as radius >4 R_⊕) are excluded.

As seen in previous studies, planetary systems that underwent a dynamical instability closely resemble the observed transiting planet population. The high mutual inclinations from unstable systems leads to a large number of systems with only a single planet in transit, explaining the so-called Kepler dichotomy between single and multi-transiting systems (Johansen et al. 2012). Furthermore, the observed period ratio distribution depends strongly on the total mass, but matches the observed period ratio distribution.

The strength of the produced dichotomy depends strongly on the initial mass. Higher mass planets systems generate larger mutual inclinations during the instability, which directly strengthens the dichotomy (Izidoro et al. 2017; Millholland et al. 2021). Furthermore, in our simulations, low mass systems have a smaller instability fraction and so the final system catalog is dominated by high multiplicity resonant chains. We find that the 40steep and 40straight simulations most closely reproduce the observed Kepler multiplicity distribution of all of the tested configurations. The lower mass 20flat and 20straight simulations still produce a dichotomy, although it is too weak to match the nearly 80% of single transiting systems observed by Kepler. Finally, the lowest mass 10straight simulation, with its low inclinations and low rate of instability, does not show a distinct dichotomy.

The orbital spacing of planetary systems post-instability is much broader for more massive systems. The observed Kepler period ratio distribution is most closely matched by the unstable systems that started with 20 M_⊕ (orange curve, Figs. 5 and 8). However, the overall period ratio distribution is strongly influenced by the stable systems because, despite their rarity, they are frequently observed as multiples. In the end, a very close match to the observed distribution is obtained in the 40straight run, although 20flat is a satisfactory match.

Now that we have confirmed that the final states of these simulations are an acceptable match to observed mature systems, we consider the new constraint of the time evolution of the fraction of systems near mean-motion resonance. In the absence of robust dynamical information for most systems, Dai et al. (2024) defined “near-resonant” systems as having a planet pair with −0.015 ≤ Δ ≤ 0.03 for first-order resonances and −0.015 ≤ Δ ≤ 0.015 for second-order resonances (Δ is defined in Section 4.2). This period ratio-based definition is different than our stable/unstable distinction because some instabilities can occur without breaking all resonances, and because unstable systems occasionally produce pairs that are serendipitously near MMR in period ratio. For the latter reason, even a catalog of entirely non-resonant systems will be measured as having a 10−20% incidence of resonance by the empirical criterion.

Figure 10 shows this empirical criterion applied to our five simulation suites as a function of time including the extrapolated instability fractions. The results follow similar trends to Fig. 4. The 10straight and 20straight simulations have too few instabilities to reproduce the 38% near-resonance incidence for systems between 100 Myr and 1 Gyr. On the other hand, 40straight and 40steep have very early instabilities that result in few detected resonant systems at 50 Myr, in conflict with the data. The best result is obtained in the 20flat simulations, which remain mostly resonant at 50 Myr but steadily continue breaking out of resonance through the entire integration. Furthermore, the incidence of resonance in this model extrapolated to 5 Gyr is very close to the measured value.

To summarize briefly, the 10straight and 20straight are a poor fit to the data, being consistently too low mass, too resonant, and too compact. In contrast, 20flat, 40steep, and 40straight qualitatively reproduce many aspects of the data well. However, each of them fails at least one constraint. 20flat has under-excited inclinations, leading to too few singletransiting systems and too compact observed systems. 40steep and 40straight are a better match to the Kepler statistics but have too many early instabilities. It is likely that a mixture of these models would be an even better fit (Goldberg & Batygin 2022). Generically, there is a wide variety of protoplanetary disks and the resulting planetary systems inherit this diversity (Emsenhuber et al. 2021b).

4 Discussion

4.1 Why some systems go unstable

Figure 11, showing a typical example of the the 20flat simulation suite, gives a clue as to how these initially resonant systems can slowly unravel. Upon exiting the disk phase, this system contains two resonant chains separated by a large gap of 3.45 in period ratio. The inner chain is relatively widely spaced, with pairs in the 2:1, 3:2, and 4:3 resonances. The outer chain is far more compact, containing pairs in 5:4, 7:6, 8:7, 1:1 (coorbital), 6:5, and 5:4 resonances. All resonances are librating upon exiting the disk phase, although the outermost pairs have fairly high libration amplitude. From then until t=30 Myr, it is clear that pure N-body dynamics are causing the system to chaotically diffuse in action space, increasing the amplitude of libration of some pairs while decreasing that of others. Around t=30 Myr, the outermost resonance rapidly grows in libration amplitude and then begins sporadically circulating, which starts to gradually excite the libration of the remaining resonances in the chain. As the eccentricities grow, new resonant modes appear for the co-orbital pair (near t=80 Myr) and the 6:5 resonance enters circulation. At around t=110 Myr, the growth in libration amplitudes suddenly accelerates and the entire outer chain breaks simultaneously. Without the protection against close conjunctions provided by the resonant libration, the outer planets excite each other onto crossing orbits and collide, leaving just two survivors.

During nearly this entire process, the inner four-planet chain remains effectively unchanged. All resonances librate at constant amplitudes and chaos is not apparent. However, at the moment that the outer chain becomes unstable and its planets reach e ≳ 0.1, the libration amplitudes in the inner chain jump to a large amplitude. From there, they quickly grow until they reach circulation, the inner chain becomes unstable, and the middle two planets merge. The final system is non-resonant and excited in eccentricity and inclination.

A particularly remarkable fact is that the inner system appears to be stable on its own. While it receives extra secular precession from the outer system, the separation of timescales between secular and resonant motion ensures that it precesses as a solid object without affecting the mean motion resonances (Petit et al. 2017; Tamayo et al. 2021). The outer cluster of planets is unlikely to be observed in transit, especially after gaining significant mutual inclination during the instability, yet plays the key role in the evolution of the entire system. We refer the interested reader to the recent work by Ogihara & Kunitomo (2026) for details on the importance of secular perturbations from an outer system to instabilities of inner resonant chains.

This two-phase, outer to inner mode is the most common way that resonant chains break in our simulations, but other pathways are possible. Some outer chains go unstable but the inner chain is unaffected (these are counted as stable in our analysis). In a few cases, the instability clearly arose from the inner planets and spread to the outer system. Finally, in other cases it is ambiguous which set of planets is the main source of diffusion.

In general, it is not clear why some resonant chains go unstable, although one instability in the three-planet case is understood (Pichierri & Morbidelli 2020) and heuristics for longer chains are known (Goldberg et al. 2022). We hypothesize that the higher incidence of instability in the simulations with a pressure bump is because low masses and slow migration in the 1 au region consistently create a very compact resonant chain. Although this outer chain is stable to chaotic overlap of adjacent first-order resonances (Deck et al. 2013), a large network of higher-order resonances can drive diffusion away from the resonant center and, once the instability is eventually triggered, an impulsive kick can be passed to the inner planets. Instabilities rarely occur in the inner chain because planets massive enough to saturate the corotation torques and reach the inner system typically capture into wider resonances which are less prone to resonant overlap and chaos. In contrast, in purely power-law disks (e.g., Fig. 12), there is no clear inner/outer system separation. The outermost planets have very low masses and the entire chain has very low resonant libration amplitude.

4.2 Remaining resonant systems

A fundamental prediction of the breaking-the-chains model is that the surviving stable systems should closely resemble the observed resonant chains (Bitsch & Izidoro 2024). The simplest comparison is to the orbital period architectures; we focus on period ratios of adjacent planets and proximity to two- and three-body resonances. For each adjacent planet pair with periods P_in<P_out, we describe proximity to two-body resonance by the parameter $$\Delta =\frac{P_{\text{out}}}{P_{\text{in}}}\frac{p}{p+q}-1$$(6) where we consider only the closest first- or second-order p+q: p resonance. To estimate proximity to zeroth-order three-body resonance, we followed the technique of Goldberg & Batygin (2021). For each planet triplet in the system with mean motions n₁=2 π/P₁>n₂=2 π/P₂>n₃=2 π/P₃, we computed the value $$B=p{n}_1-(p+q)n_2+q{n}_3$$(7) for all positive coprime integer values of p and q with p+q< 14 and normalize it by 〈 n〉=(n₁+n₂+n₃)/3. Then, only the smallest value of |B|/〈n〉 was kept for each system.

Figure 13 shows these statistics for each simulation suite after applying the synthetic observations filter from Section 3.4. We also plot these statistics for the observed 3+ planet resonant chain sample identified by Kajtazi et al. (2023) with the addition of TOI-1136 (Dai et al. 2023) and again removing systems with giant planets. Simulations suites 20flat and 40steep closely match the observed period ratio distribution, in particular reproducing the strong abundance of 3:2 pairs. In contrast, the simulations with pure power law disks produced too many systems with 5:4 and 4:3 pairs compared to the observed population. However, our simulated planet pairs are much closer to exact resonant commensurability, as expressed by Δ (top right panel), than observed. Nevertheless, the proximity to three-body commensurability is very close the observations, especially for 20flat and 40steep. In other words, many systems show a ∼ 1% positive deviation from two-body resonance while being clearly in three-body resonance. The inability of Type I migration to reproduce this feature of the observations is well known (e.g., Charalambous et al. 2023, and references therein). Tidal dissipation, which we do not include in our simulations for simplicity, can increase Δ while preserving B and potentially explain this mismatch (Papaloizou 2015; Goldberg & Batygin 2021; Charalambous et al. 2023).

Finally, we note that many of these chains have more detailed mass and eccentricity constraints measured by transit timing variations (TTVs), but a detailed comparison to these inhomogeneous analyses is left to further work.

4.3 Neglected effects

Our scheme neglects physics that may be important in planet formation. Real disks are turbulent and these stochastic eddies can suppress capture into resonance during migration (Batygin & Adams 2017; Chen et al. 2025). The true impact will depend on the location of the magneto-rotational instability turn-off (Cevallos Soto & Zhu 2025). Nevertheless, this choice appears to be justified. Izidoro et al. (2017) considered the effect of stochastic kicks caused by turbulence in their model with a much larger α than we use and found that the results were nearly indistinguishable from non-turbulent runs.

Disks are also not static or uniformly depleting in time. Our surface density profile is inspired by the results of Suzuki et al. (2016), who consistently found a flat or increasing surface density in the inner disk. However, the position of the peak in surface density evolves outwards in time as the disk dissipates. Ogihara et al. (2024) performed N-body simulations with Type I migration directly using these disk models and found similar super-Earth systems to ours. It is likely that the details of the disk profile are most important when planets are ∼ 1 M_⊕ and the corotation torque can act as a mass filter. This occurs in our simulations around 100 kyr to 1 Myr, when our disk profiles closely match the results of Suzuki et al. (2016).

We also ignored interior structure and atmospheric evolution of the forming planets by always assuming their density to be 3 g cm⁻³. In reality, young planets are expected to have puffy, extended envelopes that slowly cool and contract. The connection between this and widespread giant impacts after the gas disk dissipation is not clear. Although SPH simulations of super-Earth and sub-Neptune collisions imply relatively little atmospheric loss from the mechanical shock (Denman et al. 2020), Biersteker & Schlichting (2019) proposed that post-collision heating would thermally expand the envelope and lead to complete loss via Parker winds. Izidoro et al. (2022) showed that this effect, along with expanded radii of water-rich planets that accreted outside of the ice line, can even produce a radius gap between super-Earths and sub-Neptunes.

The radii themselves also have dynamical importance. Large planets are more likely to collide, terminating instabilities before too much dynamical excitation. But, it is not clear what the proper radius to use is, because overlap of atmospheres may not be sufficient for collision and merger (see discussion in Kimura et al. 2025). We therefore chose 3 g cm⁻³ as a compromise between core-determined and envelope-determined collision and merger.

4.4 Possible missing components

We return to the question of whether pure N-body dynamics of planetary systems are sufficient to break almost all resonant chains and generate the observed small transiting planet sample. It appears that many features of the exoplanetary census are indeed reproduced by our simplified model, at least for certain choices of parameters. Specifically, the observed period ratio and resonant fraction lie between our 20flat and 40steep models, suggesting that a range of available solids to form rocky planets is necessary to match the data.

Interestingly, the strength of the observed Kepler dichotomy is not reached even in our highest mass simulations. In other words, violent instabilities that produce enough single transiting planets also tend to make too widely spaced systems (e.g., Fig. 6) and vice versa (Fig. 5). Our 40straight suite was a good match to both the observed multiplicity and period ratio distribution (Fig. 9) but predominately makes planets of mass ∼ 10 M_⊕ (Fig. A.4), seemingly too high compared to the typical Kepler super-Earth or sub-Neptune of ∼ 3 M_⊕ (He et al. 2020; Rogers & Owen 2021). This problem also appears in Izidoro et al. (2021), who tested multiple pebble accretion models with different mass flux. Their sole model that matches the observed multiplicity distribution (Model-I, t_start=0.5 Myr in their notation), has a typical planet mass of 10−20 M_⊕.

One potential way to bridge this gap is an additional source of inclination excitation. Spalding & Batygin (2016) showed that mutual inclinations can be excited in multiplanet systems through the quadrupolar precession of a primordially misaligned host star. During spindown, a secular resonance can be excited that tilts the outer planet to a non-transiting configuration. Although this process can itself trigger a dynamical instability (Spalding et al. 2018), the quiescent version appears a promising mechanism to raise mutual inclinations in some systems without requiring large masses or orbital spacings between planets.

Orbital eccentricities are more difficult to compare to models because they are not directly measurable from transit observations. Transit timing variations provide more information on the eccentricity through the mean-motion resonant dynamics. Yee et al. (2021) show that TTV systems have eccentricities several times lower than that required for stability, suggesting some amount of dissipation that damped eccentricity following a giant impact phase. However, TTV systems represent a biased sample of near-resonant planets, which, in the context of the breaking-the-chains model, are less likely to have undergone an instability that excites eccentricity (Goldberg & Batygin 2022).

Transit duration measurements, on the other hand, can constrain the average eccentricity of an entire population of planets and have indicated that Kepler multiplanet systems have typical eccentricities of ∼ 0.05 (Van Eylen et al. 2019; Mills et al. 2019). Gilbert et al. (2025) refined this strategy by dividing the population in planet radius bins and find that small planets (≲ 4 R_⊕) have mean eccentricity of 0.06 in single transiting systems and 0.03 in multitransiting systems. After applying the synthetic observational bias, we find mean eccentricities of 0.06 and 0.05 in single and multitransiting systems in the 20flat simulations and 0.09 and 0.07 in the 40steep simulations. This is consistent with a small degree of eccentricity damping over the lifetime of the system after a dynamical instability, such as that provided by tidal dissipation (Delisle & Laskar 2014).

This model predicts that the youngest observed systems should be predominantly chains of planetary deep in resonance. Although Dai et al. (2024) identify a surplus of systems near resonance according to an orbital period based criterion, most of these systems do not have sufficient dynamical constraints to determine whether they are actually resonant according to a dynamical definition that requires libration of the resonant angle within a separatrix-bounded region. Hu et al. (2025) studied several of these systems in detail and find several planet pairs that, while classified as near-resonant in terms of period ratio, are sufficiently distant from exact commensurability that their resonant angles are likely to be circulating. As near-resonance in our synthetic systems typically implies resonant libration, these intriguing observational findings imply additional physics not captured in our simulations during the disk phase or shortly after disk dissipation that can modify the resonant state.

Finally, we note that the very steep decline in the observed incidence of resonance between the first and second bin of Fig. 10 is not well reproduced in our simulations. While Dai et al. (2024) present a compelling argument for a population trend, there is reason to believe that the quantitative details are not set in stone. For one, Dai et al. (2024) did not attempt to debias the sample. Young transiting planets present numerous difficulties in detection and confirmation due to stellar activity and inhomogeneous distribution on the sky. Furthermore, near-resonant planets experience transit timing variations that smear transits after phase-folding and impede detection (Leleu et al. 2021b). Neither of these effects nor their combination has been carefully treated in the literature. We expect that the PLATO mission, which is intended to provide a robust statistical sample as well as accurately measured ages for many planet hosts, will greatly refine the observational constraints and potentially clarify whether there is in fact a preferred time of instability.

5 Conclusion

In this work, we performed a large suite of very long duration N-body integrations incorporating planet growth in a ring, Type I migration, and dynamical instability. Planetary systems are almost entirely resonant chains at the end of the disk phase, which gradually break over the course of our 500 Myr simulations. We compared the final systems to observed transiting planet data after applying a synthetic observation filter.

We present the following conclusions:

• We find that certain combinations of disk surface density profiles and initial solid mass can reproduce the Kepler sample demographics as well as the observed multiplanet system resonant fraction as a function of age. This scenario requires that 90−95% of initially resonant systems undergo an instability. The remaining resonant systems closely resemble known resonant chains in orbital architecture with deviations being explained by tidal dissipation.

• Pressure bumps and their resulting migration traps are effective ways to increase the rate of instabilities compared to purely power-law disks. There is frequently an outer resonant chain that is unstable on long timescales and passes its instability to the inner chain, which would have been stable on its own. These outer companions are not easily observable with current transiting or radial velocity techniques, but may be detectable by microlensing. The role of eccentric outer planets in triggering instabilities in an inner chain has also been studied simultaneously and independently by Ogihara & Kunitomo (2026).

• Compact planetary systems continuously evolve over hundreds of millions of years, or billions of orbits. Comparing planet formation models that end at 100 Myr or less to mature systems with ages 1−10 Gyr can be deeply misleading. The common issue that many population synthesis models produce too many resonant systems may simply be an artifact of ending N-body integrations too early. If the high fraction of mean-motion resonance in young systems reported by Dai et al. (2024) is confirmed by further observations, this may be a correct prediction of these models.

Our model is highly simplified. We intentionally neglect complex physics in order to highlight the celestial mechanics of multiplanet systems. While many population synthesis models have implemented additional physical effects such as disk evolution, dust drift and pebble accretion, or atmospheric evolution, these processes are highly uncertain. We wish to emphasize that better-understood mechanisms like planet migration, and especially Newtonian gravity, should be tested to their limits to determine where they fall short. More complex models cannot neglect dynamical evolution when comparing to known systems.

As a whole, our findings strongly support the “breaking-thechains” model as the dominant pathway for the formation of systems of small planets within 100 d. It has become increasingly clear that the architectures of planetary systems are not set at the moment of disk dissipation (Berger et al. 2020; Yang et al. 2023). Current and upcoming observations with TESS and PLATO, especially those of young systems and with well-measured host star ages, will help us map the evolution of planetary systems over their lifetimes. The onus then lies on theorists to produce models of dynamical evolution extending beyond the first 1−10% of the lifetimes of typical systems that can be used to reliably infer details of their initial conditions and formation. This will be a formidable challenge for traditional N-body methods, but continuing advancements in accurate semi-analytical models may be a promising path forward (Kimura et al. 2025; Kokubo et al. 2025).

Acknowledgements

We are grateful to referee Dan Tamayo for useful comments that improved the paper. We thank Alessandro Morbidelli, Masahiro Ogihara, Masanobu Kunitomo, and Konstantin Batygin for insightful discussions. This work was funded by the ERC project N. 101019380 “HolyEarth.”. A.P. is supported by the French government, through the UCA^J.E.D.I. Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.

References

---

Appendix A Intrinsic system statistics

Figures A.1–A.4 show population statistics for the 40steep, 10straight, 20straight, and 40straight runs.

Appendix B Numerical tests

We performed a series of numerical tests to ensure that the instability mechanism and timescale were not numerical artifacts caused by our choice of N-body integrator. We chose a system from the 20flat sample that underwent an instability in the original simulation after approximately 110 Myr. This particular system had a dynamically separated inner and outer resonant chain, and the instability arose after a gradual increase in the libration amplitudes of the resonances in the outer chain. To test the robustness of this instability and the quantitative accuracy of the timescale, we removed all of the planets in the inner chain and re-integrated the outer system starting from the end of disk dissipation (t=5.1 Myr).

First, our re-integrations varied the integrator timestep as dt₀/2ⁱ, where i ranged from 0, 1, …, 6 and dt₀ is 1/20 of the period of the innermost planet. Figure B.1 shows time to instability in these runs (solid circles). In each case, an instability of the same nature arose in the system, with the time to instability within a factor of two of the original value of 110 Myr. Such scatter represents a fundamental limitation in instability time predictions for chaotic systems (Hussain & Tamayo 2020). The absence of a trend, especially the lack of an increase in the instability timescale for small timesteps, suggests that the slow growth in libration amplitude is not caused by truncation of the N-body Hamiltonian. If this effect were purely due to roundoff error or truncation error, there would be a clear increasing or decreasing trend with the timestep, respectively.

Second, we re-integrated the outer system with the SABA(10, 6, 4) symplectic integrator implemented in rebound (Laskar & Robutel 2001; Rein et al. 2019b), which has a much higher order than WHFast. SABA indeed conserves angular momentum and energy much better than WHFast for the same timestep (thick lines in Fig. B.2), at least until a close encounter that it was not designed to handle. However, the actual time of instability is nearly the same with SABA (crosses in Fig. B.1) and TRACE with WHFast (circles). This suggests that that the instability is robust to choice of integrators and likely not a numerical artifact.

All Tables

All Figures