Orbital Architectures of Kepler Multis From Dynamical Instabilities

The high-multiplicity exoplanet systems are generally more tightly packed when compared to the solar system. Such compact multi-planet systems are often susceptible to dynamical instability. We investigate the impact of dynamical instability on the final orbital architectures of multi-planet systems using N-body simulations. Our models initially consist of six to ten planets placed randomly according to a power-law distribution of mutual Hill separations. We find that almost all of our model planetary systems go through at least one phase of dynamical instability, losing at least one planet. The orbital architecture, including the distributions of mutual Hill separations, planetary masses, orbital periods, and period ratios, of the transit-detectable model planetary systems closely resemble those for the multi-planet systems detected by Kepler. We find that without any formation-dependent input, a dynamically active past can naturally reproduce important observed trends including multiplicity-dependent eccentricity distribution, smaller eccentricities for larger planets, and intra-system uniformity. On the other hand, our transit-detectable planet populations lack the observed sub-population of eccentric single-transiting planets, pointing towards the Kepler dichotomy. These findings indicate that dynamical instabilities may have played a vital role in the final assembly of sub-Jovian planets.


INTRODUCTION
Numerous successful space-based transit surveys such as NASA's Kepler (Borucki et al. 2010;Borucki 2016), K2 (Howell et al. 2014;Cleve et al. 2016), and TESS (Ricker et al. 2015), numerous radial velocity surveys such as HIRES (Vogt et al. 1994), HARPS (Mayor et al. 2003), HARPS-N (Cosentino et al. 2012), SOPHIE (Bouchy & Sophie Team 2006), and ESPRESSO (Pepe et al. 2021) and follow-up missions like LAMOST-Kepler survey (Dong et al. 2014;Luo et al. 2015) and California-Kepler Survey (Petigura et al. 2017), have led to the discovery and characterization of more than 5000 exoplanets (NASA Exoplanet Archive, 1  Akeson et al. 2013) in more than 3500 planetary systems with diverse orbital architectures.Among these, NASA's Kepler mission stands out with over 2700 confirmed planet discoveries, more than half of all planets discovered to date.The majority of these planets are smaller than Neptune and have relatively short orbital periods ( orb /yr ≲ 1).Apart from the sheer number of exoplanet discoveries, the Kepler mission also found an abundance of multi-planet systems (multis in short).These multi-transiting systems not only inform us about the orbital architecture of planetary systems at present, but also provide key constraints on their formation, assembly, and dynamical evolution and they are also crucial for inferring the intrinsic exoplanet populations not directly observable due to selection biases (e.g., Lissauer et al. 2011;Fang & Margot 2012;Fabrycky et al. 2014;Pu & Wu 2015;Weiss et al. 2018;He et al. 2019He et al. , 2020)).
Multi-planet systems are often tightly packed and susceptible to dynamical instabilities.Indeed, several studies have suggested that many observed multi-planet systems are very close to their stability limit (e.g., Deck et al. 2012;Fang & Margot 2013;Lissauer et al. 2013;Volk & Gladman 2015;Hwang et al. 2017;Volk & Malhotra 2020).The timescale for the onset of dynamical instabilities in a planetary system is dependent on the spacing between the planets in multiples of the so-called mutual Hill radius where   ,  p  are the semi-major axis and mass of the j-th planet and  * is the mass of the host star.The stability of two planet systems can be analytically determined from this parameter (e.g., Gladman 1993;Deck et al. 2013).However, for systems with higher multiplicities, we rely on numerical simulations to estimate the stability timescale normalised by the innermost planet's  orb , .Extensive numerical studies have revealed an empirical relation between  and , where  and  are constants that depend on the system's multiplicity and planetary masses (e.g., Chambers et al. 1996;Funk et al. 2010).Usually, to limit the parameter space, numerical investigations of  often adopt equal mass planets separated at equal .Even with these simplifications,  shows substantial statistical variations, sometimes with ranges spanning orders of magnitude, and deviates from this relation near mean motion resonances (e.g., Chatterjee et al. 2008;Funk et al. 2010).How  behaves for systems hosting planets of different masses and  is not yet properly mapped.Qualitatively though, it is well understood that over time, the low- orbits interact and  increases through collisions, ejections, and scattering.The systems, as a result, evolve towards longer .Thus, the -distribution of a population of multis essentially indicates a statistical measure of the dynamical state of the population.
Getting the spacing distribution for the observed multis is straightforward (Equation 1) if we know the semi-major axes (SMA) and planet masses ( p ).However,  p measurements are unavailable for a large fraction of Kepler planets.If we circumvent this challenge by using a mass-radius relationship (e.g., Chen & Kipping 2016) to estimate  p from the measured radii ( p ), the -distribution for Kepler multis peaks around  ∼ 12 − 15, expected for a stability timescale of ≳billion years (Pu & Wu 2015).There is a sharp drop in the number for  ≲ 12 and a power-law-like drop-off for higher .Based on this, it has been suggested that the present-day  distribution may be a result of past dynamical sculpting (e.g., Pu & Wu 2015;Volk & Gladman 2015).Past dynamical sculpting not only affects the distribution of  but also leaves tell-tale signatures in the orbital properties including eccentricities () and mutual inclinations ( m ) (e.g., Rasio & Ford 1996;Chatterjee et al. 2008;Jurić & Tremaine 2008;Nagasawa & Ida 2011).Several studies investigated the effects of planetary dynamics in planetary systems with a variety of initial conditions, usually motivated by a preferred formation scenario.For example, Izidoro et al. (2017Izidoro et al. ( , 2021)); Goldberg & Batygin (2022) studied the effects of dynamical instability on systems of super-earths, initially in compact resonant chains.Wu et al. (2019) investigated the influence of dynamical instability on an initially flat period ratio distribution.
In contrast, in this study we are intentionally agnostic towards any particular formation scenario and make simple assumptions about the initial conditions of planetary systems as they start to evolve freely without any gas disk.We want to investigate, if we do not provide any input from formation scenarios or inject any intra-system trends, how close we can get to the observed properties and correlations for the Kepler multis simply as a result of dynamical instabilities.Since the majority of the Kepler planets are smaller and the larger giant planets are expected to form differently, we focus on planetary systems containing sub-Jovian ( p ≲ 6 ⊕ ) planets only in this paper.In section 2, we describe the numerical setup and the adopted initial conditions.In section 3, we present the outcomes from our suite of simulations and compare our models with Kepler's observations after correcting for transit geometry and Kepler's detection efficiency.Finally, we outline the key results and conclude in section 4.

NUMERICAL SETUP
We construct -body models of planetary systems with planets randomly assigned according to simple assumptions on the initial conditions without any input on intra-system features or correlations.Our fiducial ensemble initially consists of 8 planets in each system following Kepler-90, the system with the highest number of confirmed exoplanets (Shallue & Vanderburg 2018).However, we create other ensembles by varying the initial number of planets ( p,ini ) to investigate the effects of  p,ini on our results.Furthermore, we explore the effects of initial  and inclination () distributions on our results.Below we describe the details of how we construct our initial planetary systems in different ensembles.

Star and Planet Properties
The central star in each model planetary system is chosen randomly from the planet-hosting stars discovered by Kepler.We draw  p randomly from a power law distribution  ( p ) ∝  p   with   = −1.08 based on the estimated intrinsic  p distribution (He et al. 2019) in the range 0.5 ≤  p / ⊕ ≤ 3. The occurrence rate of sub-jovian planets drops significantly for  p / ⊕ ≳ 3 (e.g., Petigura et al. 2013).On the other hand, the detection probability of planets smaller than  p / ⊕ ≲ 0.5) is low.We estimate  p from  p using the probabilistic relationship given in Chen & Kipping (2016).We randomly draw  p within 1 of the prescribed normal distribution for a given  p .These choices ensure that all planets in our ensembles are initially sub-Neptunes and have  p / ⊕ ≲ 16.

Interplanetary separations & initial planetary orbits
Depending on the focus, past studies adopted various schemes for spacing the initial planets, such as a Gaussian distribution centering a fixed value (Pu & Wu 2015) and a flat distribution in period ratios (Wu et al. 2019).In contrast, we space our planets using a  distribution motivated by that for the Kepler multis and the following considerations.
The observed  distribution shows a sharp decline for  ≲ 12 and a power-law like decline for  ≳ 15.Pu & Wu (2015) suggested that the latter is genuine and not a result of observational biases.Planetary orbits with low separations ( ≲ 12) are expected to become unstable within the typical age of the observed host stars and hence, the sharp decline in low  may be due to past dynamical instabilities.On the other hand, orbits separated by high  ≳ 15 likely remained stable.Thus, the high end of the present-day  distribution is expected to be more representative of the initial.Hence, we adopt a power-law distribution for the initial  of the form  () ∝    in the range  = 4 to 40.The lower limit is slightly higher than the analytic limit for the stability of two-planet systems (Gladman 1993).The upper limit comes from practical considerations.While the tail can extend to very high , we truncate at  = 40 to limit the physical size of our initial model systems.Planets in such high-SMA orbits are usually significantly beyond the detection capability of the Kepler mission, the focus of our study.Moreover, the fraction of observed planet pairs with  > 40 is low (< 14%).
Determining   , is not straightforward.One naive choice could be to fit the observed  distribution for large- to estimate   .However, there is a significant caveat.A planetary system with low initial  pairs may dynamically evolve to create high- pairs and contribute to the present-day  ≳ 15.Hence, simply finding the bestfit power-law from the high end of the present-day  distribution may not be justified.We want to find out the initial  distribution that would result in an observed present-day distribution after dynamical evolution considering detection biases.Hence, we run several test simulations with different   .Based on these tests we adopt   = −1.1 throughout this study.
For each planetary system, we assign the SMA of the innermost planet ( 1 ) by randomly choosing (with replacement) from the observed  1 in Kepler multis.We place the remaining  p,ini − 1 planets in orbit one by one from inside out by randomly drawing  rand from the adopted  distribution.The SMA of the (  + 1)th planet is then given by where, We repeat this process until the SMA of the outermost planet is calculated.Note that our initial assembly does not specifically put any planet pair in or out of resonance.In effect, without any dissipation in the system, they are all initially non-resonant.For each orbit, orbital  and  are drawn randomly from Rayleigh distributions with adopted mean values ē, ī.We create several ensembles by varying  p,ini , ē, and ī (summarised in Table 1).Other orbital elements, such as the longitude of ascending node, the argument of pericenter, and true anomaly, are chosen uniformly in their respective full ranges.

Evolving the Planetary Orbits
We use the mercurius hybrid integration scheme (Rein et al. 2019) implemented in the rebound simulation package (Rein & Liu 2012) to evolve the planetary systems with initial time step   = 1/25 of the innermost planet's  orb .We resolve close encounters by switching from the WHFast to the IAS15 integrator when the planet-planet or planet-star separations are 3 times the sum of their Hill radii.
The typical energy error in our simulations is / ∼ 10 −6 , after correcting for energy losses during collisions and ejections.For each ensemble, we simulate 360 planetary systems.We find that energy was not conserved well (/ > 10 −3 ) in ≲ 2% of these simulations, so we discard them from our analysis.Table 1 lists the number of systems that satisfy the minimum accuracy for energy conservation and are considered for further analysis.We treat collisions adopting the "sticky-sphere" approximation-whenever two planets physically touch each other, we merge them conserving mass, linear momentum, and physical volume.We remove planets if the SMA is > 100 times the initial outermost planet's initial  to reduce computational cost and consider it an ejection (happened in case of ≲ 0.7% simulated planets).

Transit Probability and Kepler's Detection Bias
The -body simulations provide us with an intrinsic set of planetary systems sculpted through dynamical instability.To compare these planetary systems with those observed, we account for the geometric transit probability and Kepler's detection biases.We use the CORBITS algorithm (Brakensiek & Ragozzine 2016) to identify the planets that would transit when viewed from a random line of sight (LOS).For each planetary system, we collect 1000 LOSs for which at least one transiting planet is found.2Even if a planet transits, detection is not guaranteed due to sensitivity biases inherent in Kepler's pipeline.We adopt a global detection efficiency model for Kepler (Burke et al. 2015;Christiansen et al. 2015Christiansen et al. , 2016;;Christiansen 2017) to estimate the detection probability ( det ) in the Kepler pipeline based on the signal to noise ratio (SNR) as a proxy to the so-called multi-event statistic (MES).We estimate SNR based on the combined differential photometric precision (CDPP 4.5hr where  * is the stellar radius,  obs is the time span of Kepler's observation,  is the fraction of  obs that contain valid data, and  dur is the transit duration of the planet.We estimate  dur using equation 1 of Burke et al. (2015).
We estimate  det in the Kepler pipeline using the gamma cumulative probability distribution given in Christiansen (2017).We multiply  det with an analytic window function (Burke et al. 2015, their Equation 6) which estimates the probability of detecting at least three transits, the minimum number adopted by Kepler to be considered detected, accounting for the limits of data coverage.We consider a transiting planet to be detected by Kepler if a random number from a uniform distribution U [0, 1] is less than the planet's  det .Finally, from the intrinsic population, we find a population that would be detected by Kepler and call this population 'model-detected' henceforth.Since we explicitly use Kepler's detection pipeline, we restrict ourselves to confirmed multis detected by Kepler for a fair comparison.

Integration Stopping Criteria
It is usually computationally impractical to continue -body integration till the physical age of a typical system, neither can  > 2 systems be proven stable.Traditionally, the adopted stopping time in -body simulations are somewhat ad hoc and guided by practical considerations.Any dynamically active system becomes more and more stable through repeated dynamical encounters by decreasing  p or by increasing  or both.Hence, the integration is usually stopped after some predefined number of orbits when scattering encounters become sufficiently rare.The challenge is that since the dynamically active system, over time moves towards increased stability and the stability timescale increases exponentially with increasing , it is not clear what can be considered as sufficiently rare.
We adopt a more sophisticated and physically meaningful approach to determine the integration stopping time.This is above and beyond the usually adopted criteria of 'sufficiently' rare subsequent interactions (see section A for more details).Since our testing hypothesis is that exoplanet systems formed more dynamically packed and we observe them after long dynamical sculpting, we want to bring our model population to a similar dynamical state as Kepler's multis.Since for any system, dynamical stability is determined by , the observed  distribution should represent the present-day dynamical state of the observed population.3However, the observed  distribution is not intrinsic and subject to geometric transit probability and pipeline detection efficiency.Hence, we must compare the synthetic  distribution after taking into account these observational biases.We achieve this in the following way.
We simulate all ensembles to at least 10 8  1 and examine the model-detected  distribution.If the model-detected  distribution peaks at a higher value than the observed, the ensemble is too dynamically evolved compared to the observed population.If so, we look back in time to find a snapshot where these two distributions are statistically identical.Similarly, if the model-detected  distribution peaks at a lower value, the ensemble has not yet reached the dynamical state of the observed population.We integrate these ensembles further until agreement is achieved or we reach our hard stop at 10 9  1 (to limit computational cost).The time taken to achieve the desired population-level dynamical state  stop depends on the distributions of initial orbital properties of the ensemble (Table 1).Only ensemble n6-e040-i024, with low initial  p,ini does not quite reach the observed dynamical state at 10 9  1 .

RESULTS
We first discuss the results from our fiducial ensemble (n8-e040-i024) and investigate its similarity with Kepler's multis.Afterward, we discuss the effects of different initial conditions on the final outcomes of our simulations.

Fiducial Ensemble
We find that about 97% of our synthetic planetary systems are affected by dynamical instabilities that result in the loss of at least one planet via collision with another planet or the host star or ejection from the system.Since our ensemble is dominated by close-in planetary systems with Safronov number Θ ≪ 1, strong scattering ultimately results in collisions rather than ejections in most cases.4

Orbital and Planetary Properties
The final systems emerging from this chaotic dynamical evolution have little memory of their initial conditions.Hence, the orbital architectures of these systems are very different from their initial configuration.By design, our models initially have more planets at closer separations, and systems with such small interplanetary spacing inevitably face dynamical instabilities.Dynamical instability leads to close encounters between planets, increasing their dynamical excitation.The degree of dynamical excitation of a planetary system can be quantified by the normalised angular momentum deficit (Chambers 2001), We show the evolution of NAMD for three different representative systems in Figure 1.In general, dynamical instabilities raise the NAMD of the systems while physical collisions decrease NAMD (also see, Chambers 2001;Laskar & Petit 2017).Competition between these two processes shapes the final interplanetary spacing and orbital properties in a system (Dawson et al. 2016).On the other hand, NAMD remains conserved for systems that do not experience significant instability.
Figure 2 shows the distributions for  from our fiducial ensemble.Clearly,  ≲ 12 is heavily sculpted by dynamical instabilities.The final model-detected  distribution shows excellent agreement with that for the observed.A two-sample Kolmogorov-Smirnov (KS) test suggests that the null hypothesis that the observed and modeldetected  are drawn from the same underlying distribution can not be rejected with  ≳ 0.3.This indicates that this ensemble of model planetary systems has achieved a similar dynamical state as Kepler's observed systems.
Although, the initial distributions for  orb ratios and  p are significantly different from those observed, once a similar dynamical state is achieved, the key observable properties including the distribution of adjacent-planet  orb ratios,  orb , and  p show excellent agreement between the model-detected and observed populations (Figure 3).While the  orb distribution remains relatively unchanged via dynamical sculpting, it is severely affected by Kepler's detection biases, making the  orb distribution for the model-detected population resemble the Kepler data.A closer inspection shows that our model-detected  orb distribution peaks at a marginally lower value compared to the observed (middle panel, Figure 3).This can be attributed to our initial conditions.
The simultaneous agreement in  and  orb ratios suggests that  p in our models are consistent with the observed.This is also evident when we convert  p into  p and compare the model-detected population with the observed (right panel, Figure 3).The simulated systems initially only contained planets with  p / ⊕ < 3. Thus, all planets with  p / ⊕ ≳ 3 in our models are created by planet-planet collisions.Interestingly the model-detected  p distribution shows excellent agreement with the Kepler data even for  p / ⊕ ≳ 3.This indicates the possibility that in nature, smaller planets form more abundantly, and most larger planets are formed via mergers of the smaller planets during dynamical instabilities.On the other hand, we find a turnover in the final intrinsic PDF at  p / ⊕ ≈ 1, about 2× higher than the smallest planet we have considered initially (subsection 2.1).The turnover in the model-detected populations is at  p / ⊕ ≳ 1.5, which illustrates Kepler's bias towards detecting larger planets.Since the turnovers in the final  p distributions are at much higher values compared to the lower cutoff we used in our initial power-law  p distribution, the effects of this initial choice are expected to be small.
It is also worth noting that in the sub-Neptune regime,  p is strongly dependent on the atmosphere properties, and the ability for late-stage atmosphere retention by the planet can significantly affect the observed  p distribution.Our idealised simulations do not incorporate additional physics for atmosphere loss, e.g., via collisions loss (e.g., Hwang et al. 2017;Denman et al. 2020), photoevaporation (e.g., Owen & Wu 2017), or core-powered evaporation (e.g., Ginzburg et al. 2018).Hence, a comparison of the  p distributions with more granular details (e.g., the bimodal distribution of  p , Fulton et al. 2017) is beyond the scope of this work.
We show the final ,  m , and  p in Figure 4.The final  distribution of our model-detected population approximately follows a Rayleigh distribution with ē ≈ 0.042 which agrees well with the population level observed  estimates for Kepler's multis (Xie et al. 2016;Mills et al. 2019).The model-detected single-planet systems exhibit a marginally higher ē ≈ 0.048 than the model-detected multis.However, the ē of the apparently single-planet systems are considerably lower compared to the estimated ē ∼ 0.2 − 0.3 for Kepler's singles (e.g., Xie et al. 2016;Mills et al. 2019;Van Eylen et al. 2019) (more discussion on this later).
In contrast to giant planet scattering (e.g., Chatterjee et al. 2008),  m does not increase significantly for most systems we have studied here.In our fiducial ensemble,  m increases from a Rayleigh distribution with ī m = 0.024 initially to ī m = 0.026 for the final intrinsic population.Nevertheless, the final intrinsic  m distribution exhibits a prominent tail extending to large  m (Figure 4).
The intrinsic final population shows a positive correlation between  and  m with a Pearson correlation coefficient, C = 0.56, consistent with the estimate for the Kepler data from the 'maximum AMD' model of He et al. (2020).We also find that  p is anti-correlated with  with C = −0.18(bottom panel, Figure 4).This is expected because the same level of AMD would alter the orbit of a lower-mass planet more compared to a higher-mass planet.Another subtle effect may also contribute to this anti-correlation.In our setup, large planets (e.g.,  p / ⊕ > 3, equivalently  p / ⊕ ≳ 16) are collision products.Since collisions reduce AMD (Figure 1), the orbits of these planets are expected to be less eccentric.
Figure 5 shows the final multiplicity distributions of our simulated systems.Each planetary system in the fiducial ensemble has  p,ini = 8.Dynamical instabilities reduce  p predominantly due to collisions.Only ∼ 3% systems in this ensemble survive with all 8 planets.The final intrinsic multiplicity ranges from 2 to 8.However, in most cases, only one planet transits.The detected number of systems declines steeply with increasing  p .Although this trend is similar to the Kepler data, we find a lack of model-detected singleplanet systems compared to what is observed.A similar discrepancy was noted in several past studies (e.g., Lissauer et al. 2011;Johansen et al. 2012;Hansen & Murray 2013;Ballard & Johnson 2016) while attempting to explain the observed multiplicity distribution under various assumptions and is most commonly referred to as the "Kepler dichotomy".The most common explanation of this apparent excess of single transiting systems is to consider a second population of either intrinsic single-planet systems or multis with higher  m (e.g., Fang & Margot 2012;Mulders et al. 2018;He et al. 2019).Indeed, we find that the distribution of detected singles as well as the higher multiplicities can be explained simply by artificially injecting additional single transiting planets in the mix, such that the fraction of injected population,  ∼ 0.48.This is only slightly higher than the estimated fraction of observed systems with high  m (  ∼ 0.4; Mulders et al. 2018;He et al. 2019).Moreover, since our single transiting planets have lower  than the Kepler's singles, our results support the hypothesis that an additional population of dynamically hotter (or intrinsically single) planets may be needed to explain the apparent excess of observed singles.This may indicate a different dynamical origin for the dynamically hotter single-planet population or a population with higher-mass planets as perturbers (which were not included in this study).On the other hand, it has also been suggested that a carefully curated single-population model can also explain the Kepler data (e.g., Zhu et al. 2018;Sandford et al. 2019;He et al. 2020;Millholland et al. 2021).
Interestingly, the low-multiplicity systems in our models are forged from systems with  p,ini = 8 predominantly via collisions.Thus, in the intrinsic population we find an anti-correlation between  p and multiplicity.However, this trend is significantly washed-out in our model-detected population, which is in excellent agreement with the estimated  p for the observed planets (Figure 6).This trend in  p with intrinsic multiplicities may provide a test for the dynamical models if these quantities are better understood for the observed planetary systems in the future.Of course, in reality, differences in the host mass and metallicity correlations (which we intentionally did not include) may make comparisons non-trivial.

Peas In a Pod
One of the most interesting observed trends in Kepler's multis is that the planets within a particular system are more similar in size and regularly spaced compared to what would have been expected if they were drawn at random from all observed planets.This is commonly referred to as 'peas in a pod' (Weiss et al. 2018).Several studies argued that this trend likely was set early during the planet formation process (Adams 2019;Adams et al. 2020;Batygin & Morbidelli 2023).As mentioned already, we do not introduce any intra-system uniformity in spacing or in planet properties.In fact, by design, the planets in our planetary systems are drawn completely randomly from the full set of possible planets and spacings.Hence, our simulations provide a clean test for the possibility of emergence of intra-system uniformity as a result of post-formation dynamical processes.Note that, since mass is more fundamental in our simulations, we quantify the intra-system uniformity in our models using the dimensionless mass uniformity metric where   p  and m p  are the standard deviation and average  p of the th system, and  sys is the number of systems in an ensemble (Goldberg & Batygin 2022).Lower the value of D, higher the intra-system uniformity.For Kepler's observed sub-jovian planetary systems D obs = 0.49 ± 0.01.5 Initially, D ≈ 1 in our ensemble, indicating complete lack of intra-system uniformity.We find that during dynamical instabilities, planet-planet collisions (and ejections in rare occassions) redistribute masses among the remaining planets in a way that decreases D. For the final intrinsic and model-detected populations D = 0.75 and 0.49, respectively.The latter is in excellent agreement with D obs .This indicates that late-stage dynamical interactions in multi-planet systems can create the observed level of intra-system uniformity even if initially no such uniformity were present.Most recently, Lammers et al. (2023, , hereafter L23) come to the same conclusion independently despite a very different initial setup (see section B for a detailed discussion).Interestingly, Goldberg & Batygin (2022) found that when started with a higher level of intra-system uniformity than observed, dynamical interactions, in fact, push systems towards reduced uniformity.On the other hand, L23 and this work independently come to the same conclusion that planetary systems with no initial intra-system uniformity, through dynamics, become more uniform.Thus, the observed level of uniformity in the Kepler's multis may represent a natural outcome of planetary dynamics independent of the initial setup.Investigating D obs as a function of observed  p may be a potential way to discern between whether the observed systems came from more or less initial intra-system uniformity.For example, if they were born more uniform and dynamical evolution reduced the uniformity to what we observe today, then the observed higher-multiplicity systems that are expected to have been less dynamically morphed should have a higher intra-system uniformity (lower D obs ) compared to the observed lower-multiplicity systems.On the other hand, if they were born with less intra-system uniformity and dynamics increased the uniformity, then the opposite trend would be observed.We find that D obs increases with increasing  p in the Kepler multis (Figure 7).This trend favors the scenario where a less uniform initial population attains more intra-system uniformity through dynamics.Nevertheless, we caution the readers that the number of high-multiplicity observed systems decreases rapidly with increasing  p , so a more careful investigation on this trend may be warranted in a future study.Also, note that the past studies do not take into account transit and detection biases like ours.
While D shows the overall intra-system uniformity in a population of planetary systems, it does not show the distribution of uniformity in the systems within the population.We explore this using the adjusted Gini index (G  ; Goyal & Wang 2022) for our ensemble. 6For a 6 G is routinely used for measuring the economic inequality in a population (Gini 1912) which was later adjusted for small sample sizes by Deltas (2003).system with  p planets, we can write By construct, G  = 0 for a completely uniform system and tends to 1 for nonuniform systems.Figure 8 shows the distribution of G  in our fiducial ensemble.Since we independently draw the initial planets in each system, initially G  are very large (median 0.57).Dynamical instabilities redistribute masses among the remaining planets to decrease G  in each system.Thus, the G  distribution shifts to lower values for the final intrinsic population indicating increased uniformity.Moreover, the G  distribution for the final model-detected population is in excellent agreement with that for the observed multis.Thus, starting from planetary systems with completely uncorrelated planets, through a combination of dynamical encounters and transit and detection biases, our models produce planetary systems containing the same degree of intra-system uniformity as Kepler's multis.
We investigate this more in depth.To understand if the reduced G  is due to the reduced  p , we conduct a numerical experiment.We randomly choose pairs of planets from the initial planetary systems and merge them until the number of collisions is the same as that found in the -body simulations.We find that the resulting planetary systems are not as uniform as the post-instability systems in our -body simulations (Figure 8).This shows that the level of intra-system uniformity in the final planetary systems is not simply due to a reduced  p .The reason behind the induced uniformity is intricately related to the details of how planets of different  p interact with each other as the NAMD of the system increases during dynamical instabilities.For example, orbits of lower- p planets would have higher excitation for the same NAMD which makes them more prone to collisions.Collisions between preferentially lower- p planets would tend to decrease the range in  p in a system, thus increasing uniformity.
Kepler's planets show similarities between planetary sizes within a system and are comparatively dissimilar to planets randomly drawn from other systems (Weiss et al. 2018).This trend is apparent when the G  distributions are compared for the observed population and a population created by replacing the planets in each multi with those drawn randomly from all observed multis (Figure 8).Our modeldetected population reproduce this trend qualitatively (Figure 8).However, we find that our randomly drawn final detectable systems are somewhat less diverse than the observed inter-system diversity.This is expected.By design, we do not introduce any correlations between the host and the planets.In real systems, however, several birth and evolutionary conditions may introduce additional intersystem diversity relative to our controlled experiment.
Another aspect of 'peas in a pod' is that the  orb ratios between adjacent planet pairs within a system with  p ≥ 3 tend to be more correlated than those in the overall population of multis (Weiss et al. 2018).We do not find any enhanced correlations arising from dynamics either in the intrinsic or the model-detected populations.If indeed the observed systems show significant additional intra-system regularities in spacings, it may have a different origin; for example, forming in resonant chains (Izidoro et al. 2017(Izidoro et al. , 2021;;Goldberg & Batygin 2022;Ghosh & Chatterjee 2023) or a preferred initial  orb distribution different from ours (L23).

Variations In Initial Conditions
So far, we have discussed the effects of planetary dynamics using our fiducial ensemble (n8-e040-i024).Here, we examine whether our findings depend on the initial conditions by varying the initial properties including , , and  p,ini one at a time.In each case, we compare the ensembles when the planetary systems achieve the same dynamical state as the observed (subsection 2.5).

Effects of the initial eccentricity distribution
Higher initial ē leads to earlier onset of instabilities.Hence, ensembles with higher initial ē reach the desired dynamical state earlier (Table 1).Nevertheless, when an ensemble reaches the same dynamical state as observed, independent of the initial ē, the distributions of observable properties are consistent with those for the fiducial ensemble and the observed (Figure 9).Note that, since the desired dynamical state is determined by comparing the model-detected population with the observed, and the detectability of a particular planet does depend on  via differences in the transit duration, small differences may remain in the final intrinsic distributions.
Figure 10 shows the initial and final  distributions scaled by the reduced Hill radius ℎ ≡ (  p 3 * ) 1/3 , for several ensembles differing from each other only in the initial ē.We find that for ē ≲ 0.035, the final intrinsic /ℎ distributions are approximately similar to each other.Moreover, the /ℎ distributions move towards higher values through dynamics in these ensembles with low initial ē.In contrast, for ensembles with initial ē ≳ 0.04, the final intrinsic /ℎ distributions move towards lower values compared to the initial.We find that the observed /ℎ distribution for Kepler's multis is consistent with  Kepler -------0.76± 0.02 0.16 ± 0.009 0.06 ± 0.006 0.021 ± 0.0034 0.009 ± 0.0022 0.001 ± 0.0008 a This ensemble contains intrinsic planetary systems with  p = 9 (  p ,9 = 0.05) and 10 (  p ,10 = 0.03).We combine them with  p ,8 .our ensembles with lower initial eccentricities.Interestingly, because of the contrasting nature of how /ℎ distributions evolve via dynamical interactions in initially low and high ē ensembles, it appears that independent of the initial eccentricities, the final /ℎ distribution always approaches the observed /ℎ distribution.Nevertheless, note that at present it is challenging to undertake an unbiased detailed study of this given the very small (56) number of Kepler planets with both  p and  measurements.
We find an anti-correlation between the NAMD of the planetary systems and multiplicity, except for ensembles with very high initial NAMD (Figure 11).The anti-correlation in our models is easy to understand.Planet-planet interactions increase the overall NAMD in a system, while collisions and ejections decrease it.However, the decrease is usually not sufficient to bring the NAMD back to the pre-scattering level (Figure 1).In our models, the intrinsically lower-multiplicity systems have gone through more dynamical activities, hence, they tend to exhibit higher NAMD compared to highmultiplicity systems.A similar anti-correlation has been inferred for the observed multis (He et al. 2020;Turrini et al. 2020).This trend, however, disappears in ensembles with very high initial ē (≳ 0.04).In those ensembles, the initial NAMD is so high that even the highmultiplicity systems that are less dynamically evolved, retain the initial high NAMD.So our high eccentricity ensembles, or in general, the ensembles with high initial NAMD, are not commensurate with the observed data if the inferred NAMD anti-correlation with multiplicity is real.This likely indicates that the pre-instability Kepler planets had minimal orbital excitation.On the other hand, the NAMD values for the observed high-multiplicity systems may be effectively used to put constraints on the initial NAMD since there is not much difference between the final model detected NAMD and initial.
As expected, this trend is directly translated to orbital eccentricities.The anti-correlation between  and  p in the final intrinsic population ranges from C = −0.18 in the lowest initial ē ensemble The box and whiskers represent the 25th-75th and 16th-84th percentiles.
The filled boxes represent the final intrinsic NAMD distributions in each bin, and the empty boxes represent the initial NAMD distributions of the same systems.Note that the NAMD values for the bin with  p = 2 may be unreliable due to very few data points (Table 3).
Other properties such as the multiplicity distributions of the modeldetected populations do not show any significant differences depending on the initial ē (Table 2, Table 3).

Effect of Initial Inclination
Similar to the finding in the case of ensembles with different initial ē, we find that when the dynamical states of our ensembles are consistent with that of Kepler's multis, the distributions of other observable properties are consistent with each other and the observed, independent of the initial ī (Figure 12).The only difference is that the ensembles with higher initial ī reach the target dynamical state earlier (Table 1).Of course,  m plays a significant role in determining transit probability in multis.As a result, when the  distributions for the model-detected populations are consistent with the observed, the intrinsic distributions may show subtle differences.For example, the intrinsic distributions for ,  orb ratios, and  p shift to lower values with increasing initial ī (Table 2).
The final intrinsic  m /ℎ distribution depends on the initial ī (Figure 13).Through dynamical evolution the  m /ℎ distribution shifts to lower (higher) values in the ensembles created with initial ī = 0.05 (0.01).For the ensemble with a medium initial ī = 0.024, the  m /ℎ distribution shows little change.Thus, it appears that a natural out- come of dynamical evolution (for the class of multis studied here) is a preferred intrinsic  m /ℎ distribution; planetary systems initially more aligned than this become more misaligned and vice-versa.This trend is qualitatively similar to the trend with /ℎ distributions for ensembles modeled with different initial ē (subsubsection 3.2.1).It will be very interesting to verify whether indeed the intrinsic  m and  distributions exhibit such a tendency.However, estimating intrinsic  m and  distribution for the observed multis, ideally in a model agnostic way, is challenging and is beyond the scope of this work.We find some subtle effects in the multiplicity distribution depending on the initial ī (Table 3).For example, ensemble n8-e040-i050 with initial ī = 0.05 intrinsically retains more systems with  p > 6.However, due to higher final  m , these planets have lower transit probabilities.As a result, the fraction of detected systems with  p > 4 is lower in this ensemble relative to the ensemble with ī = 0.01 (n8-e040-i010, Table 3).While the final intrinsic  m distribution has large spreads in all ensembles, the median  m shows a modest anti-correlation with  p , especially in the initially low ī ensemble (C = −0.11for n8-e040-i010).

Effect of Initial Number of Planets
As expected, systems with higher  p,ini become unstable at a shorter  and the ensemble reaches the desired dynamical state earlier (Table 1).When the model-detected  distributions agree with the observed, the other properties are also in agreement independent of  p,ini (Figure 14).Although the model-detected populations show excellent agreement, the intrinsic  orb distributions are quite different between these ensembles.The lower- p,ini ensembles produce intrinsic systems with lower  orb .This is directly related to our setup (subsection 2.2).Systems with lower  p,ini are naturally smaller in extent and contain planets restricted to lower  orb .In spite of this, the model-detected  orb distributions are consistent, independent of  p,ini (Table 2).Note that our lowest- p,ini ensemble n6-e040-i024 does not quite reach the desired dynamical state even at  = 10 9  1 .For this ensemble, model-detected distributions for , period ratios, and  orb peak at somewhat lower values when compared to the Kepler data (Table 2).We believe, had we integrated n6-e040-i024 longer, this ensemble also would have converged to similar distributions in these properties.
We also find that the final  and  m distributions peak at higher values for ensembles with higher  p,ini , despite having identical ini- tial distributions (Table 2), indicating that the orbits in the systems with a higher  p,ini get more churned.In our simulations, the intrinsic planet multiplicities can go up to  p,ini .However, the model-detected populations contain very few (< 0.07%) systems with  p > 6, irrespective of  p,ini .Ensembles with higher  p,ini have higher fractions of detectable planetary systems with  p ≥ 4 (Table 3).
Independent of all the variations in initial properties, none of the ensembles could produce the observed apparent excess of singletransiting planets.

SUMMARY AND CONCLUSIONS
In this study, we have investigated the role of dynamical instabilities in the final stages of assembly for planetary systems initially in non-resonant configurations.Starting from a power-law  distribution we simulate our planetary systems with varying initial ē, ī, and  p,ini until the model-detected systems in ensemble reach the observed dynamical state, equivalently,  distribution for the Kepler multis.Note that this stopping criteria is more stringent and physically motivated compared to the typically used criteria of stopping after a fixed number of orbits when further interactions slow down (section A).Also, note that because of our dynamically motivated stopping criteria, the match between the model-detected and observed  distributions should not be considered as validation of models.On the other hand, the fact that dynamical instabilities can produce the same  distribution as observed from a large variety of initial conditions, indicates that the final distributions of properties are rather insensitive to the details of the initial conditions, thus eliminates any need of fine tuning.By design, we do not inject any inter or intra-system correlations at the beginning to find what emerges solely from dynamical encounters.Below we summarize the main findings.
• When the dynamical states of our model-detected populations match that of the observed multis, distributions of all other key observable properties agree as well (e.g., Figure 3).This agreement does not strongly depend on the initial distributions of orbital properties or  p,ini (Figure 9,12,14).
• Without any initial input from formation scenarios, host-star dependence, or intra-system correlations, all model ensembles naturally evolve to produce planetary systems very similar in properties to the observed.
• Starting from completely uncorrelated planets, dynamical evolution naturally induces intra-system uniformity in  p in the modeldetected populations similar to the observed 'peas-in-a-pod' trend (subsubsection 3.1.2,Figure 8).The Kepler multis seem to favor an initially less correlated population attaining intra-system uniformity through dynamics (Figure 7).
• The anti-correlation between NAMD and multiplicity manifests as a natural consequence of dynamical instabilities if the initial NAMD is sufficiently low.In this scenario, present-day systems with higher  p have suffered less dynamical instability in the past.
• This NAMD-multiplicity trend, however, disappears in ensembles where the initial NAMD is so high that even the less dynamically evolved high-multiplicity systems retain the initial high NAMD.So if the inferred NAMD-multiplicity anti-correlation in Kepler multis is real, pre-instability Kepler planets likely had low orbital excitation.On the other hand, the initial NAMD may be constrained from that observed in present-day high-multiplicity systems.
• The multiplicity distributions in model-detected populations are similar to those observed for  p > 2. However, the observed systems show an excess of singles.This discrepancy cannot be explained simply by varying the initial orbital properties or  p,ini within our setup.The observed excess can be explained if  ≈ 0.48 single planets are added to the mix.
• Likely related to the previous point, although, the modeldetected ē for multis is similar to those observed, the model-detected singles have significantly lower ē compared to the observed singles.This indicates the necessity of an additional population of intrinsically singles, or systems containing higher-mass perturbers which can create dynamically hotter inner planetary systems of sub-Neptunes (e.g., Hansen 2017;Lai & Pu 2017;Pu & Lai 2018;Bitsch & Izidoro 2023).
The success of our model in reproducing several key trends of the observed multis, despite our initial setup was intentionally devoid of any such correlations, indicates that many of the famous observed trends and correlations may have resulted from dynamical processes post-formation and does not point towards specific formation scenarios.

APPENDIX A: INTEGRATION STOPPING TIME
Here, we compare our choice of simulation-stopping criteria with the common practice in the literature.Traditionally, the ensembles of planetary systems are integrated for ∼ 10   1 , (the adopted x depends on the computational resources, typical values are 7, 8, and 9 depending on the study).Then usually, if the distributions of some critical properties do not change significantly by extending the integration to a few ×10   1 , the system is considered 'settled'.Essentially, it means that the occurrence rate of instabilities has decreased sufficiently.This can be problematic because the stability timescale grows exponentially with the Hill spacing (e.g., Chambers et al. 1996;Chatterjee et al. 2008;Obertas et al. 2017).So, if a system is integrated till 10   1 , it may acquire properties consistent with a stability timescale of ∼ 10   1 .But, the next significant change is expected only at ∼ 10 +1  1 .So, checking for changes at a few ×10   1 may not be indicative.It is also not clear when to stop.Since the systems integrated to 10 +1  1 may again become unstable if integrated to say, 10 +2  1 , and so on.One natural point to stop would be when only two planets remain.Then we can use a suitable analytic stability criteria (e.g., Gladman 1993) to determine the possibility of orbit crossing.Else, a natural stopping point would be the system age.Numerical integration for a large number of systems till their actual ages is computationally impractical.
Instead, our criteria (as described in subsection 2.5) is more physically motivated.The  distribution in a population essentially indicates how stable or unstable the population is.Hence, bringing all model ensembles to the observed  distribution essentially ensures that the model-detected ensembles are dynamically as active or inactive as the observed.This demanding criteria already fulfills the traditional settling-down criteria.For example, our fiducial ensemble (n8-e040-i024) satisfies our stopping criteria earliest at ∼ 10 8  1 .We find no significant changes in the  distribution after further integration (Figure A1, left panels) to 2 × 10 8  1 or even 10 9  1 .In fact, changes in the ensemble properties slow down significantly after ∼ 10 7  1 .On the other hand, different ensembles may reach this dynamical state at different multiples of  1 based on their initial conditions (Table 1), since the initial conditions dictate the timescale for the onset of instabilities.For example, the ensemble (n8-e060-i024) satisfies our stopping criteria at 5 × 10 6  1 , and we find no significant changes in the Hill spacing distribution, even if we continue the integration for more than an order of magnitude longer (Figure A1, right panels).In contrast, the ensembles with lower initial orbital excitation require longer simulations (up to ∼ 10 9  1 performed in our models) to reach the same dynamical state (Table 1).Thus, using a fixed multiple of  1 for stopping time would not necessarily mean that the dynamical states of the ensembles, in the end, are similar.Thus, when ensembles with varied initial orbital properties are considered, we argue that all ensembles should be brought to a compatible dynamical state consistent with the observed before comparison.

APPENDIX B: INTRA-SYSTEM UNIFORMITY
While preparing our manuscript, we came to know about the results of L23.In this fantastic work, L23 independently reach the same conclusion regarding the intra-system uniformity induced by dynamical instabilities despite a very different initial setup compared to ours.This strengthens the result that dynamical instabilities in multiplanet systems with little intra-system uniformity, independent of the initial setup, can make systems more uniform over time.Here we compare our results with an additional simulated ensemble using the exact initial setup of L23.The initial setup of L23 is significantly different from ours.Moreover, unlike ours, the scope of their study is limited to investigating the origins of intra-system uniformity.Hence, we limit our comparisons only to this aspect.For completeness, L23 consider initial planetary systems with 10 planets orbiting a sun-like star.They assume a flat initial period ratio distribution within [1.10, 1.50), and draw  p from a Gaussian distribution centered around  p / ⊕ = 3.They use the  p - p relationship given in Wolfgang et al. (2016).We replicate these initial conditions in this additional ensemble.
L23 find initial (final intrinsic) D = 0.59 (0.48) after 10 9  1 .In our simulated ensemble using their initial setup, we also find a similar result, initial (final intrinsic) D = 0.60 (0.48).L23 have not considered the effects of various detection and sensitivity biases in their analysis.Our results suggest that this can make a difference (Figure 8, subsubsection 3.1.2).When we apply our planet detection algorithm (subsection 2.4) to create a Kepler-detectable population from this ensemble, we find a model-detected D = 0.35, significantly lower than observed (subsubsection 3.1.2,Goldberg & Batygin (2022), Lammers et al. (2023)).
Interestingly, L23 report intra-system uniformity in spacing measured in period ratios.They find that period ratios within a system with  p ≥ 3 in their intrinsic population are more correlated (correlation coefficient C = 0.39) than those overall.We find a similar result, C = 0.48, in our intrinsic final ensemble using their initial setup.However, this enhanced intra-system correlation in spacing disappears when transit and detection biases are taken into account.Our model-detected population using their initial setup shows C = 0.03.Interestingly, if we restrict our fiducial ensemble to 4 <  < 20, roughly equivalent to the range of period ratios considered in L23, we also find significant intra-system uniformity in spacing in the final intrinsic population, C = 0.21.However, this uniformity disappears when transit and detection biases are taken into account, C = 0.03.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Evolution of the NAMD for three different example systems.Although all systems initially had roughly the same NAMD, they follow different dynamical evolution with different final  p and NAMD.The crosses denote planet-planet collisions in each system.

Figure 2 .
Figure 2. Probability distribution function (PDF, top) and cumulative density function (CDF, bottom) for the mutual Hill separations () in the fiducial ensemble (n8-e040-i024).Black denotes Kepler's multis.The dotted, dashed, and solid blue lines denote the initial, final intrinsic, and final model-detected populations.The model-detected and observed  distributions are in excellent agreement.

Figure 3 .
Figure3.PDF (top)  and CDF (bottom) for adjacent-planet  orb ratio (left),  orb (middle), and  p (right) for the fiducial ensemble n8-e040-i024.Line colors and styles denote the same populations as in Figure2.When the dynamical state of the simulated population is similar to that observed, other observable properties show good agreement with the Kepler data (shown in black).The final radii of the synthetic planets from our model are obtained from the masses using the mass-radius relationship fromChen & Kipping (2016).

Figure 4 .
Figure 4.  vs  m (top) and  p (bottom) for the final intrinsic planetary systems in the fiducial ensemble.Dots denote individual planets and the contours confine 11.8%, 39.3%, 67.5%, and 86.5% of all systems.The corresponding marginalised PDFs are shown on the top and side panels where dotted, dashed hatched, and filled histograms denote the initial, final intrinsic, and final model-detected populations, respectively.Black dots (with errorbars) denote Kepler planets with measured  p and .

FractionFigure 5 .Figure 6 .
Figure 5. Multiplicity distribution.Dashed blue, filled blue, and grey bars denote the final intrinsic, final model-detected, and observed populations, respectively.Red bars denote the multiplicity distribution after artificially adding single-planet systems (  = 0.48) to model-detected population.Errorbars for model-detected denotes 1 estimated using bootstrap with a sample size equal to the number of Kepler systems.Poisson errors are used for the Kepler data.There is an excess of observed single-transiting systems relative to our model-detected population.Injection of additional single-planet systems significantly improves the agreement.

Figure 7 .Figure 8 .
Figure 7. D vs  p for our model-detected population (blue) and Kepler multis (grey).The error bars in D for the Kepler multis in each bin represent 25th and 75th percentiles resulting from the intrinsic scatter in the  p - p relationship from Chen & Kipping (2016) within 1.

Figure 10 .Figure 11 .
Figure 10.CDF for /ℎ for ensembles created with different initial ē.Solid (dotted) lines denote the final intrinsic (initial) population.The black line represents Kepler's planets with measured  p and .

Figure A1 .
Figure A1.CDF for  at different snapshots in time for our fiducial ensemble n8-e040-i024 (left) and n8-e060-i024, an ensemble with high initial  (right).Top and bottom panels show the final intrinsic and model-detected distributions.Different colors show the  distributions at different times (see legend).Black denotes Kepler's multis.The snapshot at the stopping time determined by our integration stopping criteria (subsection 2.5) is shown in blue and marked by an asterisk in the legend.Longer simulations show no significant difference.

Table 1 .
Initial properties of the simulated ensembles and the corresponding  stop (subsection 2.5).Ensemble names contain  p,ini , ē, and ī for easy understanding of the corresponding initial properties.

Table 2 .
Median, 16and 84 percentiles for final properties for the intrinsic and model-detected populations.ē and ī m denotes the mean of the best-fit Rayleigh distributions for eccentricities and mutual inclinations (measured with respect to the system invariant plane).

Table 3 .
Fractions (  p ,  ) of planetary systems with  p = .Errors represent 1 estimated using bootstrap with a sample size equal to the number of Kepler systems.The limits in the Kepler dataset represents Poisson errors. p ,3  p ,4  p ,5  p ,6  p ,7  p ,8