{'fraction_non_alphanumeric': 0.030810062312485575, 'fraction_numerical': 0.019655358104469573, 'mean_word_length': 4.874152733845459, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Criterion for the Stability of Planets in Chains of Resonances

Uncovering the formation process that reproduces the distinct properties of compact super-Earth exoplanet systems is a major goal of planet formation theory. The most successful model argues that non-resonant systems begin as resonant chains of planets that later experience a dynamical instability. However, both the boundary of stability in resonant chains and the mechanism of the instability itself are poorly understood. Previous work postulated that a secondary resonance between the fastest libration frequency and a difference in synodic frequencies destabilizes the system. Here, we use that hypothesis to produce a simple and general criterion for resonant chain stability that depends only on planet orbital periods and masses. We show that the criterion accurately predicts the maximum mass of planets in synthetic resonant chains up to six planets. More complicated resonant chains produced in population synthesis simulations are found to be less stable than expected, although our criterion remains useful and superior to machine learning models.


Introduction
Although compact systems of sub-Neptune planets are abundant, a detailed understanding of their formation remains incomplete. According to most theories of planet formation, planets form in gaseous protoplanetary disks where interactions between the planets and gas are inevitable. These interactions cause inward migration of the planets towards the disk's inner edge and capture them into chains of mean-motion resonances (Terquem and Papaloizou, 2007;Cresswell and Nelson, 2008;Lin, 2008, 2010;Cossou et al., 2014;Hands et al., 2014). Indeed, we expect that resonant chain systems such as TRAPPIST-1 Kepler-80, Kepler-223, and GJ 876 formed in this way (Mills et al., 2016;Luger et al., 2017).
Yet, population studies of exoplanet systems have revealed that resonant chains are in fact rare and that systems of multiple sub-Neptune planets are typically not in resonance (Fabrycky et al., 2014). Thus, either some process prevents the formation of resonances in the first instance, or primordial resonant chains are disrupted after the gaseous nebula dissipates. Recent work (Izidoro et al., 2017(Izidoro et al., , 2021Goldberg and Batygin, 2022) argues for the latter scenario, hypothesizing that widespread dynamical instabilities break the resonances and then a phase of giant impacts sculpts the system. Detailed simulations of such a process produce results matching the observed period ratio distribution, transit multiplicities, and peas-in-a-pod patterns of intrasystem uniformity. However, the mechanism of the instability itself is not well understood from a fundamental level, nor is there a practical way to predict which resonant chain systems are unstable.
The stability of planetary systems has been a topic of research for centuries since the development of celestial mechanics (Laplace, 1799;Le Verrier, 1840;Poincaré, 1899). With the introduction of numerical integration, the Solar System was recognized to be chaotic (Wisdom, 1983;Roy et al., 1988;Laskar, 1989) and hence unpredictable on gigayear timescales, at least on a quantitative level (Batygin and Laughlin, 2008;Laskar, 2008). Now, the rapidly growing population of exoplanetary systems, and their exotic architectures, has spurred a renewed interest in fully understanding the stability of general systems of planets (Deck et al., 2012;Batygin et al., 2015).
Previous studies, while extensive, have generally focused on two-planet systems (Gladman, 1993;Deck et al., 2013;Hadden and Lithwick, 2018;Petit et al., 2018) and the nonresonant 3+ planet regime (Chambers et al., 1996;Quillen, 2011;Petit et al., 2020;Tamayo et al., 2021;Rath et al., 2021). On the other hand, the stability of planetary systems in chains of resonances has received only limited attention. Early work was primarily empirical: Matsumoto et al. (2012) performed numerical integrations of equal mass planets locked into : − 1 resonance and found that the maximum number of planets that could be captured into the chain decreased with increasing and planet mass. Later work (Matsumoto and Ogihara, 2020) confirmed these conclusions and uncovered the unexpected result that a nominally stable resonant chain could be made unstable by a decrease in the mass of either the planets or the star.
On the analytical side, Pichierri and Morbidelli (2020) considered an equal-mass three-planet system as the simplest instance of a first-order resonant chain. Through involved perturbation theory, they showed that a secondary resonance between the fastest resonant libration frequency and a difference of the synodic frequencies can drive an instability. Rather than continuing their analytical approach, in this work we simplify their results and generalize to unequal masses and an arbitrary number of planets. Our analysis of the Pichierri and Morbidelli (2020) mechanism naturally leads to a criterion for the stability of a resonant chain, and a limit on the planet mass-or alternatively, multiplicity-in a resonant chain. We verify these results numerically on a suite of synthetic planetary systems.

Analytical Estimate of Stability
We define a resonant chain as a system of three or more planets in which each adjacent pair of planets is locked into mean-motion resonance. One can construct a wide variety of oscillation frequencies from the orbital elements, but important frequencies can be broadly separated into three categories: synodic, resonant, and secular. Synodic frequencies are linear combinations of the mean motions and do not depend on planet masses. Resonant frequencies describe the oscillations of critical resonant angles which, for twobody first-order resonances, take the form Here, is the resonant index, is the mean longitude, and is the longitude of pericenter of the -th or + 1-st planet. Finally, secular frequencies, which arise from orbit-averaged perturbations, are typically much slower than synodic and resonant frequencies and thus are not considered in this work. Pichierri and Morbidelli (2020) hypothesized that the onset of dynamical instability in compact resonant chain systems is triggered by the commensurability, or near equality, between a resonant libration frequency and a difference of synodic frequencies. Modulation of the resonant angles by synodic perturbations allows the resonant locks to break, leading to chaotic behavior. Our goal is to extend that work to more than three planets of unequal mass.
Consider a resonant chain of planets with masses 1 , ..., and in pairwise first-order resonances of 1 : 1 − 1, … , −1 : −1 − 1 so that the period ratios are ∕ +1 ≈ ( − 1)∕ . Studying the secondary resonance of Pichierri and Morbidelli (2020) would require writing the Hamiltonian of the entire system. However, we can take a simpler approach by comparing the libration frequencies of the individual resonances to the differences in synodic frequencies throughout the system.
For the purposes of computing libration frequency, we will ignore the contributions of planets that are not in the pair being considered. We verify this assumption below in our n-body simulations. In that case, the angular frequency of libration for the angle , +1 = +1 − ( − 1) − is approximately (Batygin, 2015) , +1 = 3 2 1 + 2 * 2∕3 (( − 1) 10 ∕ 4 ) 1∕9 (3( res ) 2 ) −1∕3 ( 2) where * is the stellar mass and res ≈ −0.8 + 0.34 is a constant derived from Laplace coefficients (e.g. Deck et al., 2013). As discussed in Pichierri and Morbidelli (2020), the ( ∕ * ) 2∕3 scaling is appropriate only at low eccentricities where a shift in the equilibrium point induces a forced eccentricity. At higher , the scaling is ( ∕ * ) 1∕2 . 1 Synodic frequencies are straightforward to compute. Following Pichierri and Morbidelli (2020), we have as the angular frequency of conjunctions of planets and +1. However, the analysis of Pichierri and Morbidelli (2020) specifically identifies the difference in synodic frequencies as the slower and more relevant frequency. This is Synodic frequencies are typically faster than libration frequencies. Thus, overlap is most likely to occur when the slowest difference of synodic frequencies is commensurate with the fastest libration frequency. We define the characteristic quantity for resonant chain stability to be where the minimum and maximum are taken over all synodic and libration frequencies in the chain, respectively. So far we have maintained generality, but for simplicity we will now assume that all the resonances have the same index and the mass of each planet is . Now, the slowest difference of synodic frequencies is Δ ̇ −2, −1, and the fastest libration frequency is 1,2 . Setting = 1 leads to a maximum mass of planets in the chain of This is an explicit computation of the critical mass identified by Pichierri and Morbidelli (2020).

Numerical Tests
We ran a suite of numerical experiments to test the validity of Equation 6 for different values of and . We place planets of mass ∕ * = 3 × 10 −6 on orbits 1 − 2% wide of the : −1 resonance. The semi-major axis of the inner planet is fixed and eccentricity damping is applied to all planets. To ensure sequential capture into resonance, migration is turned on for every planet except the innermost one, using a ratio of migration to eccentricity damping timescales of ∕ = 3 × 10 2 . Once the two-body resonant angles begin to librate, we remove migration and eccentricity damping exponentially so that the system settles to its stable multiresonant state. Typical orbital eccentricities at this point are ∼ 0.01. Then, we begin exponentially increasing the mass of Figure 1: The evolution of resonant libration and synodic frequencies as planet masses are increased until the instability occurs, for two initial planet configurations shown as cartoons above each grid. The left grid corresponds to a system with 3 planets started in 5:4 resonances. Each panel represents one of the four resonant angles; the red colormap is a spectrogram, or the amplitude of the Fourier transform over time, of that angle. Each resonant angle has multiple libration modes, the frequencies of which increase with mass. The bright red lines plot the analytical estimate of the libration frequency from Equation 2. Horizontal colored lines indicate synodic frequencies: blue lines are the synodic frequencies themselves and the green line is the difference of synodic frequencies (Equation 4). The right grid is the same as the left, but with 6 planets in a chain of 4:3 resonance. Only the libration frequencies for the innermost (top) and outermost (bottom) pairs of planets are plotted. each planet adiabatically (i.e. with ≫ 1∕ 1,2 ) in order to increase the libration frequency. Once an instability occurs (defined as any planet orbit becoming hyperbolic) we stop the simulation and record the planet masses. We attempted this process for integer values of between 2 and 8, and between 3 and 9. For each pair of and , we repeated the simulations 10 times with slightly different initial conditions to smooth over the chaotic behavior, although in all cases the scatter was very small. Our simulations use the whfast n-body integrator from the rebound software package and a maximum timestep of 1∕20 of the inner orbital period (Rein and Tamayo, 2015). Migration and eccentricity damping are included from the reboundx extension (Tamayo et al., 2020b).
An example of such a simulation is shown in the left panel of Figure 1, which has = 3 and = 5. As the planetary masses are increased, the libration frequencies increase, but the synodic frequencies remain constant. After 2 × 10 6 orbits of the inner planet, the resonant angles begin to circulate and the orbital eccentricities grow rapidly until there is a close encounter. Within a few orbits, at least one planet orbit becomes hyperbolic and the simulation ends. The onset of instability happens almost precisely when the highest-frequency mode of the innermost resonant angle (involving 1 and 2 ), as estimated by Fourier transform of the libration angle, intersects the difference of synodic frequencies. The bottom panels show that the libration of the resonant angle involving 2 and 3 is slower and a resonance with the synodic frequencies does not occur within the simulation timeframe. Figure 1 also demonstrates the accuracy of the analytical estimate for libration frequency, which remains within a factor of 2 of the true value throughout the simulation. The libration frequency approximation predicts that the instability will arise at max ∕ * = 1.0 × 10 −4 , whereas in the simulation the instability comes slightly later, at max ∕ * = 1.4 × 10 −4 . Nevertheless, the numericallyestimated libration frequencies grow more steeply with mass than the analytical estimate, suggesting that the low eccentricity assumption in Equation 2 has been violated.
A more complicated example is shown in the right panel of Figure 1 in which = 6 and = 4. Here, the frequency structure is more complex and the outermost synodic frequencies are slower. At the first crossing of libration and synodic frequencies, there is a resonant kick and the libration amplitudes increase instantaneously (visible as the blue synodic frequency lines becoming thicker). Upon the equality of the fastest libration frequency and the secondslowest synodic frequency, the resonant angles begin to circulate and the instability is triggered. Because the instability happens after the libration frequency of the inner planet pair has 'overshot' the difference in synodic frequencies of the outer triplet, the analytical maximum mass prediction is an underestimate of the simulation results by a factor of ∼ 2.
The full set of simulations is summarized in Figures 2  and 3. Figure 2 explores how the maximum planet mass varies with resonant index for constant multiplicity. Our analytical estimate is an excellent fit to the numerical results over a broad range of parameter space. In particular, Equation 6 maintains accuracy for values of between 3 and 7, correctly reproducing the downward trend with . This trend is in fact analogous to the Hill spacing stability criterion in non-resonant systems because the semi-major axis ratios are smaller for higher . However, instability in non-resonant systems can be fully accounted for by averaging over synodic terms and considering two-body resonance overlap and three-body resonance diffusion (Petit et al., 2020;Rath et al., 2021). Indeed, Figure 2 demonstrates that the non-resonant stability boundary from Petit et al. (2020), including the > 4 planets correction, consistently predicts a smaller maximum planet mass than is actually seen in resonant chains. Resonant chains can be stable at separations for which nonresonant systems are unstable because resonant chains reside at a fixed point in the phase space. However, interactions between synodic and resonant frequencies can excite the system away from this fixed point and into the surrounding chaotic region. Figure 3 contains the same data but shows how the maximum mass varies with multiplicity for a constant resonant index. Our analytical estimate predicts an exponential decrease in max with . While this is true for small , the dependence on multiplicity seems to flatten out near ∼ 6 planets. This may be because the inner resonance is less able to "communicate" its frequency to the outer planets for high values of , and as a result, the chain behaves like one with fewer planets. It is worth noting that a similar pattern of saturation, in which stability decreases with but flattens after ≳ 5, occurs in the non-resonant case (Chambers et al., 1996). A somewhat more accurate definition of (Equation 5) would therefore consider only adjacent subsystems of 5 − 6 planets. However, to maintain simplicity, for this work we will use the previous definition that assumes perfect coupling among all planets. Figure 3 also shows poor agreement between our prediction and simulations for the 2:1 resonance. This is likely due to the presence of indirect terms and asymmetric libration in that resonance (Beauge, 1994). That is, when eccentricities grow past ∼ 0.03, the libration centers shift away from 0 and . This is exceeded in our numerical experiments for the 2:1 resonance and our analytical estimates do not consider the asymmetric libration.

Applications to the Formation of Planetary Systems
While the results of the previous section are promising, it remains to be demonstrated whether the criterion for resonant chain stability is relevant to the more complex system architectures that are anticipated in the formation of compact super-Earth systems. Here, we apply the criterion to simulations of super-Earth system formation to show that it effectively predicts their long-term stability as well.
Our test sample is the set of synthetic planetary systems produced in the simulations of Izidoro et al. (2021). The final systems closely replicate many aspects of the observed sample of compact super-Earth systems. Beyond n-body dynamics, these simulations incorporate orbital migration and eccentricity and inclination damping due to planet-disk interactions as well as pebble accretion. In other words, the simulation suite of Izidoro et al. (2021) constitutes a successful instance of population synthesis. Within the context of these formation simulations, the gas disk dissipates at = 5 Myr, but the integrations continue until = 50 Myr in order to allow for instabilities that were suppressed by the protoplanetary disk to arise. We consider "stable" systems to be those that do not experience an instability after = 5 Myr, and "unstable" systems to be those that did experience an instability after = 5 Myr.
We removed planets with masses below 0.3 ⊕ because they tend to interfere with analyzing the chain while not contributing significantly to the dynamics. We also removed systems with a pair of planets that have semi-major axis ratios less than 1.05 because our criterion does not account for the 1:1 resonance. After these cuts, there were 54 unstable systems and 30 stable ones. The next step is to identify the likely resonances within the chain. We do this by computing the period ratio of adjacent planets. If the ratio is within 3% of a first-order resonance : − 1, for 1 < < 11, we assume the planet pair lies in that resonance. If not, we search for second-and third-order resonances with the same method but halve the threshold distance. In the case that no candidate resonance is found, we consider the chain to end at that point. For each planetary system, this process generates a collection of resonant chains separated by secular architecture. Chains with fewer than three planets are discarded because they have no difference of synodic frequencies. 2 We then calculate the libration frequency , +1 for each first-order resonance using Equation 2 (higher-order resonances are ignored) and the difference of synodic frequencies Δ ̇ , +1, +2 for each (adjacent or non-adjacent) planet triplet using Equation 4. Finally, the stability criterion is computed using Equation 5. Because some systems contain multiple resonant chains, and an instability in any one of the chains classifies the system as unstable, the system is taken to be the smallest of any of the chains.
Our hypothesis is that the secondary resonance sets a stability threshold of crit ≈ 1, above which the chain is stable. Of the 30 stable systems, 29 have > 1. However, the unstable systems have a broader distribution, clustering around ∼ 2. To quantify the boundary, we used a logistic regression (Figure 4) to model the probability of stability given only the log 10 of the chain. The fit suggests that the threshold is crit ∼ 3. That is, stability over 10 9 orbits seems to prefer wider spacing between synodic and libration frequencies than our criterion predicts.
Such a result is surprising in the context of our experiments in Section 3, which suggest that crit is near unity for ≲ 6 and smaller for higher-multiplicity systems. While fully understanding this discrepancy is outside the scope of our work, we can speculate on possible sources. Our resonant chains formed in Section 2 are especially "clean," that is, all two-body and three-body resonant angles librate with small amplitudes. In contrast, larger libration amplitudes in the Izidoro et al. (2021) systems could render them more vulnerable to higher-order secondary resonances that appear at > 1. Another possible explanation is that modulation of the resonant frequencies and widths by secular interactions with other planets in the system causes to vary over long timescales (Tamayo et al., 2021). : Relationship between our stability criterion and the true stability of synthetic and real resonant chains. Top: the cumulative distribution of log 10 for unstable and stable synthetic resonant chains from Izidoro et al. (2021), and the half-dozen well-characterized chains with small planets. Bottom: blue and orange points mark the same synthetic systems as in the top panel. Gray bars show the fraction of systems within that bin that are stable; the red curve is a logistic regression fit of the probability of stability.
As for the observed set of resonant chains, Figure 4 shows that they generally have ∼ 30 − 100, with the exception of TRAPPIST-1, for which ≈ 3. Finally, it is important to note that these simulations only capture the first 50 Myr, but instabilities can occur after billions of orbits (Petit et al., 2020). Some systems recorded as 'stable' might actually be unstable with a longer integration time that is representative of the age of typical exoplanet systems.

Discussion
Inspired by the analytical study of resonant chains, we have identified a criterion to quantify the stability of planets locked in a chain of resonances in accordance with the Pichierri and Morbidelli (2020) mechanism. We argue that the overlap between a fast resonant libration frequency and a slow difference of synodic frequencies leads to chaotic behavior and a dynamical instability. Our criterion predicts a maximum planet mass in a chain of planets with : − 1 resonances and closely agrees with numerical simulations for > 2 and ≲ 6.
In addition, this mechanism explains the counterintuitive result found by Matsumoto and Ogihara (2020) wherein an instability in a maximally-packed resonant chain can be triggered by decreasing the planet masses by 10%. Specifically, resonant chain formation occurs in a dissipative environment that suppresses the instability. During migration, the maximum libration frequency can approach, or even exceed, the slowest synodic frequencies, but the system settles into a local island of stability with < 1. After the disk is removed and the masses are decreased, the libration frequencies change and the system enters the chaotic region between the island of stability and the = 1 boundary. While initially discovered numerically, mass loss may in fact be a plausible candidate for the trigger of dynamical instabilities in packed resonant chains. Indeed, mass loss of this magnitude is physically reasonable as a result of photoevaporation (Owen, 2019). Furthermore, the highest libration frequency typically comes from the innermost resonance and hence depends only on the masses of the inner two planets. Those planets are most susceptible to photoevaporation by virtue of their proximity to the star.
Machine learning models have been especially successful in analyzing the stability of multi-planet systems. To compare our results to previous work, we used the state-ofthe-art SPOCK model (Tamayo et al., 2020a) to predict the stability of the synthetic chains presented in Section 4. We set the probability threshold to be 0.5 and ran the model in two different ways. First, we use as input each system in its entirety from Izidoro et al. (2021), only excluding the planets below 0.1 ⊕ . SPOCK correctly predicted the stability of 13 of the 30 stable systems and 42 of the 54 unstable ones. Second, we input each of the individual resonant chains identified in the systems of Izidoro et al. (2021), as described in Section 4, and use SPOCK to compute the probability of their stability. We treat those as independent random variates, and for each system we estimated the probability of stability of the whole system by computing the probability that every chain within it is stable. In that case, SPOCK correctly predicted the stability of 20 of the 30 stable systems and 31 of the 54 unstable ones. As a comparison, our onedimensional logistic regression (Figure 4) achieves 25∕30 for stable systems and 40∕54 for unstable with the same probability threshold. Accordingly, the specific problem of resonant chain stability constitutes an instance where a careful analytical treatment is comparable to or surpasses general supervised machine learning techniques.
If it is indeed true as some have suggested that nonresonant systems of small planets are the products of instabilities, the mechanism of instability is of considerable importance. Previous suggestions include changes in the stellar 2 moment (Spalding and Batygin, 2016), stellar or planetary mass loss (Matsumoto and Ogihara, 2020), or a simple overpacking of the system during the disk phase (Izidoro et al., 2017). Our work does not rule out any of these mechanisms, but clarifies the dynamics underpinning the onset of the instability. Future work should explore the consequences of each of these instability mechanisms to determine whether they leave signatures detectable in the planet population.
insightful discussions and André Izidoro for providing simulation results. K. B. is grateful to Caltech, the Caltech Center for Comparative Planetary Evolution, the David and Lucile Packard Foundation, and the Alfred P. Sloan Foundation for their generous support. A. M. acknowledges support from the ERC advanced grant HolyEarth N. 101019380.