© The Authors 2023

1 Introduction

Simulations of the mechanical motion of many bodies are generally computationally expensive. Consider the N-body problem. Here, each of the N particles moves under the influence of the gravitation of all other particles and each particle can collide with any other particle. Therefore, there are N² interactions to account for. Various methods have been invented to speed up the simulation or increase the particle count: (i) direct N-body simulations with dynamic time-steps (see Dehnen & Read 2011, for an overview); (ii) the octree code for collision detection (Bentley 1975; Meagher 1980); (iii) the Barnes-Hut algorithm (Barnes & Hut 1986; Barnes 1990; Hamada et al. 2009; Burtscher & Pingali 2011) for mutual gravity, where nearby particles are grouped so that their effect on a distant particle can be combined, which requires O(N log N) computational steps; (iv) the fast multi-pole Greengard and Rokhlin method (FMM), where higher order moments of the particle groups are included (Rokhlin 1985; Greengard 1990); (v) parallelization of these methods (Warren & Salmon 1993); (vi) particle mesh methods, where the N force vectors are calculated using the Newton potential and the Poisson Equation for the potential is solved numerically with fast Fourier transforms (Bodenheimer et al. 2007); (vii) the finite-elements method (FEM); and finally, (viii) for a Keplerian system (with a large central mass), where the particles move in slowly precessing Kepler ellipses described by the Laplace-Lagrange equations for the orbital elements (Murray & Dermott 2009); because the Kepler ellipses change slowly over time, the time-step in these numerical integration propagators can be many orbital periods.

In astronomy, collision detection is the problem of finding the precise moment at which asteroids, planets, or satellites collide. Here the difficulty in predicting collisions, or calculating the collision probability, stems from the fact that the objects are very small compared to the size of their orbits. In numerical simulations the number of nearby particles that need to be considered is a function of the step size. Because the particles travel during each time-step, the volume of space around the particle that needs to be probed for collision partners has a radius of the time-step times velocity. In order to limit the number of collision partners, this volume needs to remain small. Therefore, the step size decreases as O(N^−1/3) and the total number of steps for a fixed simulation time grows as O(N^1/3). Efficient codes, such as octree codes (Meagher 1982) or spatial hashing codes, have an algorithmic efficiency of O(N log N) per time-step. If these are used in collision detection, the number of steps grows as O(N^4/3 log N). This makes the problem of collision detection in astronomy even more challenging than pure gravitational evolution without collisions (see Dehnen & Read 2011, for a comparison between codes with and without collision detection).

In this paper, we apply collision detection to Keplerian systems, such as astrophysical disks, where all particles feel a dominant gravity force from one heavy central mass. Each particle is in a Kepler orbit given by parameters a, ϵ, ϖ, I, Ω, and v (see Table 1 for the symbols). However, the advantage of implementing collision detection in an analytic propagator is that the algorithm can be very efficient. In Sect. 2, we analyze the timescales, evaluate the numerical efficiency, and compare it with algorithms based on tree code. As there is no numerical integration of the orbits, many physical effects are neglected (see Sect. 2.3).

We describe (in Sect. 3) the algorithm for an N-body code with collision detection in detail. Initially, it compares particles sorted by radial distance using equations from the seminal paper by Öpik (1951); see Fig. 1. This is effectively an implementation of the apoapsis/periapsis filter of Hoots et al. (1984) and an example of a sweep and prune method. The algorithm then uses analytic evaluation of the points of collision (near the nodal line in Fig. 2, following Hoots et al. 1984; Manley et al. 1998, as explained in Sect. 4), of the earliest crossing time (Sect. 5), and of the time of collision (derived in Sect. 6). The algorithm keeps track of pairs of particles that are on a collision course, from the earliest to the latest moment of collision. Each step of the simulation involves only the calculation of the next collision, and updating the list. In the method, time-steps increase with decreasing s, which allows long simulation times. Indeed, for the limiting case s → 0, the algorithm stops after initialization, as it finds that there are no collisions. In contrast, collision detection using a numerical integration propagator always requires a nearest-neighbor search for every particle. The time spent on this search is independent of s.

To our knowledge, the idea of bookkeeping a list of future possible collisions has not been studied elsewhere. The algorithm relies on a novel method to quickly find the exact collision times. We derived new expressions, Eqs. (11) and (12), for the difference in the eccentric anomaly between two given points on an orbit that are also accurate at small eccentricities, when the eccentric anomalies themselves are ill defined. These formulas were needed to calculate the time for a particle to get to the collision point.

Fig. 1 Two orbits separated by a sphere (dashed). Orbit 1 (blue) has apoapsis a1 + c1 and orbit 2 (purple) has periapsis a2 − c2. Filtering out such collision-avoiding apoapsis-periapsis pairs is an efficient sweep and prune method.

2 Applications and estimates of timescales

The algorithm is intended for simulation of the dynamical evolution of a large planetary ring system or a debris disk around a star. In the latter, nearby passages also occur, where one planet is inside the sphere of gravitational influence of another. The resulting gravitational scattering may be modeled with an elastic collision. These systems therefore have four different timescales: the orbit time T, the collision time T_coll, the scattering time T_scat, and the secular time T_prec.

In order to estimate how many collisions happen per unit of time, we first consider a thin disk with a completely randomized (homogeneous) distribution of particles in near-circular orbits. A particle with radius s traces a cylindrical volume of size (2πa/T)πs² per unit of time. The disk is a cylinder of radius a and height 2Ia, meaning that the particle density is N/2πIa³. Accounting for the N²/2 pairs, the rate of collisions is estimated to be

$$\frac{1}{T_{\text{coll}}}=\underset{\text{velocity}}{\underbrace{\frac{2\pi a}{T}}}\cdot \underset{\begin{array}{c}\text{effective}\\ \text{cross section}\end{array}}{\underbrace{\pi {\left(2s\right)}^2}}\cdot \underset{\text{density}}{\underbrace{\frac{N}{2\pi I{a}^3}}}\cdot \underset{\text{pairs}}{\underbrace{\frac{N}{2}}}=\frac{2\pi N^2{s}^2}{I{a}^{2}T}.$$

If we want to include close encounters, we may substitute s into the formula for the radius of the sphere of influence s = (m/M)^2/5a. The timescales are therefore

$$\begin{array}{ccc}{T}_{\text{scatt}}=\frac{I}{2\pi N^2}{\left(\frac{M}{m}\right)}^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$5$}\right.}T,& T_{\text{coll}}=\frac{I{a}^2}{2\pi N^2{s}^2}T,& T_{\text{prec}}=\frac{4M}{Nm}T,\end{array}$$

The formula for the precession is taken from Murray & Dermott (2009). If we model the early inner Solar System by N = 10⁶ planetesimals of characteristic size s = 100 km in a disk with a = 4 au and I = .1, we have

$$\begin{array}{ccc}{T}_{\text{scatt}}\approx 45\min ,& T_{\text{coll}}\approx 35\text{\hspace{0.17em}yr},& T_{\text{prec}}\approx {10}^5\text{yr},T\approx 10\text{yr},\end{array}$$

Next, we consider a ring system around a planet. If we assume its radius is only a few times that of the planet and the ring particles have the same density as the planet, the collision time is comparable to the scattering time. For the Uranus ring system, we take N = 10¹³, s = 1m, and a = 10⁵ km, which results in

$$\begin{array}{ccc}{T}_{\text{scatt}}\approx T_{\text{coll}}\approx {10}^{-9}\text{s,}& T_{\text{prec}}\approx 190\text{\hspace{0.17em}d,}& T\approx 1\text{\hspace{0.17em}d}\end{array},$$

The precession is now entirely due to planet oblateness (J₂ = 3 × 10⁻³). We now estimate the deflection angle due to scattering. When the scattering at impact parameter b is integrated over all values, for a path length of 2πa we find that

$$\frac{\text{deflection}}{\text{orbit}}=\frac{N}{\underset{\text{density}}{\underbrace{2\pi I{a}^3}}}\cdot {\displaystyle \underset{0}{\overset{s}{\int }}\frac{2am}{\underset{\text{deflection}}{\underbrace{bM{\in }^2}}}}\cdot \underset{\text{volume\hspace{0.17em}shell}}{\underbrace{2\pi a2\pi b\text{d}b}}=\frac{4\pi N}{I{\in }^2}{\left(\frac{m}{M}\right)}^{\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$5$}\right.}$$

The orbital eccentricity e accounts from the fact that the relative velocity is roughly $\sqrt{GM/aє}$, which becomes small for orbits with the same sense of rotation. Although there are many close encounters where the mutual gravity takes over the central force, the (very crude) estimate of the deflection angle is about 10⁻⁵ and 10⁻⁷ per orbit for the inner Solar System and the Uranus ring system, respectively, mainly due to high relative velocities.

The algorithm maintains a list of all particle pairs that are on a collision course. As any pair can only collide near the line of intersection (in Fig. 2) of the two orbital planes, any random pair has the probability ≈2s/a of being on a collision trajectory. An estimate of the number of pairs on a collision course, or “collision pairs”, is therefore N² s/a. The results of simulations shown in Figs. 3 and 4 validate these estimates.

As the algorithm steps from one collision to the next, the (average) time-step is equal to the collision time T_coll. Although exact precision is already lost in one orbit if corrections for secular motion are not included, we expect collision detection using the Kepler orbits to be able to give reliable statistical results for T/N < T_coll < T_prec.

Fig. 2 Orbits of two planetoids m1, m2 in their orbital planes. Because the bodies are much smaller than the orbits, sj << aj, collisions happen near the mutual nodal line of intersection of the orbital planes, even for small inclinations. The tangent vectors (gray) indicate the linear approximation that may be used to find the collision point.

Fig. 3 Pairs vs. particle radius s for a homogeneous disk with a ≤ amax = 2au, I ≤ 10−3, ϵ ≤ 10−3, and with N = 104, from Aliberti (2022a,b). Orange: twice the number of pairs that needed to be checked in the apoapsis/periapsis filter. Blue: actual number of collision possibilities. The guideline with slope 1 (dashed) shows the approximate linear dependence of the collision pairs on s. The influence of the planet size can be seen in the apoapsis/periapsis filter for large s. Because there are two collision possibilities for each pair, the blue dots lie below or on the orange dots.

2.1 Comparing time-steps

We are interested in comparing this approach with numerical integration propagators with collision detection. In order to detect a collision in numerical integration, one must find nearest neighbors. In the tree code, one uses boxes of volume dx³ of a sufficiently small size to contain only one or a few particles. In one time-step, dt, the change in position should not move the particle too many boxes away from its original position; otherwise, it becomes impossible to select neighbors. Alternatively, in an algorithm that uses a sorted list of the coordinates (so-called spatial hashing codes), a particle coordinate, say x, can overtake the values of other particles when its position changes by dx. In this latter case, the number of particles that one particle overtakes in one step should also remain small in order to limit the number of neighboring particles that need to be inspected. As a result, numerical integration with collision detection requires not only several time-steps for one orbital period (dt ⪅ T) but also spatial steps of the order of the inter-particle distance $\left(\text{d}x\underset{˜}{\underset{˜}{<}}\overline{r_{ij}}\right)$. We may estimate the average distance by assuming a homogeneous distribution of the particles. We then find for a disk with a = 4au, I = .1:

$$\begin{array}{l}\begin{array}{cc}\overline{r_{ij}}=a\Gamma \left(\frac{4}{3}\right){\left(\frac{3I}{2N}\right)}^{1/3}=\frac{2\text{au}}{N^{1/3}},& \text{d}t\underset{˜}{\underset{˜}{<}}\frac{T\Gamma \left(\frac{4}{3}\right)}{2\pi }{\left(\frac{3I}{2N}\right)}^{1/3}=\frac{220\text{d}}{N^{1/3}}.\end{array}\hfill \end{array}$$

Clearly, the time-steps for numerical integration actually need to be quite small.

Fig. 4 Pairs vs. particle count N, as in Fig. 3, but with particle radius s = 2 × 10−3au. Orange: number of pairs that needed to be checked in the apoapsis/periapsis filter. Blue: actual number of collision pairs (that could ultimately collide). These determine the runtime and memory for the creation of the collision list; see Table 2. Solid lines: N(N − 1)/2. Dashed line: guideline with a slope of 2.

2.2 Efficiency of the algorithm

Apart from the advantage of having large time-steps, another benefit of our method is that only the collision products need to be tested for possible future collisions with the set of existing particles. This involves ∝ N computations per collision. However, a list of collision possibilities needs to be maintained. The length of this list is expected to be of the order N² s/a. This list takes up memory and therefore requires careful manipulation. At the creation or the removal of a particle, an average number of 2Ns/a new collision possibilities is added or removed, respectively. Hence, after N/2 collisions the list has mostly changed. This is understandable, because then most particles are replaced by new particles. It also means that for Ns/a ≫ 1, most possibilities do not actually happen. Figure 4 shows the initial list length and Fig. 5 shows how this length changes during the simulation.

Table 2 sums up the efficiencies in the various steps of the algorithm. Figure 6 shows the measured runtime for a large set of simulations with different particle numbers. As we did not include defragmentation but only mergers, Fig. 7 shows the values of particle radius s and particle count N for which the Kepler collision detection is more efficient than the algorithmic efficiency (T_sim/dt)N log N ∝ N^4/3 log N of numerical integration with nearest neighbor search using a tree code.

In order to make another comparison between the methods, consider a Solar System with a fixed amount of material volume. We take a disk with a = 4au, I = .1, and a mass N_m = 3M_Earth and with particles of the density of Earth. We thus have Ns³/a³ = 10⁻¹⁵ The resulting average s is also plotted in Fig. 7. The efficiencies for collision detection for the cases without and with precession are then (N² s/a)(T/T_coll) ∝ N³ and N²T/T_prec ∝ N², respectively. These are shown in Fig. 8, and are also compared with the algorithmic efficiency numerical integration with a tree code. For a Solar System model with these parameter values, the octree code is faster for N > 10⁸.

Fig. 5 Length of the list of collision pairs vs. particle count (same parameters as in Fig. 4). Only mergers are simulated, and consequently the N value decreases by one at every step. The graph shows 64 runs.

Fig. 6 Initialisation- (orange) and simulation runtime (blue) vs. particle count N, as in Fig. 4. Dashed lines are guidelines with slopes of 1, 2, and 3. The initialization time is respectively linear and quadratic in N for N < 1/є and N > 1/є, in accordance with the estimates in Table 2. However, because only mergers were simulated, there were at most N collisions, which results in a runtime of O(N3) instead of O(N4).

Fig. 7 Particle radius s vs. particle count N for a disk system with I ≤.1 and a ≤ 4au. Above the black line, the numerical integration propagator with tree code for collision detection beats the analytic propagator using Kepler orbits (in terms of time efficiency). Precession due to a Jupiter has a small effect above the orange line, where Tcoll = Tprec, and has a negligible effect above the orange dashed line, aTcoll = sTprec. If the mass of three Earths is distributed over equal-sized planetoids (of Earth density), one is constrained to the blue line. On this line, precession is small for N > 103 but only attains the much smaller error of s per orbit for N > 109. Clearly, tree-code is better for these high number densities.

2.3 Applications and neglected effects

We can think of the following applications: (i) gravity assists at planetary flybys of a probe traversing the Solar System: the collision detection scheme calculates the times of passage using the sphere of influence. (ii) The tracking of a comet as its orbit is perturbed at close encounters by the planets. (iii) The rings of Saturn, with contact collisions and/or the gravitational influence of shepherd moons, (iv) Growth by merging planetesimals into protoplanets, in young planetary systems, (v) A simple model to study the Kessler syndrome, where artificial satellites break apart due to impacting space debris.

Because of the approximations, collision detection with Kepler orbits is often inaccurate. However, sometimes accuracy is not the aim, or not even possible due to the chaotic nature of the problem. Instead, we may simply want to find out what could happen, and calculate the probabilities of the various outcomes. The speedup allows sampling of initial states, either by adding many small perturbations to one initial state or by adding one perturbing body with many initial states from a large phase-space volume. Analytic propagation with collision detection based on the Kepler orbits neglects the following effects: (i) orbital precession due to mutual gravity or oblateness of the central body (as discussed in Sect. 6.3), (ii) three-body gravitational scattering, (iii) planetary migration due to interaction with the gas in the protoplanetary disk, (iv) capture of planets in mean-motion resonances, (v) the Kozai mechanism, (vi) moons and binaries, (vii) atmospheric and (viii) tidal drag, (ix) solar wind, and (x) the Poynting-Robertson effect.

Fig. 8 Theoretical algorithmic efficiency: number of computational steps (big O) for one orbital period vs. particle count N. Blue: intersecting Kepler orbits without secular dynamics. Orange: with secular dynamics. The blue dotted line indicates where precession due to secular dynamics cannot be neglected. The black line shows the estimated number of steps for a tree code and/or spatial hashing. The slopes of the lines are 3, 2, and 4/3, respectively.

3 The algorithm

3.1 Initialization

We have a system of planets, or particles orbiting a central mass. The particles are numbered j = 1,…,N. For each particle, we store the following set of variables:

$$\left\{\begin{array}{lllllllll}{t}_j^0,\hfill & a_j,\hfill & c_j,\hfill & s_j,\hfill & m_j,\hfill & r_j^0,\hfill & L_j,\hfill & {\in }_j\hfill & {\omega }_j\hfill \end{array}\right\},$$(1)

Here, t⁰ is the time of creation, α, c are the orbital radius and focal distance, s is the particle radius, m is the particle mass, r⁰ is its position, L its angular momentum vector, є is the eccentricity vector, and ω is the frequency.

An initial state would consist of many particles in nearly circular orbits and therefore with small є. If one draws random numbers for the mean anomaly from the interval [0,2π], the resulting (smoothed) phase-space distribution will become stationary; if the values of ϖ Ω are sampled from [0,2π], the distribution will become axisymmetric; if also cos I is drawn from [0,1], it will become spherically symmetric (see Savransky et al. 2011). To simulate a thin disk, one takes small values for I. Equations (B.6) and (B.7) give the vectors L and є in terms of α and e and the angles I, ϖ , and Ω.

We then sort the particle list by increasing value of periapsis. This will allow the implementation of the apoapsis/periapsis filter (Baraff & Witkin 1992). The next step is to consider all particle pairs, and list the pairs that can collide. For these, we also store the calculated collision time $t_{\left(i,j\right)}^1$. If sufficient memory is available, it is possible to store the parameters in Eq. (1) for the new particle that would be formed after the collision. The soonest collision is at the top of the list.

3.2 Main loop

1. If the collision list is empty, end the simulation.

2. Take the pair (i, j) with the soonest collision: the first in the list.

3. Update the time t to the time $t_{\left(i,j\right)}^1$of the collision.

4. Remove any pair containing i and any pair containing j from the pairs list.

5. Remove the particles and j from the particle list.

6. If the orbit of the new particle intersects the central mass or is unbound, go to the next collision on the list.

7. Create new particle(s) defined by $\left\{t_{\left(i,j\right)}^1,a,c,s,m,r^0,L,є,\omega \right\}$

8. For any new particle, consider the other particles and decide if the pair is on a collision course. If this is the case, calculate the time of the earliest collision.

9. Make a sorted list of the new collision possibilities with a record of the collision time and the pair, soonest collision first.

10. Merge this sorted list with the existing sorted list of collision possibilities into a full list of pairs, sorted by time of the collision event, soonest collision first.

3.3 Determining if a pair is on a collision course

At the initialization, pairs of particles need to be considered for a possible future collision. Also, during the simulation, each time a new particle is created, all existing particles need to be paired with the new particle and considered for a possible future collision. However, as we implement the sweep and prune method, only pairs need to be considered with an overlap in the range of radial motion. The radial coordinate for each particle i ranges over the interval from the periapsis to the apoapsis:

$$\left[a_i-c_i-s_i,a_i+c_i+s_i\right],$$

including an extension s_i of the size of the particle, or with the substitution of its gravitational reach s_i = (m /M)^2/5a_i . It is sufficient to compare each particle i with the particles j = i + 1, j = i + 2,… Because the list is sorted, we have a_i −- c_i − s_i ≤ a_j − c_j − s_j. As long as a_i + c_i + s_i ≥ a_j − c_j − s_j, the intervals for i and j overlap and the pair is a candidate for collision. Once we encounter the first j where a_i + c_i + s_i < a_j − c_j − s_j, there are no more particles that can interact with i and we can go to particle i + 1. If є is the average eccentricity, only a fraction of 2є of the total number of pairs N²/2 need to be checked. The resulting reduced number of checks in our numerical simulations is shown in Figs. 4 and 6.

For a particle i created during the simulation, the selection of pairs is slightly different. Again, it is sufficient to consider only particles j with periapsis smaller than the apoapsis of i. However, this time we have to start at j= 1.

1. Consider a pair, say (i, j) = (1, 2); we assume that a₁ > a₂. First we need to find the minimal orbit intersection distance to decide whether or not there can be a collision. For this, retrieve m₁ L₁ є₁ and m₂, L₂, є₂. Then calculate the direction of the nodal line (see Fig. 2), the semi-latus recti, the intersection points, and the velocities at these points:

$$\begin{array}{l}K=L_1\times L_2,\hfill \\ K=\sqrt{K•K,}\hfill \\ \begin{array}{cc}{\mathcal{l}}_1=\frac{L_1•L_1}{GM{m}_1^2},& {\mathcal{l}}_2=\frac{L_2•L_2}{GM{m}_2^2},\end{array}\hfill \end{array}$$

$$\begin{array}{cc}{r}_1=\frac{K{\mathcal{l}}_1}{\pm K+{\in }_1•K},& r_2=\frac{K{\mathcal{l}}_2}{\pm K+{\in }_2•K},\\ r_1=\sqrt{r_1•r_1,}& r_2=\sqrt{r_2•r_2,}\end{array}$$(2)

$$\begin{array}{cc}{\upsilon }_1=\frac{L_1}{m_1{\mathcal{l}}_1}\times \left({\in }_1+\frac{r_1}{r_1}\right),& {\upsilon }_2=\frac{L_2}{m_2{\mathcal{l}}_2}\times \left({\in }_2+\frac{r_2}{r_2}\right),\end{array}$$(3)

We write G for Newton’s constant. Equation (2) for the intersection points is derived in Sect. 4.1. Equation (3) is Eq. (B.1). The two pairs (r₁, r₂) at the opposite sides are indicated with the plus/minus symbol. The steps that now follow must be performed on both of the pairs.

2. Next, calculate for both particles the (approximate) points on the orbits where the distance is minimal:

$$\begin{array}{c}d=r_2-r_1,\\ w={\upsilon }_1\times {\upsilon }_2,\\ \left(w^2\right)=w•w,\end{array}$$

$$\begin{array}{ll}{r^{\prime }}_1=r_1+\left(d•\frac{{\upsilon }_2\times w}{\left(w^2\right)}\right){\upsilon }_1,\hfill & {r^{\prime }}_2=r_2+\left(d•\frac{{\upsilon }_1\times w}{\left(w^2\right)}\right){\upsilon }_2,\hfill \\ r_1={r^{\prime }}_1,\hfill & r_2={r^{\prime }}_2,\hfill \\ d=r_2-r_1,\hfill & \hfill \\ d=\sqrt{d•d},\hfill & \hfill \end{array}$$(4)

Equation (4) for the collision points r′₁ and r′₂ is derived in Sect. 4.2.

3. Retrieve a₁, s₁ and a₂, s₂. Decide whether or not an interaction can take place.

d < s₁ + s₂   ⇒ 1 and 2 collide.

d < (m₁/M)^2/5a₁ ⇒ 1 can perturb orbit 2.

d < (m₂/M)^2/5a₂ ⇒ 2 can perturb orbit 1.

If there is no contact or perturbation: go to the next pair.

3.4 Deterministic collision time

We now give the steps in the calculation of the exact collision moment.

1. Retrieve $t_1^0$, ω₁, $r_1^0$ and $t_2^0$, ω₂, $r_2^0$ of the particles involved (with ω₁ < ω₂). Calculate the first time a particle passes the crossing point in Fig. 2. These times are denoted by $t_1^1$ ,$t_2^1$.

$$\begin{array}{c}\left({\in }_1^2\right)={\in }_1•{\in }_1,\\ \left({\in }_2^2\right)={\in }_2•{\in }_2,\\ r_1=\sqrt{r_1•r_1,}\\ r_2=\sqrt{r_2•r_2,}\end{array}$$

$$\begin{array}{l}z=\left(\frac{r_1}{a_1}-\frac{\text{i}{r}_1{\upsilon }_1}{a_1^2{\omega }_1}\right)•\left(\frac{r_1^0-{\in }_1{\in }_1•r_1^0}{a_1-\left({\in }_1^2\right)a_1}+{\in }_1\right)+\frac{r_1^0•{\in }_1}{a_1}+\left({\in }_1^2\right),\hfill \\ \begin{array}{cc}\Delta E_1=\arg z,& 0\le \Delta E_1<2\pi ,\end{array}\hfill \\ z=\left(\frac{r_2}{a_2}-\frac{\text{i}{r}_2{\upsilon }_2}{a_2^2{\omega }_2}\right)•\left(\frac{r_2^0-{\in }_2{\in }_2•r_2^0}{a_2-\left({\in }_2^2\right)a_2}+{\in }_2\right)+\frac{r_2^0•{\in }_2}{a_2}+\left({\in }_2^2\right),\hfill \\ \begin{array}{cc}\Delta E_2=\arg z,& 0\le \Delta E_2<2\pi ,\end{array}\hfill \end{array}$$(5)

$$\begin{array}{c}{t}_1^1=t_1^0+\frac{\Delta E_1}{{\omega }_1}-\frac{{\in }_1\times \left(r_1-r_1^0\right)}{1-\left({\in }_1^2\right)}•\frac{L_1}{GM{m}_1},\\ t_2^1=t_2^0+\frac{\Delta E_2}{{\omega }_2}-\frac{{\in }_2\times \left(r_2-r_2^0\right)}{1-\left({\in }_2^2\right)}•\frac{L_2}{GM{m}_2},\end{array}$$(6)

These equations are derived in Sect. 5. Equation (5) is obtained by combining Eqs. (11) and (12). The complex number has unit modulus. Equation (6) for the passage times follows from Eq. (10).

2. Evaluate the small dimensionless parameter.

$$\begin{array}{l}u={\upsilon }_2-{\upsilon }_1,\hfill \\ \delta =\frac{1}{\left|t_1^1-t_2^1\right|}\sqrt{u•u\frac{{\left(s_1+s_2\right)}^2-d^2}{\left(w^2\right)}}\hfill \end{array}$$(7)

3. Next calculate the exact number of periods k that planet 1 makes before colliding with planet 2. The method uses the continued fraction of the ω₂/ω₁. Initialize the loop:

$$\begin{array}{cc}{q}_0=\frac{2\pi }{{\omega }_1\left|t_1^1-t_2^1\right|},& k_0=1,\\ q_1=\frac{2\pi }{{\omega }_2\left|t_1^1-t_2^1\right|},& k_1=0,\end{array}$$

4. Start loop counter at n = 0. The loop creates integer sequences α_n, k_n and the positive sequence q_n. These are the digits, the denominators of the convergents, and the remainders of the continued fraction (we do not need the numerators, denoted l_n). The loop performs the iterations

$$\begin{array}{ll}{\alpha }_{2n}=\left|\frac{q_{2n}}{q_{2n+1}}\right|,\hfill & q_{2n+2}=q_{2n}-{\alpha }_{2n}{q}_{2n+1},\hfill \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}if}\hfill & q_{2n+2}=0\Rightarrow \text{\hspace{1em}next\hspace{0.17em}pair}\hfill \\ \hfill & k_{2n+2}=k_{2n}-{\alpha }_{2n}{k}_{2n+1},\hfill \\ {\alpha }_{2n+1}=\left|\frac{q_{2n+1}}{q_{2n+2}}\right|,\hfill & q_{2n+3}=q_{2n+1}-{\alpha }_{2n+1}{q}_{2n+2},\hfill \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} if}\hfill & q_{2n+3}=0\Rightarrow \text{\hspace{1em}next\hspace{0.17em}pair}\hfill \\ \hfill & k_{2n+3}=k_{2n+1}-{\alpha }_{2n}{k}_{2n+2},\hfill \end{array}$$

We use the notation $\lfloor \cdot \rfloor$ and $\lceil \cdot \rceil$ for the floor and the ceiling function. The time T_sim to be simulated sets an upper bound for the solution:

$$x_{\max }=\frac{\left(T_{\text{sim}}-t_1^1\right)\frac{{\omega }_1}{2\pi }{q}_{2n+1}-\left(1+\delta \right)k_{2n+1}}{k_{2n}{q}_{2n+1}-q_{2n}{k}_{2n+1}}$$

We then test the points with coordinates $x=\lceil \frac{1-\delta }{q_{2n}}\rceil ,\dots ,\lfloor \frac{1+\delta }{q_{2n+2}}\rfloor$

if x>x_max⇒ next pair $y=\max \left(0,\lceil \frac{q_{2n}x-1-\delta }{q_{2n}{}_{+1}}\rceil \right)$

If for any of these points 1 − δ < xq_2n − yq_2n+1, then we have found the solution. If not, we increase n by 1 and check this next.

5. The solution is

$$\begin{array}{cc}k=x{k}_{2n}-y{k}_{2n+1},& t^0=t_1^1+\frac{2\pi k}{{\omega }_1},\end{array}$$(8)

Equation (8) for the deterministic collision time is derived in Sect. 6.1.

The solution k = 0 means that the pair (1, 2) is about to collide within the present period. If a gravitational scattering between the same pair (1 , 2) has happened at the previous time-step, the solution k = 0 corresponds to the scattering that was just simulated, and therefore is invalid. However, for a different pair where one particle participated in this last scattering, the solution k = 0 is valid, as it implies an immediate collision with a third particle. For three (or more) bodies inside each others sphere of influence, the multi-body gravitational scattering will therefore be treated as three (or more) successive two-body interactions.

3.5 Stochastic collision time

Although the time of collision between two bodies is uniquely determined by the initial conditions, it is highly sensitive to the precise values of the creation times and the orbital periods. Therefore, if the physical or numerical error in the (initial) values is bigger than s/a, the collision time becomes unpredictable.

If the collision process is assumed to be stochastic, one may adopt the following Monte Carlo method. First, a random number ξ is drawn from the interval [0, 1]. The moment of collision is calculated with

$$t^0=t_1^1+\frac{\left(-\log \xi \right){\left(2\pi \right)}^2}{2GM}\sqrt{\frac{\left(w^2\right)}{u•u}\frac{a_1^3{a}_2^3}{{\left(s_1+s_2\right)}^2-d^2},}$$(9)

Equation (9) is derived in Sect. 6.2.

3.6 The collision

In order to simulate the collision or scattering event, a physical model of the merger or break-up of the particles needs to be implemented. In the case of pure gravitational scattering (close passage), one can use the formulas from Appendix C for the momentum exchange. This elastic collision is depicted in Fig. 9. We now outline how to find the orbit of the new particle in case of a merger between particles 1 and 2.

1. For a simple merger, the new particle has radius, mass, position, velocity, and angular momentum¹ calculated from the basic conservation laws:

$s={(s_1^3+s_2^3)}^{\frac{1}{3}}$ (material volume),

m = m₁ + m₂ (mass),

$r=\frac{m_1{r}_1+m_2{r}_2}{m}$ (center-of-mass motion),

$\upsilon =\frac{m_1{\upsilon }_1+m_2{\upsilon }_2}{m}$ (momentum),

L = mr × υ (angular momentum).

2. Decide whether or not the new particle collides with the central body. We suppose that it is a sphere of radius S.

$$\begin{array}{c}\mathcal{l}=\frac{L•L}{GM{m}^2},\\ r=\sqrt{r•r},\\ \in =\frac{\upsilon \times L}{GMm}-\frac{r}{r},\\ \in =\sqrt{\in •\in },\end{array}$$

In the general case that several collision fragments are created in the collision, there are five cases, with three outcomes (see Fig. 10) for a fragment. First, consider the cases where there is no crossing:

ℓ> (1 + є)S and є < 1⇒ new particle stays.

ℓ> (1 + є)S and є > 1⇒ particle escapes.

In the remaining cases, ℓ ≤ (1 + є)S and the orbit crosses the central body.

є ≥1 and r•υ > 0 ⇒particle escapes.

є ≥ 1 and r • υ ≤ 0⇒collides with M.

є<1     ⇒collides with M.

Only in the first case does the particle stay in a bound orbit, and we continue; otherwise, the particle is removed.

For a merger, there is only one collision product, and the newly formed particle will always have є < 1 because (total) energy can only decrease.

3. Calculate the required orbital parameters of the new particle:

$$\begin{array}{ccc}a=\frac{\mathcal{l}}{1-{\in }^2},& c=a\in ,& \omega =\sqrt{\frac{GM}{a^3}.}\end{array}$$

4. Register {t⁰, a, c, s, m, r, L, є, ω} of the new particle.

5. Sort the list of new particles by increasing a − c − s.

6. Merge these new-particle lists with the existing-particle list.

Go to the next time-step.

Fig. 9 Gravitational scattering between particles 1 and 2 in the center-of-mass frame. The initial velocities υ1, υ2 are scattered in the directions ${{\upsilon }^{\prime }}_1,{{\upsilon }^{\prime }}_2$, ${{\upsilon }^{\prime }}_1,{{\upsilon }^{\prime }}_2$ along the asymptotic lines (dotted). The actual orbits are hyperbolas (blue, purple). The focal points (black dots), with O being the common focal point of the orbits, lie on the orange circles. The four points of intersection of a circle with the asymptotes form rectangles with dimensions of the transverse axis 2aj by the conjugate axis 2bj (Adams & Essex 2021). The impact parameter is the sum b = b1 + b2.

Fig. 10 Five cases for the obit of a collision fragment. (a) The particle is bound. (b), (c) The particle collides with the central mass. (d), (e) The particle leaves the system. For long intervals between collisions (Tcoll ≫ T), cases ((b)≫(e)) can be dealt with by removing the particle. Otherwise, the fragments need to be stored and paired in a similar manner to the bound particles.

4 Calculating points of closest approach

In this section, we derive an approximation for the points on two Kepler orbits with minimal separation. This distance is called the minimum orbit intersection distance (MOID). As we are considering possible collisions between planets, we are interested in the case where the MOID for the orbits of planet 1 and planet 2 is less than s₁ + s₂. We assume that the planet radii s_j are many orders of magnitude smaller than the orbital radii a_j. Then, for an angle I between the orbital planes of larger than (s₁ s₂)/(a₁ a₂), the MOID will be close to the line of intersection of the orbital planes (Hoots et al. 1984; Manley et al. 1998). This principle is illustrated in Fig. 2: outside a cylinder with radius (s₁ + s₂)/$\left(s_1+s_2\right)/\left(2\sin \frac{I}{2}\right)$ about the “mutual nodal line”, all points in orbit 1 are separated by more than s₁ + s₂ from points in orbit 2. Any collision must therefore happen inside the cylinder. The cylinder is only large for very small inclinations. If the system is a disk with an average inclination angle $\overline{I}$, these small inclinations are rare for $N<\overline{I}{a}^2/s^2$.

Because the range for gravitational scattering can be larger, our approach will only work for low-mass planets. The iterative scheme converging to the MOID that projects the points onto the orbit followed by linearization is described in Hoots et al. (1984). Various other methods to obtain the MOID have been found (see e.g. Gronchi 2005; Milisavljević 2010; Segan et al. 2011; Wiźniowski & Rickman 2013; Hedo et al. 2018).

Fig. 11 Elliptical Kepler orbit (a) and orbit squeezed into a circle (b) from scaling the x-axis by b/a. The travel time from r0 to r1 is equal to the ratio of the swept (cyan) area in (a) to πab times the period. In (b), this ratio is the (yellow) circle segment plus the (gray) triangle minus the image of the (gray) triangle in (a) over πb2. The circle segment has area ΔEb2/2, with ΔE = E1 − E0 being the difference in eccentric anomaly.

4.1 Intersecting orbit 1 with orbital plane 2

The Kepler orbit of a planet is entirely determined by its angular momentum L and eccentricity vector є (Goldstein 1964). The angular momentum is normal to the place of the orbit and the vector aє points from the center of the ellipse to the focal point where the central mass is located (see Fig. 11). Now let us consider two orbits, for planet 1 and planet 2, specified by L₁, є₁ and L₂, є₂, respectively. The line of intersection of the two orbital planes, or the nodal line, can be found as follows. Because the angular momenta L₁ and L₂ are both perpendicular to the nodal line, a direction vector of the nodal line is

$$K=K\widehat{K}=\pm L_1\times L_2.$$

The plus/minus symbol indicates the two opposite directions in which the intersection points with an orbit are found. Because the eccentricity vector є of an orbit points from the central mass towards the periapsis, the true anomaly v of the intersection point in the direction K is given by

$$\cos v=\frac{\in •\widehat{K}}{\in }.$$

The point of intersection can now be found from the formula for the orbit Eq. (B.3). We find

$$r=r\widehat{K}=\frac{\left(1-{\in }^2\right)a}{1+\in \cos v}\widehat{K}=\frac{\left(1-{\in }^2\right)a}{1+\in •\widehat{K}}\widehat{K}=\frac{\left(1-{\in }^2\right)a}{K+\in •K}K.$$

Therefore, for the two pairs of intersection points, we obtain

$$\begin{array}{c}{r}_{1,\pm }=\frac{\left(1-{\in }_1^2\right)a_1{L}_1\times L_2}{\pm \left|L_1\times L_2\right|+{\in }_1•\left(L_1\times L_2\right)},\\ r_{2,\pm }=\frac{\left(1-{\in }_2^2\right)a_2{L}_1\times L_2}{\pm \left|L_1\times L_2\right|+{\in }_2•\left(L_1\times L_2\right)},\end{array}$$

4.2 Pair of closest points between two orbits

Next, we approximate the points where the MOID is found. In order to do so, we consider the tangent lines of the orbits at the points r₁ and r₂ of intersection with the nodal line. The tangent lines point in the direction of the velocities υ₁ and υ₂ at r₁ and r₂. These can be found using Eq. (B.1). The distance between the two lines is given by

$$d=\left|\left(r_2-r_1\right)•\frac{\left({\upsilon }_1\times {\upsilon }_2\right)}{\left|{\upsilon }_1\times {\upsilon }_2\right|}\right|.$$

This is the projection of the difference vector onto the direction of shortest distance. We refer to the two points on the lines where the distance is minimal as ${r^{\prime }}_1$ and ${r^{\prime }}_2$. These positions are given by

$$\begin{array}{c}{r^{\prime }}_1=r_1+\left[\left(r_2-r_1\right)•\frac{{\left|{\upsilon }_2\right|}^2{\upsilon }_1-\left({\upsilon }_1•{\upsilon }_2\right){\upsilon }_2}{{\left|{\upsilon }_1\right|}^2{\left|{\upsilon }_2\right|}^2-{\left({\upsilon }_1•{\upsilon }_2\right)}^2}\right]{\upsilon }_1,\\ =r_1+\left[\left(r_2-r_1\right)•\frac{{\upsilon }_2\times \left({\upsilon }_1\times {\upsilon }_2\right)}{{\left|{\upsilon }_1\times {\upsilon }_2\right|}^2}\right]{\upsilon }_1,\\ {r^{\prime }}_2=r_2+\left[\left(r_2-r_1\right)•\frac{{\upsilon }_1\times \left({\upsilon }_1\times {\upsilon }_2\right)}{{\left|{\upsilon }_1\times {\upsilon }_2\right|}^2}\right]{\upsilon }_2,\end{array}$$

One may verify that $\left({r^{\prime }}_2-{r^{\prime }}_1\right)•{\upsilon }_1=\left({r^{\prime }}_2-{r^{\prime }}_1\right)•{\upsilon }_2=0.$. This proves that the minimal distance between the lines is realized at the points ${r^{\prime }}_1$ and ${r^{\prime }}_2$. We also have that $\left({r^{\prime }}_2-{r^{\prime }}_1\right)•\left(r_2-r_1\right)=d^2$, implying that $\left|{r^{\prime }}_2-{r^{\prime }}_1\right|=d<\left|r_2-r_1\right|.$

If the inclination between the orbital planes, which is given by I ≈ K/L₁L₂, is not much larger than (s₁ + s₂)/(a₁ + a₂), the linear approximation is inaccurate. One can improve this approximation by reducing the lengths of ${r^{\prime }}_1$, so that the points lie on the respective orbits, and then finding the shortest distance between the tangent lines to the orbits at these new points. Here, one may iterate as in the method of Newton Raphson.

5 Calculating the earliest passage of the crossing

In this section, we derive expressions for the time it takes a particle on a Kepler orbit to go from r⁰ to r¹ . In the algorithm, r⁰ is the particle’s creation point and r¹ is the collision point. A standard approach is to use the eccentric anomaly E values at these points. However, E becomes ill defined for small eccentricities². Here, we derive the exact expression, Eq. (6), that does not depend on the value of E; only the difference between the values of E enters the derivation.

Kepler’s second law states that the position vector sweeps out equal areas in equal times. The swept area can be decomposed into a segment of the ellipse (with an angle determined by the eccentric anomaly) plus the triangle between −aє, 0, and r⁰ minus the triangle between −aє, 0, and r¹ in Fig. 11a. Therefore,

$$\text{area}=\underset{\text{ellipse\hspace{0.17em}segment}}{\underbrace{\frac{\left(E^1-E^0\right)ab}{2}}}+\underset{\text{added\hspace{0.17em}triangle}}{\underbrace{\frac{a\in \times r^0}{2}•\widehat{L}}}-\underset{\text{substracted\hspace{0.17em}triangle}}{\underbrace{\frac{a\in \times r^1}{2}•\widehat{L}}}.$$

Kepler’s second law therefore implies that the time it takes a body to move from r⁰ to r¹ is equal to

$$t^1-t^0=\frac{\text{area}}{\pi ab}T=\frac{2\text{\hspace{0.17em}area}}{wab}=\frac{E^1-E^0}{w}-\frac{є\times \left(r^1-r^0\right)}{wb}•\widehat{L}.$$(10)

As the ellipse has its major axis pointing in the direction ofє, we can transform the elliptic orbit into the circle in Fig. 11b with the linear transformation

$$r\stackrel{\text{Squeeze}}{\to }\left(1-\frac{\in {\in }^{\text{T}}}{{\in }^2}+\frac{b}{a}\frac{\in {\in }^{\text{T}}}{{\in }^2}\right)r=r-\frac{a\in •r}{a+b}\in .$$

This squeezes the semimajor axes by a factor b/a, and leaves the semiminor axes intact. The eccentric anomaly is defined with respect to the center of the circle. We therefore require the vector pointing from the center to the point on the squeezed r. This vector is found by

$$r\stackrel{\text{trranslate}}{\to }r+\left(a\in \right)\stackrel{\text{squeeze}}{\to }\left(r-\frac{a\in •r}{a+b}\in \right)+\left(b\in \right).$$

When r is on the orbit, the transformed vector has a length of b. The cosine of the difference in eccentric anomalies is the dot product between the directions of the squeezed vectors:

$$\cos \left(E^1-E^0\right)=\frac{1}{b^2}\left(r^1-a\in \frac{\in •r^1}{a+b}+b\in \right)•\left(r^0-a\in \frac{\in •r^0}{a+b}+b\in \right).$$

This simplifies to

$$\cos \left(E^1-E^0\right)=\frac{r^1•r^0}{b^2}+\frac{\left(r^1+r^0\right)•\in }{a}-\frac{r^1•\in \in •r^0}{b^2}+{\in }^2.$$(11)

By noting that the cross product between the two direction vectors gives us the sine, we find in terms of the position vectors:

$$\sin \left(E^1-E^0\right)=\frac{r^0\times r^1}{ab}•\widehat{L}+\frac{\in \times \left(r^1-r^0\right)}{b}•\widehat{L}.$$

If we combine this with Eq. (10), we obtain Eq. (2.69) in Murray & Dermott (2009) for the so-called g-function:

$$t^1-t^0=\frac{E^1-E^0-\sin \left(E^1-E^0\right)}{\omega }+\frac{m\left(r^0\times r^1\right)•L}{\underset{\text{g-function}}{\underbrace{L^2}}}.$$

Although this Equation is remarkable because it does not contain є or the values a, b, we will not need it.

Differentiating Eq. (11) with respect to time t¹ in the end-point gives another equation:

$$\frac{-\omega a}{r^1}\sin \left(E^1-E^0\right)=\frac{{\upsilon }^1•r^0}{b^2}+\frac{{\upsilon }^1•\in }{a}-\frac{{\upsilon }^1•\in \in •r^0}{b^2}.$$(12)

These results can be verified by direct substitution of Eqs. (B.4), (B.5), and (B.7) into the right-hand side of Eq. (12). Equation (10) with the smallest non-negative value for E¹ − E⁰ that satisfies Eqs. (11) and (12) gives the time to get from r⁰ to r¹.

6 Calculating the time to collision

We consider two planets 1 and 2, with a MOID of less than s₁ + s₂, and we want to determine the time at which the planets collide. To this end, let r₁ be the point on orbit 1 with minimal distance to orbit 2, and r₂ the corresponding point on orbit 2 that is closest to r₁. Let υ₁ and υ₂ be the respective velocities of the planets if they pass these points. Now, assume that there is a possible collision:

$$|r_2-r_1|<s_1+s_2.$$

Let ${r^{\prime }}_i$ be the first time that planet 1 passes r₁, and $t_1^1$ be the first time planet 2 passes r₂. A collision occurs at time t(k, l) when both planets are near the points where the distance to the other ellipse is minimal. At that time, planet 1 then passes the point for the kth time, and planet 2 for the lth time. The collision is therefore at

$$t\left(k,l\right)=t_1^1+T_{1}k+\text{d}{t}_1=t_2^1+T_{2}l+\text{d}{t}_2.$$

Here k and l are integers and dt₁ and dt₂ are small shifts that allow for the fact that the planets only need to be close to the point where the distance between the orbits is minimal. Because the algorithm moves forward in time,

$$\begin{array}{cc}k\ge 0,& l\ge 0.\end{array}$$

The shifts in time from the point of closest approach are therefore

$$\begin{array}{cc}\text{d}{t}_1=t\left(k,l\right)-t_1^1-T_{1}k,& \text{d}{t}_2=t\left(k,l\right)-t_2^1-T_{2}l,\end{array}$$(13)

and these need to be small. We linearize the motion about the collision time t(k, l), as

$$\begin{array}{cc}{r}_1+{\upsilon }_1\text{d}{t}_1,& r_2+{\upsilon }_2\text{d}{t}_2.\end{array}$$

By solving for the closest approach between the two particles (in contrast to the MOID, the smallness of the differentials will be a consequence of the fact that the minimal distance for colliding particles is smaller than s₁ + s₂. For this, we introduce the difference vectors

$$\begin{array}{cc}{r}_{12}\left(t\right)=r_2+{\upsilon }_2\text{d}{t}_2-r_1-{\upsilon }_1\text{d}{t}_1,& {\upsilon }_{12}={\upsilon }_2-{\upsilon }_1.\end{array}$$

We note that dt₁ and dt₂ need to be considered as functions of the collision time t. At this time t, the distance is minimal, which is at (see Eq. (A.1))

$$t\left(k,l\right)=\frac{-{\upsilon }_{12}•r_{12}\left(0\right)}{{\upsilon }_{12}^2},$$

with υ₁₂ = |υ₁₂|. The value of the distance must be smaller than the sum of the planet radii (see JeongAhn & Malhotra 2017):

$$\frac{|{\upsilon }_{12}\times r_{12}\left(0\right)|}{{\upsilon }_{12}}<s_1+s_2.$$

When we expand this equation, we find

$$|\left(r_2-r_1-{\upsilon }_2{t}_2^1-{\upsilon }_2{T}_{2}l+{\upsilon }_1{t}_1^1+{\upsilon }_1{T}_{1}k\right)\times {\upsilon }_{12}|<\left(s_1+s_2\right){\upsilon }_{12}.$$

Because (r₂ − r₁) • υ₁ = (r₂ − r₁) • υ₂ = 0, this is equivalent to

$$|\left({\upsilon }_1{t}_1^1\text{+}{\upsilon }_1{T}_{1}k-{\upsilon }_2{t}_2^1-{\upsilon }_2{T}_{2}l\right)\times {\upsilon }_{12}|<\sqrt{{\left(s_1+s_2\right)}^2-|{r_2-r_1|}^2}{\upsilon }_{12}.$$

Using υ₁₂ = υ₂ − υ_ı, this can be further simplified to

$$|t_1^1+T_{1}k-t_2^1-T_{2}l|<\frac{\sqrt{{\left(s_1+s_2\right)}^2-|{r_2-r_1|}^2}{\upsilon }_{12}}{|{\upsilon }_1\times {\upsilon }_2|}.$$

We need to find the smallest non-negative integers k, l for which this inequality is satisfied. We can now recast the problem of the time to collision as finding the smallest integers k, l, so that

$$1-\delta <kp-lq<1+\delta .$$(14)

In this inequality, we use dimensionless parameters p, q, δ, which are defined as

$$\begin{array}{ccc}p=\frac{T_1}{|t_1^1-t_2^1|},& q=\frac{T_2}{|t_1^1-t_2^1|},& \delta =\frac{\sqrt{{\left(s_1+s_2\right)}^2-|{r_2-r_1|}^2}{\upsilon }_{12}}{|t_1^1-t_2^1||{\upsilon }_1\times {\upsilon }_2|}\end{array}.$$

The linearization of the motion around the crossing points of the nodal line translates the collision problem into integer linear programming (in two dimensions). We can find the exact solution in a few steps, even if k and l turn out to be very large numbers.

We assume, without loss of generality, that T₁ > T₂. Consequently, p > q > 0 and p > 1. Equation (14) says that we need an integer linear combination of the irrationals p and q that is within a distance δ from 1. We therefore need to find the point in the grid ℕ² closest to the origin that lies in between the lines

$$\begin{array}{cc}x=\left(\frac{q}{p}\right)y+\frac{1+\delta }{P}\text{\hspace{0.17em}\hspace{0.17em}and}& x=\left(\frac{q}{p}\right)y+\frac{1+\delta }{P}\end{array}$$

in the first quadrant of the xy/-plane (ℝ²). This is shown in Fig. 12. The (horizontal) width of the narrow band is δ/p, which will be of order of magnitude of s/a. For the following solution method to give the correct values of k and l, the numerical accuracy needs to be smaller than this number.

Fig. 12 Search space. The problem of finding the collision time for two planets is equivalent to finding the grid point (k, I) in the small blue band that is closest to the origin. The horizontal and vertical axes are given by $x=\left(t-t_1^1\right)/T_1$ and $y=\left(t-t_2^1\right)/T_2$ in Eq. (13), respectively. The slope is the ratio p/q = T1/T2 of orbital periods. The center of the intersection with the horizontal axis is $T_2/|t_1^1-t_2^1|$, where $|t_1^1-t_2^1|$ is the difference between the time-oſ-passage of the collision point for the two planets. The width of the band depends on the minimal distance of the orbits and the planet radii. The solution can be rapidly found using the bases {b2n,b2n+1}

6.1 Deterministic collision time

Here we present an iteration scheme that finds (k, l) in log(q/δ) steps. We will need the continued-fraction representation of p/q.

This is found from the recursion relation that calculates the successive remainders:

$$\begin{array}{lllll}{q}_0=p,\hfill & q_1=q,\hfill & q_{n+2}=\left(q_n\text{mod\hspace{0.17em}}{q}_{\text{n+1}}\right),\hfill & \text{for}\hfill & n=0,1,2,\mathrm{...}\hfill \end{array}$$

This is similar to the Euclidian algorithm for finding the greatest common divisor. For rational p/q, the sequence becomes zero in finite steps, and for irrational p/q the sequence decreases to zero (Khinchin 1964; Rockett & Szüsz 1992):

$$\begin{array}{cc}0<q_{n+1}<q_n\to 0,\text{\hspace{0.17em}\hspace{0.17em}as}& n\to \infty .\end{array}$$

The q_n are all integer linear combinations of p, q, and are therefore elements of ℤ-span{p, q}. When we write q_n = k_np − l_nq, then the rationals l_n/k_n are precisely the successive convergents of the continued-fraction representation.

Next, we define the bases {b_n, b_n+1} for ℤ², with slopes equal to the fractions:

$$\begin{array}{cc}{b}_n={\left(-1\right)}^n\left(\frac{k_n}{l_n}\right),& \frac{l_{\text{2n}}}{k_{\text{2n}}}\end{array}<\frac{p}{q}<\frac{l_{2n+1}}{k_{2n+1}}.$$

The slopes of the even base vectors increase and the slopes of the odd base vectors decrease to p/q. This is shown in Fig. 12. The sequence of remainders may be found from

$$\left(\begin{array}{c}{q}_{n+2}\\ q_{n+1}\end{array}\right)=\left(\begin{array}{c}-\perp q_{\text{n}}/q_{n+1}\perp \\ 1\end{array}\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{c}{q}_{n+1}\\ q_n\end{array}\right),$$

and base vectors can be found from the recursion relation:

$$\begin{array}{ccc}{b}_0=\left(\begin{array}{c}1\\ 0\end{array}\right),& b_1=\left(\begin{array}{c}0\\ 1\end{array}\right),& b_{\text{n+2}}=b_n+\lfloor \frac{q_{\mathrm{n}}}{q_{n+1}}\rfloor b_{n+1}.\end{array}$$

Therefore, each basis is related to the preceding basis by the transformation

$$\left(b_{n+1}|b_{n+2}\right)=\left(b_n|b_{n+1}\right)\left(\begin{array}{c}0\\ 1\end{array}\begin{array}{c}1\\ \lfloor q_n/q_n{}_{+1}\rfloor \end{array}\right).$$

The transformation matrix is unimodular, which implies that the transformation is a bijection between the points of ℤ². The integer coordinates (x_n, y_n) in each base are defined by

$$\left(\begin{array}{c}x\\ y\end{array}\right)=x_n{b}_n+y_n{b}_{n+1},$$

which means that x₀ = x and y₀ = y. Substitution in Eq. (14) gives us the inequalities

$$1-\delta <x_{0}p-y_{0}q={\left(-1\right)}^n\left(x_n{q}_n-y_n{q}_{n+1}\right)<1+\delta .$$

Consider the following bases, which are composed of the successive even- and odd-numbered vectors:

$$\begin{array}{llll}\left\{b_0,b_2\right\},\hfill & \left\{b_1,b_3\right\},\hfill & \left\{b_2,b_4\right\},\hfill & \left\{b_3,b_5\right\},\mathrm{....}\hfill \end{array}$$

The positive spans (linear combinations) of these pairs partition the first quadrant of ℝ² into segments (see again Fig. 12). As we may assume that (k, l) = (0,0) is not a solution, the band intersects the y-axis at negative y and therefore lies below the segments spanned by the odd bases. However, in any even basis, the intersection of the band with a segment is a trapezoid. Because the union of all segments is the entire first quadrant, the solution of Eq. (14) must lie in one of the trapezoids. However, the pairs of even and odd base vectors do not form a unimodular matrix, and therefore they do not necessarily span ℤ². For this reason, we must instead search in the original bases {b_n,b_n+1} inside the parallelogram of the strip between y′ = 0 and the y′-value where the bottom of the strip intersects span{b_n+2}. This is the top-right point and the bottom-right point of the subsequent parallelogram. The parallelogram has corner points:

$$\begin{array}{ll}\left(x^{\prime },y^{\prime }\right)=\left(\frac{1-\delta }{q_{\text{n}}},0\right),\hfill & \left(\frac{1-\delta }{q_{\text{n}}}+\frac{\left(1+\delta \right)\left(q_{\text{n}}-q_{\text{n+2}}\right)}{q_{\text{n}}{q}_{\text{n+2}}},\frac{\left(1+\delta \right)\left(q_{\text{n}}-q_{\text{n+2}}\right)}{q_{\text{n}}+1q_{\text{n+2}}}\right),\hfill \\ \left(\frac{1+\delta }{q_{\text{n}}},0\right),\hfill & \left(\frac{1+\delta }{q_{\text{n+2}}},\frac{\left(1+\delta \right),\left(q_{\text{n}}-q_{\text{n+2}}\right)}{q_{\text{n+1}}{q}_{\text{n+2}}}\right),\hfill \end{array}$$

where n is even and the coordinates are defined by

$$\left(\begin{array}{c}x\\ y\end{array}\right)=x^{\prime }{b}_{\text{n}}+y^{\prime }{b}_{\text{n+1}}=\left(\begin{array}{c}{x}^{\prime }{k}_{\text{n}}-y^{\prime }{k}_{\text{n+1}}\\ x^{\prime }{l}_{\text{n}}-y^{\prime }{l}_{\text{n+1}}\end{array}\right).$$

Successive parallelograms overlap. One needs to search one entire parallelogram first before going to the next, because the integer point closest to the origin will be found at the earliest occasion. Also, for each y′-value, only the lowest integer x′-value to the right of the left line segment needs to be checked³.

When the calculated orbital periods T₁ and T₂ have a ratio close to that of two small integers, the planets could be in mean-motion resonance and may never collide (as is the case with Neptune and Pluto). The method neglects these cases, and will erroneously find a collision at a high k number. If the continued fraction is actually finite (because p/q is rational), the strip has an integer slope at the final step. Figure 13 shows the results of a numerical test with random collision pairs. We find that the total number of checks grows as k^1/2, while the number of iterations grows as log k with the solution k.

Fig. 13 Numerical checks needed to calculate the exact collision time for random collision pairs. Here we show the number of checks (orange), the number of iterations (black), and the average number of checks per iteration (blue) performed by the algorithm. On the horizontal axis is the number of orbits k before the collision, which is the solution calculated by the iteration scheme. The gray line is #checks = k1/2 and shows the trend.

6.2 Stochastic collision times

For small δ and irrational p/q, a generic solution (k, l) will form a pair of large integers, with

$$\frac{l}{k}\approx \frac{p}{q}.$$

The precise value depends very sensitively on q and p. If we assume that we cannot obtain the required numerical accuracy to find the exact solution, we may use a statistical approach. The integer points (k, l) are uniformly distributed over the plane. If we assume that the points are statistically independent (this it clearly an approximation valid for δ << 1), the distribution is that of a Poisson point process. The probability that there is a grid point (k, l) between the lines with k in the interval [x, x + dx] is equal to the area of the small parallelogram. This area is (2δ/q)dx. Now, the area between the lines below k = x is given by

$$\begin{array}{cc}\frac{2\delta }{q}x-\frac{2\delta -4{\delta }^2}{pq}\approx \frac{2\delta }{q}x,\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for}& x\gg 1.\end{array}$$

Therefore, the respective probabilities for the solution k to be found above x and below x are

$$\begin{array}{ccc}\text{Prob}\left(k>|x\right)=& \begin{array}{c}\lim \\ \text{dx}\to 0\end{array}{\left(1-\frac{2\delta }{q}\text{dx}\right)}^{\text{x}/\text{dx}}& =e^{-2x\delta /q}\end{array},$$

$$\begin{array}{cc}\text{Prob}\left(k\le x\right)=& 1-e^{-2x\delta /q}.\end{array}$$

The latter formula is the cumulative distribution function for k. We obtain a realization by generating a random real number ξ inside the interval [0,1] and use

$$\begin{array}{cc}k=\frac{\left(-\log \xi \right)q}{2\delta },& l=\frac{\left(-\log \xi \right)p}{2\delta }-\frac{1}{q}\end{array}.$$

Because these values are large, rounding off to the nearest integer is not important. The time of the collision is then given by

$$t=t_1^1+k{T}^1=t_1^1+\frac{\left(-\log \xi \right)T_1{T}_2|{\upsilon }_1\times {\upsilon }_2|}{2\sqrt{{\left(s_1+s_2\right)}^2-|{r_2-r_1|}^2{\upsilon }_{12}}}.$$

The formula for the average waiting time $\overline{k}{T}^1$ is precisely the reciprocal of the collision probability for one orbit. The fact that this reproduces the formula Eq. (23) in Öpik (1951), Eq. (29) in JeongAhn & Malhotra (2017), and Eqs. (2) and (3) in Diserens et al. (2020) for this probability P validates our Eq. (14).

The Monte Carlo method for finding the collision time also follows from assuming homogeneous distributions of the mean anomalies of two fixed Kepler orbits, as in Öpik’s scheme. This method assumes large k, implying that the initial crossing time t₁ cannot be precisely known (due to numerical inaccuracy or neglected physics effects). However, when the system contains one or more large planets, the collision or nearby passage could happen after a few revolutions, that is, for small k.

6.3 Including orbital precession

Our method for calculating the collision time outlined in the previous subsection assumes perfect Kepler orbits. However, if one intends to make accurate predictions over long timescales, the slow precession of the periapsis and of the orbital plane cannot be neglected. For example, the perihelion shift for the planet Jupiter in one revolution is about twice the planet’s diameter (see Fig. 14 for the precession rates for the Solar System planets). For a planet ring system, precession is mainly due to the oblateness of the planet. Hence, if one is not just interested in statistical averages, then the secular dynamics must be included on the orbital timescale.

A method to remedy this problem is to numerically propagate the Laplace-Lagrange equations. The time-steps are now set by the secular timescale dt ⪅ T_prec. There are N² terms in the system of differential equations and there are N²ϵ/2 collision possibilities, which need to be evaluated. The theoretical algorithmic efficiency is shown by the orange curve in Fig. 8.

It may be possible to include the collision detection in the following way. At each time-step (now shorter than the collision time), one calculates the instantaneous change is the orbital elements є, ϖ, I, Ω. One expresses the resulting linear change in the points r_j and velocities υ_j near the points of closest approach, as linear functions in the passage numbers k and l. The step where the minimal distance d is compared with the sum s_i + s_j is skipped. Instead one directly uses the modified inequality Eq. (14). This would then lead to a problem from integer linear programming, as before. For this modified case, we would expect the two lines bounding our search domain to be nonparallel. The solution for the collision problem has a spatial accuracy of less than a planet radius on the longer timescale where the perihelion shift can be approximated as linear motion. A thorough development of this idea is a possible direction for follow-up research.

Fig. 14 Timescale vs. particle count N. The horizontal lines are the orbital periods Tj and the precession periods Tprecj. of the Solar System planets (gray, brown, blue, red, coral, golden, turquoise, azure are Mercury through Neptune) from Murray & Dermott (2009). The black curve shows the collision time Tcoll estimated for a disk with a = 4au, I = .1, and a mass of three Earths containing N particles of Earth density. For N < 10, the collision time is comparable to the precession time (105yr). The radius of gravitational influence is roughly 103 times the planet radius s; this determines Tscatt, the dashed curve.

7 Conclusions

We describe an algorithm that simulates collisional Keplerian systems: N bodies in the Coulomb potential of a central mass. The method uses the orbital elements and has three basis ingredients, of which the third is novel: (i) for a new particle i, a small set of possible collision candidates j is selected using the apoap-sis/periapsis filter. (ii) The MOIDs between the particle pairs (i, j) can be approximated numerically. (iii) For the pairs (i, j) on a collision trajectory, one can obtain the collision time t_(i,j) with integer linear programming. During the simulation, sorted lists of the particles and the collision pairs are maintained. The algorithm steps from one collision to the next as it updates the particle orbits and collision possibilities.

We show that the problem of finding the collision time is mathematically equivalent to the problem in integer linear programming of finding the grid point (k, l) in ℕ² between two parallel lines that is closest to the origin. The exact solution uses the continued-fraction representation of the ratio T_i/T_j of the orbital periods.

Because at most N new collision possibilities have to be added to the list, less than N interactions need to be considered at each step. The length of the collision list is O(N²s/a) and the total number of collisions is O(N²s²/a²), resulting in an algorithmic efficiency of O(N⁴s³/a³). This may be compared to the efficiency O(N^4/3 log N) of a numerical integration propagator with collision detection (tree code or spatial hashing), which is independent of the particle radius s. In the astronomical applications, the radii are usually small compared to the orbits. The collisions are therefore rare, and the proposed collision-detection method can be fast. However, the perturbations we neglect become increasingly important, and, at the same time, the result becomes progressively sensitive to the initial state. Needless to say, including collisions is important, even if they are rare. Studying statistics of outcomes of the dynamics requires many simulations with near-identical initial states. For relatively small particle numbers, say for N < 10⁶, the individual realizations can be fast.

Acknowledgements

The author would like to thank John Chambers for acting as referee and for improving the algorithm, and Dylan Aliberti, Aron Schouten, and Philip Soliman for their meticulous checking and verification of the formulas and algorithms in this paper.

Appendix A Closest approach

Consider two non-interacting particles in linear motion:

$$r_1\left(t\right)=r_1+{\upsilon }_{1}t,r_2\left(t\right)=r_2+{\upsilon }_{2}t.$$

We call the intitial distance vector and the relative velocity:

$$d=r_2-r_1,u={\upsilon }_2-{\upsilon }_1.$$