Issue 
A&A
Volume 669, January 2023



Article Number  A97  
Number of page(s)  15  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202243754  
Published online  18 January 2023 
Collision detection for Nbody Kepler systems
Delft Institute of Applied Mathematics, Delft University of Technology,
Mekelweg 4,
2628 CD
Delft, The Netherlands
email: p.m.visser@tudelft.nl
Received:
11
April
2022
Accepted:
29
October
2022
Context. In a Keplerian system, a large number of bodies orbit a central mass. Accretion disks, protoplanetary disks, asteroid belts, and planetary rings are examples. Simulations of these systems require algorithms that are computationally efficient. The inclusion of collisions in the simulations is challenging but important.
Aims. We intend to calculate the time of collision of two astronomical bodies in intersecting Kepler orbits as a function of the orbital elements. The aim is to use the solution in an analytic propagator (Nbody simulation) that jumps from one collision event to the next.
Methods. We outline an algorithm that maintains a list of possible collision pairs ordered chronologically. At each step (the soonest event on the list), only the particles created in the collision can cause new collision possibilities. We estimate the collision rate, the length of the list, and the average change in this length at an event, and study the efficiency of the method used.
Results. We find that the collisiontime problem is equivalent to finding the grid point between two parallel lines that is closest to the origin. The solution is based on the continued fraction of the ratio of orbital periods.
Conclusions. Due to the large jumps in time, the algorithm can beat tree codes (octree and kd tree codes can efficiently detect collisions) for specific systems such as the Solar System with N < 10^{8}. However, the gravitational interactions between particles can only be treated as gravitational scattering or as a secular perturbation, at the cost of reducing the timestep or at the cost of accuracy. While simulations of this size with highfidelity propagators can already span vast timescales, the high efficiency of the collision detection allows many runs from one initial state or a large sample set, so that one can study statistics.
Key words: gravitation / methods: analytical / methods: statistical / celestial mechanics / planets and satellites: formation / protoplanetary disks
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
Simulations of the mechanical motion of many bodies are generally computationally expensive. Consider the Nbody 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^{2} interactions to account for. Various methods have been invented to speed up the simulation or increase the particle count: (i) direct Nbody simulations with dynamic timesteps (see Dehnen & Read 2011, for an overview); (ii) the octree code for collision detection (Bentley 1975; Meagher 1980); (iii) the BarnesHut 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 multipole 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 finiteelements 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 LaplaceLagrange equations for the orbital elements (Murray & Dermott 2009); because the Kepler ellipses change slowly over time, the timestep 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 timestep, the volume of space around the particle that needs to be probed for collision partners has a radius of the timestep 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 timestep. 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 Nbody 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, timesteps 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 nearestneighbor 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.
List of symbols and notation.
Fig. 1 Two orbits separated by a sphere (dashed). Orbit 1 (blue) has apoapsis a_{1} + c_{1} and orbit 2 (purple) has periapsis a_{2} − c_{2}. Filtering out such collisionavoiding apoapsisperiapsis 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 nearcircular orbits. A particle with radius s traces a cylindrical volume of size (2πa/T)πs^{2} per unit of time. The disk is a cylinder of radius a and height 2Ia, meaning that the particle density is N/2πIa^{3}. Accounting for the N^{2}/2 pairs, the rate of collisions is estimated to be
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/5}a. The timescales are therefore
The formula for the precession is taken from Murray & Dermott (2009). If we model the early inner Solar System by N = 10^{6} planetesimals of characteristic size s = 100 km in a disk with a = 4 au and I = .1, we have
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^{13}, s = 1m, and a = 10^{5} km, which results in
The precession is now entirely due to planet oblateness (J_{2} = 3 × 10^{−3}). 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
The orbital eccentricity e accounts from the fact that the relative velocity is roughly , 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^{−5} and 10^{−7} 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^{2} 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) timestep 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 m_{1}, m_{2} in their orbital planes. Because the bodies are much smaller than the orbits, s_{j} << a_{j}, 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 ≤ a_{max} = 2au, I ≤ 10^{−3}, ϵ ≤ 10^{−3}, and with N = 10^{4}, 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 timesteps
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^{3} of a sufficiently small size to contain only one or a few particles. In one timestep, 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 (socalled 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 timesteps for one orbital period (dt ⪅ T) but also spatial steps of the order of the interparticle distance . 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:
Clearly, the timesteps 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^{−3}au. 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 timesteps, 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^{2} 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^{3}/a^{3} = 10^{−15} 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^{2} s/a)(T/T_{coll}) ∝ N^{3} and N^{2}T/T_{prec} ∝ N^{2}, 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^{8}.
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. 
Algorithmic efficiency.
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(N^{3}) instead of O(N^{4}). 
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 T_{coll} = T_{prec}, and has a negligible effect above the orange dashed line, aT_{coll} = sT_{prec}. If the mass of three Earths is distributed over equalsized planetoids (of Earth density), one is constrained to the blue line. On this line, precession is small for N > 10^{3} but only attains the much smaller error of s per orbit for N > 10^{9}. Clearly, treecode 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 phasespace 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) threebody gravitational scattering, (iii) planetary migration due to interaction with the gas in the protoplanetary disk, (iv) capture of planets in meanmotion resonances, (v) the Kozai mechanism, (vi) moons and binaries, (vii) atmospheric and (viii) tidal drag, (ix) solar wind, and (x) the PoyntingRobertson 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:
Here, t^{0} is the time of creation, α, c are the orbital radius and focal distance, s is the particle radius, m is the particle mass, r^{0} 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) phasespace 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 . 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
If the collision list is empty, end the simulation.
Take the pair (i, j) with the soonest collision: the first in the list.
Update the time t to the time of the collision.
Remove any pair containing i and any pair containing j from the pairs list.
Remove the particles and j from the particle list.
If the orbit of the new particle intersects the central mass or is unbound, go to the next collision on the list.
Create new particle(s) defined by
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.
Make a sorted list of the new collision possibilities with a record of the collision time and the pair, soonest collision first.
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:
including an extension s_{i} of the size of the particle, or with the substitution of its gravitational reach s_{i} = (m /M)^{2/5}a_{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}/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_{1} > a_{2}. First we need to find the minimal orbit intersection distance to decide whether or not there can be a collision. For this, retrieve m_{1} L_{1} є_{1} and m_{2}, L_{2}, є_{2}. Then calculate the direction of the nodal line (see Fig. 2), the semilatus recti, the intersection points, and the velocities at these points:
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_{1}, r_{2}) 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:
Equation (4) for the collision points r′_{1} and r′_{2} is derived in Sect. 4.2.
3. Retrieve a_{1}, s_{1} and a_{2}, s_{2}. Decide whether or not an interaction can take place.
d < s_{1} + s_{2} ⇒ 1 and 2 collide.
d < (m_{1}/M)^{2/5}a_{1} ⇒ 1 can perturb orbit 2.
d < (m_{2}/M)^{2/5}a_{2} ⇒ 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 , ω_{1}, and , ω_{2}, of the particles involved (with ω_{1} < ω_{2}). Calculate the first time a particle passes the crossing point in Fig. 2. These times are denoted by ,.
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.
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 ω_{2}/ω_{1}. Initialize the loop:
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
We use the notation and for the floor and the ceiling function. The time T_{sim} to be simulated sets an upper bound for the solution:
We then test the points with coordinates
if x>x_{max}⇒ next pair
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
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 timestep, 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 multibody gravitational scattering will therefore be treated as three (or more) successive twobody 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
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 breakup 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^{1} calculated from the basic conservation laws:
(material volume),
m = m_{1} + m_{2} (mass),
(centerofmass motion),
(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.
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:
4. Register {t^{0}, 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 newparticle lists with the existingparticle list.
Go to the next timestep.
Fig. 9 Gravitational scattering between particles 1 and 2 in the centerofmass frame. The initial velocities υ_{1}, υ_{2} are scattered in the directions , 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 2a_{j} by the conjugate axis 2b_{j} (Adams & Essex 2021). The impact parameter is the sum b = b_{1} + b_{2}. 
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 (T_{coll} ≫ 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_{1} + s_{2}. 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_{1} s_{2})/(a_{1} a_{2}), 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_{1} + s_{2})/ about the “mutual nodal line”, all points in orbit 1 are separated by more than s_{1} + s_{2} 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 , these small inclinations are rare for .
Because the range for gravitational scattering can be larger, our approach will only work for lowmass 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 xaxis by b/a. The travel time from r^{0} to r^{1} 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 πb^{2}. The circle segment has area ΔEb^{2}/2, with ΔE = E^{1} − E^{0} 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_{1}, є_{1} and L_{2}, є_{2}, respectively. The line of intersection of the two orbital planes, or the nodal line, can be found as follows. Because the angular momenta L_{1} and L_{2} are both perpendicular to the nodal line, a direction vector of the nodal line is
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
The point of intersection can now be found from the formula for the orbit Eq. (B.3). We find
Therefore, for the two pairs of intersection points, we obtain
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_{1} and r_{2} of intersection with the nodal line. The tangent lines point in the direction of the velocities υ_{1} and υ_{2} at r_{1} and r_{2}. These can be found using Eq. (B.1). The distance between the two lines is given by
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 and . These positions are given by
One may verify that . This proves that the minimal distance between the lines is realized at the points and . We also have that , implying that
If the inclination between the orbital planes, which is given by I ≈ K/L_{1}L_{2}, is not much larger than (s_{1} + s_{2})/(a_{1} + a_{2}), the linear approximation is inaccurate. One can improve this approximation by reducing the lengths of , 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^{0} to r^{1} . In the algorithm, r^{0} is the particle’s creation point and r^{1} 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^{2}. 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^{0} minus the triangle between −aє, 0, and r^{1} in Fig. 11a. Therefore,
Kepler’s second law therefore implies that the time it takes a body to move from r^{0} to r^{1} is equal to
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
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
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:
This simplifies to
By noting that the cross product between the two direction vectors gives us the sine, we find in terms of the position vectors:
If we combine this with Eq. (10), we obtain Eq. (2.69) in Murray & Dermott (2009) for the socalled gfunction:
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^{1} in the endpoint gives another equation:
These results can be verified by direct substitution of Eqs. (B.4), (B.5), and (B.7) into the righthand side of Eq. (12). Equation (10) with the smallest nonnegative value for E^{1} − E^{0} that satisfies Eqs. (11) and (12) gives the time to get from r^{0} to r^{1}.
6 Calculating the time to collision
We consider two planets 1 and 2, with a MOID of less than s_{1} + s_{2}, and we want to determine the time at which the planets collide. To this end, let r_{1} be the point on orbit 1 with minimal distance to orbit 2, and r_{2} the corresponding point on orbit 2 that is closest to r_{1}. Let υ_{1} and υ_{2} be the respective velocities of the planets if they pass these points. Now, assume that there is a possible collision:
Let be the first time that planet 1 passes r_{1}, and be the first time planet 2 passes r_{2}. 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
Here k and l are integers and dt_{1} and dt_{2} 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,
The shifts in time from the point of closest approach are therefore
and these need to be small. We linearize the motion about the collision time t(k, l), as
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_{1} + s_{2}. For this, we introduce the difference vectors
We note that dt_{1} and dt_{2} 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))
with υ_{12} = υ_{12}. The value of the distance must be smaller than the sum of the planet radii (see JeongAhn & Malhotra 2017):
When we expand this equation, we find
Because (r_{2} − r_{1}) • υ_{1} = (r_{2} − r_{1}) • υ_{2} = 0, this is equivalent to
Using υ_{12} = υ_{2} − υ_{ı}, this can be further simplified to
We need to find the smallest nonnegative 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
In this inequality, we use dimensionless parameters p, q, δ, which are defined as
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_{1} > T_{2}. 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 ℕ^{2} closest to the origin that lies in between the lines
in the first quadrant of the xy/plane (ℝ^{2}). 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 and in Eq. (13), respectively. The slope is the ratio p/q = T_{1}/T_{2} of orbital periods. The center of the intersection with the horizontal axis is , where is the difference between the timeoſ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 {b_{2n},b_{2n+1}} 
6.1 Deterministic collision time
Here we present an iteration scheme that finds (k, l) in log(q/δ) steps. We will need the continuedfraction representation of p/q.
This is found from the recursion relation that calculates the successive remainders:
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):
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_{n}p − l_{n}q, then the rationals l_{n}/k_{n} are precisely the successive convergents of the continuedfraction representation.
Next, we define the bases {b_{n}, b_{n+1}} for ℤ^{2}, with slopes equal to the fractions:
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
and base vectors can be found from the recursion relation:
Therefore, each basis is related to the preceding basis by the transformation
The transformation matrix is unimodular, which implies that the transformation is a bijection between the points of ℤ^{2}. The integer coordinates (x_{n}, y_{n}) in each base are defined by
which means that x_{0} = x and y_{0} = y. Substitution in Eq. (14) gives us the inequalities
Consider the following bases, which are composed of the successive even and oddnumbered vectors:
The positive spans (linear combinations) of these pairs partition the first quadrant of ℝ^{2} into segments (see again Fig. 12). As we may assume that (k, l) = (0,0) is not a solution, the band intersects the yaxis 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 ℤ^{2}. 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 topright point and the bottomright point of the subsequent parallelogram. The parallelogram has corner points:
where n is even and the coordinates are defined by
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^{3}.
When the calculated orbital periods T_{1} and T_{2} have a ratio close to that of two small integers, the planets could be in meanmotion 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 = k^{1/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
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
Therefore, the respective probabilities for the solution k to be found above x and below x are
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
Because these values are large, rounding off to the nearest integer is not important. The time of the collision is then given by
The formula for the average waiting time 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_{1} 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 LaplaceLagrange equations. The timesteps are now set by the secular timescale dt ⪅ T_{prec}. There are N^{2} terms in the system of differential equations and there are N^{2}ϵ/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 timestep (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 followup research.
Fig. 14 Timescale vs. particle count N. The horizontal lines are the orbital periods T_{j} and the precession periods T_{precj}. 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 T_{coll} 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 (10^{5}yr). The radius of gravitational influence is roughly 10^{3} times the planet radius s; this determines T_{scatt}, 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 apoapsis/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 ℕ^{2} between two parallel lines that is closest to the origin. The exact solution uses the continuedfraction 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^{2}s/a) and the total number of collisions is O(N^{2}s^{2}/a^{2}), resulting in an algorithmic efficiency of O(N^{4}s^{3}/a^{3}). 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 collisiondetection 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 nearidentical initial states. For relatively small particle numbers, say for N < 10^{6}, 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 noninteracting particles in linear motion:
We call the intitial distance vector and the relative velocity:
We want to decide if there is a collision in the interval [0, dt]. Therefore, we calculate the distance between the particles at t = 0 and at t = dt to see if there is a collision at the endpoints:
If not, the only possibility for a collision on the interval is that the distance obtains a minimum on the (interior) of the interval. This means that the relative velocity in the direction between the particles is first decreasing and then increasing:
The time of minimal distance is at
which is then indeed between 0 and dt. We then decide if the distance at this time t is smaller than the sum of the radii. When we substitute t back into the Equation for the distance, we find that it is equal to the component of d perpendicular to the direction of u. Hence, the condition for a collision is equivalent to:
Appendix B Orbital elements
Consider a single particle of mass m in a Kepler orbit about the central mass M. The orbit is an ellipse in a fixed plane. The angular momentum vector is defined by
The LaplaceRungeLenz vector is proportional to the dimensionless eccentricity vector:
The orbit is fixed by the orthogonal pair of vectors L and є.
For the problem of solving the MOID, we need a formula that expresses the velocity υ along the orbit as a function of the position r and the orbital elements. Given L, є, r, we equate
and therefore the velocity can be expressed as:
The Equation for the energy is called the visviva equation
With this Equation and Kepler’s third law,
the orbital elements a, ϵ and mean motion ω can be calculated from position and velocity:
The orbit is parametrized by the true anomaly v or the eccentric anomaly E. If we know the eccentric anomaly, we can calculate the time since periapsis from the Kepler equation
The formula for the radial distance in terms of the parameters is
The semimajor axes, a and b, the distance from the center to a focus, c, and eccentricity ϵ are related by
The position vector and the velocity vector can now be expressed as:
and, using Eqs. (B.2) and (B.3),
In the expressions Eqs. (B.4) and (B.5), ℛ is a constant rotation matrix, which can be expressed as a product of three elementary rotations
Here, ϖ is the argument of perihelion, I is inclination, and Ω is the ascending node. We have for the angular momentum and the eccenticity vectors,
and
These five orbital elements a, ϵ, ϖ, Ω, I represent five independent conserved quantities. They can be found from the position and velocity vectors as follows. The equations at the beginning of this section give us L and є in terms of r, υ. Then, the component l = L_{3} can give us the angle I. The component L_{1} can give us Ω, and ϵ_{3} can give us the angle ϖ.
Appendix C Elastic collisions
When planets approach each other so closely that the gravity they exert on each other is of the same strength as the gravity from the star, the mutual gravitational interaction cannot be neglected. The close encounter can be approximated as a scattering event. Here, for the short duration of the encounter, we neglect the gravity from the central star (and of all other particles in the system). As no energy is dissipated, the process is an elastic collision. It is described by a hyperbolic orbit for the relative coordinate r_{2} – r_{1}. The hyperbola has parameters a, b, c, and is given by the following equations in Cartesian and in polar coordinates
The parameter a is the semitransverse axis of the hyperbola. It is a characteristic distance where the gravity between the two bodies becomes noticeable given the relative velocity. The semiconjugate axis b of the hyperbolic orbit is equal to the impact parameter. The semifocal separation is c (see Fig. 9). We have
When the scattering between two planets is significant, the velocities of the two planets will change. The centerofmass moves with velocity
In the center of mass frame, both velocities rotate over an angle π – 2 arctan(b/a). This is described by:
Appendix D Mathematica code
The following Mathematica code tests the deterministic collisiontime algorithm. It first finds the solution by considering all integers k, and then uses the continuedfraction representation. Both procedures give the same result for random orbital periods.
tmax = 100; T1 = RandomReal [tmax ]; T2 = RandomReal [tmax ]; If[T1 < T2, T1 = tmax − T1; T2 = tmax − T2]; t1 = RandomReal [ T1 ] ; t2 = RandomReal [ T2 ] ; de l ta t = Abs [ t1 − t2 ] ; p = T1 / deltat q = T2 / deltat del = deltat/10^5 qlist = {p, q}; n = 1; While[ qlist [[n + 1]] > del , qnm = q list [ [ n ] ] ; n++; qn = qlist [ [ n ] ] ; AppendTo [ qlist, qnm − Floor [ qnm/qn ] qn ]] qlist For [k = 1, k < 10^6, k++, l = Ceiling[p k/q − (1 + del)/q]; If[l < p k/q − (1 − del)/q, Break[]]] Plot[{ p x /q − (1 − del ) /q , p x /q − (1 + del ) / q}, { x, k − del , k + del } , PlotRange → {l − del, l + del }, GridLines → { { k}, { l } } , AspectRatio → 1, Axes → None , Frame → None] k l q0 = p ; q1 = q ; k0 = 1 ; k1 = 0; l0 = 0; l1 = −1; plot = { } ; found = Fals e ; For [ n = 0 , Not [ found ] , n++ , a0 = Floor [ q0 /q1 ] ; q2 = q0 − a0 q1 ; I f [ q2 == 0, Break [ ] ] ; k2 = k0 − a0 k1 ; 12 = l0 − a0 l1 ; a1 = Floor [ q1 /q2 ] ; For [ x = C eili ng [ ( 1 − del ) / q0 ] , x < (1 + del)/q2, x++, y = Max[ 0, Ceiling [ ( q0 x − 1 − del ) / q1 ] ] ; If [ 1 − del < x q0 − y q1, k = x k0 − y k1; l = x l0 − y l1 ; found = True ; Break[]] ]; q3 = q1 − a1 q2 ; If[q3 == 0, Break []] ; k3 = k1 − a1 k2 ; 13 = l1 − a1 l2 ; q0 = q2 ; q1 = q3 ; k0 = k2 ; k1 = k3 ; l0 = l2 ; l1 = l3 ;] k l Plot[{p x/q − (1 − del)/q, p x/q − (1 + del)/q}, { x , k − del , k + del }, PlotRange → {l − del, l + del }, GridLines → {{k}, {l}}, AspectRatio → 1, Axes → None , Frame → None]
References
 Adams, R. A., & Essex, C. 2021, Calculus A Complete Course, 10th edn. (North York, Ontario: Pearson Education) [Google Scholar]
 Aliberti, G. D. 2022a, Astrophysics Source Code Library [record ascl:2211.002] [Google Scholar]
 Aliberti, G. D. 2022b, Bachelor thesis, TU Delft, The Netherlands [Google Scholar]
 Baraff, D., & Witkin, A. 1992, SIGGRAPH Comput. Graph., 26, 303 [CrossRef] [Google Scholar]
 Barnes, J., & Hut, P. 1986, Nature, 324, 446 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, J. E. 1990, J. Comput. Phys., 87, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Bentley, J. L. 1975, Commun. ACM, 18, 509 [CrossRef] [Google Scholar]
 Bodenheimer, P. H., Laughlin, G., Różyczka, M., & Yorke, H. W. 2007, Numerical Methods in Astrophysics: an Introduction, Series in Astronomy and Astrophysics (New York London: Taylor & Francis) [Google Scholar]
 Burtscher, M., & Pingali, K. 2011, GPU Computing Gems Emerald Edition (Amsterdam: Elsevier) [Google Scholar]
 Dehnen, W., & Read, J. 2011, Euro. Phys. J. Plus, 126, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Diserens, S., Lewis, H. G., & Fliege, J. 2020, J. Space Safety Eng., 7, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Goldstein, H. 1964, Classical Mechanics, Ninth dover printing, Tenth gpo printing edn. (New York: Dover) [Google Scholar]
 Greengard, L. 1990, Comput. Phys., 4, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Gronchi, G. F. 2005, Celest. Mech. Dyn. Astron., 93, 295 [Google Scholar]
 Hamada, T., Nitadori, K., Benkrid, K., et al. 2009, Comput. Sci., 24, 21 [Google Scholar]
 Hedo, J., Ruiz, M., & Pelaez, J. 2018, MNRAS, 479, 3288 [NASA ADS] [CrossRef] [Google Scholar]
 Hoots, F. R., Crawford, L. L., & Roehrich, R. L. 1984, Celest. Mech., 33, 143 [NASA ADS] [CrossRef] [Google Scholar]
 JeongAhn, Y., & Malhotra, R. 2017, AJ, 153, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Khinchin, A. Y. 1964, Continued Fractions (USA: University of Chicago Press) [Google Scholar]
 Manley, S. P., Migliorini, F., & Bailey, M. E. 1998, A&AS, 133, 437 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meagher, D.J.R. 1980, Octree Encoding: A New Technique for the Representation, Manipulation and Display of Arbitrary 3D Objects by Computer, IPLTR80111 [Google Scholar]
 Meagher, D.J.R. 1982, Comput. Graph. Image Process., 19, 129 [CrossRef] [Google Scholar]
 Milisavljevic, S. 2010, Serbian Astron. J., 180, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C., & Dermott, S. 2009, Solar System Dynamics (New York: Cambridge University Press) [Google Scholar]
 Öpik, E. J. 1951, Proc. R. Irish Acad. Sect. A Math. Phys. Sci., 54, 165 [Google Scholar]
 Rockett, A., & Szüsz, P. 1992, Continued Fractions (Singapore: World Scientific) [CrossRef] [Google Scholar]
 Rokhlin, V. 1985, J. Comput. Phys., 60, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Savransky, D., Cady, E., & Kasdin, N. J. 2011, ApJ, 728, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Schouten, A. 2022, Bachelor thesis, TU Delft, The Netherlands [Google Scholar]
 Segan, S., Milisavljevic, S., & Marceta, D. 2011, Acta Astron., 61, 275 [NASA ADS] [Google Scholar]
 Soliman, P. 2022, Bachelor thesis, TU Delft, The The Netherlands [Google Scholar]
 Warren, M., & Salmon, J. 1993, in A parallel Hashed OctTree ABody Algorithm (USA: ACM), 12 [Google Scholar]
 Wizniowski, T., & Rickman, H. 2013, Acta Astron., 63, 293 [NASA ADS] [Google Scholar]
The need for the following calculations was pointed out by Soliman (2022).
It was pointed out by Schouten (2022) that by looping through the integer x values instead of the integer y values the number of checks is significantly lower.
All Tables
All Figures
Fig. 1 Two orbits separated by a sphere (dashed). Orbit 1 (blue) has apoapsis a_{1} + c_{1} and orbit 2 (purple) has periapsis a_{2} − c_{2}. Filtering out such collisionavoiding apoapsisperiapsis pairs is an efficient sweep and prune method. 

In the text 
Fig. 2 Orbits of two planetoids m_{1}, m_{2} in their orbital planes. Because the bodies are much smaller than the orbits, s_{j} << a_{j}, 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. 

In the text 
Fig. 3 Pairs vs. particle radius s for a homogeneous disk with a ≤ a_{max} = 2au, I ≤ 10^{−3}, ϵ ≤ 10^{−3}, and with N = 10^{4}, 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. 

In the text 
Fig. 4 Pairs vs. particle count N, as in Fig. 3, but with particle radius s = 2 × 10^{−3}au. 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. 

In the text 
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. 

In the text 
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(N^{3}) instead of O(N^{4}). 

In the text 
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 T_{coll} = T_{prec}, and has a negligible effect above the orange dashed line, aT_{coll} = sT_{prec}. If the mass of three Earths is distributed over equalsized planetoids (of Earth density), one is constrained to the blue line. On this line, precession is small for N > 10^{3} but only attains the much smaller error of s per orbit for N > 10^{9}. Clearly, treecode is better for these high number densities. 

In the text 
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. 

In the text 
Fig. 9 Gravitational scattering between particles 1 and 2 in the centerofmass frame. The initial velocities υ_{1}, υ_{2} are scattered in the directions , 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 2a_{j} by the conjugate axis 2b_{j} (Adams & Essex 2021). The impact parameter is the sum b = b_{1} + b_{2}. 

In the text 
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 (T_{coll} ≫ 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. 

In the text 
Fig. 11 Elliptical Kepler orbit (a) and orbit squeezed into a circle (b) from scaling the xaxis by b/a. The travel time from r^{0} to r^{1} 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 πb^{2}. The circle segment has area ΔEb^{2}/2, with ΔE = E^{1} − E^{0} being the difference in eccentric anomaly. 

In the text 
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 and in Eq. (13), respectively. The slope is the ratio p/q = T_{1}/T_{2} of orbital periods. The center of the intersection with the horizontal axis is , where is the difference between the timeoſ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 {b_{2n},b_{2n+1}} 

In the text 
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 = k^{1/2} and shows the trend. 

In the text 
Fig. 14 Timescale vs. particle count N. The horizontal lines are the orbital periods T_{j} and the precession periods T_{precj}. 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 T_{coll} 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 (10^{5}yr). The radius of gravitational influence is roughly 10^{3} times the planet radius s; this determines T_{scatt}, the dashed curve. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.