Carving Out the Inner Edge of the Period Ratio Distribution through Giant Impacts

We consider five characteristics of exoplanet systems when constructing our simulations: mass, eccentricity, period ratio, correlation of period ratios, and number of exoplanets in the system.

Perhaps most importantly, dynamical instabilities are qualitatively different for isolated pairs of planets vs. systems with three or more planets. In two-planet systems, there exists a critical "Hill-limit" orbital separation beyond which the pair is guaranteed to never undergo close encounters (C. Marchal & G. Bozis 1982; B. Gladman 1993). This boundary separates pairs that are too close and typically destabilize within  ∼ 100 orbits, from pairs that are stable beyond systems' Gyr lifetimes. By contrast, systems with three or more planets exhibit instabilities for significantly wider interplanetary separations, and the time to such instabilities spans the full dynamic range from orbital timescales to many billions of orbits (e.g., J. Chambers et al. 1996; A. W. Smith & J. J. Lissauer 2009; A. Obertas et al. 2017). This is thought to be due to the effects of three-body mean motion resonances (MMRs; A. C. Quillen 2011; A. C. Petit et al. 2020; J. Rath et al. 2022; C. Lammers et al. 2024), as well as long-term oscillations in the eccentricities that cause MMRs to adiabatically expand and contract, thereby sweeping out regions of chaos (D. Tamayo et al. 2021; Q. Yang & D. Tamayo 2024). Given that, on average, observed systems have three planets with orbital periods within 400 days (W. Zhu et al. 2018) and such systems are the lowest multiplicity to exhibit the general dynamical behavior of compact multi-planet systems, we choose to generate synthetic populations of three-planet systems.

We draw our synthetic populations from the set of observed systems taken from the NASA exoplanet archive on 2024 August 4 (NASA Exoplanet Science Institute 2020). For simplicity, we begin by filtering for all planets in 3+ planet systems that have an inner period ratio (P₂/P₁) below 1.5.³ We then calculate planet–star mass ratios and remove any planets where these are undefined (due to either a missing planet, or stellar mass), or giant planets with mass ratios >10⁻⁴ ≈ 2 Neptune masses. After filtering, we are left with 73 compact pairs.

Given that Newtonian gravity is scale invariant, we rescale all our synthetic systems to have central stars of 1M_⊙ and an inner planet with an orbital period of 1 yr. We then draw the planet–star mass ratios and orbital eccentricity for each planet randomly and independently from the observed distributions described above (with replacement).

Graphs plotting both the observed and synthetic distributions for the eccentricities and mass ratios can be found in Figure 2.

Zoom In Zoom Out Reset image size

We randomly assign the longitude of pericenter for each orbit uniformly between [0, 2π). For simplicity, we also make all systems coplanar since typical mutual inclinations of transiting exoplanets of less than a few degrees were found by D. Tamayo et al. (2021) to have a negligible effect on dynamical stability on timescales of less than a billion orbits.

We choose the period ratio for the inner planet pair randomly between our range of interest of 1.1–1.5. For the outer pair, we try to capture the observed preference toward uniform orbital spacings (L. M. Weiss et al. 2018), since this also should affect dynamical stability (A. C. Petit et al. 2020; J. Rath et al. 2022). We do this by randomly drawing from a Gaussian distribution centered at the inner period ratio (corresponding to uniform spacing) with a standard deviation that we empirically determine in order to reproduce the observed spread. We quantify this deviation from uniform spacings by calculating a normalized dispersion D_i for each system

$$\begin{eqnarray}&&\begin{array}{r}{D}_{i}=\displaystyle \frac{{\sigma }_{i}}{{\mu }_{i}},\end{array}\end{eqnarray} \tag{ 1 }$$

where i is an index labeling the multi-planet system, and μ_i and σ_i are the corresponding mean and standard deviation of the ${\mathrm{log}}_{10}$ of the period ratios between adjacent planets in that system. From the NASA Exoplanet Archive, we calculated a mean dispersion for our observed sample (filtered as described above) of D = 0.30.

Traditionally, one would simulate dynamical evolution via direct numerical integrations; however, this method is computationally expensive. Utilizing the WHFast integrator (H. Rein & D. Tamayo 2015) in the REBOUND N-body package (H. Rein & S.-F. Liu 2012), a 10⁹-orbit integration takes  ≈ 6 CPU hr using a 2.1 GHz Intel Xeon Silver 4116. For the numerical experiments we perform in this paper, we estimate N-body methods would require almost  ≈ 3 million CPU hr.

We therefore instead use the Stability of Planetary Orbital Configurations Klassifier (SPOCK), an open-access machine learning model that predicts the stability of compact systems of three or more planets (D. Tamayo et al. 2016, 2020). SPOCK runs short 10⁴ orbit integrations, and calculates dynamically informed features from which the gradient-boosted decision-tree model is trained to predict stability over 10⁹ orbits. SPOCK is roughly 10⁵ times faster than direct integration, and outputs a probability of survival between 0 and 1. This reduces the computational cost to ≈50 CPU hr.

For each experiment, we follow the procedure outlined in Section 2.1 to repeatedly generate synthetic three-planet systems. We build up our simulated population by retaining each sampled system with a probability given by its corresponding likelihood of stability as estimated by SPOCK. We stop once we reach a population that contains the same number of compact pairs as the observations (73). This allows us to naturally weight systems with higher probabilities of stability more heavily (D. Tamayo et al. 2021).

We now assess the sensitivity of our results to our various population parameters, which are both uncertain and potentially observationally biased.

We first perform a sequence of experiments, where we reproduce 100 synthetic trials like in Figure 3, except after drawing the mass ratio for each planet, we multiply all the masses by a common scaling factor. We vary the scaling factor logarithmically from 0.5 to 2 in the left panel of Figure 5, where each blue point corresponds to the mean (and error bars to the standard deviation) of the p-values across the 100 trials of that experiment.

Zoom In Zoom Out Reset image size

We see that we retain qualitatively similar results over the fairly broad range of parameters probed, showing that the result is robust against these choices. As a simple test of how the biases inherent to different observational methods affect our results, we also considered only the subset of planet pairs with masses measured through radial velocities 71, dropping the remainder (largely derived from mass–radius relations). We find this yields an even larger p-value of 0.66.

The middle panel does the equivalent test, scaling the orbital eccentricities. We see that it does not significantly affect the p-values, mostly because a large fraction of our observed sample has circular orbits. Of course, eccentricities of zero in the Exoplanet Archive do not necessarily imply circular orbits, but typically reflect that the eccentricities were too small to measure observationally. As an opposite limiting case, if we instead use our nominal scaling factor of 1, but remove all zero-eccentricity planets, we obtain a p-value of 0.25. We interpret this as a lower bound, given that this sample should instead be biased toward high eccentricities (since high eccentricities are easier to measure). Indeed, we find that the average non-zero-eccentricity in this sample is 0.067, somewhat higher than population-level constraints from transit durations in multi-planet systems that fit a Rayleigh distribution with a Rayleigh parameter of σ_e = 0.049 (V. Van Eylen & S. Albrecht 2015).

We similarly find that varying the dispersion multiplier to vary the uniformity of our period ratios does not significantly affect our results.