Statistical Distribution Function of Orbital Spacings in Planetary Systems

The minimum orbital separation of planets in long-stable planetary systems is often modeled as a step function, parameterized with a single value Δmin (measured in mutual Hill radius of the two neighboring planets). Systems with smaller separations are considered unstable, and planet pairs with greater separations are considered stable. Here we report that a log-normal distribution function for Δmin , rather than a single threshold value, provides a more accurate model. From our suite of simulated planetary systems, the parameters of the best-fit log-normal distribution are μ = 1.97 ± 0.02 and σ = 0.40 ± 0.02, such that the mean, median, and mode of Δmin are 7.77, 7.17, and 6.11, respectively. This result is consistent with previous estimates for Δmin threshold values in the range ∼5–8. We find a modest dependence of the distribution of Δmin on multiplicity within the system, as well as on planetary mass ratios of the closest planet pair. The overall distribution of nearest-neighbor planetary orbital spacings (measured in the mutual Hill radii and denoted simply as Δ) in long-term stable systems is also well fit with a log-normal distribution, with parameters μ = 3.14 ± 0.03 and σ = 0.76 ± 0.02. In simulations of sets of many planets initially packed very close together, we find that the orbital spacings of long-term stable systems is statistically similar to that in the observed Kepler sample of exoplanetary systems, indicating a strong role of sculpting of planetary architectures by dynamical instabilities.


INTRODUCTION
NASA's Kepler mission (Borucki et al. 2010) has discovered a very large number of exoplanetary systems that can be studied in a consistent manner for planetary system demographics.One of the properties of particular interest for planetary system demographics is the orbital spacings (equivalently, orbit period ratios) of nearest-neighbor planets.Studies of the Kepler population statistics have shown that the distribution of orbital period ratios of adjacent planets has an overall broad peak and several narrow peaks associated with low-integer orbital period ratios (Fabrycky et al. 2014;Malhotra 2015;Mulders et al. 2018).The distribution of orbital spacings observed in the Kepler sample peaks near 10-20 mutual Hill radii, with significant dependence on multiplicity (Pu & Wu 2015;Gilbert & Fabrycky 2020).There is also a long tail of wide spacings, indicative of large gaps between planets where ad-ditional undetected planets may exist (e.g., Gilbert & Fabrycky 2020;Dietrich & Apai 2020).He et al. (2019) generated sets of simulated "Kepler-like" planetary systems, simulated observing them with the Kepler detection pipeline, and reported that the best-fit models to the actual observational data requires minimum orbital spacings of about 8 mutual Hill radii.In the present study, our goal is to obtain a better characterization of the statistical distribution of the minimum orbital spacing of planets in long term stable systems, including its dependence on the mass ratios of the planets and on planet multiplicity in the system.
We will use the parameter ∆ defined as the orbital spacing in units of the mutual Hill radius: ∆ = 2(a 2 (1 − e 2 ) − a 1 (1 + e 1 )) where a is the semi-major axis, e is the eccentricity, i is the mutual inclination, and m is the mass of each planet.The subscripts 1 and 2 indicate the inner and outer planet, respectively, and M * is the mass of the host star.(Note that the inclination does not appear in this definition.)Other studies also denote this dynamical spacing parameter as K (see e.g., Pu & Wu 2015;Malhotra 2015).We will specifically use the notation "∆ min " for the minimum value of ∆ taken over all nearest-neighbor planet pairs in a planetary system; unless otherwise indicated, this will be the value at an initial epoch.For a system of two planets of initially circular and coplanar orbits, the theoretical minimum value of ∆ for Hill stability is given by (Birn 1973;Gladman 1993;Tremaine 2023) In the more general case, there is no known simple, analytic stability criterion for multi-planet systems, therefore insights have been sought with numerical simulations.
Numerical studies of the dynamical stability of multiplanet systems do not seek a definitive statement over all time, but only over a finite time range, usually limited by computational resources.In modern times, such investigations have extended to timescales of 10 7 − 10 9 orbital periods of the innermost planet in the system (e.g., Chambers et al. 1996;Smith & Lissauer 2009).Numerical techniques that reduce computational costs and reliably predict long-term stability based on the short-term behaviors of orbital parameters (inclination, eccentricity, and angular momentum) are also being developed (see e.g., Tamayo et al. 2020;Volk & Malhotra 2020).
In reviewing the literature on numerical simulations of the stability of multi-planet systems, Malhotra (2015) reported that for systems of 3-20 equal-mass planets in initially coplanar and nearly circular orbits, stability times exceeding 10 8 orbital periods require ∆ min ≥ 8 for low mass planets (planet-to-star mass ratio in the range 10 −9 -10 −5 ), whereas systems of higher planet masses (planet-to-star mass ratio in the range 10 −3.4 -10 −2.4 ) are stable with somewhat smaller minimum separations, ∆ min ≥ 5.
In almost all previous modeling studies, the minimum value of ∆ required for long-term stability has been treated as a simple cutoff, where systems with ∆ min below some threshold are considered unstable and systems with ∆ min above that threshold considered stable.However, even a cursory examination of previous numerical results shows that a single cutoff value is merely an approximation, and that this threshold is not a single valued parameter but has a range of values that may depend upon many characteristics of a planetary system, such as the stellar mass, the number of planets, their masses, orbital eccentricities, mutual inclinations, their ordering, and the mass ratio of the planet pair with the minimum orbital separation.Because there are potentially many dependencies on parameters that are themselves random variates, we make the reasonable hypothesis that the threshold value of ∆ min can be better characterized as a random variate whose distribution can be illuminated with N-body numerical simulations.
Our approach is to carry out N-body numerical simulations of a large set of simulated planetary systems with close-packed planets, measure the minimum dynamical separation, ∆ min , in the long-term stable systems, and thereby quantitatively determine a more accurate characterization of the distribution of the threshold value of this parameter.This paper is organized as follows.We describe the methods we used to create the synthetic planetary system populations in Section 2. The results of the N-body integrations to assess the dynamical stability of these systems are described in Section 3. In Section 4 we discuss our results for the distribution of ∆ min and its dependence on planet mass and multiplicity.We summarize our findings in Section 5.

METHODS
Our general process for determining the distribution of the minimum stable value ∆ min with simulated planetary systems is the following: 1. We generate our initial planetary systems based on prior model distributions for planet multiplicity, orbital period, and planet mass.
2. We use N-body integrations to evolve each system in time for 1 million orbits of the innermost planet and measure the parameters at evenly spaced time intervals.
3. Systems that become unstable (either via ejection, collision, or final orbits that lead to instability at later times) are separated.
4. The systems (both stable and unstable) are ordered based on the minimum initial ∆ value ∆ min from the original system.
5. We also search for correlations between ∆ min and the period and mass ratios of the corresponding planet pair, as well as the planet multiplicity of the system.
To generate our simulated planetary systems, we needed to initialize the number of planets in a system as well as the orbital periods, masses, inclination, and eccentricity for each of the planets.The distribution of orbital periods between planet pairs is motivated by the Kepler population of exoplanets (Mulders et al. 2018).The Kepler sample from Mulders et al. (2018) includes planets with orbital periods in the range [2, 400] days and planet radii in the range [0.5, 6] R ⊕ , with cuts on the stellar population to remove giant and sub-giant stars via the effective temperature dependent surface gravity criterion (Huber et al. 2016).This led to 1201 total planets in 487 multi-planet systems, with an additional 1840 observed single systems.The other parameters were initialized in a wide range so as to account for as many different planetary system architectures as possible.Specifically, we chose the following parameters and ranges for each simulated system: • Number of planets: randomly chosen with a uniform prior distribution between 3 and 10 planets in the system.
• Orbital periods: the first planet's orbital period randomly chosen with a uniform prior in log period space between 1 to 100 days.Additional planets had their period ratio with the previous planet drawn from the distribution of observed pairwise period ratios in the Kepler multi-planet systems, which are taken from the sample in Mulders et al. (2018).
• Masses: randomly chosen for each planet with a uniform prior in log mass between 1/3 Earth mass and 3 Jupiter masses.
• Inclination and Eccentricity: We set each planet to have a coplanar and circular orbit, but these parameters are allowed to evolve due to interactions during the N-body integrations.
• Angular Variables: The longitude of ascending node, argument of periastron, and mean anomaly were initialized for each planet randomly from a uniform distribution across the range [0, 2π].
Based on provisional results for the initial distributions of ∆ min , we generated two different sets of simulated planetary systems.The most significant difference between the two sets is how closely packed the systems are at the initial time.The first set, identified as the Nominal Kepler Analog Set, uses the period ratio distribution from Kepler as described above, whereas the second set -the Very Closely Packed Set -increases the number of neighboring planet pairs with ∆ ≲ 8.

Nominal Kepler Analog Set
We generated a set of 10,000 simulated planetary systems using the procedure described above and calculated ∆ between each of the nearest-neighbor planet pairs (see Figure 1).This set of draws for the planetary parameters yields an approximately log-normal distribution of ∆ with parameters µ = 3.13 ± 0.03 and σ = 0.79 ± 0.02.The log-normal distribution has the following probability density function, such that the natural logarithm of x is normally distributed with mean µ and standard deviation σ.Note that we only use "log" as the natural logarithm within the name of this distribution; for all other uses "log" refers to the base-10 logarithm.This functional form was previously introduced by Malhotra (2015) in an empirical way for the distribution of dynamical spacings of planetary systems.We also tested the beta, exponential, gamma, Lorentzian, and skewed normal distributions, but found that the log-normal function was the best fit as measured by residual sum of squares and/or reduced chi-square metrics.The best fit model parameters and their uncertainties are derived using the Python module lmfit to fit the distribution, with up to 600 repetitions (200 × n for n = 3 parameters in the fit: amplitude, µ, and σ) using the Levenberg-Marquardt algorithm (Levenberg 1944;Marquardt 1963).The modest but visible mismatch between the fit and the distribution is a result of the hard lower limit we placed on the initial ∆ values (see Eq. 2).If we remove that limitation, a small fraction (∼2%) of additional systems are generated with ∆ < ∆ min , mostly removing the mismatch.

Very Closely Packed Set
Our initial setup (see previous paragraph) yielded a low number of closely separated nearest neighbor planet pairs (i.e., planet pairs with ∆ < 8) due to the distribution of Kepler period ratios peaking at ∼1.9.This is an understandable result of a survival bias: The Kepler period distribution is derived from the observational data of systems that are presumably long-term stable if we observe them now, and thus few are marginally unstable.Therefore, to mitigate the bias against low values of ∆ in our Nominal Kepler Analog Set, we also generated additional systems to fill in the low ∆ bins to obtain a nearly uniform distribution in log ∆ at ∆ ≲ 10; the latter is the value of ∆ min,Kepler (Mulders et al. 2018).While this overly tuned our simulated planetary systems to be very closely packed, it allowed us to determine if the ∆ min distribution of the dynamically stable systems that survive until the end of our N-body integrations are determined more by dynamical stability itself or by the bias in the prior distribution of planets.Our procedure for generating the additional systems was as follows: • We generated 10,000 additional systems based on the same criteria as above, except the period ratios were drawn from a uniform distribution in the range [1, 3], with the imposed condition that all resulting ∆ values, for every nearest-neighbor planet pair in each system, were less than 10.
• We accepted a system only if it did not push the number of systems for the ∆ bins that were much lower than the peak of the initial distribution (at ∆ ≈ 10) above that peak value in the histogram.This filled in much of the left side of the initial distribution, but still produced very uneven results in some of the bins in log ∆ due to low number statistics from the low probability to generate certain planet pair parameter combinations.
• To achieve a more even range across all the bins, we again generated additional systems with the same planet multiplicity and planet mass criteria, except that we further limited the period ratio draws to only produce ∆ values within the bins with fewer systems.
This method of generating simulated planetary systems yielded a relatively uniform initial distribution (see Figure 2) of log ∆ min in the range ∆ = 2 √ 3 (the theoretical minimum separation for Hill stability of two planets) to ∆ ∼ 10 (the peak of ∆ min of the Nominal Kepler Analog Set).Of the other variables we track (multiplicity, period ratios, mass ratios), this procedure mostly affected only the period ratio distribution by increasing the number of low period ratios between the nearestneighbor planet pairs.The multiplicity distribution and mass ratio distribution between nearest-neighbor planets remained similar to those from the Nominal Kepler Analog set.

N-body Integrations
We ran each system for both sets of simulated planetary systems through N-body integrations using the mercurius integrator in the REBOUND package (Rein & Liu 2012;Rein et al. 2019).We integrated each system over 5 million orbits of the innermost planet, recording the orbital elements at 3,000 evenly spaced time intervals along the way.At the end of those integrations, we identified each system as unstable or stable based on whether collisions (two planets' physical radii overlapping) and/or ejections occurred or did not occur during the time span of the integration.We note that stability for up to 5 million orbits is not necessarily indicative of long-term stability on the timescale of billions of years, but this timescale of integration allows to run thousands of N-body simulations over a reasonable computational time period.
For the systems from both sets that were deemed stable (i.e., experiencing no collisions or ejections), we also ran the spectral fraction analysis designed by Volk & Malhotra (2020) to predict longer-term stability from short-term numerical integrations.This analysis computes the power spectrum of the orbital parameters in the exoplanet system propagated for timescales of ∼ 10 6 orbits of the innermost planet and looks for empirical indicators of long-term stability in the power spectrum, specifically the fraction of the spectral frequencies having power above a certain threshold fraction of the power of the peak frequency.The best empirical measure of stability was found to be as follows: with a threshold value of 5% of the peak power, a spectral fraction below 1% is indicative of long-term stability (on at least (10 9 ) orbital periods), whereas a larger spectral fraction is very likely long-term unstable.
We also checked the stability criterion as measured by SPOCK (Tamayo et al. 2020), which uses shorter N-body integrations (sim10 4 orbits) and machine learning algorithms trained on typical stable systems to predict long-term stability.In general, we found that the SPOCK analysis found more systems unstable than the spectral fraction analysis or the ∼ 10 6 orbit N-body survival rate, similar to results on a single system found by Dietrich et al. (2022).Unsurprisingly, SPOCK found a much higher fraction of the Very Closely Packed sample to be unstable compared with the Nominal Kepler Analog sample; this is because the former includes a much larger population of low ∆ systems.

RESULTS
For both the Nominal Kepler Analog Set and the Very Closely Packed Set, the overall fraction of unstable systems is less than 50%, but significantly larger in the latter set.The Nominal Kepler Analog Set had an unstable system fraction of ∼11%, whereas the Very Closely Packed Set had an unstable system fraction of ∼37%.The unstable cases for both sets of simulated planetary systems are concentrated around a ∆ min ≲5 within the system.The spectral fraction analysis found a small subset (≤5%) of originally stable systems to be classified as unstable in both the Nominal Kepler Analog Set and the Very Closely Packed Set.

Nominal Kepler Analog Set
We calculated the best-fit initial distributions of ∆ and ∆ min for both the stable and unstable planetary systems, as seen in Figure 3.The stable systems' initial ∆ distribution essentially matched the distribution of all systems in this set, as over 90% of systems were found stable.The ∆ min distribution most closely matched a log-normal distribution with parameters µ = 1.97 ± 0.02 and σ = 0.40 ± 0.02.These values correspond to a median value of ∆ min of 7.17 ± 0.15, a root-mean-square value ⟨∆ 2 min ⟩ 1 2 = 8.4 ± 0.3, and the range (4.8, 10.7) contains two-thirds of ∆ min values.The distribution for the unstable systems' initial ∆ min was also fit to a log-normal distribution, with parameters µ = 1.39 and σ = 0.13.

Initial vs final ∆
The final ∆ values for the unstable planetary systems showed a wide dispersion, ranging from ∆ < 0.001 for neighboring planet pairs about to collide, to ∆ ≳ 1000 for an ejected planet in originally-neighboring planet pairs, as seen in the top left panel of Figure 4.For the long-term stable systems, the ∆ distribution changes very little from the initial to the final state.This is an expected outcome since long-term stability is an indicator of few interactions and/or little orbital evolution that would change the planet pair spacings.In the unstable systems, a majority of the surviving planet pairs have final values of ∆ near 10, similar to those of the stable systems.

Planet multiplicity
The stable systems with higher planet multiplicity had consistently smaller initial ∆ min values, as seen in the top right panel of Figure 4 by the green triangles denoting the median of the ∆ min distributions for each planet multiplicity.The stable fraction also decreases monotonically with planet multiplicity, with a large fraction (> 80%) maintaining stability even at high multiplicities (N ≳ 8), as seen in the bottom left panel of Figure 4.

Planet pair mass ratio
The bottom right panel of Figure 4 shows a scatter plot of the fraction of stable systems versus the planet pair mass ratio.The stable fraction is found to be in the relatively narrow range of 0.78-0.94,with the highest stable fractions around | log 10 m 2 /m 1 | ∼ 1.There is a weak trend of larger stable fraction in cases where the outer of the planet pair is more massive.There is also a hint of lower stable fraction in cases of nearly equal mass planet pairs and cases where the planet paris are very disparate in mass, | log 10 m 2 /m 1 | ≳ 2. There is large scatter in this plot so these trends are not conclusive.

Very Closely Packed Set
For the Very Closely Packed Set (generated with a near-uniform distribution of log ∆ in the range log 3.5 − log 10), the initial ∆ min for stable systems also follows a roughly log-normal distribution, with parameters µ = 1.98 ± 0.02 and σ = 0.40 ± 0.01.The initial ∆ min for unstable systems follow a log-normal distribution with parameters µ = 1.30 ± 0.02 and σ = 0.30 ± 0.02 (see Figure 5).The stable distribution is very similar to that of the stable subset in the Nominal Kepler Analog Set, but the unstable distribution has a significantly smaller µ and larger σ.There are a few other notable differences that we discuss below.

Initial vs final ∆
The comparison of all initial to all final ∆ values shows the systems with a larger number of closely packed planets mostly tend to have planetary collisions (final values of ∆ ≪ 1 for any neighboring planet pair) or ejections (final values of ∆ ≫ 10 for any neighboring planet pair).There is an even wider dispersion of the final ∆ values of the unstable systems in the Very Closely Packed Set, as can be seen in the top left panel of Figure 6, as the excess of close-separation planet pairs caused a higher number of encounters resulting in collisions or ejections.

Planet multiplicity
The trends in ∆ min with planet multiplicity have some similarities with those found in the Nominal Kepler Analog set, namely, lower initial ∆ min values for all planet pairs in systems with more planets (green triangles in the top right panel of Figure 6).However, in some contrast with the Nominal Kepler Analog set, the fraction of stable systems does not decrease monotonically in the Very Closely Packed Set (as seen in the bottom left panel of Figure 6); rather, it is flat for N = 5 − 7, followed by an increasing trend above N = 9.

Planet pair mass ratio
In the scatter plot of the stable fraction versus planet pair mass ratio (bottom right panel of Figure 6, we observe several features in the Very Closely Packed Set that contrast with what we found for the Nominal Kepler Analog Set.The fraction of stable systems in the Very Closely Packed Set has a range of 0.60-0.85,which is lower and has little overlap with that found for the Nominal Kepler Analog Set; this is not surprising as a larger fraction of the closer packed planetary systems of the Very Closely Packed Set is unstable.Another interesting difference is that we observe a strong pattern in this scatter plot for the Very Closely Packed Set, in contrast with the significant scatter and only weak trends found in the Nominal Kepler Analog Set.The stable fraction is predominantly high for | log 10 m 2 /m 1 | ≈ 2 − 2.5, prominently low in nearly equal mass planet pairs, and trends lower also in cases where the planet pairs are very disparate in mass, | log 10 m 2 /m 1 | ≳ 2.5.There is an evident nearly linearly increasing trend of the stable fraction as | log 10 m 2 /m 1 | increases from ∼ 0 to ∼ 2.  The primary goal of this analysis was to find a more accurate model for the threshold of dynamical stability than a simple cutoff value of ∆ min somewhere greater than the minimum theoretical limit of ∆ * min = 2 √ 3 for the stability of a two planet system.In our procedure for generating initial states, we excluded systems with initial values of ∆ min < 2 √ 3, as these would not have formed (or survived) in a stable region of the parameter space.We also classified as unstable any system which survived with ∆ min < 2 √ 3 at the end of the integrations (which occurred in less than 0.1% of all integrations), or if it did not pass the spectral fraction test for long term stability (which occurred in less than 5% of all integrations).The resulting distributions of initial ∆ min are very similar between the Nominal Kepler Analog Set and the Very Closely Packed Set, and this similarity is true for both the stable systems and the unstable systems.The distributions of ∆ min for the two tested initial planetary spacing methods are lognormals of (1.97, 0.40) and (1.98, 0.40) for the stable systems and log-normals of (1.39, 0.13) and (1.30, 0.30) for the unstable systems.The unstable systems have a significant decrease in µ that defines the center of the distribution, as well as a significantly higher σ denoting increased scatter.
We tested a factor of ∼1.5 more systems with closely separated planets in the second set, such that the initial distribution of ∆ across all the systems was roughly uniform in log space from ∼3.5 − 10 (see Figure 5).Thus, the main difference in the number and fraction of stable vs. unstable systems comes from a higher number of these closer systems being found unstable.The nearly identical distributions of initial ∆ values of the stable cases in the Nominal Kepler Analog Set and the Very Stability in the Very Closely Packed Set Closely Packed Set provides evidence that the orbital spacings in the Kepler sample ∆ Kepler most likely arise from dynamical stability, rather than formation or migration, consistent with the results obtained by Pu & Wu (2015).
The stable systems have a mean ∆ min of ∼7.81 and a median of ∼7.03, which is similar to but slightly smaller than the limiting value of ∆ min = 8 adopted by the SysSim forward modelling framework (He et al. 2019, who also found ∆ min ∼ 10 in additional models).Malhotra ( 2015) noted a value of ∆ min ≳ 8 for equal-mass planets in the low-planetary-mass regime ( mp m * < 10 −5 ) but ∆ min ≈ 5 for planets of higher mass ( mp m * ∼ 10 −3 ), based on results from studies by Chambers et al. (1996) and Zhou et al. (2007).The values reported by both Malhotra (2015) and He et al. (2019) are consistent with our results across the entire set of mass ranges mentioned above, as our distribution of ∆ min peaks at ∼5.69.However, there is no evidence for a correlation between ∆ min and m p -only an increase in instability -when ∆ min in a system is between a neighboring pair containing a higher-mass planet.
Another value for ∆ min from the literature is ∆ min ≈ 10 − 12 for certain assumptions using the Kepler distributions (Pu & Wu 2015).We note that ∆ = 10 is at the 80th percentile of our ∆ min distribution and ∆ = 12 is at the 89th percentile; the result found by Pu & Wu (2015) is consistent with our distribution but significantly near the higher end.In particular, their study includes non-circular and non-coplanar starting orbits with a measured dispersion in eccentricity and mutual inclination, as well as a factor to measure stability over a set timescale.The dispersion in mutual inclinations and eccentricities induced an increase in the dynamical spacings measured (Equation 14 of Pu & Wu 2015), and the dynamical spacings for stability timescales at the length of our N-body integrations, we find a median stable spacing of ∼7 mutual Hill radii (Equation 12of Pu & Wu 2015).
We find no significant difference (i.e., no significant evolution) between the initial and final ∆ min values for both of our sets of stable systems.We allowed en- counters to freely change the orbital inclination and eccentricity of the planets, but since we did not use an initial seed of mutual inclinations or eccentric orbits, all of the planetary systems remained coplanar while a very large majority of the stable systems remained circular.Our results were not significantly different using a mutual inclination distribution (∼ 1 − 5 • ) and an eccentricity distribution (∼ 0.01 − 0.07) that have been empirically found in known stable systems (e.g.He et al. 2020).Therefore, we conclude there must be an additional factor that could describe why the ∆ Kepler distribution as measured by Pu & Wu (2015) has ∆ min ∼ 10 − 12 whereas our Nominal Kepler Analog Set has ∆ min ∼ 6 − 7. One such factor could be the similarity of planet masses within individual Kepler systems (see e.g., Millholland et al. 2017;Weiss et al. 2018;Gilbert & Fabrycky 2020), a correlation which is not well-preserved in either our Nominal Kepler Analog set or our Very Closely Packed Set.

Multiplicity dependence
We noticed in both sets of simulated planetary systems a modest trend of lower ∆ min with increasing planet multiplicity.In generating the initial suite of planetary systems as detailed in Section 2.1, we did not place an upper limit on the outermost planetary period.Therefore, systems with more planets correspondingly had planets with longer orbital periods, and we did not place a restriction on the spacing that would artificially cause a higher number of low-separation pairs in systems with higher multiplicity.Thus, we are able to exclude as a possibility the potential observed degeneracy between high multiplicity and low planetary spacing, which would also have increased the strength of planetary perturbations and led to likely decreased stability (see e.g., Pu & Wu 2015;Rice & Steffen 2023).Other studies have also shown a correlation between higher multiplicity and closer and more ordered spacings of planets (see e.g., Wu et al. 2019).
For stable planetary systems with 10 planets, the initial distribution of ∆ min had a median of ∼5.5 and a 95th-percentile value of ∼12.5, whereas for stable planetary systems with 3 planets the median was ∼9.2 and the 95th percentile was ∼40.This is shown in the middleleft panels of Figures 4 (for the Nominal Kepler Analog Set) and 6 (for the Very Closely Packed Set).The upper end of the range of initial ∆ min values for 3-planet stable systems almost reaches 100, while the upper end for 10-planet stable systems is less than 30.The medians of the initial ∆ min distribution for each planet multiplicity are denoted by the green triangles, which show a significant trend of decreasing initial ∆ min values with increasing planet multiplicity (note that the x-axis is in log scale for this figure).
Therefore, we infer the difference in ∆ min between low-and high multiplicity systems is due to the higher chance of having Jovian-mass planets in the system, creating lower ∆ values.Since the planet masses were drawn log-uniformly from 10 −6 to 10 −3.5 , higher multiplicity systems were more likely to contain one or more Jovian-mass planets.Any Jovian-mass planet in a system with similar period ratios will decrease the ∆ value between it and any neighboring planets, so each Jovianmass planet likely decreases two ∆ values in a system.This makes the chance of every planet pair having a high ∆ value much less likely for high planet multiplicity than for low planet multiplicity.
Even with the difference in ∆ min corresponding to multiplicity for the stable planetary systems, the upper bound of ∆ min for the unstable planetary systems in both subsets of data was relatively evenly spaced in multiplicity.The 95th-percentile value for the unstable cases in the Nominal Kepler Analog Set was ∼5, and this stayed relatively constant across all multiplicities.Similarly, the 95th-percentile value for the unstable systems in the Very Closely Packed Set was ∼7 for each planet multiplicity.
For the Nominal Kepler Analog Set, there was an increase in the fraction of unstable systems with increasing planet multiplicity.The unstable fraction for 3-planet systems was ∼4% but for 9-and 10-planet systems was ∼15%.This decrease in stability is consistent with the expectation that higher-multiplicity systems with close spacings are more likely to go unstable at an earlier time (e.g., Pu & Wu 2015;Petit et al. 2020).However, for the Very Closely Packed Set, the rejection fraction seems at most weakly dependent on the number of planets in the system, as 3-and 10-planet systems had roughly the same unstable fractions that were the lowest amongst the different planet multiplicities.

Mass ratio dependence
We also wanted to determine if there were any correlations between the system stability and the mass ratios within the system.The random initialization of the systems with initial masses between sub-Earths and super-Jupiters allowed for a large absolute mass ratio between nearest-neighbor planets.Positive values of log m2 m1 mean that the outer planet is more massive and negative values mean that the inner planet is more massive.This is intended to cover a wide variety of planetary systems.For both our samples of simulated planetary systems, we found a pattern between system stability and the mass ratios in the system.Systems with equal-mass planets are slightly less likely to be stable than those with unequal planetary mass ratios; however, as the mass ratios become more extreme (below ∼ 0.01 or above ∼ 100), the systems tend to become less stable (see the bottom right panels of Figures 4 and 6).
This pattern in stable fraction with planet mass ratio may have its origin in dynamical evolutionary pathways.Goldberg & Batygin (2022) found that interactions in systems with planets of almost-equal mass tend to decrease the uniformity of the planets in the system by either removing planets from the system or in some cases merging together in a collision, increasing the mass inequality of initially similar-mass planets.In addition, Lammers et al. (2023) showed that a level of intra-system mass uniformity, but specifically not rising to mass equality, is a likely result of dynamical interactions causing ejections and collisions during the planetary systems' evolution.
In general, the strength of the destabilizing effect scales proportionally with the mass of the perturbing bodies (in this case the neighboring planet pair; see e.g., Volk & Malhotra 2020).For individual pairs of similarmass planets, the destabilizing effect of a close interaction between the two planets is lessened if they have a slight but significant mutual inclination (e.g., Rice et al. 2018).In addition, unequal mass planets (especially with unequal spacings) tend to reach a larger part of the parameter space and to have a longer stability time than equal mass planets on equal spacings (Pu & Wu 2015).

SUMMARY
In this work, we carried out numerical simulations to study the distribution of the orbital separations in dy-namically stable planetary systems.The critical parameter of interest is the dynamical separation (Eq.1), that is, the orbital separation in units of the mutual Hill radius of nearest-neighbor planets.We studied the overall distribution of ∆ as well as that of its minimum value, ∆ min , across nearest-neighbor pairs in multi-planet systems.We also studied how planet multiplicity as well as orbital period ratios and mass ratios affect the minimum value of this parameter.We generated two sets of simulated planetary systems: the Nominal Kepler Analog Set following the observed Kepler period distribution, and the Very Closely Packed Set containing more planet pairs with short separations.The key findings of our study are as follows: • The ∆ min distribution is not well fit by a single threshold value, a common assumption in previous studies.Instead, a log-normal function with parameters µ ≈ 1.97±0.02and σ ≈ 0.40±0.02best fits the distribution of ∆ min .This corresponds to a median value of ∆ min of 7.17 ± 0.15, a rootmean-square value ⟨∆ 2 min ⟩ 1 2 = 8.4 ± 0.3, and the range (4.8, 10.7) contains two-thirds of ∆ min values.The Nominal Kepler Analog Set and the Very Closely Spaced Set differ insignificantly from each other in the ∆ min distributions of their long-term stable systems.
• The overall distribution of ∆ in long-term stable systems is well fit with a log-normal function with best-fit parameters µ = 3.14 ± 0.03 and σ = 0.76 ± 0.02 (the Nominal Kepler Analog Set), and µ = 3.08 ± 0.04 and σ = 0.76 ± 0.02 (the Very Closely Packed Set).
• The distribution of planetary spacings of the Kepler sample is consistent with sculpting by dynamical stability.Specifically, there is little difference between the overall ∆ distributions found in the final stable architectures of planetary systems of the Nominal Kepler Analog Set and in the Very Closely Packed Set.
• There is a significant trend that long-term stable systems with high planet multiplicity have smaller initial minimum ∆ values, which is attributed to the higher probability of a Jovian-mass planet in the system given our random planetary mass draws.
• Planetary systems with large numbers of planets are in general less likely to be stable, but for the Very Closely Packed Set the long-term stable fraction is non-monotonic with the planet multiplicity, increasing again at the highest planet multiplicities (N > 8).
• Systems with nearest-neighbor planet pairs having similar mass ratios tend to be less stable than those with somewhat dissimilar mass ratios, but the stable fraction tends to decrease again for mass ratios ≳ 10 2 .
Our results on the statistical distribution function of the critical minimum value ∆ min of nearest-neighbor planet pairs in long-term stable planetary systems provide more precise and accurate parameters for planetary system evolution models, as well as improvements in statistically predictive tools, such as Dynamite (Dietrich & Apai 2020;Dietrich et al. 2022).Our results on the overall distribution of planetary spacings support the hypothesis that the ∆ distribution of the thousands of exo-planetary systems discovered by the Kepler space telescope likely evolved from more tightly packed systems that underwent orbital evolution through dynamical instabilities (including ejections and collisions) to reach their current states.This suggests that dynamical interactions amongst planets decisively shape the architectures of compact planetary systems.
The results reported herein benefited from collaborations and/or information exchange within the program "Alien Earths" (supported by the National Aeronautics and Space Administration under Agreement No. 80NSSC21K0593) for NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate.RM additionally acknowledges funding from NASA grant 80NSSC18K0397.We acknowledge use of the software packages NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), and REBOUND (Rein & Liu 2012;Rein et al. 2019).An allocation of computer time from the UA Research Computing High Performance Computing (HPC) at the University of Arizona is gratefully acknowledged.The citations in this paper have made use of NASA's Astrophysics Data System Bibliographic Services.

Figure 1 .
Figure 1.The distribution of initial ∆ (nearest-neighbor orbital separations in units of mutual Hill radius) of the Nominal Kepler Analog Set of simulated planetary systems.The random draws we performed generated a roughly log-normal distribution of ∆ in a system, as shown by the overplotted fit.

Figure 2 .
Figure 2. The distribution of initial ∆ (nearest-neighbor orbital separations in units of mutual Hill radius) of the Very Closely Packed Set of simulated planetary systems.Note that this distribution is roughly uniform in the range ∼3.5 -10 (compare with Figure 1) .

Figure 3 .
Figure 3. Results for the Nominal Kepler Analog Set.Top: Histogram of all the ∆ values for each nearest-neighbor planet pair in stable systems (left) and unstable systems (right).Bottom: Histogram of the initial ∆min distribution and the best-fit log-normal function for stable systems (left) and unstable systems (right).
The threshold of stability Stability in the Nominal Kepler Analog Set

Figure 4 .
Figure 4.More results for the Nominal Kepler Analog Set.Top left: Scatter plot and histograms of all initial ∆ and final ∆ for stable systems (in blue) and unstable systems (in orange).Top right: Scatter plot and histograms of the planet multiplicity versus the initial ∆min distribution stable systems (in blue) and unstable systems (in orange).The green triangles indicate the median of the initial ∆min for stable systems.Bottom left: the fraction of stable systems versus planet multiplicity.Bottom right: the fraction of stable systems versus the mass ratio between the closest nearest-neighbor planet pair.

Figure 5 .
Figure 5. Results for the Very Closely Packed Set.Top: Histogram of all the ∆ values for each nearest-neighbor planet pair in stable systems (left) and unstable systems (right).Bottom: Histogram of the initial ∆min distribution and the best-fit log-normal function for stable systems (left) and unstable systems (right).

Figure 6 .
Figure 6.More results for the Very Closely Packed Set.Top left: Scatter plot and histograms of all initial ∆ and final ∆ for stable systems (in blue) and unstable systems (in orange).Top right: Scatter plot and histograms of the planet multiplicity versus the initial ∆min distribution stable systems (in blue) and unstable systems (in orange).The green triangles indicate the median of the initial ∆min for stable systems.Bottom left: the fraction of stable systems versus planet multiplicity.Bottom right: the fraction of stable systems versus the mass ratio between the closest nearest-neighbor planet pair.