Differences between Stable and Unstable Architectures of Compact Planetary Systems

We present a stability analysis of a large set of simulated planetary systems of three or more planets based on architectures of multiplanet systems discovered by Kepler and K2. We propagated 21,400 simulated planetary systems up to 5 billion orbits of the innermost planet; approximately 13% of these simulations ended in a planet–planet collision within that time span. We examined trends in dynamical stability based on dynamical spacings, orbital period ratios, and mass ratios of nearest-neighbor planets as well as the system-wide planet mass distribution and the spectral fraction describing the system’s short-term evolution. We find that instability is more likely in planetary systems with adjacent planet pairs that have period ratios less than 2 and in systems of greater variance of planet masses. Systems with planet pairs at very small dynamical spacings (less than ∼10–12 mutual Hill radii) are also prone to instabilities, but instabilities also occur at much larger planetary separations. We find that a large spectral fraction (calculated from short integrations) is a reasonable predictor of longer-term dynamical instability; systems that have a large number of Fourier components in their eccentricity vectors are prone to secular chaos and subsequent eccentricity growth and instabilities.


INTRODUCTION
Exoplanet surveys have revealed a large number and wide variety of multi-planet systems in the solar neighborhood.The close orbital spacings of many of the observed multi-planet systems seemingly place many of them on the edge of dynamical stability, despite the gigayear ages of many of their host stars (e.g., as previously discussed in the context of the surprisingly packed Kepler-11 system; Lissauer et al. 2011).The study of planetary system stability is a long-standing area of celestial mechanics research (see, e.g., recent discussions in Petit et al. 2020;Tamayo et al. 2021;Weiss et al. 2023).The stark contrast between the observed Kepler multiplanet systems and our own solar system has also sparked many investigations into different models of planetary system formation and migration (see recent reviews by Weiss et al. 2023;Lissauer et al. 2023).The statistical distribution of planet masses and orbits based on dynamical stability have also been investigated (e.g., Malhotra 2015;Tremaine 2015).In this paper, we use Corresponding author: Kathryn Volk kat.volk@gmail.comthe observational data of multi-planet system architectures as the basis for a statistical study of planetary system architectures with a focus on two questions: What are the distinguishable features of a stable planetary system architecture?What are the drivers of dynamical instabilities in systems whose architectures are broadly similar to those observed?
To answer these questions, we first identify a few parameters of planetary systems that broadly quantify their architectures, and then we examine how long term stability depends upon small variations of these parameters (within the uncertainties of the observational data of the Kepler sample).The parameters we identify are the following: the dynamical spacings, orbital period ratios, and mass ratios of nearest-neighbor planets; the coefficient of variation in system-wide planet masses; and the spectral fraction describing the system's secular orbital evolution.These parameters have been investigated in a few recent studies, though mostly in isolation.For example, planetary dynamical spacings have been explored in randomly-generated simulated systems by Pu & Wu (2015); Obertas et al. (2017); Rice & Steffen (2023), period ratios have been discussed by Steffen & Hwang (2015), planet mass distributions have been discussed by Mishra et al. (2023), and the spectral in-dex of secular orbital evolution has been discussed by Volk & Malhotra (2020).In the present work, we examine all of these parameters for a large set of planetary architectures modeled on observed systems.In addition to expanding our understanding of the features of longterm stable planetary system architectures in general, our study can help identify observed exoplanet systems where dynamical considerations can provide useful constraints on planetary properties or even identify cases where some observational constraints may be spurious.
In Section 2, we describe our starting set of observed planetary architectures and the large suite of numerical simulations we performed to investigate system stability.Section 3 describes our statistical exploration of stability in these simulations as a function of various system parameters.Section 4 provides a brief discussion of a number of notable observed multi-planet systems that warrant additional future dynamical investigations to improve constraints on planetary masses.Section 5 provides a discussion of our results, focusing on the system-wide properties that are most predictive of stability/instability.We summarize our findings in Section 6.

A LARGE SUITE OF PLANETARY SYSTEM SIMULATIONS
We initiated this study in 2019 with the then-known set of Kepler and K2 planetary systems hosting three or more confirmed planets and having a published estimated stellar mass in the NASA Exoplanet Archive 1 composite data table (NASA Exoplanet Archive 2019).Subsequently we added all the additional reported systems with four or more planets as of January 2023 (NASA Exoplanet Archive 2023).We also added new simulations of any of our original systems with additional planets reported between 2019 and January 2023.The resulting 165 multi-planet Kepler and K2 systems are listed in Table 2 in Appendix A, including citations to published planet radius, planet mass, planet eccentricity, and stellar mass estimates we make use of.We do not include circumbinary planets in our sample, but we do include 8 planetary systems that orbit one star within likely or confirmed multi-star systems.These systems are noted in Table 2.For these systems, we do not account for the stellar companions in our modeling as we are not attempting high precision orbital propagation, but investigating only broad dynamical stability properties of their orbital architectures.
Our sample of planetary systems includes 165 systems with three or more transiting planets; a majority of these 1 https://exoplanetarchive.ipac.caltech.eduplanets have radii smaller than 4R ⊕ and only a minority have measured masses.For 162 of these systems, we run a set of simulations with planetary masses assigned from a prescribed (statistical) mass-radius relationship.Henceforth, we will refer to this simulation set as the "MR relationship cases".Our sample includes 47 systems with published mass and/or eccentricity constraints from either radial velocity (RV) or transit-timing variations (TTV) for at least one planet that provide mass estimates more stringent than a statistical MR relationship.For these, we run a set of simulations using those mass/eccentricity constraints.Henceforth, we will refer to this simulation set as the "TTV/RV cases"; note that 42 of the 45 systems in this set are also in the MR-relationship set.Full details of the initialization of each of these cases are given below.We emphasize that none of our simulations represent the exact configuration of the observed planetary systems.Our goal is to probe general trends in planetary system stability based on system-wide properties such as period ratio distributions and typical planetary mass distributions.

MR relationship cases
For the 162 systems with transiting planets smaller than 4R ⊕ , we initialize 100 simulation instances as follows using the same procedures as in Volk & Malhotra (2020) (which presented results from a subset of the wider range of simulations in this work).We assign individual planetary masses by randomly sampling each planet's radius uniformly within the reported 1σ uncertainties and then use that radius to sample a mass from the statistical MR relationship from Wolfgang et al. (2016) 2 .That mass-radius relationship can be sampled based on either RV or TTV priors, so we evenly split our simulations between those two cases.For the small number of systems that have some planets whose radii exceed 4R ⊕ (the size range over which Wolfgang et al. (2016)'s MR relationship is valid), we revert to TTV or RV mass constraints if they exist or to the Exoplanet Archive's listed mass estimate based on the MR relationship from Chen & Kipping (2017) 3 .These special cases are noted in Table 2.We then assign stellar masses and planet orbital periods chosen randomly within the 1σ error estimates.We choose to sample all the parameter uncertainties uniformly and only within 1-σ because available computational and time resources limit us to a relatively small number of simulations.This allows us to probe a fair range of parameters with reasonable density of initial conditions; we are also not attempting to reproduce the exact observed systems.We assign small initial orbital eccentricities e and inclinations i from Rayleigh distributions of widths σ e = 0.02 and σ i = 1.4 • , which are consistent with estimates of the intrinsic eccentricity and inclination distributions for Kepler multiplanet systems (e.g., Hadden & Lithwick 2014;Xie et al. 2016;Van Eylen & Albrecht 2015).We randomize the planets' longitudes of ascending node Ω, arguments of perihelion ω, and mean anomalies M .This resulted in 16,500 simulated systems representing 162 unique observed multiplanet systems.Three observed systems were included twice: K2-138 due to including additional planets detected since our initial exploration of that system architecture, Kepler-90 for an additional choice of planetary mass assumptions for the R > 4R ⊕ planets, and Kepler-1542 for significantly updated planetary radius measurements.Each instance was integrated forward using the mercurius integrator within rebound (Rein & Liu 2012), which uses the whfast (Rein & Tamayo 2015) integrator for most timesteps (using a timestep equal to 1/40th of the shortest orbital period) and switches to the adaptive step-size ias15 (Rein & Spiegel 2015) integrator to resolve close-encounters between massive bodies.We include a first-order post-Newtonian general relativity correction as a user-defined force based on a simplified potential (the implementation is described in Tamayo et al. 2020).We integrated each simulated system until a planetary collision was detected or the integration reached 5 × 10 9 orbital periods of the innermost planet.As discussed in Volk & Gladman (2015), ejection is highly unlikely for planets orbiting as close to their stars as the observed transiting multiplanet systems and stellar collisions are also expected to be rare; while our simulations did have inner and outer boundary conditions, these were never triggered, and all simulations either survived the maximum integration length or ended in planetary collisions.

TTV/RV cases
For the 47 observed systems with at least one planet with a reasonably well-constrained TTV or RV mass (i.e., a mass measurement ranges smaller than those provided by the above MR relationship) or eccentricity measurements for multiple planets, we generated 100 simulations with planetary masses (and eccentricities, where reported) assigned from within 1σ uncertainties of those measurements.For planets in these systems without TTV/RV measurements, we assigned masses according to the same MR relationship as described above.We sampled the same Rayleigh distribution in e as above for planets without eccentricity constraints.We assigned all planets inclinations from the Rayleigh distribution described above and randomized all other orbital angles.For planets with multiple reported masses, we defaulted to the masses reported in the NASA Exoplanet Archive composite data table (references for our chosen masses are given in Table 2).For two systems (Kepler-107, Kepler-11) that proved to be particularly prone to collisions, we ran simulations both with and without the reported eccentricity constraints.We generated 4900 simulations representing the 47 observed systems and integrated them as described above.

Summary of simulation outcomes
Overall, we ran simulations of 21300 cases of Keplerlike planetary architectures.Of these, 2764 simulations (∼ 13%) ended with planet collisions occurring within 5 × 10 9 orbital periods of the innermost planet.Of the 165 observed system architectures we consider, 60 experience at least one simulation that ends in a collision.When only considering the MR-relationship cases, 30% of the observed systems (49/162) have at least one unstable simulation.In some contrast, 55% of the system architectures for which we used reported planet masses (26/47) have at least one simulation that ends in collision.
The distributions of instability times we found in our simulations are summarized in Figure 1.We show the results for the entire simulation set as well as the distributions for the MR relationship cases and the TTV/RV cases.Additionally, we separate systems that contain near-resonant (or potentially actually resonant) planets and those that do not have planets in near-resonant configurations.To classify the planets within systems as either potentially near or not near strong mean motion resonances (MMRs), we use analytical resonance width estimates from Murray & Dermott (1999) for a subset of low-order MMRs, including all internal and external first order MMRs from the 2:1 through the 8:74 .These analytical resonance width estimates are far from perfect, especially at very low eccentricities (see, e.g., Malhotra & Chen 2023), but they afford a quick way to estimate whether the observed period ratio of a pair of planets likely falls inside, near, or well outside a mutual MMR for some assumption about their masses.For each system we simulated, we include our assessment of their resonant status for their initially low-eccentricity orbits in Table 2 as either likely non-resonant (no planets are near or in strong MMRs), likely resonant (at least one pair of planets is in or quite near a low order MMR), or ambiguous (apparently non-resonant, but modest eccentricity enhancement could bring a planet pair near a low order resonance).For each curve in Figure 1, the simulations are weighted such that all observed systems contribute equally; for example, in the 'all simulations' curve, systems that were simulated 200 times because they were included in both the MR relationship and TTV/RV cases were given weights of 50% compared to systems that appeared in only one of those simulation sets and thus were simulated 100 times.
From Figure 1, it is clear that the near-resonant and TTV/RV cases have higher instability rates (∼ 30% of the simulations) than the non-resonant and MR relationship cases (∼ 7% of simulations).There is significant overlap between the observed systems that contribute to the curves for the near-resonant and TTV/RV cases as close proximity to strong MMRs increases the likelihood of detectable TTVs.The near-resonant systems account for approximately one-third of all the collisions in our simulations.In addition to the higher rate of instabilities, the TTV/RV and near-resonant cases also have instabilities that occur on shorter timescales.This is consistent with some of these instabilities being directly related to the MMRs, which have shorter dynamical timescales than the secular dynamics that can also drive instabilities; this is discussed further in Section 3. Overall, the rate of instabilities in the simulations is linear in log(t), which has has been found in previous studies of long term stability in ensembles of dynamical systems (e.g.Holman & Wisdom 1993;Chambers 2007;Volk & Gladman 2015).We note that the shape of the near-resonant systems curve is distinct from the other curves; it displays a shallow slope similar to that of the non-resonant curve up to 10 7 orbits before changing to a steeper slope.However, this could be the effect of a few observed system architectures as the number of observed systems contributing to the near-resonant curve is smaller than for any of the other curves.We discuss the distribution of key parameters for the unstable cases and the implications of these overall instability rates in more detail below.

STABILITY AS A FUNCTION OF DIFFERENT SYSTEM PARAMETERS
Here we examine the characteristics of all our simulated planetary architectures as well as individual simulations to identify the parameters most indicative of stability and instability.We separate our entire suite of simulations into stable and unstable categories based 10 3 10 4 10 5 10 6 10 7 10 8 10 9 on whether a collision occurred within the simulated 5 × 10 9 orbital periods of the innermost planet.We also label system architectures (i.e., the combination of observed orbital periods and planetary radii in the MRrelationship cases or observed orbital periods and constrained masses in the TTV/RV cases) as either stable or unstable based on whether any of the 100 simulations of that architecture experienced a collision.For these two sets of our simulation suite, we then examine the similarities and differences in nearest-neighbor planet dynamical spacings, period ratios, and mass ratios.We note that only ∼10% of the collisions detected in our simulations occur between non-adjacent planets; in these cases, a planet's orbit becomes eccentric enough to cross the orbits of more than its nearest neighbor(s) and happens to collide with an initially non-neighboring planet.We also examine stability as a function of the system-wide coefficient of variation in planet masses, which was recently suggested by Mishra et al. (2023) as useful for characterizing planetary architectures.Finally, we examine stability as a function of short-term secular evolution as characterized by the spectral fraction parameter, which was found by Volk & Malhotra (2020) to be predictive of long-term stability in planetary systems with four or more planets.The details of all of the considered parameters are given below.In all of our analyses below, individual simulations are weighted such that the sum of the number of all simulations of an observed adjacent planet pair is one; thus in our histograms there can be fractions (< 1) of simulated planet pairs represented because for most observed pairs, each of the 100 simulations containing that pair has a weight of 0.01 (some observed pairs have lower weight if, for example, they have 200 simulations in the MR-relationship cases due to significantly differing ra-dius measurements; see Sections 2.1 and 2.2).We examine stability as a function of our dynamical parameters in two ways: 1) examining the distribution of parameters when observed architectures are labeled as entirely stable (i.e., none of the 100 simulations of that architecture experienced a collision) vs architectures with at least one collision, and 2) on a simulation-by-simulation basis, separating simulations with and without collisions.We also examine separately our simulations that use a mass-radius relationship for assigning planet masses and those that use RV or TTV mass constraints.
We begin with examining the dynamical separations, ∆, of neighboring planets with semimajor axes a 1 and a 2 and masses m 1 and m 2 in units of their mutual hill radius R mH : where M * is the mass of the planets' host star.Figure 2 shows the distribution of ∆ in our MR-relationship cases; Figure 3 shows the same distributions for our TTV/RV cases.In Figures 2 and 3, system architectures in which no collisions occurred (i.e., all 100 simulations reached 5×10 9 orbits of the innermost planet) are shown in grey while system architectures with at least one collision (out of the 100 simulations) are shown in pink; note that these histograms are stacked.
Figure 2 shows that a large fraction of the system architectures with planet pairs at ∆ ≲ 10 are prone to instabilities; this is consistent with previous investigations (e.g., Gladman 1993;Chambers et al. 1996;Smith & Lissauer 2009;Funk et al. 2010;Malhotra 2015;Pu & Wu 2015;Obertas et al. 2017;Lissauer & Gavino 2021).The peak in the pink histogram (architectures with unstable cases) occurs at a smaller value of ∆ than in the grey histogram (stable architectures); the bottom panel of Figure 2 confirms that the two distributions are distinct.There are two additional notable features in the ∆ distribution of the unstable architectures: (i) the peak is relatively flat in the range ∼ 12 − 18, and (ii) the distribution extends to fairly large values (consistent with findings in Volk & Gladman 2015), with ≳ 10% of planet pairs at separations larger than 25 mutual hill radii.We note that the ∆ distributions of entirely stable vs unstable architectures for the systems with TTV/RV mass measurements (Figure 3) are less distinct from each other than in the MR relationship cases (Figure 2), but the unstable architectures span a similar range of ∆ as in the MR relationship simulations.It is also notable that a larger fraction of the stable TTV/RV system architectures have planet pairs with ∆ < 10 (∼ 30% of In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.1).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
planet pairs in the TTV/RV cases compared to < 5% in the MR cases).
Figure 4 shows the ∆ distribution of the actually colliding pairs of planets compared to both the noncolliding pairs in unstable simulations and planet pairs in stable simulations.This figure combines all simulations regardless of the source of the planet masses, and we use a log scale on the histogram to better highlight the colliding planet pairs.We see that collisions do indeed occur between planets at initially quite large dynamical separations, with ∼ 30% of collisions occurring between planets initially separated by ∆ > 20.
We note that ∆, which scales with the planet masses to the one-third power (Equation 1), is not the only possible formulation for dynamical separations.Resonance overlap considerations yield a dynamical separation parameter that instead scales with planet mass to the onefourth power (e.g.Hadden & Lithwick 2018;Tamayo et al. 2021).For the simulations presented here, both In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.2).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
variations of the dynamical separation parameter yield very similar results.We now examine the period ratio distributions of neighboring planets in the systems simulated; note that we define the period ratio to be larger than one (i.e., for each pair we take the period of the outer planet divided by the inner planet).Figure 5 shows this distribution in our simulations based on MR relationship cases and Figure 6 shows the same distributions for TTV/RV cases.It is immediately notable that the period ratio distribution is much less smooth than the ∆ distributions.This is partly due to the nature of the dataset on which we base our sample.Due to the large uncertainty in each planet's mass, the ∆ value for each observed planet pair is, in effect, spread out over a wide range in the 100 simulation instances; in both the MR and TTV/RV cases, these uncertainties smooth out the contributions to the distribution for each observed planet pair.In contrast, each planet's orbital period is much better constrained by the observations and has a very small uncertainty, so every instance of each simulated planet pair has a nearly identical period ratio.However, there is also an expectation on grounds of dynamical stability that the ∆ distribution is fairly smooth (e.g.Dietrich et al. 2023) whereas the period ratio distribution has strong variability near first order resonances.In particular, there are notable features in the period ratio distribution near the strongest mean motion resonances (MMRs), e.g., the peak at the 3:2 MMR and the trough near the 2:1 MMR, as discussed in previous studies of the Kepler sample (e.g., Fabrycky et al. 2014;Steffen & Hwang 2015).
A very striking result of our suite of simulations is that ∼ 90% of planet pairs in system architectures that are prone to collisions are at period ratios ≲ 2 (see the bottom panels of Figures 5 and 6).This can also be seen in Figure 7, which shows the distribution of period ratios for colliding pairs of planets across both sets of our simulations compared to the period ratios of other planet pairs in those simulations as well as planet pairs in stable simulations.Some of this tendency can likely be ex- In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.1).In pink are period ratios in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
plained by proximity; planets with smaller period ratios are on orbits that are spatially closer together than those with larger period ratios, and thus require only modest eccentricities to achieve crossing orbits.There is also a higher density of low-order mean motion resonances at period ratios less than two, which can enhance dynamical chaos and eccentricity excitation.We will discuss the implications of this period ratio trend below as well as in Section 5. We also examine the mass ratios of neighboring planets in our systems; note that we define the planet mass ratio to be the smaller of m inner /m outer or m outer /m inner , so it is always less than one.Figures 8  and 9 show the distribution of mass ratios for neighboring planets in the MR relationship cases and the TTV/RV cases, respectively.There are no strong apparent trends in these plots to differentiate the stable architectures from those that are prone to instabilities.It is notable, however, that the TTV/RV cases have a higher abundance of neighboring planets with smaller mass ratios (more unequal masses). .TTV/RV cases: Distribution of period ratios between neighboring planets in 47 observed planetary system architectures with 3 or more planets and TTV or RV mass constraints on at least one planet (top: stacked histogram that sums to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.2).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
shows the distribution of mass ratios for the colliding pairs of planets in all unstable simulations compared to the non-colliding planets in those simulations and the planet pairs in the stable simulations.We note a weak trend that the unstable simulations have slightly smaller mass ratios than the stable simulations, and that the colliding planet pairs in the unstable simulations have slightly smaller mass ratios than the stable pairs in those simulations.
Following Mishra et al. (2023), we also examine the mass distributions in our individual simulations using the coefficient of variation in planet masses; this coefficient, C v is defined as the standard deviation of the planet masses in the system divided by the average planet mass.Figure 11 shows this characterization of system-wide planet masses for all of our MR relationship simulations.The vast majority of our simulated systems have reasonably small coefficients of variation for the planet masses (C v ≲ 1, meaning the planet masses vary by less than the average planet mass).Across this 1.0 1.5 2.0 2.5 3.0 period ratio of adjacent planets  well-sampled range, the fraction of stable simulations decreases roughly linearly with increasing C v , from ∼ 98% in the smallest C v bin to ∼ 90% at C v ≈ 1.The trends at larger C v values are noisier due to the smaller number of simulated systems, but the stable fraction drops off for C v ≳ 1 and is very low at the largest C v values.This trend in C v is consistent with the weak trend noted above that unstable systems are more common for more unequal planetary mass ratios.The modest ∼ 10% decrease in the stable fraction may contribute to the high abundance of systems with nearly equal size planets (the so-called 'peas-in-a-pod' trend, discussed in, e.g.Weiss & Petigura 2020;Weiss et al. 2023).We discuss this further in Section 5. Lastly, we examine stability as a function of the secular dynamics of each simulated system.Briefly, classical linear secular theory describes the orbit-averaged perturbations of the planets on each other, which cause their orbital eccentricities and inclinations to vary over time in a quasi-periodic way with a superposition of several secular modes.Beyond the linear secular approximation, the secular evolution of multi-planet systems In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.1).In pink are period ratios in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
is more complex than the simple superposition of a few secular modes.We characterize the secular dynamics of our simulated planetary systems using the spectral fraction parameter from Volk & Malhotra (2020), which is a modification of the spectral number described by Michtchenko et al. (2002) for asteroid dynamics.To calculate a spectral fraction, we perform a short integration (5 × 10 6 orbits of the innermost planet) of each of our simulated systems with high-cadence output of the planets' orbital parameters (3000 outputs evenly spaced in time).We then perform a Fast Fourier Transform (FFT) of the orbital parameter of interest for a planet and identify the frequency with the highest power (the square of the amplitude) in the FFT.We determine how many of the frequencies sampled in the FFT have a power ≥ 5% of that maximum power, and this number divided by the total number of frequencies in the FFT is then the spectral fraction for that orbital parameter for that planet.
We considered the spectral fraction of many dynamically meaningful parameters, and we determined empir- .TTV/RV cases: Distribution of mass ratios between neighboring planets in 47 observed planetary system architectures with 3 or more planets and TTV or RV mass constraints on at least one planet (top: stacked histogram that sums to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.2).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
ically that the results for the angular momentum deficit, amd, of each planet were the most relevant to predicting long-term stability.The amd of an individual planet in a system is given by: where the index j indicates the j-th planet's parameters.The amd spectral fraction captures the same stability trends we see in the e and i spectral fractions, allowing us to rely on just one parameter rather than two; we note that the semimajor axis spectral fraction is not strongly correlated with stability in our simulation set.Following Volk & Malhotra (2020), we choose the largest spectral fraction of the amd of any planet in each system to characterize the whole system.We note that some of our simulations ended in a collision within the short 5 × 10 6 orbit simulations used for the spectral fractions.We only calculate spectral fractions for simulations that survived at least 10 6 orbits (the minimum time to compute the FFT of the amd over secular  timescales), so our analysis here excludes shorter-lived systems.The limited length of our spectral fraction integrations does limit how well the low-frequency evolution of the planets is captured by the spectral fraction parameter.The integration length of 5 × 10 6 orbits was chosen with both efficiency and expected secular timescales in mind.However there may be individual systems, e.g.ones with outer planets at orbital periods particularly larger than the innermost planet, where longer integration timescales are necessary to adequately determine the nature of their secular evolution.The general usefulness of the spectral fraction calculated from our chosen integration length (see below) indicates that timescales important to instabilities are reasonably well-captured for the systems we investigate here, but this limitation should be considered when applying the spectral fraction calculation in future investigations.
Figure 12 shows the stable fraction of our simulations versus the maximum amd spectral fraction for systems that appear not to be near any strong mean motion resonances and systems that are potentially in or near strong resonances (see Section 2.3).It is useful to separate sys- tems in this way because our definition of the spectral fraction is motivated by considering secular rather than resonant dynamics.In both the near-resonant and not near-resonant simulations, the highest stable fractions are at the lowest values of the amd spectral fractions.
The trend of decreasing stability with increasing spectral fraction is clear in the non-resonant systems; this trend is also clear in the near-resonant systems at small spectral fractions, though stability flattens out and becomes noisier at amd spectral fractions larger than ∼ 0.01 (Figure 12).This supports the idea that secular chaos is driving the onset of instability in the non-resonant systems; secular chaos can play a role in the resonant and near-resonant systems, but the spectral fraction is not as clearly predictive of stability or instability as in the non-resonant cases.Figure 13 shows the average simulation length for the non-resonant systems as a function of the maximum amd spectral fraction.Over the wellsampled range of spectral fractions (values ≲ 0.1), there is a clear trend toward shorter stability times with increasing spectral fraction.We discuss the implications of this further in Section 5.
While none of our simulations are meant to perfectly reproduce the observed Kepler or K2 systems they are  2).These plots include all simulations performed (whether the planets masses were assigned from the mass-radius relationship or from TTV/RV measurements) that survived for at least 10 6 orbits, the minimum timespan to perform a useful FFT.
based on, they still offer insights into the likely dynamics and properties of the real systems.Here we highlight interesting results for a subset of our systems and compare them to other studies that have examined these systems in detail.We discuss three systems that are potentially in resonant chains (Kepler-223, Kepler-60, and Kepler-80) and the role those resonances play in their stability.Next, we discuss the K2-138 system, comparing our simulations of its system architecture before and after confirmation of a sixth planet.Finally, we discuss a number of systems with particularly unstable architectures: Kepler-23, Kepler-402, Kepler-444, Kepler-1542, and K2-384.
Kepler-223 : The four planets in Kepler-223 are observed to be in a resonant chain (Mills et al. 2016) with period ratios of 4/3 (e/d), 3/2 (d/c), and 4/3 (c/b).When adopting the nominal TTV-derived masses and eccentricities for all four planets (  orbits of planet b (101 Myr).When adopting MR derived masses for the inner two planets (the outer two are too large for the assumed MR relationship), and eccentricities chosen from a Rayleigh distribution of width σ e = 0.02 for all planets, only 21% of the simulations become unstable on the same timescale.There is not a strong correlation between stability and the assumed masses of the inner two planets (our assumed MR distribution overlaps the TTV mass constraints), but the TTV-derived eccentricities for planets b and c are large enough to place them initially on orbits that are not far from crossing orbits, hence close to instability unless the planets are deep in resonance libration.It thus seems likely that libration in resonance is important to maintaining long-term stability at larger eccentricities.Interestingly, however, we find a similar rate of stability on 10 Myr timescales for our simulations based on the TTV-constrained masses as Mills et al. ( 2016) (25% of our simulations vs 30% of theirs) despite our simulations randomizing all the planets' orbital angles; the Mills et al. (2016) simulation initial conditions were based on photodynamical modeling, matching all dynamical parameters to the observed transits.A statistically equivalent percentage of both sets of the 10 Myr stable simulations are also stable to 100 Myr; from their 10 Myr stable simulations they randomly selected to extend 25 of them, 16 of which were stable compared to 20 of our 25 simulations.Mills et al. (2016) also ran a set of simulations based on orbital parameters that ensure libration of the Laplace angles associated with the resonant chain, finding 100% of their resulting 100 Myr simulations were stable; they thus argue that libration in the resonant chain is favorable to long-term stability.This simulation set also had smaller planetary eccentricities than their nominal TTV solution.It is interesting to note that in our low-eccentricity, MR mass simulation set, with no constraints on planetary angles, 79% of the simulations are still stable at 100 Myr.Our MR simulations have slightly lower planetary eccentricities than the Mills et al. (2016) resonant solution set, particularly for planet b (see their Extended Data Table 3), but no constraints on resonant angles.It is indeed likely, as Mills et al. (2016) argue, that libration in the resonant chain promotes long-term stability given that our stability rate is smaller than theirs.However, our simulations show that smaller eccentricities can account for some portion of enhanced stability, so libration within the resonant chain might not be strictly required for long-term stability.
Kepler-60 : The three observed planets in Kepler-60 are potentially in a resonant chain, with period ratios of 4/3 (d/c) and 5/4 (c/b).When adopting MR masses for the planets, 61% of our simulations are unstable over the 97 Myr simulations.However, none of our simulations are unstable when we adopt the TTV planet masses from Jontof-Hutter et al. ( 2016) despite not imposing any constraints on the relative orbital phases of the three planets in the system.Jontof-Hutter et al. (2016) did complete fits to the transit data and found that both resonant and non-resonant solutions are consistent with the data.Goździewski et al. (2016) argue that longterm stability in this system is more likely for planets in the resonant configuration.Amongst our MR set of 200 simulations of Kepler-60 analogs, we found (by visual examination of the resonant angles) only one case of resonant libration in the resonance chain, and this case was long-term stable.Given that all of our TTVmass-based simulations were stable to nearly 100 Myr (and the vast majority of these do not have resonant librations), it is clear that stability on that timescale is not reliant on the resonant configuration, though it is possible that gigayear stability does depend on resonant libration.We note that the Exoplanet Archive lists an additional, much longer period, unconfirmed planet candidate in this system; future studies might consider whether a long-period planet could help stabilize the architecture of Kepler-60's inner planets.
Kepler-80 : Kepler-80 is a five-planet system with multiple librating three-body resonances; the lowest-order librating angles identified by MacDonald et al. ( 2016) are ϕ = 3λ b − 5λ e + 2λ d and ϕ = 2λ c − 3λ b + λ e .In our simulations of this system architecture, only 2 of the 100 TTV-mass based simulations were stable to 5 billion inner planet orbits (∼13 Myr), and just 10 of the MR relationship simulations were stable.From visual inspection of the short integrations for our TTVbased initial conditions, ∼25% showed libration of the 3-body angles above for the resonances between planets c-b-e and planets d-e-b; we note that both stable simulations showed libration.For the MR relationship simulations, ∼30% showed libration of those 3-body angles, including 5 of the 10 stable simulations.None of our simulations showed the potential three-body resonance between planets g-c-b identified by MacDonald et al. (2021).Given the low stability rates in our simulations, and the higher rate of stability in the subset of our simulations that showed libration of the resonances for this system, it seems likely that resonant libration can enhance the stability in this system.K2-138 : The first five transiting planets in the K2-138 system were reported by Christiansen et al. (2018) with a sixth planet later reported by Hardegree-Ullman et al. (2021).Our initial set of simulations for this system was based on MR relationship masses for the initial five planets; this set was almost entirely unstable as only 3 of the 100 simulations survived the 32 Myr simulations.We ran a new set of MR simulations including the sixth planet, which has an orbital period of 42 days, significantly longer than the 2-13 day orbital periods of the inner five planets.In this set, the stability rate jumped to just over half of the simulations (54 out of 100 survived to 32 Myr).This highlights how the presence of an additional, longer period planet can help stabilize a tightly packed inner planetary system.When we examine the amd spectral fractions of the planets in the 5-planet version of the K2-138 system, they are systematically larger than those in the 6 planet system.In other words, the presence of the additional longer-period planet tends to cause the secular evolution of all the planets to be more 'orderly' with fewer secular frequencies controlling the system's evolution.Figure 14 shows the maximum amd spectral fraction for the three versions of the K2-138 system we considered, including the 6-planet version with masses constrained from RV measurements (Lopez et al. 2019).The simulations using RV masses had slightly larger spectral fractions than the MR mass simulations despite having a higher survival rate (77 of the 100 RV simulations survive to 32 Myr); however, simulations of the K2-138 system with 6 planets have more cases with significantly smaller spectral fractions than the 5-planet version.It is possible that other observed architectures that are particularly prone to instability could be similarly stabilized by additional not-yet-observed planets.Though we note that any additional outer planets would need to not induce significant inclination evolution that would make transits much less likely (see, e.g., Becker & Adams 2017).
Kepler-23 : For the three observed planets in Kepler-23, we note a significant difference in the stability rate when using TTV and MR derived masses.When adopting the TTV masses from Hadden & Lithwick (2014), the majority (87%) of our simulations become unstable within the 97 Myr simulations compared to none of the MR simulations.The planets in this system are similar sizes, with radii in the range 2-3R ⊕ , so the MR relationship from Wolfgang et al. (2016) yields masses for all planets of less than ∼ 25M ⊕ .The TTV solution from Hadden & Lithwick (2014) yields mass estimates for planets b and d similar to the MR estimates, but the mass estimate for planet c is 50 − 70M ⊕ .This large mass for planet c does not seem compatible with longterm stability, suggesting the mass estimate is too large.This is consistent with the note in Table 2 of Hadden & Lithwick (2014) that planet c is likely to have an eccentricity large enough to cause the nominal TTV mass to be an overestimate.
Kepler-402 : There are four confirmed transiting planets in the Kepler-402 system, with one additional candidate planet.In our simulations of the four confirmed planets using MR masses, 50 of the 100 simulations were unstable within 55 Myr.This is another system that warrants future investigation to see if the longer-period planet candidate enhances the stability of the observed architecture.
Kepler-444 : The five observed sub-Earth radius planets in Kepler-444 are in orbit around one star in a triple star system (see Zhang et al. 2023 for the latest stellar system properties).Two planets (d and e) are very close to a 5:4 period ratio (e.g.Mills & Fabrycky 2017), but none of the other planets are particularly close to loworder resonant ratios.Our simulations of this system's architecture using a MR relationship for the planets has a high rate of instability (81 of 100 simulations ended in a collision).We find no significant statistical differences of the key architectural parameters (planetary masses, mass-ratios, or initial Hill radius separations) between the stable and unstable instances of our Kepler-444 architectures.The TTV-derived masses of planets d and e (Mills & Fabrycky 2017) are on the small end of our MR mass ranges, but we find stable simulations spanning the full MR range for these planets, which includes masses 2-3 times more massive than the TTV-derived limits.We note that the Mills & Fabrycky (2017) TTV measurements are not listed in the exoplanet archive and thus were not included in our simulation initial conditions.The MR masses for the other planets span the mass ranges suggested by Stalport et al. (2022) to be consistent with stability based on a frequency analysis of their shorter-term semimajor axis evolution; Stalport et al. (2022) also include the effects of the Kepler-444BC binary star around Kepler-444A in their analysis.Like most Kepler multi-planet systems, the Kepler-444 planetary system is estimated to be quite old (11 Gyr;Campante et al. 2015), implying the real planets are stable on timescales much longer than our ∼ 50Myr integrations despite most of those integrations ending in instability.Possibly unseen planets or the influence of the companion stars play a role in stabilizing the real system.
Kepler-1542 : There are four confirmed transiting planets in the Kepler-1542 system.We ran simulations of this system using the original reported planetary radii (Morton et al. 2016) for the MR masses.We found that 77 of the 100 simulations were unstable within 40 Myr.Berger et al. (2018) reported planetary radii for this system that are consistently ∼ 20% larger than those reported earlier, so we ran a second set of simulations with MR masses based on the updated radii which resulted in 70 of the 100 simulations being unstable; the stability rate is not meaningfully improved with the updated radii.The Exoplanet Archive lists one additional, slightly larger period unconfirmed planet candidate for this system; a future investigation might consider whether such an additional planet would improve system stability.K2-384 : There are five confirmed transiting planets in the K2-384 system (Christiansen et al. 2022).Our simulations assuming MR masses for the planets show an instability rate similar to the 5-planet version of K2-138 above (75 of the 100 simulations end in collisions).
The presently known planets have orbital periods ranging from 2-14 days.We highlight this system as one where the dynamical architecture could point to an additional, stabilizing planet that is presently undetected and could be detected in future observations.

RESULTS AND DISCUSSION
From our analysis above, we can identify a few factors that appear to be important in determining long-term stability: (1) In most of the unstable simulations, the majority of neighboring planets have period ratios below 2 (see, e.g., Figure 7), which increases the chances that planetary mean motion resonances can overlap at moderate to large eccentricities.( 2) The unstable simulations have distributions of dynamical separations of neighboring planets that skew toward smaller separations (though ∼ 30% of collisions in these systems occur between planets at initial separations larger than 20 Hill radii).( 3) Most of the unstable simulations have larger spectral fractions, meaning that the system's secular evolution is not dominated by a small number of secular modes; this increases the chances that secular chaos can increase eccentricities over time.( 4) The unstable simulations tend to have slightly larger coefficients of variation in planetary masses, meaning the planets are less likely to have very similar masses; this trend is weak, but large inter-planet mass variations might contribute to larger secular eccentricity increases.
We can explore the relative influence of these planetary system features for predicting long-term stability using machine learning.We trained and tested a simple gradient boosting classifier 5 described at https: //scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html#sklearn.ensemble.GradientBoostingClassifier.For simplicity, we left all input parameters set to their default values.using the entire ensemble of simulations that survived long enough for spectral fractions to be calculated (∼ 20000 simulations, of which ∼ 2000 were unstable on less than 5 × 10 9 orbits of the innermost planet).We split our simulation set into a training set consisting of 60% of the simulations and a testing set of the remaining 40%, and trained the classifier to predict stability based on different combinations of data features from the simulation set; the training and testing sets remained fixed for all combinations of data features.To characterize the performance of the classifier, we use three parameters: 1) the positive predictive value (PPV) of the classifier predicting a system is unstable (i.e., the fraction of systems predicted to be unstable that are actually unstable); 2) the classifier's true positive rate (TPR) for identifying unstable systems (i.e., the fraction of all unstable systems correctly classified as unstable); and 3) the area under the receiver operating characteristic curve (ROC-AUC), a summary statistic based on a binary classifier's true positive and false positive rates that is 1 for a perfect classifier and 0.5 for random chance.We first tested the predictive power of several data features individually: the largest amd spectral fraction amongst the planets in a simulation, the coefficient of variation for the masses of the planets, the average period ratio for neighboring planets, and the average dynamical separation of neighboring planets.The PPV, TPR, and ROC-AUC for each of these single-feature classifiers is given in Table 1.
The amd spectral fraction and average period ratio both performed much better than average ∆ or mass C v .Training a classifier using the combination of amd spectral fraction and average period ratio yields a PPV of 0.79, TPR of 0.67, and ROC-AUC of 0.82; adding mass C v and average ∆ results in very small improvements.When we add in a variety of other easily calculated data features (system-wide minimum ∆, minimum and maximum period ratios, the coefficients of similarity and variation for mass and radius, the maximum spectral fractions for a, e, and i) we get modest improvements in the classifier's performance: a PPV of 0.81, TPR of 0.78, and ROC-AUC of 0.88.Including our assessment of whether planets in these systems are in, near, or not near MMRs does not change the classifier's performance in these metrics.This supports the idea that MMRs are not always directly driving the instabilities.As noted in Section 2.3, the near-resonant systems do have a higher instability rate, but do not account for the majority of the instabilities.However, MMRs may play a role in destabilizing systems that experience eccentricity growth due to secular interactions, especially in the systems with smaller period ratios.
Overall, our machine learning analysis shows that the amd spectral fraction is a particularly useful and easily calculated instability indicator alongside the period ratios and dynamical separations of nearest-neighbor planets.The importance of the amd spectral fraction in predicting (in)stability supports the finding in Volk & Malhotra (2020) that many instabilities seen in simulations can be attributed to secular chaos.There are additional potential drivers of instability in multi-planet systems that warrant additional investigations, and there remain many insights to be gained from detailed studies of stability as a function of planetary spacings, masses, and orbital period distributions (see, e.g., some recent ex-Table 1.The positive predictive value (PPV), true positive rate (TPR) and area under the receiver operating characteristic curve (ROC-AUC) for machine learning classifiers trained using the listed planetary system parameters to make a binary prediction for whether a system is stable or unstable.

CONCLUSIONS
We have analysed a large set of simulated planetary systems with three or more planets based on architectures of observed exoplanet systems.Systems that experience instabilities within 5 × 10 9 orbits of the innermost planet tend to share the following properties: • The secular evolution is not dominated by just a small number of secular frequencies, but instead has many frequencies in the evolution of the planets' eccentricities and inclinations.
• The system includes one or more adjacent planet pairs of period ratio below two.
• Smaller planetary separations (measured in units of mutual Hill radius) and more unequal planetary masses are more common, though they are not necessary for instability.
Spectral fractions that characterize the secular evolution of a system based on very short integrations can provide a prediction of long-term stability, especially when combined with other slightly weaker predictive features such as period ratios, planet dynamical spacings, and characterizations of the mass distribution.These stability predictors might allow dynamical refinement of mass estimates for some observed planets or identify systems whose stability could be enhanced by as-yet-undetected additional planets.
Facilities: Exoplanet Archive, ADS Software: rebound,numpy (Harris et al. 2020),scikitlearn APPENDIX A. SYSTEM DETAILS Table 2 lists details for every set of 100 simulations we ran (described in Section 2).Note that many systems appear more than once in the table; for example, systems with TTV/RV planet mass constraints appear once for the simulation set based on those constraints and once for the simulations set based on MR relationship planet masses.Some systems also appear with differing numbers of planets.These are systems for which we ran an initial simulation set based on data available in 2019 but have subsequently had additional planets detected in the system.Note-If the value of n pl for a system is marked with an asterisk, at least one additional planet was detected in that system since its data was queried for the initial set of simulations; in these cases the system has a second entry in the table/simulation set with the additional planet included.In all cases except one (Kepler-271; see notes for that system), the orbital periods for all planets were taken to be the default value and uncertainties in the Exoplanet Archive composite planet table 6 ; in almost all cases, the orbital periods for each planet from different lightcurve analyses agree within uncertainties.Unless otherwise indicated in the notes column, inclinations and eccentricities for all planets were assigned from the Rayleigh distributions described in Section 2. MR masses are generated from the statistical mass-radius relationship from Wolfgang et al. (2016) for planets with radii of 4R⊕ or smaller.For planets exceeding that size limit, we either use a TTV/RV mass measurement (indicated in the notes column) or we use the mass and uncertainty given in the Exoplanet Archive Composite data

Figure 2 .
Figure2.MR-relationship cases: Distribution of mutual hill radius separation (∆) between neighboring planets in 162 observed planetary system architectures with 3 or more planets (top: stacked histogram with simulations weighted to sum to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.1).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.

Figure 3 .
Figure3.TTV/RV cases: Distribution of mutual hill radius separation (∆) between neighboring planets in 47 observed planetary system architectures with 3 or more planets and TTV or RV mass constraints on at least one planet (top: stacked histogram weighted to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.2).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.

Figure 4 .
Figure 4.All simulations: Distribution of mutual hill radius separation (∆) between neighboring planets (top: stacked histograms that sum to the observed number of planet pairs; bottom: normalized cumulative distributions).Note the log scale in the top panel used to better highlight the distribution of the colliding pairs.Colliding planet pairs are shown in orange with the other adjacent planet pairs in those simulations shown in black.Adjacent planet pairs in simulations with no collisions are shown in blue.Note that a few percent of the colliding pairs of planets were not adjacent at the start of the simulations.

Figure 5 .
Figure5.MR-relationship cases: Distribution of period ratios between neighboring planets in 162 observed planetary system architectures with 3 or more planets (top: stacked histogram that sums to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.1).In pink are period ratios in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
Figure6.TTV/RV cases: Distribution of period ratios between neighboring planets in 47 observed planetary system architectures with 3 or more planets and TTV or RV mass constraints on at least one planet (top: stacked histogram that sums to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.2).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.

Figure 7 .
Figure 7.All simulations: Distribution of period ratios between neighboring planets (top: stacked histogram that sums to the observed number of planet pairs; bottom: normalized cumulative distributions).Note the log scale in the top panel used to better highlight the distribution of the colliding pairs.Colliding planet pairs are shown in orange with the other adjacent planet pairs in those simulations shown in black.Planet pairs in simulations with no collisions are shown in blue.Note that a few percent of the colliding pairs of planets were not adjacent at the start of the simulations.

Figure 8 .
Figure8.MR-relationship cases: Distribution of mass ratios between neighboring planets in 162 observed planetary system architectures with 3 or more planets (top: stacked histogram that sums to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.1).In pink are period ratios in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.
Figure9.TTV/RV cases: Distribution of mass ratios between neighboring planets in 47 observed planetary system architectures with 3 or more planets and TTV or RV mass constraints on at least one planet (top: stacked histogram that sums to the observed number of planet pairs; bottom: normalized cumulative distributions).In grey are pairs of planets in system architectures where no collisions occurred within 5 × 10 9 orbits of the shortest period planet across 100 instances of that architecture (initial conditions described in Section 2.2).In pink are planet pairs in system architectures for which at least one of the 100 instances experienced a collision within 5 × 10 9 inner planet orbits.

Figure 10 .
Figure 10.All simulations: Distribution of mass ratios between neighboring planets (top: stacked histogram weighted to the observed number of planet pairs; bottom: normalized cumulative distributions).Note the log scale in the top panel used to better highlight the distribution of the colliding pairs.Colliding planet pairs are shown in orange with the other adjacent planet pairs in those simulations shown in black.Planet pairs in simulations with no collisions are shown in blue.Note that a few percent of the colliding pairs of planets were not adjacent at the start of the simulations.

Figure 11 .
Figure 11.MR-relationship cases: Fraction of stable simulations as a function of the coefficient of variation for the masses of the planets in each simulation (top) and number of simulations per bin in that coefficient (bottom).

Figure 12 .
Figure 12.All simulations: Fraction of stable simulations as a function of the largest amd spectral fraction for any planet in that simulation (top) and number of simulations per spectral fraction bin (bottom).Planetary systems without planet pairs starting obviously near strong mean motion resonances (systems with MMRs labeled as 0 or 2 in Table 2) are shown as blue dots, while those that are likely in or near resonant are shown as red squares (systems with MMRs labeled as 1 in Table2).These plots include all simulations performed (whether the planets masses were assigned from the mass-radius relationship or from TTV/RV measurements) that survived for at least 10 6 orbits, the minimum timespan to perform a useful FFT.

Figure 13 .
Figure13.All simulations: Average simulation stopping time (in units of inner-most planet orbital periods) vs maximum amd spectral fraction for systems not obviously near mean motion resonances (blue dots) and those that are likely resonant or near-resonant (red squares).The standard deviation of the simulation lengths in each spectral fraction bin are shown with black error bars, except for the right-most three bins, which only have one simulation per bin.Note that the two sets of data points are artificially offset on the horizontal axis to better differentiate the points.Simulations ran until 5 × 10 9 orbits of the innermost planet (grey horizontal line) or until two planets collided.

Figure 14 .
Figure 14.Cumulative distributions of maximum amd spectral fractions.
and eccentricities for b-f, planet radii, and stellar mass fromLissauer et al. (2013); MR mass for planet g except for planet e; planet radii, TTV mass for planet e, and stellar mass from Lissauer et al. for all planets except g, TTV eccentricities for planets d, e, f, planet radii and stellar mass fromLissauer et al. (2013); MR mass for planet planets b,c, d, g, eccentricities for planets b, c, g, and radii for planets b, c, d, e, f fromBuchhave et al. (2016); MR masses for planets e and f; planet g does not transit and its radius is taken from the Exoplanet Archive's estimate based on a mass-radius; relationship; stellar mass fromMorton et al. (2016); unall planets except g; planet radii and RV mass for planet g fromBuchhave et al. (2016); planet g does not transit and its radius is taken from the Exoplanet Archive's estimate based on a mass-radius; stellar mass fromMorton et al. (2016); unfor b, c, d, e from MacDonald et al. (2016); MR planet f mass; planet radii and stellar mass from MacDonald et al. masses; planet radii and stellar mass fromChristiansen et al. (2022) (mass for planet c from alternate MR relationship); planet radii from Rowe et al. (2014); stellar mass from Morton et al. ; planet radii for b, c from Rowe et al. (2014); planet radii for d, e, f and stellar mass from Morton et al. masses; planet b, c, d, e radii from Rowe et al. (2014); planet f radius and eccentricity from Torres et al. (2015); stellar mass from Morton et al. (2016) masses; planet b, c, d, e radii from Rowe et al. (2014); planet f radius from Torres et al. (2015); stellar mass from Morton et al. (2016) and radii for e, f from Xie (2014); MR planet masses for b, c, d; radii for b, c, d from Rowe et al. (2014); stellar mass from Morton et alfor all planetx except e; planet e mass from Xie (2014); radii for b, c, d from Rowe et al. (2014); radii for e,f from Xie (2014); stellar mass from Morton et al. ; planet radii and stellar mass from Fabrycky et al. for b, c from Hadden & Lithwick (2014); MR planet masses for d, e, f; planet radii and stellar mass from Fabrycky et al. all planets except b from Hadden & Lithwick (2016); MR mass for planet b; planet radii and stellar mass from Lissauer et al. ; planet eccentricities, radii, and stellar mass from Campante et al. (2015); there are two unmodelled stellar ; planets d, e, f radii from Rowe et al. (2014); planets b, c radii from Steffen et al. (2013); stellar mass from Morton et al. , c masses from Hadden & Lithwick (2014); MR planet masses for d, e, f; planets d, e, f radii from Rowe et al. (2014); planets b, c radii from Steffen et al. (2013); stellar mass from Morton et al. and radii for b, c, f from Freudenthal et al. (2019); MR planet masses for d, e; planet d, e, radii from Berger et al. (2018); l stellar mass fromFreudenthal et al. (2019) (2014); stellar mass from Morton et al. ; planets b, c radii from Xie (2014); planets d, e, f radii from Rowe et al. (2014); stellar mass from Morton et al. for b, c from Hadden & Lithwick (2014); MR planet d, e, f masses; planets b, c radii from Xie (2014); planets d, e, f radii from Rowe et al. (2014); stellar mass from Morton et al. for d, e from Rodriguez et al. (2018); MR planet masses for b, c; planet radii and stellar mass from Rodriguez et al. (2018); unmasses; planet radii and stellar mass from Rodriguez et al. (2018); unmasses; planet radii and stellar mass from Palle et al. (2019) for b, c, planet radii, and stellar mass from Palle et al. (2019); MR masses for dfor all planets and eccentricities for planets d, e, f, planet radii, and stellar mass from Lillo-Box et al. (2020); eccentricities from Van Eylen & Albrecht (2015); planet radii from Rowe et al. (2014); stellar mass from Morton et al. for b, c, e from Bonomo et al. (2019); MR mass for d; planet radii fromRowe et al. (2014); stellar mass from Morton et al. ; planet radii and stellar mass from Morton et al. ; planet radii and stellar mass from Morton et al. (mass for planet e from alternate MR relationship); planet radii and stellar mass from Kipping et al. (2016); un-modeled stel-; planet b, c, d radii from Rowe et al. (2014); planet e radius and stellar mass from Morton et al. for c from Hadden & Lithwick (2014); MR masses for planets b, d, e; eccentricity for planet e from Van Eylen & Albrecht (2015); planet radii from Rowe et al. (2014); setllar mass from Morton et al. (2016); un-modeled stellar companion Ford et al. (2012); radii for planets d, e fromRowe et al. (2014); stellar mass from Morton et al. for b, c from Hadden & Lithwick (2014); MR planet masses for d, e; planet b, d, c radii from Berger et al. (2018); planet e radius and stellar mass from Morton et al. ; planet b, c, d radii from Rowe et al. (2014); planet e radius and stellar mass from Morton et al. for b, c from Jontof-Hutter et al. (2016); MR planet masses for d, e; radii for d, e, from Rowe et al. (2014); radii for b, c and stellar mass from Jontof-Hutter et al. ; radii for d, e, from Rowe et al. (2014); radii for b, c and stellar mass from Jontof-Hutter et al. ; planets b, c radii from Rowe et al. (2014); planets d, e radii from Xie (2014); stellar mass from Morton et al. and radii for d, e from Xie (2014); MR planet masses for b, c; radii for b, c from Rowe et al. (2014); stellar mass from Morton et al. ; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016) masses; planet b, c, d radii fromRowe et al. (2014); planet e radius and stellar mass fromMorton et al. (2016) ; planet radii for c, b and stellar mass from Crossfield et al. (2016); planet d radius from Mayo et al. for c, d; planet radii, RV mass for planet b, and stellar mass from Marcy et al. ; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016); un-modeled stellar companion for c, d from Xie (2014); MR planet mass for b; planet radii from Rowe et al. (2014); stellar mass from Morton et al. ; planet b radius form Rowe et al. (2014), planets c and d radii from Xie (2014); stellar mass from Morton et al. , planet eccentricities, and stellar mass from Jontof-Hutter et al. (2015) ; planet radii for b, c from Rowe et al. (2014); planet radius for d and stellar mass from Morton et al. ; planet radii for b,c from Rowe et al. (2014); planet radius for d and stellar mass from Morton et al. masses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016) masses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016) masses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. planet masses and planet c and d eccentricities from Malavolta et al. (2017); planet b radius and stellar mass from Ballard et al. (2011); planets c and d do not transit, so radii are taken from the Exoplanet Archive's estimate based on a mass-; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. ; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. ; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. ; planet radii from Rowe et al. (2014); stellar mass from Morton et al. ; planet radii from Rowe et al. (2014); stellar mass from Morton et al. ; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. ; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. from Hadden & Lithwick (2014); planet b, c, radii from Ford et al. (2012); planet d radius from Rowe et al. (2014); stellar mass from Morton et al. (2016) masses; planet b, c, radii from Ford et al. (2012); planet d radius from Rowe et al. (2014); stellar mass fromMorton et al. (2016) , radii for planets b, c and stellar mass fromMarcy et al. (2014); planet d does not transit, so its radius is taken from the Exoplanet Archive's estimate based on a mass-radius relationship; for b, c; radii for planets b, c, RV mass for planet d, and stellar mass fromMarcy et al. (2014); planet d does not transit, so its radius is taken from the Exoplanet Archive's estimate based on a mass-radius relationshipmasses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016) masses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016) ; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. ; planet radii, planet orbital periods, and stellar mass from Morton et al. (2016); note that conflicting solutions for planet b's orbital period have been published, and there are two additional unconfirmned and unmodeled planet candidates in this system ; planet b radius from Rowe et al. (2014); planets c, d radii from Xie (2014); stellar mass from Morton et al. ; planet radii from Rowe et al. (2014); stellar mass from Morton et al. planets b, d; TTV mass for planet c and planet radii from Schmitt et al. (2014); stellar mass from Morton et al. and eccentricities and planet radii from Sanchis-Ojeda et al. (2012); stellar mass from Fabrycky et al. planet b; TTV planet masses for c, d and planet radii from Sanchis-Ojeda et al. (2012); stellar mass fromFabrycky et al. (2012) ; planets b, c radii from Xie (2014); planet d radius from Rowe et al. (2014); stellar mass from Morton et al. masses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016) masses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. (2016) masses; planet radii from Rowe et al. (2014); stellar mass from Morton et al. and planet radii for c, d from Xie (2014); MR planet mass for b; planet b radius from Rowe et al. (2014); stellar mass from Morton et al. ; planet radii for c, d from Xie (2014); planet b radius from Rowe et al. (2014); stellar mass from Morton et al. ; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. ; planet radii from Rowe et al. (2014); stellar mass from Morton et al. ; planet b, c radii from Rowe et al. (2014); planet d radius and stellar mass from Morton et al. ; planet radii from Rowe et al. (2014); stellar mass fromMorton et al. (2016)

Table 2 .
List of system architectures used in this work.

Table 2 continued
Table 2 (continued) ; planets b, c radii from Xie (2014); planet d radius from Van Eylen & Albrecht (2015); stellar mass from Silva Aguirre et al. (2015) table, which is based on the massradius relationship from Chen & Kipping (2017), see https://exoplanetarchive.ipac.caltech.edu/docs/composite_calc.html; for systems where we used this mass estimate, the specific planet is indicated in the notes column as coming from an alternate MR relationship.a 0 indicates no planets are near strong MMRs (mean motion resonances); 1 indicates one or more pair of planets is in or very close to a strong MMR; 2 indicates one or more pair of planets is close enough to a strong MMR that resonant interactions cannot entirely be ruled out (especially if simulated planets evolve to large eccentricities).