Open Access
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/0004-6361/202243754
Published online 18 January 2023

© The Authors 2023

Licence Creative CommonsOpen 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 Subscribe-to-Open 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 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 N2 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(N1/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(N4/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.

Table 1

List of symbols and notation.

thumbnail Fig. 1

Two orbits separated by a sphere (dashed). Orbit 1 (blue) has apoapsis a1 + c1 and orbit 2 (purple) has periapsis a2c2. 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 Tcoll, the scattering time Tscat, and the secular time Tprec.

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)πs2 per unit of time. The disk is a cylinder of radius a and height 2Ia, meaning that the particle density is N/2πIa3. Accounting for the N2/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/5a. 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 = 106 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 = 1013, s = 1m, and a = 105 km, which results in

The precession is now entirely due to planet oblateness (J2 = 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 N2 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 Tcoll. 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 < Tcoll < Tprec.

thumbnail 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.

thumbnail Fig. 3

Pairs vs. particle radius s for a homogeneous disk with aamax = 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 dx3 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 . 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 time-steps for numerical integration actually need to be quite small.

thumbnail 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 N2 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 (Tsim/dt)N log N ∝ N4/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 Nm = 3MEarth and with particles of the density of Earth. We thus have Ns3/a3 = 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 (N2 s/a)(T/Tcoll) ∝ N3 and N2T/TprecN2, 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 > 108.

thumbnail 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.

Table 2

Algorithmic efficiency.

thumbnail 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).

thumbnail 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.

thumbnail 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:

(1)

Here, t0 is the time of creation, α, c are the orbital radius and focal distance, s is the particle radius, m is the particle mass, r0 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 . 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 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

  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:

including an extension si of the size of the particle, or with the substitution of its gravitational reach si = (m /M)2/5ai . It is sufficient to compare each particle i with the particles j = i + 1, j = i + 2,… Because the list is sorted, we have ai- cisiajcjsj. As long as ai + ci + siajcjsj, the intervals for i and j overlap and the pair is a candidate for collision. Once we encounter the first j where ai + ci + si < ajcjsj, 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 N2/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 a1 > a2. First we need to find the minimal orbit intersection distance to decide whether or not there can be a collision. For this, retrieve m1 L1 є1 and m2, L2, є2. Then calculate the direction of the nodal line (see Fig. 2), the semi-latus recti, the intersection points, and the velocities at these points:

(2)

(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 (r1, r2) 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:

(4)

Equation (4) for the collision points r1 and r2 is derived in Sect. 4.2.

3. Retrieve a1, s1 and a2, s2. Decide whether or not an interaction can take place.

d < s1 + s2   ⇒ 1 and 2 collide.

d < (m1/M)2/5a1 ⇒ 1 can perturb orbit 2.

d < (m2/M)2/5a2 ⇒ 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 ,.

(5)

(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.

(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 ω2/ω1. Initialize the loop:

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

We use the notation and for the floor and the ceiling function. The time Tsim to be simulated sets an upper bound for the solution:

We then test the points with coordinates

if x>xmax⇒ next pair

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

5. The solution is

(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

(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 momentum1 calculated from the basic conservation laws:

 (material volume),

m = m1 + m2 (mass),

 (center-of-mass 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 {t0, a, c, s, m, r, L, є, ω} of the new particle.

5. Sort the list of new particles by increasing acs.

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

Go to the next time-step.

thumbnail 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 , 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.

thumbnail 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 (TcollT), 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 s1 + s2. We assume that the planet radii sj are many orders of magnitude smaller than the orbital radii aj. Then, for an angle I between the orbital planes of larger than (s1 s2)/(a1 a2), 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 (s1 + s2)/ about the “mutual nodal line”, all points in orbit 1 are separated by more than s1 + s2 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 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).

thumbnail 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 L1, є1 and L2, є2, respectively. The line of intersection of the two orbital planes, or the nodal line, can be found as follows. Because the angular momenta L1 and L2 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 r1 and r2 of intersection with the nodal line. The tangent lines point in the direction of the velocities υ1 and υ2 at r1 and r2. 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 IK/L1L2, is not much larger than (s1 + s2)/(a1 + a2), 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 r0 to r1 . In the algorithm, r0 is the particle’s creation point and r1 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 eccentricities2. 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 r0 minus the triangle between −aє, 0, and r1 in Fig. 11a. Therefore,

Kepler’s second law therefore implies that the time it takes a body to move from r0 to r1 is equal to

(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

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

(11)

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 so-called g-function:

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 t1 in the end-point gives another equation:

(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 E1E0 that satisfies Eqs. (11) and (12) gives the time to get from r0 to r1.

6 Calculating the time to collision

We consider two planets 1 and 2, with a MOID of less than s1 + s2, and we want to determine the time at which the planets collide. To this end, let r1 be the point on orbit 1 with minimal distance to orbit 2, and r2 the corresponding point on orbit 2 that is closest to r1. 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 r1, and be the first time planet 2 passes r2. 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 dt1 and dt2 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

(13)

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 s1 + s2. For this, we introduce the difference vectors

We note that dt1 and dt2 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 (r2r1) • υ1 = (r2r1) • υ2 = 0, this is equivalent to

Using υ12 = υ2υı, this can be further simplified to

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

(14)

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 T1 > T2. 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.

thumbnail 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 = T1/T2 of orbital periods. The center of the intersection with the horizontal axis is , where 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:

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 qn are all integer linear combinations of p, q, and are therefore elements of ℤ-span{p, q}. When we write qn = knplnq, then the rationals ln/kn are precisely the successive convergents of the continued-fraction representation.

Next, we define the bases {bn, bn+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 (xn, yn) in each base are defined by

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

Consider the following bases, which are composed of the successive even- and odd-numbered 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 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 ℤ2. For this reason, we must instead search in the original bases {bn,bn+1} inside the parallelogram of the strip between y′ = 0 and the y′-value where the bottom of the strip intersects span{bn+2}. This is the top-right point and the bottom-right 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 checked3.

When the calculated orbital periods T1 and T2 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 k1/2, while the number of iterations grows as log k with the solution k.

thumbnail 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

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 t1 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 dtTprec. There are N2 terms in the system of differential equations and there are N2ϵ/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 rj 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 si + sj 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.

thumbnail 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 ℕ2 between two parallel lines that is closest to the origin. The exact solution uses the continued-fraction representation of the ratio Ti/Tj 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(N2s/a) and the total number of collisions is O(N2s2/a2), resulting in an algorithmic efficiency of O(N4s3/a3). This may be compared to the efficiency O(N4/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 < 106, 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:

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

(A.1)

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 Laplace-Runge-Lenz vector is proportional to the dimen-sionless 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:

(B.1)

The Equation for the energy is called the vis-viva 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

(B.2)

The formula for the radial distance in terms of the parameters is

(B.3)

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:

(B.4)

and, using Eqs. (B.2) and (B.3),

(B.5)

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,

(B.6)

and

(B.7)

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 = L3 can give us the angle I. The component L1 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 r2r1. 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 semi-transverse axis of the hyperbola. It is a characteristic distance where the gravity between the two bodies becomes noticeable given the relative velocity. The semi-conjugate axis b of the hyperbolic orbit is equal to the impact parameter. The semi-focal 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 center-of-mass 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 collision-time algorithm. It first finds the solution by considering all integers k, and then uses the continued-fraction 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

  1. Adams, R. A., & Essex, C. 2021, Calculus A Complete Course, 10th edn. (North York, Ontario: Pearson Education) [Google Scholar]
  2. Aliberti, G. D. 2022a, Astrophysics Source Code Library [record ascl:2211.002] [Google Scholar]
  3. Aliberti, G. D. 2022b, Bachelor thesis, TU Delft, The Netherlands [Google Scholar]
  4. Baraff, D., & Witkin, A. 1992, SIGGRAPH Comput. Graph., 26, 303 [CrossRef] [Google Scholar]
  5. Barnes, J., & Hut, P. 1986, Nature, 324, 446 [NASA ADS] [CrossRef] [Google Scholar]
  6. Barnes, J. E. 1990, J. Comput. Phys., 87, 161 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bentley, J. L. 1975, Commun. ACM, 18, 509 [CrossRef] [Google Scholar]
  8. 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]
  9. Burtscher, M., & Pingali, K. 2011, GPU Computing Gems Emerald Edition (Amsterdam: Elsevier) [Google Scholar]
  10. Dehnen, W., & Read, J. 2011, Euro. Phys. J. Plus, 126, 1 [NASA ADS] [CrossRef] [Google Scholar]
  11. Diserens, S., Lewis, H. G., & Fliege, J. 2020, J. Space Safety Eng., 7, 274 [NASA ADS] [CrossRef] [Google Scholar]
  12. Goldstein, H. 1964, Classical Mechanics, Ninth dover printing, Tenth gpo printing edn. (New York: Dover) [Google Scholar]
  13. Greengard, L. 1990, Comput. Phys., 4, 142 [NASA ADS] [CrossRef] [Google Scholar]
  14. Gronchi, G. F. 2005, Celest. Mech. Dyn. Astron., 93, 295 [Google Scholar]
  15. Hamada, T., Nitadori, K., Benkrid, K., et al. 2009, Comput. Sci., 24, 21 [Google Scholar]
  16. Hedo, J., Ruiz, M., & Pelaez, J. 2018, MNRAS, 479, 3288 [NASA ADS] [CrossRef] [Google Scholar]
  17. Hoots, F. R., Crawford, L. L., & Roehrich, R. L. 1984, Celest. Mech., 33, 143 [NASA ADS] [CrossRef] [Google Scholar]
  18. JeongAhn, Y., & Malhotra, R. 2017, AJ, 153, 235 [NASA ADS] [CrossRef] [Google Scholar]
  19. Khinchin, A. Y. 1964, Continued Fractions (USA: University of Chicago Press) [Google Scholar]
  20. Manley, S. P., Migliorini, F., & Bailey, M. E. 1998, A&AS, 133, 437 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Meagher, D.J.R. 1980, Octree Encoding: A New Technique for the Representation, Manipulation and Display of Arbitrary 3-D Objects by Computer, IPL-TR-80-111 [Google Scholar]
  22. Meagher, D.J.R. 1982, Comput. Graph. Image Process., 19, 129 [CrossRef] [Google Scholar]
  23. Milisavljevic, S. 2010, Serbian Astron. J., 180, 91 [NASA ADS] [CrossRef] [Google Scholar]
  24. Murray, C., & Dermott, S. 2009, Solar System Dynamics (New York: Cambridge University Press) [Google Scholar]
  25. Öpik, E. J. 1951, Proc. R. Irish Acad. Sect. A Math. Phys. Sci., 54, 165 [Google Scholar]
  26. Rockett, A., & Szüsz, P. 1992, Continued Fractions (Singapore: World Scientific) [CrossRef] [Google Scholar]
  27. Rokhlin, V. 1985, J. Comput. Phys., 60, 187 [NASA ADS] [CrossRef] [Google Scholar]
  28. Savransky, D., Cady, E., & Kasdin, N. J. 2011, ApJ, 728, 66 [NASA ADS] [CrossRef] [Google Scholar]
  29. Schouten, A. 2022, Bachelor thesis, TU Delft, The Netherlands [Google Scholar]
  30. Segan, S., Milisavljevic, S., & Marceta, D. 2011, Acta Astron., 61, 275 [NASA ADS] [Google Scholar]
  31. Soliman, P. 2022, Bachelor thesis, TU Delft, The The Netherlands [Google Scholar]
  32. Warren, M., & Salmon, J. 1993, in A parallel Hashed Oct-Tree A-Body Algorithm (USA: ACM), 12 [Google Scholar]
  33. Wizniowski, T., & Rickman, H. 2013, Acta Astron., 63, 293 [NASA ADS] [Google Scholar]

1

As pointed out by the referee, orbital angular momentum is not strictly conserved, because it is transferred into spin for oblique collisions: this spin becomes S = S1 + S2 + (m1m2/m)d × u, provided

2

The need for the following calculations was pointed out by Soliman (2022).

3

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

Table 1

List of symbols and notation.

Table 2

Algorithmic efficiency.

All Figures

thumbnail Fig. 1

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

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

In the text
thumbnail Fig. 3

Pairs vs. particle radius s for a homogeneous disk with aamax = 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.

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

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

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

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

In the text
thumbnail 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 (TcollT), 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
thumbnail 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.

In the text
thumbnail 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 = T1/T2 of orbital periods. The center of the intersection with the horizontal axis is , where 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}

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

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

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.