Simplified derivation of the collision probability of two objects in independent Keplerian orbits

TL;DR

This paper introduces a simplified, symmetric derivation of collision probabilities for two objects in Keplerian orbits, including tangential collisions, improving upon classic linearized methods.

Contribution

A new, more physically motivated derivation that handles tangential collisions and ensures symmetry between colliding bodies, extending previous linearized approaches.

Findings

Provides formulas for collision probabilities in non-tangential and tangential cases.

Regularizes singularities in tangential encounter calculations.

Offers a more straightforward derivation compared to classic methods.

Abstract

Many topics in planetary studies demand an estimate of the collision probability of two objects moving on nearly Keplerian orbits. In the classic works of \"Opik (1951) and Wetherill (1967), the collision probability was derived by linearizing the motion near the collision points, and there is now a vast literature using their method. We present here a simpler and more physically motivated derivation for non-tangential collisions in Keplerian orbits, as well as for tangential collisions that were not previously considered. Our formulas have the added advantage of being manifestly symmetric in the parameters of the two colliding bodies. In common with the \"Opik-Wetherill treatments, we linearize the motion of the bodies in the vicinity of the point of orbit intersection (or near the points of minimum distance between the two orbits) and assume a uniform distribution of impact parameter within the collision radius. We point out that the linear approximation leads to singular results for the case of tangential encounters. We regularize this singularity by use of a parabolic approximation of the motion in the vicinity of a tangential encounter.

Equations

$$\mathbf{\rho }(t)=\mathbf{r}+t\mathbf{v},$$

$${\mathbf{r}}_1+t_1{\mathbf{v}}_1={\mathbf{r}}_2+t_2{\mathbf{v}}_2,$$

$$\mathbf{U}={\mathbf{v}}_1-{\mathbf{v}}_2.$$

$$({\mathbf{r}}_1-{\mathbf{r}}_2)\times \mathbf{U}=(t_1-t_2)({\mathbf{v}}_1\times {\mathbf{v}}_2).$$

$$\frac{d}{dt}|{\mathbf{\rho }}_1(t)-{\mathbf{\rho }}_2(t){|}^2=0,$$

$$\Delta t=\frac{D_{\min }U}{\left|{\mathbf{v}}_1\times {\mathbf{v}}_2\right|},$$

$$\Delta t_{\mathrm{col}}=\frac{\tau U}{\left|{\mathbf{v}}_1\times {\mathbf{v}}_2\right|}.$$

$$P_1=\frac{2\Delta t_{\mathrm{col}}}{T_2}.$$

$$P=\frac{2\Delta t_{\mathrm{col}}}{T_1{T}_2}=\frac{2\tau U}{\left|{\mathbf{v}}_1\times {\mathbf{v}}_2\right|T_1{T}_2}.$$

$${\mathbf{r}}_1+t_1{\mathbf{v}}_1={\mathbf{r}}_2-\mathbf{s}+t_2{\mathbf{v}}_2,$$

$$({\mathbf{r}}_1-{\mathbf{r}}_2+\mathbf{s})\times \mathbf{U}=(t_1-t_2)({\mathbf{v}}_1\times {\mathbf{v}}_2).$$

$$\left|{\mathbf{r}}_1-{\mathbf{r}}_2+\mathbf{s}\right|=\sqrt{{D_{\min }}^2-s^2}.$$

$$\Delta t_{\mathrm{col}}=\frac{\tau U}{\left|{\mathbf{v}}_1\times {\mathbf{v}}_2\right|}\sqrt{1-\frac{s^2}{{\tau }^2}}.$$

$$P=\frac{\pi \tau U}{2\left|{\mathbf{v}}_1\times {\mathbf{v}}_2\right|T_1{T}_2}.$$

$$P=\frac{\pi {\eta }_1}{2\left|{\mathbf{v}}_1\right|T_1{T}_2},$$

$${\eta }_1=\frac{\tau U}{\sqrt{U^2-{\left(U_x\cos {\alpha }_1+U_y\sin {\alpha }_1\right)}^2}}.$$

$${\eta }_1=\frac{\tau U\left|{\mathbf{v}}_1\right|}{\left|\mathbf{U}\times {\mathbf{v}}_1\right|}=\frac{\tau U\left|{\mathbf{v}}_1\right|}{\left|{\mathbf{v}}_1\times {\mathbf{v}}_2\right|},$$

$${\mathbf{\rho }}_i(t)={\mathbf{r}}_i+t{\mathbf{v}}_i+\frac{1}{2}(t-t_i{)}^2\mathbf{g}-\frac{1}{2}{t_i}^2\mathbf{g},$$

$${\mathbf{U}}_0={\mathbf{v}}_1-{\mathbf{v}}_2-(t_1-t_2)\mathbf{g}.$$

$$k=v_2/v_1,-1<k<1,$$

$$t_2=\frac{t_1}{k}-\delta t,$$

$$\delta t\simeq \frac{D_{\min }}{v_2\tan \alpha }.$$

$${\mathbf{r}}_1+t_1{\mathbf{v}}_1-\frac{1}{2}{t_1}^2\mathbf{g}={\mathbf{r}}_2+t_2{\mathbf{v}}_2-\frac{1}{2}{t_2}^2\mathbf{g}.$$

$$({\mathbf{r}}_1-{\mathbf{r}}_2)\times {\mathbf{U}}_0=\frac{1}{2}({t_1}^2-{t_2}^2)\mathbf{g}\times \mathbf{U}+(t_1-t_2)\mathbf{g}\times (t_2{\mathbf{v}}_2-t_1{\mathbf{v}}_1).$$

$$D_{\min }{U}_0=\left[\frac{1}{2}({t_2}^2-{t_1}^2)U+(t_2-t_1)\delta t{v}_1\right]g\sin \alpha ,$$

$$D_{\min }=\frac{1}{2}({t_2}^2-{t_1}^2)g\sin \alpha .$$

$$t_2=\mp \sqrt{\frac{2D_{\min }}{(1-k^2)g\sin \alpha }},$$

$$\Delta t_{\mathrm{col}}=\sqrt{\frac{2(1-k)\tau }{(1+k)g\sin \alpha }}.$$

$$\Delta t_{\mathrm{col}}\simeq \sqrt{\frac{2(1-k)\tau }{(1+k)g\sin \alpha }}\pm \frac{k\delta t}{1+k},$$

$$P=\frac{2\Delta t_{\mathrm{col}}}{T_1{T}_2}=\frac{1}{T_1{T}_2}\sqrt{\frac{8(1-k)\tau }{(1+k)g\sin \alpha }}.$$

Full text

1 Introduction

An accurate estimate of the intrinsic collision probability between two objects moving on independent Keplerian orbits is essential in many topics in planetary system studies: the impact flux of interplanetary dust particles and minor planets on the Earth and other planets (Öpik, 1951; Moses et al., 1999; Ivanov, 2001), evolution of the orbits of a swarm of planetesimals (Greenberg, 1982), planet formation (Wetherill, 1990), dynamical lifetimes of small bodies (Dones et al., 1999), collisional evolution of asteroids (Bottke et al., 2005), the impact hazard of near-Earth asteroids (Harris & D’Abramo, 2015), and collisions amongst artificial Earth-orbiting satellites (Liou, 2006).

In such problems, one usually wishes to quantify the probability of collision within some volume of space which is small compared to the uncertainties of the orbital parameters of the objects that pass through that volume; the objects are assumed to be moving on independent Keplerian orbits about a central body. Öpik (1951) derived an equation for this collision probability for the case when one of the objects is assumed to be in a circular orbit. Wetherill (1967) generalized this solution to two eccentric orbits. Öpik’s and Wetherill’s approaches have two steps in the calculation of collision probability. First, the collision probability for two intersecting Keplerian orbits is calculated. We call this probability $P$, which is a function of the collision radius and the orbital elements defining the shape and the orientation of the two orbits. The two orbits are assumed to be fixed in space, and the mean anomalies are assumed to be independent (i.e., there is no mean motion resonance between the two bodies). Then, over a long period of time, the pair of objects has a well-defined probability of impact near the location where the two orbits intersect or where the distance between the two orbits is small enough for collision to be possible. (The collision radius is inflated from the sum of the physical radii of the bodies, to account for the gravitational focusing effect from their mutual interaction.) Secondly, for ensembles of bodies, the equation for $P$ works as a back-bone to calculate the average collision probability marginalized over all values of the mutual argument of pericenter, $\omega$. For most values of $\omega$, the minimum distance between two orbits is much larger than the distance that allows collision. Thus, the specific ranges of $\omega$ that allow the collision condition is calculated and $P$ over those intervals is integrated over the entire range.

In most previous works in the context of collision rates of asteroids, it has been assumed that apsidal and nodal precession rates are uniform so that $\omega$ is uniformly distributed over its range $(0,2\pi)$ (Wetherill, 1967; Greenberg, 1982; Bottke & Greenberg, 1993). In two recent studies (Vokrouhlický et al., 2012; Pokorný & Vokrouhlický, 2013), the assumption of uniform precession is discarded and the secular evolution along the Kozai-Lidov cycle is adopted to integrate $P$ for high inclination orbits. Rickman et al. (2014) showed that a Monte-Carlo method can be used to integrate $P$ over the precession cycle. Although all the cases described in Rickman et al. (2014) have uniformly distributed angular parameters, their method can be straightforwardly extended to applications with non-uniformly distributed angular parameters. JeongAhn & Malhotra (2015) implemented the Öpik-Wetherill method for the case of non-uniform precession by generating a large number of clones which follow the non-uniform angular distributions. The $P$ values of clones are then summed up to yield the total collision probability on a target object without assuming uniform precession or uniform distribution of $\omega$. This method was used to calculate the seasonal variation of the asteroid impact flux on Mars (JeongAhn & Malhotra, 2015).

We stress that collision probability, $P$, should be interpreted only from the statistical standpoint of collisions of a large population of small bodies in nearly Keplerian orbits. It is not appropriate for predicting a specific impact event in the near future, nor for estimating the long-time impact probability, for a specific pair of objects. For the former purpose, we need the information of passage time which is regarded as a random value in the calculation of $P$. For the latter purpose, the premise of fixed Keplerian orbits is invalid as orbits evolve over time. With the steady state condition for the orbital distribution of numerous colliders, however, we can statistically calculate the impact probability for a given object by integrating $P$ over the distribution of colliders.

In the present work, we first present, in Section 2, a new and simplified derivation of Wetherill (1967)’s collision probability, $P$. We then show that $P$ diverges when the two objects are moving in the same direction, i.e., when the two bodies have a tangential encounter. Several authors (Greenberg, 1982; Vokrouhlický et al., 2012; Rickman et al., 2014) have pointed out this singularity in Öpik (1951) and Wetherill (1967)’s approach, but addressed it only in the averaging of $P$ over the precession cycle. Greenberg et al. (1988) discussed outcomes of tangential encounters but did not calculate the modification of $P$ in such cases. Thus, the singularity problem in $P$ itself has not been solved. Even though near-tangential encounters are not very common, the singularly high collision probability of even just a few such cases can cause non-negligible errors in estimates of collision rates. We examine this singularity problem carefully, and we derive an improved equation for $P$ for tangential encounters which regularizes this singularity (Section 3). We comment on the practical implementation of the formulas for tangential and non-tangential encounters in a general purpose code for collision rates (Section 4), and we describe a case study to illustrate the importance of the correct treatment of tangential encounters (Section 5). We summarize and conclude in Section 6.

2 Collision Probability for Non-tangential Encounters

2.1 Derivation for Intersecting Orbits

Consider the motion of two bodies, body 1 and body 2, whose fixed Keplerian orbits are intersecting each other. Their motions are approximated to be linear near the orbit intersection where collision is possible, as illustrated in Figure 1. Each body’s position can be written as

[TABLE]

where ${\bf r}$ is an arbitrary position vector along the line at time $t=0$ and ${\bf v}$ is the constant velocity in the neighborhood of the orbit intersection point. We can set up a position relation between the two bodies using the point of intersection:

[TABLE]

where subscripts denote the object number (Figure 1). In this configuration each body passes the orbit intersection point at time $t_{1}$ and time $t_{2}$, respectively. The encounter velocity is

[TABLE]

By taking the cross product with the encounter velocity on both sides of Equation 2 and re-arranging the terms, we get

[TABLE]

Any pair of bodies moving with non-parallel constant velocities has a unique time when the distance between the two bodies, $\left|{\bf r_{1}}+t{\bf v_{1}}-{\bf r_{2}}-t{\bf v_{2}}\right|$, becomes minimum. By calling this specific time $t=0$, we get the minimum distance ${D_{\rm min}}=\left|{\bf r_{1}}-{\bf r_{2}}\right|$.

At the minimum distance, the encounter velocity vector is normal to the relative position vector. This is a trivial consequence of the fact that at the minimum mutual distance,

[TABLE]

and the left hand side is proportional to $({\bf r}_{1}-{\bf r}_{2})\cdot\bf{U}$. As this concept is a core part in further derivations, we write it down as a mathematical theorem below for frequent reference.

Theorem 1

Two moving points have their local minimum distance when their relative position vector is normal to their encounter velocity vector.

This is true for any two moving points on regular curves, except when the relative distance of the two moving points vanishes and the direction of their relative position becomes meaningless. We mention that Theorem 1 encapsulates the same concept as in the definition of the so-called “$b$-plane” adopted by Greenberg et al. (1988) in the context of Öpik-Wetherill formulas. The $b$-plane is defined as the plane passing through the target body and being normal to the inbound asymptotic (unperturbed) relative velocity. The impact parameter $b$, which is the magnitude of the projection of the inbound asymptote on the $b$-plane, is related to the minimum encounter distance by a scaling factor that accounts for the gravitational interaction between the two bodies. In fact, the gravitational interaction between the bodies makes the two bodies approach closer than the distance $b$ (this is the so-called “gravitational focusing” effect); we will include this effect when we discuss the collision radius later. We note that gravitational focusing also causes a change in the relative velocity of the two bodies, but this is usually insignificant and is neglected in the calculation of $P$ in the Öpik-Wetherill approach. Throughout this paper, we use “encounter velocity” to refer to the relative velocity which is unperturbed by the mutual gravitational interaction of body 1 and body 2.

From Theorem 1 and Equation 4, we can find that the time interval between body 1 and body 2 passing the orbit intersection point, $\Delta t=\left|t_{1}-t_{2}\right|$, is given by

[TABLE]

where $U$ is the scalar magnitude of ${\bf U}$.

For convenience we choose body 1 as the body passing the intersection point ahead of the other body ($t_{1}<t_{2}$), and illustrate this in Figure 2. Note that both $t_{1}$ and $t_{2}$ can have negative values as one or both bodies may have already passed the orbit intersection point when the two bodies approach each other with the minimum distance, i.e. when $t=0$. In the middle panel of Figure 2, body 2 is $\Delta t$ ahead of the orbit intersection point when body 1 is passing the intersection ($t=t_{1}$). If body 2 were farther away from the intersection at $t=t_{1}$, the pair would have larger minimum distance than ${D_{\rm min}}$. Likewise, if body 2 were located near the orbit intersection point at $t=t_{1}$, the minimum distance would be smaller than ${D_{\rm min}}$. Collision is possible if ${D_{\rm min}}$ is smaller than a specified collision radius $\tau$, i.e. $\Delta t<\Delta t_{\rm col}$, where $\Delta t_{\rm col}$ is given by

[TABLE]

The collision radius $\tau$ would be the sum of the physical radii of the two bodies, if we neglect their mutual gravity. If we wish to take account of their mutual gravity, we can multiply with the gravitational focusing factor.

Thus, if body 2 passes the orbit intersection point later than body 1, collision would occur if and only if, at time $t=t_{1}$, body 2 is within $\Delta t_{\rm col}$ of reaching the intersection point. Conversely, if body 2 passes the orbit intersection point ahead of body 1, collision would occur if and only if body 2 already has passed within the time interval $\Delta t_{\rm col}$ before body 1 passes the intersection point. Therefore, whenever body 1 passes the intersection point, collision would occur if body 2 is within the time interval $2\Delta t_{\rm col}$ of reaching the intersection point.

Considering that body 1 passes the intersection point once per its orbital revolution, the probability $P_{1}$ that it collides with body 2 in any single orbital period is

[TABLE]

Thus the collision probability per unit time, $P$, for body 1 to collide with body 2 is $P_{1}$ divided by the orbital period of body 1,

[TABLE]

Note that $P$ is the collision probability per unit time when two orbits exactly intersect. The general case of non-intersecting orbits will be considered in the next section.

2.2 Extension to Non-intersecting Orbits

Even when two orbits do not intersect, collision is still possible if the distance between the closest points of the two orbits, so-called “Minimum orbit intersection distance” or MOID, is less than the collision radius $\tau$. We note that, in general there exist multiple local minima of the distance between two orbits (Gronchi, 2005), and it is possible that more than one of these could be less than $\tau$. The collision probabilities introduced by those local minima can be easily integrated, if necessary, as they have the same functional form as that of the MOID. In the case study that we describe in Section 5, we include contributions by all the multiple local minima.

The calculation of $P$ for the non-intersecting case and its averaging for MOID$<\tau$ is described well in Greenberg (1982). Here we derive the same equation with an alternative approach. First we note that two skewed lines, i.e. a pair of lines not intersecting nor being parallel, have a unique minimum distance and the vector along the minimum distance is normal to both the lines. This can be generalized as the following theorem.

Theorem 2

Two smooth curves have the local minimum distance along the line orthogonal to the tangents of both curves.

Stated in other words, the normal planes of two smooth curves at a certain point on each curve should coincide with each other if there exists a local minimum distance at the given points.

Then, as before, we linearize the motion of the two bodies in the neighborhood of the MOID. Figure 3 illustrates the situation when the minimum non-zero distance between two orbits, denoted by ${\bf s}$, is along the Z axis. By Theorem 2, both bodies move normal to the Z axis in the vicinity of the MOID location. Without loss of generality, we choose the Y axis to be parallel to the motion of body 2. The path of body 2, shifted by $-{\bf s=-(0,0,s)}$, intersects the path of body 1 at the origin, as in Section 2.1. Therefore, from Equation 2 and 4,

[TABLE]

[TABLE]

As before, we consider body 1 and body 2 to have a minimum distance ${D_{\rm min}}$ at ${\bf r_{1}}$ and ${\bf r_{2}}$ at $t=0$. In Equation 10, the resultant direction of the right-hand side is parallel to the Z axis. As ${\bf U}$ is on the XY plane, ${\bf r_{1}}-{\bf r_{2}}+{\bf s}$ should also lie on the XY plane and its magnitude is

[TABLE]

from Pythagorean theorem (Figure 3). Interestingly, ${\bf U}$ is normal to both ${\bf s}$ and ${\bf r_{1}}-{\bf r_{2}}$; the latter is normal to ${\bf U}$ according to Theorem 1. Therefore, ${\bf r_{1}}-{\bf r_{2}}+{\bf s}$ is also normal to ${\bf U}$, which gives

[TABLE]

for collision radius $\tau$.

For two fixed orbits having $0<s\leq\tau$, we can use Equation 12 without any modification. However, in Monte-Carlo type numerical simulations, the averaged value of $P$, within $0<s\leq\tau$, is more useful as we can get a better estimate of impact flux from the same number of test projectiles.

In the small vicinity where collisions are allowed, $s$ can be assumed to be uniformly distributed. (This follows from the argument that, for a random distribution of lines with fixed directions in space, the fraction of cases with minimum distance less than $s$ is proportional to $s$; this is a good assumption for $\tau\ll$ heliocentric distance.) Averaging $\Delta t_{\rm col}$ within $0<s\leq\tau$, Equation 12 gives $\pi/4$ times $\Delta t_{\rm col}$ of Equation 6. Therefore, the averaged collision probability per unit time when two bodies have $s$ random in the range $(0,\tau)$ becomes

[TABLE]

2.3 Equivalence with Wetherill’s Expression

In Wetherill (1967)’s derivation, his Equation 7 gives the probability of collision per unit time for two bodies on intersecting orbits. Expressing his equation with our notation gives

[TABLE]

where

[TABLE]

In above equation $U_{x}$ and $U_{y}$ are the X and Y components of ${\bf U}$ and $\alpha_{1}$ is the angle between the X axis and the velocity vector of body 1. The Sun is located in the –X direction. The physical meaning of $\eta_{1}$ is the maximum distance from the orbit intersection point where body 1 should be present to allow collision to occur when body 2 is at the orbit intersection point. Derivation of Equation 14 and 15 is lengthy and complicated in Wetherill (1967) and his equation does not look commutative between body 1 and body 2 at first glance. Below we prove that his equation is the same as our simplified derivation of the collision probability per unit time, Equation 13, which is clearly symmetric in body 1 and body 2.

The term, $U_{x}\cos{\alpha_{1}}+U_{y}\sin{\alpha_{1}}$, in the denominator of Equation 15 is the projection of the encounter velocity ${\bf U}$ along the velocity direction of body 1. Therefore, the denominator is equal to the component of ${\bf U}$ normal to ${\bf v_{1}}$. Thus,

[TABLE]

which gives Equation 13 when it is substituted to Equation 14.

3 Collision Probability for Tangential Encounters

3.1 Derivation for Intersecting Orbits

It is evident that the right hand side of Equation 13 is singular when ${\bf v_{1}}$ and ${\bf v_{2}}$ are parallel. This situation is illustrated in Figure 4. The origin of this singularity lies in our linear approximation for the motion of bodies near the collision point. To resolve this singularity, we can use a better approximation of the true motion of the two bodies, namely parabolic motion instead of linear motion,

[TABLE]

where ${\bf g}$ is the gravitational acceleration due to the Sun, assumed to be constant in the vicinity of the encounter. The constant vector ${\bf v}_{i}$ is the velocity of body $i$ at the orbit intersection point, and body $i$ passes the orbit intersection point at time $t_{i}$. At this time, the direction of the velocity vector of both bodies is the same. The encounter geometry is illustrated in Figure 4 in a coordinate system with the origin at the orbit intersection point, the Sun is located on the negative X axis, the orbital poles of the two bodies are aligned to the Z direction, and the motion of the two bodies is approximated as parabolic paths. (Note that the time epochs, $t_{1}$ and $t_{2}$, when body 1 and body 2 pass the origin, have negative values in the case illustrated in Figure 4.) We denote by $\alpha$ the angle between the X-axis and the common direction of the velocity vectors; the range of $\alpha$ is [math] to $\pi$. The time $t=0$ is defined as the epoch when the two bodies approach the minimum distance, ${D_{\rm min}}$; at this time their position vectors are ${\bf r}_{1}$ and ${\bf r}_{2}$ and the encounter velocity is given by the relative velocity at $t=0$,

[TABLE]

We define the outer orbit, the one having smaller curvature, as body 1 and the inner orbit as body 2. Because the outer body should have higher velocity than the inner body at $t=0$, we introduce a velocity ratio

[TABLE]

where the negative values occur when body 2 orbits in the opposite direction of body 1. For simplicity, in the following we consider only the case of positive $k$ for the derivation of $P$. The derivation for the negative $k$ is similar, and we note that the final results (Equation 27, 29, 36, and 37) hold true for both positive and negative $k$. The general derivation with an alternative approach is also provided in the Appendix.

The path of body 2 from the origin to point B is slightly longer than the path of body 1 from the origin to point A in the case illustrated in Figure 4. Because the Y components of the velocities remain constant, the ratio of the travel time for $\overline{OB^{\prime}}$ for body 2 and $\overline{OA}$ for body 1 is equal to the inverse of the ratio of initial velocity, $1/k>1$. (Here $B^{\prime}$ is the location of body 2 having the same Y component as the point A.) The relation between the times $t_{1}$ and $t_{2}$ (when body 1 and body 2 pass the origin) can be written as

[TABLE]

where the small time interval, $\delta t$, for body 2 to go from $B^{\prime}$ to $B$ is given by

[TABLE]

Note that $\delta t$ is positive when $0<\alpha<\pi/2$ and negative for $\pi/2<\alpha<\pi$.

To solve for $t_{1}$ and $t_{2}$, we take a similar approach as in Section 2. Thus, similar to Equation 2, we can set up the following relation for the positions of the two bodies using the orbit intersection point:

[TABLE]

By taking the vector product with the encounter velocity, ${\bf U}_{0}$ (Eq. 18), on both sides of Equation 22, and using ${\bf v_{1}}\times{\bf v_{2}}=0$ (for tangential encounters), we get

[TABLE]

where ${\bf U}={\bf v_{1}}-{\bf v_{2}}$. With the use of Theorem 1, the above equation simplifies to the following,

[TABLE]

Because $\delta t$ is a small value, we neglect the term with $\delta t$ for the moment (but we return to it below). Also, considering that $\left|t_{1}-t_{2}\right|g$ is much smaller than $\left|{\bf v_{1}}-{\bf v_{2}}\right|$, ${\bf U}$ can be approximated as ${\bf U_{0}}$. Then, we find the minimum distance of approach of the two bodies is related to their times of passage at the origin as follows,

[TABLE]

By using $t_{1}\simeq kt_{2}$ and rearranging Equation 25, we obtain

[TABLE]

where the minus sign is for the case when the bodies have already passed the origin when the minimum distance is achieved (as in Figure 4), whereas the plus sign is for the opposite case.

Analogous to the derivation in Section 2, we can define the time interval $\Delta t_{\rm col}=\left|t_{2}-t_{1}\right|$ for collision to occur for a given collision radius $\tau$:

[TABLE]

Thus, if body 2 has passed the origin more than $\Delta t_{\rm col}$ earlier than body 1, body 1 cannot catch up to collide with body 2. On the other hand, unlike Figure 4, in order for the collision to occur before the bodies reach the orbit intersection point, the expected $t_{1}$ should not be more than $\Delta t_{\rm col}$ earlier than the expected $t_{2}$.

If we wish to be more accurate, Equations 20 and 25 give the following better estimate of $\Delta t_{\rm col}$,

[TABLE]

where the sign order is the same as in Equation 26.

In other words, regardless of the choice between the above two equations, if body 2 passes the origin within a time interval $\sim 2\Delta t_{\rm col}$ of when body 1 passes the origin, the collision would occur. (Note that the $\delta t$ term in Equation 28 cancels when we add the time intervals with the plus and minus signs.) Thus, the collision probability per unit time of these two bodies whose velocity vectors are parallel near the collision point is given by

[TABLE]

This equation may appear counter-intuitive at first glance because $P$ decreases as the velocity ratio $k$ approaches 1. It is true that the path length OA and OB increases as $k$ approaches 1. However, as the two velocities become identical, it takes longer for body 1 to catch up with body 2. Consequently, when body 1 is at the origin, body 2 should be located within a smaller distance to have a chance to collide with body 1.

3.2 Extension to Non-intersecting Orbits

Now we consider a more general collision geometry where the orbits do not intersect but their minimum approach distance, $s$, is small enough to allow collision. We define a coordinate system that has the origin as the point on body 1’s orbit where the MOID, $0<s\leq\tau$, occurs. Again, the Sun is located along the –X direction and the orbit of body 1 is on the XY plane. We assume that the variation of the vertical velocity component of body 2 is negligible in the small region where collisions are possible.

According to Theorem 2, at the location where MOID occurs, the position vector of body 2 should be normal to the velocity direction of both bodies, $(\cos{\alpha},\sin{\alpha},0)$. Thus, the position of body 2 at $t=t_{2}$ can be expressed as follows:

[TABLE]

where $\beta$ is the angle between the position vector of body 2 and the XY plane.

Figure 5 shows the orbit of body 1 on the XY plane with the projection of the orbit of body 2 on the XY plane. As before, when body 1 passes the point A, it can barely touch body 2 located at B. Because the path of body 2 shifted by $-{\bf s}$ intersects the path of body 1 at the origin, from Equations 22 and 23, we get

[TABLE]

[TABLE]

As in Section 3.1, we approximate ${\bf U_{0}}\simeq{\bf U}$ and $t_{2}\simeq kt_{1}$. As ${\bf r_{1}}-{\bf r_{2}}\perp{\bf U_{0}}$ (Theorem 1) and ${\bf s}\perp{\bf U}$ (Theorem 2), ${\bf r_{1}}-{\bf r_{2}}+{\bf s}$ is also close to orthogonal to ${\bf U_{0}}$ and ${\bf U}$, therefore we get

[TABLE]

The Z-directional motion of body 2 is negligible in the small vicinity of the MOID location, so body 2 is moving in the plane parallel to the orbital plane of body 1. Thus, as illustrated in Figure 6, ${\bf r_{1}}-{\bf r_{2}}+{\bf s}$ lies in the XY plane and has length

[TABLE]

where $\gamma$ is the angle between the line corresponding to the MOID and the XY plane. By the law of sines, the angle $\gamma$ is given by

[TABLE]

As is the case in Section 3.1, $t_{2}$ is not exactly the same as $t_{1}/k$, but the difference almost cancels with the time difference on the opposite side, i.e. the case where the minimum distance is achieved before the bodies reach the origin. Thus, from Equation 33 through 35, we obtain, for collision radius of $\tau$, the time interval $\Delta t_{\rm col}$ to be

[TABLE]

Assuming a uniform distribution within the parameter space $0<s\leq\tau$ and $-\pi/2\leq\beta\leq\pi/2$, we can obtain the average value of $\Delta t_{\rm col}$. (Care must be taken over the range of $\beta$, as the orbit of the faster body should always be outside of the slower body in the vicinity of the MOID.) With numerical integration we found the average value of the factor $(\sqrt{1-\frac{s^{2}}{\tau^{2}}\sin^{2}{\beta}}-\frac{s}{\tau}\cos{\beta})^{1/2}$ is 0.61. Thus, we obtain the collision probability per unit time,

[TABLE]

In the Appendix, we present fully analytic derivations which prove that Equation 28 and 36 are valid for $\tau/r\ll 1-k$. This condition is well satisfied in solar system applications of the impact rate of asteroids on planets, since planet sizes are much smaller than their heliocentric distances. An example where this condition may be violated is for large planets orbiting close to their host star. For example, in the case of hot jupiter WASP-18b orbiting with semimajor axis of 0.02 AU, the physical radius is 2.7% of the semimajor axis (Southworth et al., 2009). Considering the gravitational focusing factor, we calculate $\tau/r$ of WASP-18b to be equal to $1-k$ when $k=0.7$. For the more common planetary targets with larger semimajor axis and smaller radius than WASP-18b, Equation 28 and 36 are valid approximations unless the encounter velocity is vanishingly small.

4 Transition between Non-tangential and Tangential Encounters

We have seen that the collision probability for non-tangential encounters, Equation 8 and 13, increases to infinity as the velocity vectors become aligned, which is not physical. On the other hand, our formulas for the collision probability for the tangential case, Equation 29 and 37, are independent of the angle, $\theta=\arccos({\bf v_{1}\cdot v_{2}}/v_{1}v_{2})$, between the two velocity vectors. This is also not a good approximation when $\theta$ is not too small. We expect that the collision probability smoothly transitions from the non-tangential formula to the tangential case as $\theta$ decreases below some transition angle, $\theta_{c}$. However, a rigorous calculation of this transition is hard to obtain, for two reasons. First, the collision time interval $\Delta t_{\rm col}$ (Equation 12 and 36) is a function of $s$ and $\beta$, therefore, the transition will be a function of both parameters. Furthermore, for tangential collisions, the body with higher velocity should orbit outside of the body of slower velocity near the MOID location, whereas this restriction is removed for the non-tangential case.

Here, we provide an approximate condition for the transition between the non-tangential and tangential collisions by equating Equation 13 and 37. When the velocity vectors are nearly parallel, we have ${\bf v}_{2}\simeq k{\bf v}_{1}$ and $U\simeq(1-k)v_{1}$. Then, with the approximation $\sin{\theta}\simeq\theta$, we obtain the transition value $\theta_{c}$,

[TABLE]

It is useful to comment on some properties of $\theta_{c}$.

The $\sin{\alpha}$ term in the numerator can be expressed as a function of true anomaly $f$ and eccentricity $e$ of an elliptical orbit,

[TABLE]

The minimum value of $\sin{\alpha}$ is $\sqrt{1-e^{2}}$, and it occurs when $\cos{f}=-e$. For nearly circular orbits, $\sin{\alpha}$ is close to unity. For an eccentric orbit of $e=0.5$, the minimum $\sin{\alpha}$ is 0.866 when $f=\pm 120^{\circ}$, and $\sin\alpha$ can be as small as 0.5 only when the eccentricity is as large as 0.866. 2. 2.

The factor $\sqrt{1-k^{2}}/k$ is monotonically decreasing to zero as $k$ approaches unity. This means that $\theta_{c}$ is smaller for smaller encounter velocities. 3. 3.

The transition angle, $\theta_{c}$, is proportional to the square root of the collision radius divided by the heliocentric distance, $\sqrt{\tau/r}$. Therefore, the transition angle would be significantly larger for close-in large planets, such as hot Jupiters in exo-planetary systems.

As an example, consider the case of the Earth with a circular orbit and an impactor with velocity ratio of $k=0.8$. We adopt the physical radius of the Earth for $\tau$ (neglecting gravitational focusing). Figure 7 shows the collision probability $P$ for this example, calculated with Equations 13 and 37 (solid and dashed lines, respectively). (We plot the product of $P$ and the orbital period of the impactor to provide the dimensionless result, collision probability per orbital revolution of Earth.) For this case, we find $\theta_{c}\simeq 0.26^{\circ}$.

In practical Monte-Carlo type numerical simulations with a moderate number of particles, the number of cases of MOID $<\tau$ is often statistically too small to accurately calculate the integrated impact flux of a projectile population on a given target. For computational efficiency, an artificially enhanced collision radius, $\tau^{\prime}=p\tau$ with $p\gg 1$, can be adopted and the results scaled to the real collision radius (e.g. JeongAhn & Malhotra, 2015). For the case of non-tangential encounters, the number of random orbits having MOID $<\tau$ on the target orbit increases linearly with $\tau$. As the collision probability itself also linearly increases with $\tau$ (Equation 13), the total impact flux is proportional to $\tau^{2}$. Thus, the numerically computed impact flux with an adopted collision radius $\tau^{\prime}=p\tau$ can be rescaled by a factor $p^{-2}$ to obtain the real impact flux.

On the other hand, for the case of tangential encounters, the number of cases of MOID $<\tau$ is proportional to $\tau^{2}$. In practice, however, the exactly tangential cases are of measure zero, so two velocity vectors will not be strictly parallel near the MOID location, and there would be a small, non-zero angle between them. We can identify the near-tangential cases, i.e., those with $\theta<\theta_{c}$, and in these cases the minimum encounter distance should be measured along the ${\bf v}_{1}\times{\bf v}_{2}$ direction. The number of such cases having MOID $<\tau$ is proportional to $\tau$. As the collision probability in these encounters is proportional to $\sqrt{\tau}$ (Equation 36), the integrated impact flux of these cases is proportional to $\tau^{3/2}$. Therefore, in Monte-Carlo type numerical simulations, care should be taken with the artificially enhanced collision radius: the transition angle $\theta_{c}$ should be calculated with the physical radius $\tau$ (including a gravitational focusing factor), not the artificially enhanced collision radius; and the impact flux of the near-tangential collisions should be rescaled by a factor $p^{-3/2}$.

5 Case Study

In this section, we provide a case study of the impact flux of a synthetic population of NEOs on Earth to illustrate that the correct treatment of tangential encounters is crucial to determine the correct impact frequency. The orbital parameters, the mass and the physical radius of Earth are adopted and the impact frequency from $N=5\times 10^{6}$ test particles orbiting in Earth-like orbits is calculated. The test particles’ semimajor axes are chosen randomly from the range 1.1–1.2 AU; this range exceeds ten times the Hill radius of Earth, and avoids the co-orbital region of Earth. The eccentricities are chosen randomly from zero to 0.3, and inclinations from zero to 5 degrees relative to the ecliptic; these ranges are selected for more frequent near-tangential encounters. The arguments of perihelia and the longitudes of ascending node are randomly generated within the range zero to $2\pi$. For statistical analysis, we repeat the same simulation 100 times (making different realizations of the random orbital parameters in each case), and we report mean values and standard deviations of the impact counts below.

Using the code developed by Gronchi (2005), we determine that the total number of orbit crossings having local minimum distances smaller than the collision radius is $39019\pm 220$. For the collision radius, we multiply the physical radius of Earth by the gravitational focusing factor of each test particle. For this synthetic population of impactors with Earth-like orbits, the mean of the gravitational focusing factor is a significantly large value of 2.96.

The total impact frequency calculated from Equation 13 varies greatly from simulation to simulation, from a minimum of 1.6 impacts per year to the maximum of 456 impacts per year; the results are shown in the gray shaded histogram in Figure 8. The mean value and the standard deviation are 9.8 and 47, respectively. This highly variable impact frequency is owed to the divergence of Equation 13 when test particles sporadically encounter Earth at nearly tangential velocities with $\theta\simeq 0$. A careful examination of the distribution of the close encounters reveals that, among the $39019\pm 220$ local minimum distances smaller than the collision radius, only $50\pm 8$ of them have smaller $\theta$ than the critical value (Equation 38), but the smallest $\theta$ case alone contributes more than 50 percent of the total impact frequency in 18 simulations out of 100. However, the corrected total impact frequency, when Equation 13 is replaced by Equation 37 for the near-tangential encounters, is found to be a very stable quantity, $1.39\pm 0.01$ impacts per year (shown in the narrow black histogram on the left of Figure 8). These results show that the incorrect treatment of near-tangential encounters leads to systematically higher impact rate estimates and also greater scatter of the estimates.

In order to demonstrate the validity of our formulas for the impact probability in near-tangential encounters (Equation 37), we sampled 1000 test particle orbits having $\theta<\theta_{c}$, from amongst the above-described synthetic populations of NEOs. We numerically integrated these particles with SWIFT-RMVS3 (http://www.boulder.swri.edu/~hal/swift.html), in a model including the Sun and the Earth with their actual mass and radius. The mean anomalies of the test particles were randomly generated in 10 different realizations, and 10 different numerical simulations are carried out for statistical analysis. We limited the orbit integration time to 10 years to avoid orbit evolution of the particles, and we adopted an integration step size of 1.4 minutes to accurately resolve collisions. From the 10 different numerical integrations, we find an average of $10.4\pm 2.4$ impacts. We can compare this result of the numerical integrations with the expected number of impacts calculated with the unmodified Öpik-Wetherill formulas and with our corrected formulas. For the former, the expected value calculated from Equation 13 is 1802 impacts over ten years; this large value, which is unsurprising due to the nature of the singularity in Equation 13 for near-tangential encounters, is strongly in disagreement with the numerical integration result. For the latter, the expected value calculated from Equation 37 is 8 impacts over ten years, which is within 1–$\sigma$ of the numerical integration result.

6 Summary

The classical method of Öpik (1951) and Wetherill (1967) for calculating collision probabilities of pairs of objects in Keplerian orbits has been widely used in many problems in planetary dynamics. In this paper, we have given a simplified derivation of the backbone of these calculations. Our derivation is easier to understand and to relate to the underlying geometry of collisions of Keplerian orbits. Additionally, our formula for the collision probability per unit time, $P$ (Equation 13), is explicitly commutative between the two colliding orbits (in contrast with the Öpik-Wetherill formulas).

We also derived the collision probability for tangential encounters (Equation 37); this regularizes a singularity in the Öpik-Wetherill formulas. We achieve this regularization by replacing the linear approximation in the vicinity of the collision point with a parabolic approximation of the true motions of the bodies, but otherwise the derivation is similar to the derivation for the non-tangential encounters. In the Appendix, we provide an alternate, fully analytic derivation, which additionally identifies the domain of applicability of the regularized collision probability of tangential collisions. Stated qualitatively, our formulas are valid in the regime in which the collision radius is much smaller than the heliocentric distance and the encounter velocity is not a vanishingly small fraction of the orbital velocity. The quantitative condition is described in detail in the Appendix.

The domains of the non-tangential and near-tangential collisions should be chosen based on the critical angle we derived in Section 4. The additional step needed in computing collision rates is not computationally expensive by virtue of the large increase in available computing power since the 1960s. The neglect of near-tangential encounter cases has the potential to lead to erroneous results; we demonstrated this by an exemplary, although extreme, case of the collision rates on Earth of a population of particles in Earth-like orbits.

We thank the referee, Davide Farnocchia, for helpful comments that improved this paper. This work is supported by UNAM-DGAPA-PAPIIT (grant IN107316). R. M. acknowledges funding from NASA (grant NNX14AG93G) and NSF (grant AST-1312498).

Appendix A Appendix: Alternative Derivation for Tangential Collisions

In this appendix, we present an alternative derivation for $\Delta t_{\rm col}$ for tangential collisions. The notation used is the same as in Section 3, with a few exceptions as noted.

A.1 Derivations for Intersecting Orbits

Let us set up a coordinate system with origin at the point of intersection of the two orbits, the XY plane is the common orbital plane of the two orbits, and the X direction is radially outward from the Sun as in Section 2.1. In general, the minimum distance between the center-of-figure of the two bodies is not zero and it does not occur at the point of intersection of the two orbits. Unlike Section 3.1, let us assume that body 2 passes the origin at time $t=0$ with velocity ${\bf v}_{2}$, and body 1 passes the origin a time $t=\Delta t$ with velocity ${\bf v}_{1}$. (Note that $\Delta t$ here may be negative, if body 1 passes the origin before body 2.) At some time $t=t_{*}$, the two bodies achieve a minimum mutual distance. For tangential encounters, we can write (without loss of generality)

[TABLE]

In the vicinity of the origin, we can approximate the motion of body 1 and body 2 as follows,

[TABLE]

Then the square of the distance between the two bodies can be expressed as a function of time,

[TABLE]

The minimum of the mutual distance occurs at $t=t_{*}$, which can be obtained by the condition $\partial|{\bf r}_{1}-{\bf r}_{2}|^{2}/\partial t=0$.

[TABLE]

where $\varepsilon$ is defined as

[TABLE]

Then the minimum distance occurs at

[TABLE]

where

[TABLE]

It is useful to compute

[TABLE]

[TABLE]

Then, setting $t=t_{*}$ in Eq. A4, we find

[TABLE]

For collision to occur, we must have $|{\bf r}_{1}-{\bf r}_{2}|_{min}\leq\tau$, where $\tau$ is the collision radius. Thus, collision will occur provided $|\Delta t|\leq\Delta t_{\rm col}$, where $\Delta t_{\rm col}$ is given by

[TABLE]

Recall that $\varepsilon\propto\Delta t$ (Eq. A6), therefore Eq. A10 presents a quartic equation for $\Delta t_{\rm col}$, whose exact analytic solution is possible but tedious. Provided that $|\varepsilon|\ll 1$, we obtain the leading order solution,

[TABLE]

A better approximation can be achieved by plugging the first approximation, Equation A11, into Equation A10, to obtain the leading order correction; this yields

[TABLE]

These two equations above are equivalent to Equation 27 and 28.

The condition $|\varepsilon|\ll 1$ requires that the change, $g\Delta t_{\rm col}$, in the heliocentric velocity of the bodies over the time $\Delta t_{\rm col}$ is much smaller than the encounter velocity, $\Delta v=v_{1}-v_{2}$. We note that, using this approximate solution in Eq. A6 we have

[TABLE]

where $\langle v_{T}\rangle=(v_{1}+v_{2})\sin\alpha/2$ is the average transverse velocity of the two bodies at the collision point, and we used $g=GM_{\odot}/r^{2}=v_{c}^{2}/r$ ($r$ is the heliocentric distance at the collision point, and $v_{c}=\sqrt{GM_{\odot}/r}$ is the heliocentric circular velocity). Thus, the condition $|\varepsilon|\ll 1$ is equivalent to

[TABLE]

In words, we can state the condition as: the collision radius (as a fraction of the heliocentric distance) is much smaller than the velocity difference $(1-\left|k\right|)v_{1}$ (as a fraction of the orbital velocity). In physical terms, we can describe this as the condition that the collision radius $\tau$ should be much smaller than the heliocentric distance, when the velocity difference is not a vanishingly small fraction of the orbital velocity.

A.2 Derivations for Non-intersecting Orbits

For the case of non-intersecting orbits, we note that at the location of the MOID, the tangent to each orbit is normal to the relative distance vector, ${\bf s}$. We assume that the minimum mutual distance occurs in the vicinity of the MOID location. Let’s choose as origin that point on body 1’s orbit where the MOID occurs, and let’s choose the plane of body 1 as the XY plane.

Let’s assume that body 2 passes the MOID location, ${\bf s}$, at time $t=0$ with velocity ${\bf v}_{2}$ and that body 1 passes the origin at some time later, $t=\Delta t$ with velocity ${\bf v}_{1}$. Then, for near-tangential encounters we can write,

[TABLE]

We can also express the time-dependent positions of body 1 and body 2 as follows.

[TABLE]

Then the distance between the two bodies is

[TABLE]

The minimum of the mutual distance occurs at $t=t_{*}$, which can be obtained by the condition $\partial|{\bf r}_{1}-{\bf r}_{2}|^{2}/\partial t=0$. First, we derive

[TABLE]

where $\varepsilon$ is given by Eq. A6. Then, we find

[TABLE]

It is useful to compute

[TABLE]

[TABLE]

We use Eq. A23 in Eq. A20 to find that the minimum distance:

[TABLE]

where

[TABLE]

In a similar way to Equation A13, the small parameter $\lambda$ can be understood as

[TABLE]

where $\langle v\rangle=(v_{1}+v_{2})/2$ is the average velocity of the two bodies. Because $s_{x}\leq\tau$, we see that $\lambda$ is of order $\varepsilon^{2}$.

Setting $|{\bf r}_{1}-{\bf r}_{2}|_{min}^{2}=\tau^{2}$, we obtain a polynomial equation for $\Delta t_{\rm col}$:

[TABLE]

This is a polynomial of the 6th degree in $\Delta t$. Keeping only the leading order terms, we need only solve a quadratic in $(\Delta t)^{2}$:

[TABLE]

Only one of the two solutions is physical; we find

[TABLE]

which is the same as Equation 36.