Orbital stability of compact three-planet systems, II: Post-instability impact behaviour

Recent observational missions have uncovered a significant number of compact multi-exoplanet systems. The tight orbital spacing of these systems has led to much effort being applied to the understanding of their stability; however, a key limitation of the majority of these studies is the termination of simulations as soon as the orbits of two planets cross. In this work we explore the stability of compact, three-planet systems and continue our simulations all the way to the first collision of planets to yield a better understanding of the lifetime of these systems. We perform over $25,000$ integrations of a Sun-like star orbited by three Earth-like secondaries for up to a billion orbits to explore a wide parameter space of initial conditions in both the co-planar and inclined cases, with a focus on the initial orbital spacing. We calculate the probability of collision over time and determine the probability of collision between specific pairs of planets. We find systems that persist for over $10^8$ orbits after an orbital crossing and show how the post-instability survival time of systems depends upon the initial orbital separation, mutual inclination, planetary radius, and the closest encounter experienced. Additionally, we examine the effects of very small changes in the initial positions of the planets upon the time to collision and show the effect that the choice of integrator can have upon simulation results. We generalise our results throughout to show both the behaviour of systems with an inner planet initially located at $1$ AU and $0.25$ AU.


INTRODUCTION
The now retired NASA Kepler Space Telescope is responsible for observations leading to the confirmation of hundreds of multi-planet systems (Lissauer et al. 2014;Rowe et al. 2014). Of these systems, as many as six percent are thought to be compact (Wu et al. 2019), containing planets that are much more closely spaced than the inner planets of our own Solar System. These discoveries have naturally led to many questions being asked about the long-term stability of compact exoplanet systems. Indeed, it is even possible that compact planetary embryos existed interior to Venus's current orbit that have subsequently been expelled from this region due to orbital instabilities (Volk & Gladman 2015). Within the class of observed compact systems, a large population of planets have been observed with a mass (Mayor et al. 2011) and radius (Petigura et al. 2013) between that of Earth and Neptune. Moreover, the observed orbital architecture is such that mutual inclinations are small, typically in the region of 1 • to 2 • (Fabrycky et al. 2014), while eccentricities are also found to be small, on average¯≈ 0.04 (Xie et al. 2016). An archetypal example of these systems, albeit containing six planets, is Kepler-11 (Lissauer et al. 2011). Exoplanet systems with orbital spacings much greater than that required for stability are also present in the Kepler data set. It is a favoured hypothesis that this orbital architecture is a result of dynamical instabilities in much more compact systems leading to ★ E-mail: p.bartram@soton.ac.uk close encounters and orbital reconfiguration (Pu & Wu 2015). Understanding of the stability and evolution of compact exoplanet systems is therefore not only important for making sense of observations but also for understanding the planetary formation process as a whole.
Characterisation of the stability of three or more planet systems can be approached in several ways. Analytical models have been built that can predict the lifetime of three planet systems based upon resonance overlap (Wisdom 1980;Quillen 2011;Petit et al. 2020). Recently, machine learning approaches have also been developed that, after being guided by a training set of 10 9 year integrations, can use far shorter integrations to predict with surprisingly high accuracy which given exoplanet systems will remain stable for a billion orbital periods (Tamayo et al. 2016. However, the most common approach to the problem, and the one employed in this paper, is the use of n-body simulation (Chambers 1999;Smith & Lissauer 2009;Obertas et al. 2017;Hussain & Tamayo 2020;Lissauer & Gavino 2021).
The majority of studies performed take a subset of the possible input parameter space for a compact, near-circular, near co-planar system of a given number of planets and then evolve this system forward in time checking for either the first close approach, typically specified as one Hill radius, , or waiting for an orbital crossing to occur: this is then termed the instability event. Throughout this work we will use orbital crossing as our definition of an instability event and refer to the time at this point as the crossing time. Rice et al. (2018) found that systems containing four Neptune-size and Neptune-mass planets initially located at 1 AU can continue to evolve after an instability event for over ten million dynamical periods before a collision of planets, meaning that the commonly used instability metric may not capture the entire evolution of the system. Given that the manner in which these planets collide determines the final orbital architecture it is important properly to understand this phase of the exoplanet system life cycle.
Our study builds upon the work done by Rice et al. (2018) and Lissauer & Gavino (2021) by considering the post-instability evolution of compact, Earth-analogue, three-planet systems across a large range of initial orbital separations equally spaced in units of mutual Hill radii. We create three integration suites called the standard suite, perturbed suite and inclined suite, and perform 4, 835 integrations each in the first two and a further 11, 200 in the final one. We continue integrations up until the time of first collision between planets or for 10 8 or 10 9 orbits depending on the experiment.
In section 2 of this paper we describe the methodology used for our integrations including the initial conditions for each integration suite, the integration packages used and the termination criteria. Section 3 contains the results of all standard suite integrations: section 3.1 details the time scales for orbital crossing and collision between pairs of planets, and details collision probabilities over time for various initial configurations of systems; the effects of small changes in initial orbital longitudes upon these results are then examined in section 3.2; and, finally, section 3.3 examines the probabilities of particular pairs of planets colliding. Section 4 introduces the results of inclined suite integrations: in section 4.1 we explore the heating of what are initially dynamically cold systems that eventually enables orbital crossing and collision; here, we find that the three-planet Earth-mass systems behave in a similar manner to the four-planet Neptune-mass case but follow a different power law. Section 4.2 examines the time scales leading to collision in the inclined case and shows that the survival time after crossing can be a non-trivial fraction of the mainsequence lifetime of stars. In addition, this section also looks at the effects on lifetime of systems dependent on the distance from the innermost planet to the star and the initial inclination. We summarise our findings in section 5.

METHODS
We have chosen to simulate three-planet systems comprising of analogues from our own Solar System. The central body in each of our systems is a one solar mass star, 0 = 1 M . Each of the planets within the systems are Earth mass, = 1 M ⊕ where = 1, 2, 3 with a planetary radius also equal to that of Earth, R = R ⊕ . Planets are placed on initially circular orbits orbiting the star in a common direction with the innermost planet located at 1 AU. Time throughout this work is provided in units of initial orbital period of the innermost planet, this means that the crossing time is invariant to rescaling of the system so long as initial orbital period ratios between bodies are maintained along with mass-ratios of planets and star.

Initial semi-major axes
Initial semi-major axes of systems are evenly spaced in terms of mutual Hill radii. The mutual Hill radii are defined as , +1 = (1) This allows for a dimensionless value to be defined to specify the even spacing of adjacent planetary orbits in units of their mutual Hill radii as .
(2) Therefore, the initial semi-major axes of adjacent planets are chosen to be such that The innermost planet is placed such that it has a semi-major axis of 1 AU, and all other semi-major axes are chosen through Eq. (3). We refer to this configuration as a system at 1 AU. Likewise, later on, when results are generalised to include systems with an innermost planet located at 0.25 AU with other planets spaced as per Eq. (3) we refer to it as a system at 0.25 AU.

Stopping criteria and integration packages
We have opted to use the Terrestrial Exoplanet Simulator (TES) 1 package to perform our integrations (Bartram & Wittig 2021). TES is a new numerical integration package written in C++ for propagating exoplanet systems. This package combines an integrator that follows Brouwer's law (Brouwer 1937) with a new special perturbation method to allow for reduced run-times and decreased numerical error resulting in, e.g., improved energy conservation. Additionally, this tool has been designed to allow for integration all the way to collision of terrestrial mass planets to machine precision. TES can be run using C++ directly, or through a python interface allowing for ease of use and for multiple integrations to be performed in parallel. Throughout our simulations we have opted to use TES with a non-dimensional tolerance of 1 × 10 −8 which has ensured that the relative energy error in all simulations, even after collision, and for the longest lived systems, is maintained below 1 × 10 −13 . To validate our own results, we also repeated all of our standard suite integrations making use of IAS15 (Rein & Spiegel 2015) within the REBOUND package (Rein & Liu 2012). The results from this comparison can be found in Appendix A.
As mentioned before, time is measured by periods of the innermost planet in the system throughout this work, meaning that all times are specified in units of orbits or dynamical periods. Integrations run until either a collision is detected or the simulation reaches a maximum time of 10 8 or 10 9 dynamical periods, depending on the experiment.
In order to detect an orbital crossing, the orbital elements of each planet are calculated at every step within each integration. These are then compared to determine the time at which the apoapsis of a planet crosses the periapsis of the exterior adjacent planet. We define this as the crossing time and denote it . Moreover, also at each step, the mutual separations of each of the planets are calculated so that collisions can be detected. The metric of two planets coming within 2R ⊕ of one another is used for collision detection. We define the time at which this occurs as the impact time and denote it . We also define the post-crossing survival time, , of a system to be the time Simulations are run for up to 10 9 orbits in general but some are terminate at 10 8 orbits to save on computation. Orbits are specified by the initial period of the innermost planet. Impacts that take place before a crossing are highlighted by a green diamond whereas systems that did not cross within the maximum simulation time are marked with a red triangle. Models fitted to the crossing and impact times according to Eq. 4 are shown as a dashed black and a dashed red line, respectively.
All encounters closer than any experienced previously are recorded such that it is possible to generalise the collision results to systems with planetary radii greater than that of the Earth or, equivalently, initial orbital radii closer than 1 AU. We use this generalisation to consider systems at 0.25 AU and 1 AU for all integration suites. We also define the time of closest encounter prior to collision as the closest encounter time, . To ensure bit-wise identical initial conditions as in Lissauer & Gavino (2021), initial conditions are specified as orbital elements which are then entered in to the MERCURY (Chambers 1999) integration package in order to generate an initial state vector which is then provided to either TES or REBOUND. Table 1 contains a summary of all symbols related to simulations event times.

Standard integration suite
The first suite of integrations is composed of 4, 835 orbital configurations and is termed our standard suite. In this suite systems are on initially circular, co-planar orbits with an initial mean anomaly for the th planet = 2 radians where ≡ 1 2 1 + √ 5 , i.e., the golden ratio, and are merely chosen to avoid special orientations. As we wish to study the effects of the initial spacing of planets upon impact timescales we choose a high resolution in such that there are 1 × 10 3 integrations per unit over the range = [3.465, 8.3]. Generally, integrations are terminated after 10 9 orbits if a collision is not encountered. However, in certain areas we have chosen to limit integrations to 10 8 orbits to save on computation; these regions are clearly marked on any plots.

Perturbed integration suite
The second integration suite is termed our perturbed suite and is also composed of 4, 835 integrations. The only difference between the initial conditions of the standard suite and the perturbed suite is that in the latter case the innermost planet is perturbed by 100 m along its orbital arc. We strictly terminate integration at 1×10 8 orbital periods of the innermost planet in this suite. This suite is used to examine the effects of very small changes in initial conditions upon crossing and impact time.

Inclined integration suite
The final integration suite is the inclined suite and is composed of 16, 800 integrations. Of these, 15 did not complete in the available CPU time and have been excluded from the dataset. This is equivalent to 0.09% of inclined integrations, and we therefore do not believe this will have biased the dataset in any meaningful way. We choose initial conditions across a subset of the available parameter space manually rather than randomly and perform integrations for a maximum simulation time of 1×10 8 orbital periods of the innermost planet. To make best use of computational resources we limit this study to the range = [3.5, 6.3] and perform experiments uniformly spaced in with fifty values per unit . At each value of we perform one hundred Cumulative sum of integrations with a collision before orbital crossing for various initial values for semi-major axis of the innermost planet. The flat region between beta = 7 and = 8 is due to systems not experiencing an orbital crossing within the maximum simulation time in that region and integrations being terminated early (see the red triangles in Figure 1). A large fraction of integrations with an initial spacing > 6.3 were stopped at 10 8 orbits so results beyond this value cannot be considered to have been drawn from a uniform sample.
and twenty experiments where the initial values of semi-major axis, eccentricity and mean longitude are the same as in the standard suite. Planets are, however, inclined relative to each other in one of four ways: one of inner, middle or outer planet inclined above the orbital plane of the system, and also with the middle planet above and the outer planet below. For each such configuration of relative inclination fifteen initial values of inclination are logarithmically spaced between 0 = 0.06 • and 0 = 0.58 • , yielding an initial orbital height ranging from 0.10 to . The distribution of initial inclinations within this range is such that ten values are used between 0 = 0.24 • and 0 = 0.58 • and five values are used over the region 0 = 0.06 • and 0 = 0.24 • . Finally, two values are chosen for the ascending nodes Ω: either according to the golden ratio in Section 2 such that = Ω or equally spaced such that Ω = [0 • , 120 • , 240 • ].
The full state vector of each simulation is output to file once every ten thousand orbital periods; additionally, each planetary flyby closer than any other previously observed is also recorded.

STANDARD INTEGRATION SUITE
This section contains results from the standard integration suite described in section 2.3. Additionally, the results of the perturbed integration suite, described in section 2.4, are analysed in subsection 3.2.

Timescale to planet-planet collision
The crossing and impact times for the standard suite are plotted in Fig. 1. Inspection of the crossing time with respect to the initial orbital spacing shows the clear upwards trend present in other works (Smith & Lissauer 2009;Obertas et al. 2017;Hussain & Tamayo 2020;Tamayo et al. 2020;Lissauer & Gavino 2021). We also capture the large scale variations about the trend which for the most part are a result of mean motion resonances as discussed in Obertas et al. (2017). Additionally, we replicate the finding of Lissauer & Gavino (2021) in the discovery of a highly stable configuration around = 5.74 which they attribute to the distance of this configuration from any strong resonances. Throughout this work we used a linear logarithmic fit of the form in several places where = − 2 √ 3 and is used to reduce the dependency of the y-intercept upon the slope. We fit this model to three data sets such that = , or and state explicitly which at the time of use. Unless otherwise stated, we only include data points in the region = [3.465, 6.3] in the fits to avoid biasing the results due to systems that did not experience an orbital crossing within the maximum simulation time. For = , over this region, we find that = 1.352 and = 2.067 which is in strong agreement with Lissauer & Gavino (2021) and confirms the functionality of the TES tool. For impact times we find = 1.192 and = 2.42. Figure 1 highlights that the post-crossing survival time is very small compared to the crossing time for the majority of systems observed. The log scale of the plot and the relatively small magnitude of means the bulk of the impact time data points are hidden in this figure. The only exception is in the region of small where the ratio / is large due to the relatively small size of .
Finally, it can be seen that for a small subset of integrations collisions can occur before an orbital crossing has taken place. A cumulative sum showing the numbers of occurrences is shown in Fig. 2 where we believe that the increase between systems at 1 AU and 0.25 AU is not dependent purely on the physical cross-sectional area of planets but rather the enhanced cross-sectional area due to gravitational focusing (Safronov 1972). It is likely that a symplectic integrator, configured to use the standard step size of 1 /20 th of the smallest dynamical period, would miss these collisions. However, given the small number of occurrences relative to the number of integrations typically performed in stability studies, it is unlikely that these missed collisions will have biased the data-sets in any statistically meaningful way. Figure 3 shows the post-crossing survival time for all systems within the standard suite against ; Figure 5 is identical but plotted against . We find two main populations of post-crossing survival times present: those surviving for less than two orbits, and those surviving for more than ten orbits with very few outliers in between. Within the long surviving population, it can be seen that there is a clear increase in the post-crossing survival time of systems with respect to both and . We fit models of the form of Eq. (4) to both the long-lived population and the population in its entirety, we call these datasets long and all, respectively. The model coefficients and can be found in the top two rows of Table 2. Similarly, we also fit linear models to the two datasets present in Fig. 5 for log 10 ( ) against log 10 ( ). The model coefficients and can be found in the bottom two rows of Table 2. In all cases, we calculate the Pearson correlation coefficient (PCC) and also calculate the standard deviation, , of the other-neighbouring-pair extrema-pair same-pair Figure 5. Post-crossing survival time, , against orbital crossing time, , for systems initially at 1 AU. Blue dots indicate the same pair both crossed orbits and collided; orange indicates the pair that collided was not the pair that crossed; green indicates a collision between the inner and outer planets. The model (dashed black) is fitted to all data points with a survival time greater than two orbits. data minus the fitted model, e.g. (log 10 ( ) − ( + )). Clearly, there is a tendency for systems to persist for longer after an orbital crossing when the initial mutual spacing between them is greater, with a difference of a factor of three in median post-crossing survival time over the entire beta range. However, even given this increase, the post-crossing survival time for systems simulated did not ever exceed one million orbits. Given that this represents roughly one ten-thousandth of the main sequence lifetime of solar-mass stars it is possible, although very unlikely, that we could observe a compact exoplanet system that has undergone an orbital crossing but has not yet experienced a collision between planets, even if it were a truly co-planar system.
In the case of the short-lived population, there is a further subdivision of different behaviours: those systems that experience a collision almost immediately following a crossing, e.g. those bodies whereby < 10 −1 , and those which persist for longer than this but less than a couple of orbits. In the former case, we have observed that the trajectories of two planets about the star simply cross, leading to straightforward collisions, and also triggering an orbital crossing in the process. However, in the latter case, we find that the trajectories of the planets about the star are such that a very close encounter occurs which causes the two planets to become temporarily gravitationally captured. These two planets then remain within approximately a Hill radius of one another before finally experiencing a fatal collision a fraction of an orbit later. These behaviours are shown in the satellite images in Fig. 3. It can be seen here that temporary gravitational capture is not the cause of collision in the case of the outliers with a post-crossing survival time between two and ten orbits.
To give consideration to whether these results generalise to other systems of planets we have calculated and for a system with the inner planet initially placed at 0.25 AU. This is equivalent to artificially inflating the radius of all planets in systems at 1 AU by a factor of four. When thought of this way, this is akin to placing planets with a radius approximately the same size as Neptune at 1 AU; is invariant to the initial location of the inner planet. The probability, calculated as the cumulative fraction of systems that have experienced collisions over the total number of systems, of collision over time for both settings are shown in Fig.6. The separation between dashed and solid lines indicates that a given collision probability is reached sooner in systems composed of planets with a larger radius. The difference in time remains about constant over all values of , even if the log scale suggests otherwise. Figure 4 contains normalised histograms of within different regions of for systems with the inner planet initially at 1 AU and 0.25 AU. We find that the distribution of post-crossing survival times is log-skew-normal distributed across all systems; we confirmed this using a Kolmogorov-Smirnov test with a precision parameter of = 0.005. The skew-normal distribution is a generalisation of the normal distribution that allows the class to be extended to include distributions with non-zero skewness through the addition of a shape parameter (Azzalini & Capitanio 1999). Log-skew-normal probability density functions, shown in orange, are fitted to the data through a maximum likelihood estimator. We calculated the mean , standard deviation , and the skew for each distribution; we use the Fisher-Pearson coefficient of skewness throughout. We find that increases with increasing range, and also find the same pattern for in all but one case. In all cases, is negative indicating a skew towards shorter post-crossing survival times as compared to a normal distribution. This means that there is a preference for systems to collide sooner rather than later after an orbital crossing as compared to the most frequent survival times. There is a slow build up in the number of systems experiencing collisions over time after an orbital crossing but a much sharper cut-off after the peak density of collisions. This highlights the difficulty for systems to persist for long timescales after an orbital crossing in the co-planar case. Systems with a shorter mean post-crossing survival time show a skew of a smaller magnitude than those with a longer survival time, e.g. at 1 AU = −0.14 for < 4.0 whereas for >= 4.0 the smallest, in magnitude, value observed is = −0.57. We find that the distributions of post-crossing survival times at 0.25 AU are less skewed than those at 1 AU, indicating that the survival times of systems in this case are closer to a log-normal distribution.

Sensitivity to initial conditions
To examine the sensitivity to initial conditions of the results of our simulations we use our perturbed suite of integrations described in section 2.4. The crossing and collision times of each integration between the standard suite and the perturbed suite are compared to determine the effect of the perturbation. Table 3 contains the results of that comparison for the time of orbital crossing. Tables 4 and 5 contain the same comparison for but for the impact time of systems at 1 AU and 0.25 AU, respectively.
In general, the comparison between crossing times in Table 3 aligns closely with Lissauer & Gavino (2021). Percentages between the two studies rarely differ by more than a few points despite the different integration tools used: TES and MERCURY. One notable difference between the two studies is in the initially wider spaced systems. In the regions = [5.0, 5.999] and = [6.0, 6.33] we find roughly double the number of initial orbital spacings where the standard and perturbed suite integrations experience orbital crossing times within 10% of one another. Given the precise orbital evolution required in order for standard and perturbed suite systems to experience a crossing at the same time it is unlikely that numerical error would ever cause an increase in this statistic. We therefore take this as an indication that TES has maintained a higher precision than the symplectic Wisdom-Holman (Wisdom & Holman 1991) scheme within MERCURY. To further validate TES in this setting we have also repeated the standard suite integrations with IAS15 from the REBOUND package for comparison. We find very good agreement in results between the two routines. Detailed results from this experiment are included in Appendix A.
The right-most summary column for the full range of = [3.465, 6.33] in Table 4 shows there is a marked decrease in the number of collisions occurring within a factor of two, and within ten and one percent of one another; as compared to the orbital crossing times in Table 3. The largest reduction is seen in the within-a-factorof-two row where a reduction of over 15 percentage points highlights the sensitivity to close approaches in this setting. The majority of this difference in collision times is seen in the initially closely spaced systems where a reduction of almost 50 percentage points can be seen for integrations finishing within 10% of one another. However, once the crossing times exceed approximately 1 × 10 4 orbits at = 5 the effect of the perturbation disappears and values between crossing and collision times for the two data sets converge.

Which planets collide?
We find a slight discrepancy between the prevalence of orbital crossings in our standard integration suite, with the innermost pair triggering 48% of crossings compared to 52% for the outermost pair. These percentages were calculated using = 4, 835 integrations and the expected stochastic variation, about the mean, i.e. 50%, is therefore approximately 0.72% (Dobrovolskis et al. 2007).
In the following, we designate the specific pair of planets that collide as the collision pair, and analogously we refer to the pair of planets that experienced an orbital crossing as the crossing pair. We find that across all values of a collision between two planets is almost twice as likely if the same two planets were also involved in the orbital crossing. Figure 7 highlights clearly that this is the case with between 48% and 62% of collision events occurring between the crossing pair for systems at 1 AU depending on the initial orbital spacing. Moreover, these percentages appear to be invariant as to whether the inner or outer pair was involved in the orbital crossing. A clear trend can be seen with respect to , where an increase in the initial orbital spacing between planets leads to an increased probability of collision between the crossing pair. Figures 3 and 5 are coloured based on the collision and crossing pair for each system. As first crossings can only ever occur between neighbouring planets, it is possible to use only three colours for this: blue for same-pair systems whereby the same pair was involved in both the first orbit crossing and the collision, orange for otherneighbouring-pair systems to indicate that the colliding pair was the neighbouring pair that did not first cross, and green for extrema-pair systems to indicate a collision between the inner and outer planets. Across the whole range of it can be seen that for systems with a below the model fit line collisions are predominantly between the first crossing pair. Figure 8 shows how these three combinations of events are distributed over time. Collisions that take place within a single orbit of orbit crossing are almost exclusively found in samepair systems due either to simple immediate collisions or to the temporary gravitational bounding of planets discussed previously. Same-pair collisions are the most likely outcome for all systems at 1 AU, shown in the top pane, until ≈ 10 4 orbits, followed by otherneighbouring-pair systems, with extrema-pair systems being the least likely. However, after this period the probability of collision between any combination of planets becomes almost identical, indicating that the mixing of planetary orbits after crossing is sufficient to overcome the increased probability of same-pair integrations due to the initial orbital configuration. Interestingly, the peak of other-neighbouringpair and extrema-pair systems do not align, instead the former peaks first. This can be understood as the mixing process taking longer to cause the inner and outer planets orbits to overlap than to excite the middle planet enough to cross the orbits of both of its neighbours. In the bottom pane, it can be seen that at 0.25 AU the behaviour is similar; however, the number of collisions taking place within a single orbit roughly doubles.

INCLINED INTEGRATION SUITE
In the co-planar case, no systems survived for more than a million orbits after the first orbital crossing. However, Rice et al. (2018) observed a number of non co-planar systems that survived for their maximum simulation time of ten million orbits. Therefore, we now go on to examine the behaviour in the non co-planar case described by the inclined suite of initial conditions in section 2.5. As a reminder, these initial conditions include fifteen initial inclinations ranging from an initial orbital height of 0.10 to .

Dynamic heating
The systems studied in the inclined integration suite begin with modest inclinations and no eccentricities, making them dynamically cold. Figure 9 shows how the system heats up over time by plotting the root-mean-square (RMS) inclination and eccentricity over time. We calculate the mean over all runs that have experienced an orbital crossing, and fit a linear model to this mean which is shown as the solid green line. Individual integrations are shown in purple until they experience an orbital crossing and in grey thereafter. For clarity, in Fig. 9 results of individual integrations are only shown for eighty integrations in the inclined suite for = 5.98. Rice et al. (2018) found that, for four-planet Neptune-mass systems, there are two distinct growth modes of RMS eccentricity before and after an instability event: Eccentricity evolves rapidly to a quasiequilibrium at a value of 10 −2 at which point encounters begin. After a period of mixing as a result of close approaches, systems transition into a new evolutionary phase during which eccentricity growth follows a power-law form approximately ∝ 1 /6 . In the threeplanet Earth-mass case, our systems reach a quasi-equilibrium value of ≈ 10 −3 before a period of chaotic mixing and rapid growth, which finally settles into the new growth phase approximately ∝ 1 /6 . The RMS inclinations in Fig. 9, on the other hand, while similar are different to the four-planet Neptune-mass case. We also observe that the inclination of the systems remains at roughly the initial value until the first encounter, at which point they are rapidly excited before entering a new growth mode. This rapid excitation is in keeping with the findings of Matsumoto & Kokubo (2017). These behaviours can be seen by the horizontal inclination lines in the population of systems before crossing and in the power-law growth in the population afterward. Rice et al. (2018) stated that the trend towards long-lived systems depends upon only the RMS inclination being greater than the averaged ratio of Hill radius to semi-major axis, this is called the critical inclination and is marked on this plot. We also find this to be the case across all systems within our inclined suite: any systems that have experienced orbital crossing and have their RMS inclination damped below this threshold rapidly experience a collision. The key difference in results in our simulations as compared to the fourplanet, Neptune-mass case is that the power-law growth rate appears to be ∝ 1 /4 as opposed to ∝ 1 /3 . We offer two possible explanations for this: 1) our data set could be biased due to the non-random ini- . Post-crossing survival time of inclined integration suite with systems at 1 AU. Colours of data points represent the initial inclinations, with darker colours representing higher inclinations. The twenty-three systems that persisted for the full 10 8 orbits are highlighted via a red triangle, independent of their initial inclination. Note that most of these surviving systems had their initial orbital crossing in far less than 10 8 years, so they survived for almost 10 8 years post-crossing before the simulation was terminated and appear as triangles at the top of the plot; the two exceptions, which survived for < 3×10 7 years, both had initial orbital separations > 5.3. tial conditions used; or, 2) there could be an underlying dependence between either the planetary mass or the number of planets within the system and the growth rate. Further investigation is needed to distinguish between these two possibilities. Figure 10 shows the crossing time for systems within our inclined suite. We find a large variance in crossing time across the inclined suite with a difference between the maximum and minimum crossing times at each value of as large as two orders of magnitude in many cases. The spikes seen in Fig. 1 are also present in some of The upper two plots, in cyan, are for systems initially at 1 AU and the lower two plots, in grey, are for 0.25 AU. The two leftmost plots contain data for systems with the minimum initial inclination, 0 = 0.06 • , whereas the two rightmost plots contain data for systems with the maximum initial inclination, 0 = 0.58 • . Two systems survived for the full simulation time after an orbital crossing in the low inclination case at 1 AU whereas one survived in the high inclination case. No systems in the 0.25 AU case survived for the full simulation duration after an orbital crossing in any of our integrations.

Time scale to planet-planet collision
our inclined cases. A model of the type in Eq. 4 is fitted to the mean values of crossing time observed at each value of , yielding coefficients = 1.39 and = 2.18. These values are in very good agreement with those from the standard suite. This is, however, where the similarities between the co-planar and inclined cases end. Figure 11 shows the post-crossing survival time for systems within the inclined suite, the times are much higher than in the co-planar case where the longest surviving system after crossing survived for roughly one million orbits. Here, the majority of systems survive for longer than this and, in fact, there are twenty-three systems that do not experience any collision at all within the maximum simulation time (100 million orbits), equivalent to 0.14% of all integrations. Given that the post-crossing survival time is now approaching one percent of the lifetime of the Sun, it is much less unlikely that we actually could observe an inclined system between a crossing and a collision. However, at 0.25 AU no integrations survived for the full simulation duration after an integration. Figure 12 shows the probability of a collision across all integrations within the inclined suite. The probability is calculated as the cumulative fraction of systems that have experienced collisions divided by the total number of systems. Results are included for systems initially at 1 AU as well as at 0.25 AU. Decreasing the initial distance to the star by this amount is identical to having artificially inflated the planetary radius R by a factor of four, i.e., made R approxi- log 10 (minimum miss distance) (r ) Figure 17. Time between closest encounter prior to impact and impact against the time-averaged inclination range, i.e. the difference between the smallest and largest inclinations, for systems at 1 AU in the inclined integration suite. The closest encounter experienced by a system is indicated through colouring.
mately equal to that of Neptune, whilst keeping the innermost planet initially at 1 AU. It is therefore expected that the collision probability over time should increase with decreasing initial distance to the star. However, the increase is striking: for Earth analogues the probability of a collision for a given system after one million orbital periods is roughly 50%, but for a Neptune radius (1 M ⊕ ) planet at 1 AU that probability increases to over 75% across all ranges, reaching almost 90% in all but one range. Furthermore, for the 1 AU systems it can be seen that the various regions converge after roughly a million orbits. This indicates that the evolution after the first close encounter has reconfigured the system such that any prior collision probabilities due to initial orbital spacing are lost. To understand this, we can look at the collision probability in Figure 12 at one million orbits for 0.25 AU systems. These systems are equivalent to to a Neptune radius planet being placed at 1 AU and roughly 90% have experienced a collision within this timescale. We can therefore infer that the same roughly 90% of Earth radius planets at 1 AU must have experienced a log 10 (minimum miss distance) (r ) Figure 18. Time between closest encounter prior to impact and impact against the time-averaged maximum eccentricity for systems at 1 AU in the inclined integration suite. The closest encounter experienced by a system is indicated through colouring.
close encounter within 4R . The loss of prior collision probabilities due to orbital spacing after this point in time therefore appears to be driven by these particularly close encounters. Figure 13 contains the distribution of post-crossing survival times for two subsets of the inclined suite results: the subsets each contain 1120 configurations, one at the minimum initial inclination (0.06 • ) and the other at the maximum initial inclination (0.58 • ). It can be seen that the distributions are different at each initial orbital radii and inclination. Firstly, the population of collisions taking place within several orbits of an orbital crossing decreases with increasing initial inclination. In both of the most highly inclined cases there is only a single peak present in the distribution; however, this distribution is much more negatively skewed in systems initially at 1 AU. In the lowest inclination cases there are two peaks present in addition to the one caused by immediate collisions. One peak is collocated with those found in the more inclined case. The second peak is centered at approximately = 10 2.5 . In the co-planar case we have seen that the distribution of post-crossing survival times are centered at approximately 10 2.5 orbits and it is also known that if the inclination is below the critical threshold = the number of collisions occurring within a factor of three of the orbital crossing increases (Rice et al. 2018). Both of these factors combined explain the appearance of this second peak. Additionally, a larger proportion of systems at 0.25 AU experience a collision in this second peak.
The effect of increased initial inclination across the whole inclined integration suite can be seen in Fig. 14, where an increase in inclination, shown here in terms of orbital height, leads to a moderate increase in the median post-crossing survival times for systems at both 0.25 AU and 1 AU. The RMS inclination in compact three-body systems has been seen to stay approximately constant up until the time of the first close encounter, which means that observed inclinations of actual planetary systems could in fact provide information about the probable survival times of systems after an orbital crossing.
The parameter that dominates the post-crossing survival time of systems in the inclined suite is the ratio of the planetary radius to the Hill radius at 1 AU. Figure 15 shows the median of the log postcrossing survival times for all systems in the suite at 1 AU. We find almost two orders of magnitude difference in the average survival time of systems with planets where R / = 0.017 as compared to systems with planets where R / = 0.004. This outweighs the effect of initial inclination on the survival times. Interestingly, systems surviving for the full 10 8 orbits can be seen all the way down to a value of R / = 0.0157 where a rapid decrease in the lifetime of the longest lived of systems is seen. This is equivalent to a planet initially located at 1 AU with a radius 3.5 times that of Earth.
In addition to the dependence of the post-crossing survival time upon the orbital elements of the system, we also find a correlation with the distance of the closest approach. Figure 16 shows the time taken for a collision to occur after the closest encounter experienced prior to it, at time denoted , against the distance between the surfaces of the planets. Data points in the shaded grey area are excluded from our fitted models, and this area corresponds to the boundary seen in Figure 3 at approximately eight orbits. Here, we see a strong negative correlation where a least squares model fitted to the log of − and the miss distance of the encounter has a slope of −0.26 with a -intercept of 1.6. Ergo, the closer an encounter experienced by a system the longer it is likely to survive afterwards. In this plot, each point is also coloured according to the post-crossing survival time of the system. Looking vertically from top to bottom at the colouring it can also be seen that the absolute post-crossing survival time of systems depends upon the miss distance of the closest encounter. It seems that for planetary systems to survive for a long time after an orbital crossing they must risk collision.
We find that the closest encounters are responsible for driving the largest changes in both inclination and eccentricity, and we believe that it is the increase in inclination that causes the trend seen in Figure 16. Figure 17 shows the time taken for a collision to occur after the closest encounter experienced prior to it against the timeaveraged inclination range, i.e. the difference between the maximum and minimum inclinations. Systems with the largest inclination range survive for the longest after a close encounter and the minimum miss distance, indicated through colouring, is key to increasing this range. Figure 18 is identical except it shows the time-averaged maximum eccentricity in a system. Again, the minimum miss distance can be seen to be responsible for the increases in eccentricity. These increases in eccentricity will also work to increase the lifetime of systems through a reduction in the effect of gravitational focusing on the combined physical/gravitational cross-sectional area of planets (Safronov 1972). Figure 19 is the equivalent to Fig. 8 but for the inclined suite. Similarly to the co-planar case, we find that collisions within a single orbit are almost exclusively between the same pair involved in the crossing. We also find an increase in the number of collisions within this time frame in the 0.25 AU case compared to the 1 AU case. However, at a factor of approximately 3, here we see that the increase is more substantial. The distributions of survival times for systems surviving after crossing for longer than a single orbit appear very different to the co-planar case. Nonetheless, some similarities in behaviour are present: in both the co-planar and inclined case there is a peak present of same-pair collisions between 10 2 and 10 3 orbits. Adjusting for the number of systems in each suite we find that the fraction of systems colliding at this point is roughly five times smaller at 1 AU in the inclined case. The time period for mixing in the inclined case is approximately 10 4 orbits, slightly longer than in the co-planar case, after which collisions between any pair of planets become equally likely.

CONCLUSIONS
We performed more than 25, 000 integrations of compact threeplanet systems with the TES integration tool for a maximum time of 10 9 orbits of the innermost planet or until the first collision of planets. We chose to focus our attention on the effects of orbital spacing and therefore distributed system configurations across a wide range of initial values evenly spaced in . Efforts were initially focused on the co-planar case where it is easier to isolate the effects of increasing but then extended to include the inclined case as well. We find in the co-planar suite that planetary systems are doomed after an orbital crossing: they rapidly experience a collision within a maximum observed time of less than one million orbits. However, despite this prognosis, we found that systems with a wider initial spacing of planets do survive longer, exhibiting a median post-crossing survival time following a slope log 10 ( ) ∝ 0.12 . Additionally, we show that three distinct populations of post-instability impact behaviour are present, with very few outliers: (i) immediate collisions within a tenth of an orbit, (ii) prompt collisions between a tenth of an orbit and two orbits, (iii) those surviving for much longer than ten orbits.
The pathology of these different behaviours have been identified and each of them are also observed in the inclined suite.
The probabilities of a collision between specified planetary pairs were also calculated and it was found that collisions will take place between the same pair of planets that initially crossed in the majority of cases, ranging from 48% to 62% depending on the region of . These probabilities increase further depending on the radius of the planet with Neptune-radius planets experiencing probabilities as high as 76%. Despite this increase in probabilities in the co-planar case, the post-crossing survival time only weakly depends upon the planetary radius causing an increase of only 10 3 orbits. In the inclined suite, however, we observe that the planetary radius is the main driver of the post-crossing survival time. We find a decrease in median postcrossing survival time of almost two orders of magnitude between Earth and Neptune radius planets. Additionally, the initial orbital inclinations have been shown to also influence the post-crossing survival times across the full range of by as much as an order of magnitude.
Additionally, we looked at the RMS eccentricity and inclination growth of all systems within our inclined suite after an orbital crossing. Here, we replicate the eccentricity growth rate ∝ 1 /6 found in other studies. We do, however, find the growth rate of the inclination to be ∝ 1 /4 instead of the ∝ 1 /3 observed in previous work.
Finally, we have shown that systems that experience the closest encounters also survive for the longest, and planetary systems that wish to survive must therefore live dangerously.

APPENDIX A: INTEGRATOR COMPARISON
As TES is a new scheme, we have chosen to perform additional integrations making use of IAS15 such that crossing and collision times can be compared. We performed integrations using IAS15 for all runs in our standard integration suite over the range = 3.5 to 6.3 using the standard configuration where the tolerance parameter = 10 −9 . Additionally, we also perform identical comparisons with the dataset of crossing times in Lissauer & Gavino (2021) obtained using the MVS implementation within MERCURY. Throughout this section, TES and IAS15 check for an orbital crossing on every integration step whereas the MVS and hybrid schemes check for an orbital crossing once every ten years. Table A1 compares the crossing time obtained by TES and IAS15; we find very good agreement in results especially in lower ranges where the lifetime of systems is shorter. In particular, for systems where < 4.0 we find that 68% of systems experience an orbital crossing within 1% of one another, increasing to 79% if the tolerance is relaxed to being within 10% of one another. These percentages can be compared to the data presented in Table 3 comparing the performance of TES under the influence of an initial 100 m perturbation in the position of the innermost planet along its orbital arc. Comparison of the summary columns in these two tables reveals that the difference in crossing time between TES and IAS15 is smaller than the effect of the 100 m perturbation. Tables A2 and A3 compare the crossing times found in TES and IAS15, respectively, against those obtained with MERCURY. Unsurprisingly, we find that the comparison yields very similar statistics in both cases. In particular, the runs finishing within 1% and 10% of each other drop by at least 54% in the lowest region of , although the reduction is much less pronounced at higher . We also find that the difference in crossing times between TES and IAS15 at higher is very similar to that of TES and IAS15, and MVS. Finally, Table A4 compares the collision times for TES and IAS15, and this is where we find the largest differences between the two schemes. In the summary column, over the entire range, we see that the number of runs finishing within 1%, 10% and a factor of two have decreased substantially when compared to Table A1 performing the same comparison but for crossing times. The majority of the reduction in these statistics comes from the integrations where < 5 and the evolution after crossing is still a substantial fraction of the overall simulation time. These differences highlight the sensitivity of integrations to close encounters.
In addition to the quantitative comparisons between integrators thus far we have also encountered some qualitative differences in behaviour between the schemes examined so far and the hybrid integration scheme within MERCURY. Figure A1 visualises the comparison between integration schemes found in Table A1 to A3. It shows how tightly the crossing times for the three schemes are clustered to one another, and in particular it shows that in the region located around = 5.7, where a region of high stability is found, the schemes all perform identically capturing the long-lived system behaviour. The hybrid integration scheme within MERCURY combines the MVS scheme with a non-symplectic integrator to allow for close approaches to be handled, and it would therefore seem like an ideal candidate for performing all of the experiments in this article. Figure A2 shows the results of an experiment identical to that described in Section 2.3 except with initial longitudes of = [0, 10.17, 20.33] • . In this experiment, the MVS scheme has a density of one thousand integrations per unit whereas the hybrid scheme has a reduced density of one hundred per unit . Both the MVS and hybrid schemes use what is considered a conservative step size of 18 days resulting in slightly over 20 steps per orbit. For the most part, the schemes agree well with each other; however, a key difference between the schemes can be seen in the region = 6 where no integrations performed by the hybrid scheme lasted for longer than 10 8 orbits despite the majority of MVS scheme integrations lasting for 10 10 orbits. It appears that the hybrid scheme is not accurate enough properly to capture the dynamics in this region which has led to a population of short-lived systems that in fact should have been stable for a lot longer. This was a strong motivation for not using the hybrid scheme for experiments in this work instead opting to use TES which, as can be seen in Fig. A1, can properly capture the dynamics in regions of increased stability. This paper has been typeset from a T E X/L A T E X file prepared by the author.