2 The method

2.1 Assumptions

The algorithm is designed to work with two arbitrary elliptical orbits. Due to its nature some classes of problem (i.e. objects in secular, mean-motion and Kozai resonances, including extreme cases of (e, $\omega$)-coupling) may be inappropriate for this approach. The basic assumptions involved are:

1.

each orbit has fixed semi-major-axis (a), eccentricity (e) and inclination (i);

2.

the rate of variation of the argument of perihelion ($\omega$)is greater than that of any of the other elements;

3.

the variations of the respective arguments of perihelia are not correlated.

In addition, we assume that the objects are of small size compared to their orbits, and that there is no significant gravitational focusing. The algorithm is readily extended to incorporate the latter, as we show in Sect. 3.

2.2 Algorithm

Let us suppose that each orbit has semi-major axis, eccentricity and inclination (a₁, e₁, i₁) and (a₂, e₂, i₂) respectively. The longitudes of the ascending nodes ($\Omega_1, \Omega_2$ respectively) are assumed to be fixed, as are the inclinations, and hence so too is the mutual inclination $\Delta i$ of the two orbital planes. We allow the arguments of perihelion of each object, $\omega_1$ and $\omega_2$, to take all possible values in the range ($-\pi, \pi$). This defines a two-dimensional phase space in which the differential collision probability at each point, $P=P(\omega_1,\omega_2)$, can be evaluated. The total collision probability is then

$$P_{\mbox{\scriptsize tot}} = \int_{-\pi}^{\pi} \int_{-\pi}^{... ...(\omega_{1},\omega_{2}) \, {\rm d}\omega_{1} {\rm d}\omega_{2}.$$ (1)

Over most of this phase space the collision probability $P(\omega_1,\omega_2)$ is zero because interactions can only occur when the mutual nodes are sufficiently close. Without loss of generality the coordinate system can be chosen so that the origin lies at the focus of the two orbits, the (x,y)-plane coincides with the orbital plane of the first orbit, and the x-axis lies along the line of intersection of the two orbital planes (Fig. 1). In each case $\omega_1$ and $\omega_2$are measured in the orbital plane from the mutual node.

 

Figure 1: Intersection of orbital planes, illustrating the mutual inclination, $\Delta i$, and the definition of the argument of perihelion, $\omega$, of one orbit measured with respect to the mutual line of nodes

In general, the heliocentric distances r of the ascending or descending nodes are given by

$$r=\frac{a(1-e^{2})}{1\pm e \cos{\omega}}$$ (2)

and collisions may occur at the ascending node when

$$\frac{a_{1}(1-e_{1}^{2})}{1+e_{1} \cos{\omega_{1}}}=\frac{a_{2}(1-e_{2}^{2})}{1+e_{2} \cos{\omega_{2}}}$$ (3)

with a similar expression

$$\frac{a_{1}(1-e_{1}^{2})}{1-e_{1} \cos{\omega_{1}}}=\frac{a_{2}(1-e_{2}^{2})}{1-e_{2} \cos{\omega_{2}}}$$ (4)

for collisions at the descending node. Choosing the ascending node and replacing a(1-e²) with p, and solving for $\cos \omega_1$ we obtain:

$$\cos{\omega_{1}} = \frac{p_{1}-p_{2}}{p_{2} e_{1}} + \frac{p_{1} e_{2}}{p_{2} e_{1}} \cos{\omega_{2}}.$$ (5)

This relation takes the general form  

$$\cos{\omega_{1}} = A + B \cos{\omega_{2}}$$ (6)

where

$$A=\frac{p_{1}-p_{2}}{p_{2} e_{1}}$$ (7)

and  

$$B=\frac{p_{1} e_{2}}{p_{2}e_{1}}\cdot$$ (8)

This can be differentiated to give  

$$\frac{{\rm d}\omega_{1}}{{\rm d} \omega_{2}} = \frac{B \sin{\omega_{2}}}{\sqrt{1-(A+B \cos{\omega_{2}})^{2}}}\cdot$$ (9)

Relations (6) through (8) and corresponding equations for the descending node define the locus of the points in $(\omega_{1},\omega_{2})$-space where collisions may occur, illustrated in Figure 2. Equation (9) is used to determine the tangent and normal to this curve. We note that the curve is symmetrical on both axes, so that looking at one point of intersection there are at most three others with conjugate configurations.

 

Figure 2: The figure illustrates the symmetry of the solutions of Eq. (6) in ($\omega_1,\omega_2$)-space for two typical pairs of orbits. Each pair of orbits has two solutions. The first pair has solutions denoted "1" and "2"; the second has solutions denoted "3" and "4". In each case these correspond to an intersection at the ascending or descending node. The important points to note are that (i) each solution shows four-fold symmetry, and one can be obtained from the other; and (ii) $\omega_1$ is always a monotonically increasing or decreasing function of $\omega_2$ in the principal quadrant, depending on whether the intersection occurs at the the ascending or descending node respectively

We emphasize that, due to the finite size of real objects, collisions may occur not only on this curve but also in its immediate vicinity, over a range $\Delta \omega(\omega_1,\omega_2)$ measured in the $(\omega_{1},\omega_{2})$-plane perpendicular to the curve of intersection. For spherical particles $\Delta \omega$ is proportional to the combined radius of interaction R=R₁+R₂ of the two bodies. Moreover, moving perpendicular to the curve of intersection, it is straightforward to show that the differential collision probability at a distance $\vert\delta\vert<\Delta\omega(\omega_1,\omega_2)$from the collision point $(\omega_{1},\omega_{2})$ is $P(\omega_1,\omega_2)\sqrt{1-(\delta/\Delta\omega)^2}$, so the differential collision probability integrated normal to the curve of intersection is ${\pi\over2}P(\omega_1,\omega_2)\Delta\omega(\omega_1,\omega_2)$.

The total collision probability is then calculated by integrating along the curves of intersection, choosing small line elements of each curve, determining their dimensions and summing their resulting contributions to the total probability. The direction along each curve of intersection is chosen using ${\rm d}\omega_{1} / {\rm d}\omega_{2}$ from Eq. (9). The stepping can be mirrored to allow for the symmetry, but differing geometry (associated with the particles' relative velocities at collision) means that P and $\Delta \omega$may be different for different mirrored points.

Determination of $P(\omega_1,\omega_2)$ and $\Delta \omega(\omega_1,\omega_2)$ requires a model of the encounter geometry. We adopt a simple vector approach in which each body is assumed to move at a constant velocity, denoted ${\vec v}_1$ and ${ {\vec v}_2}$ respectively, during the close encounter. The respective vectors (see Fig. 3) are determined using the classical formulae for heliocentric elliptical motion (Roy 1988):

$$v^2=\mu \left(\frac{2}{r}-\frac{1}{a}\right)$$ (10)

where $\mu=G(M+m)$ and G is the gravitational constant, M is the mass of the central body and m is the mass of the orbiting particle. If we choose $G=4\pi^2$, M=1 and m=0 we have:

$$v=2\pi \sqrt{\left(\frac{2}{r}-\frac{1}{a}\right)}$$ (11)

and for the angle $\phi$ between the velocity and radius vectors at the node

$$\sin\phi= \left( \frac{1+e\cos\omega}{\sqrt{1+e^2+2e\cos\omega}} \right)\cdot$$ (12)

(Note that $\omega$ is measured from the mutual node, and that for simplicity we have dropped subscripts 1 and 2 from Eqs. (10)-(12).) For the first orbit, which lies in the (x,y)-plane, the velocity vector is (Fig. 3):

$${\vec v}_1 = v_1 (\cos{\phi_1},\sin{\phi_1},0)$$ (13)

and for the second, whose plane is inclined at an angle $\Delta i$ to the first, we have

$${\vec v}_2 = v_2(\cos{\phi_2},\cos{\Delta i}\sin{\phi_2},\sin{\Delta i}\sin{\phi_2}).$$ (14)

 

Figure 3: Diagram to illustrate the simplified vector encounter geometry described in the text

In order for a collision to take place, both particles must arrive at the point of intersection of their respective orbits at the same moment, within a small time interval $\Delta t$. We use $\Delta t$ to calculate the collision probability $P(\omega_1,\omega_2)$ at the point $(\omega_{1},\omega_{2})$ on the curve of intersection, and then use the vectors ${\vec v}_1$ and ${ {\vec v}_2}$ to evaluate the range $\Delta \omega$ associated with this point on the curve.

Moving the origin of the coordinate system to the point at which the two paths cross, the motion of the two particles near the mutual node or collision point can be described by

$$b = v_{1} \tau \sin{\theta}$$ (15)

where $\theta$ is the angle between the two position vectors ${\vec r'}_1$and ${\vec r'}_2$ in the new frame. By replacing the maximum value of impact parameter b which produces a collision, namely b=R where R is the radius of interaction for both particles (i.e. R = R₁ + R₂ for spheres of radius R₁ and R₂ respectively), we calculate the maximum value of $\tau$, and hence obtain $\Delta t$ using

$$\Delta t = R/v_1\sin\theta.$$ (16)

Since the temporal window for a collision is $2\Delta t$ per revolution, the differential collision probability per unit time for such intersecting orbits is $P(\omega_1,\omega_2)=2\Delta t/(P_{\mbox{\scriptsize orb,1}} P_{\mbox{\scriptsize orb,2}})$, where $P_{\mbox{\scriptsize orb,1}}$ and $P_{\mbox{\scriptsize orb,2}}$denote the orbital periods of each particle.

We now determine the magnitude of the range $\Delta \omega$perpendicular to the curve of intersection. We define ${\vec s}$ to be the vector that describes the relative motion of the two encounter vectors as $\omega_1$ and $\omega_2$ are varied normal to the curve of intersection, and introduce a separation parameter S, analogous to the impact parameter b, whose value equals R at the extremum $\Delta \omega$. We thus have

$$S = \frac{{\vec s}\cdot ({\vec v}_1 \times {\vec v}_2)}{\vert{\vec v}_1\vert\vert{\vec v}_2\vert}$$ (17)

where

$${\vec s} = r(0, \delta \omega_{1}+\delta \omega_{2} \cos{\Delta i}, \delta\omega_{2} \sin{\Delta i})$$ (18)

and $\delta \omega_1$ and $\delta\omega_2$ denote the respective variations of $\omega_1$ and $\omega_2$ perpendicular to the curve of intersection, with $\vert{\vec \delta \omega}\vert < \Delta \omega$. Here, r is the distance from the Sun where the interaction occurs.

Therefore $\Delta \omega$ is simply calculated from

$$\Delta \omega=\frac{R}{S}$$ (19)

and the total collision probability is evaluated by summing the values of ${\pi\over2}P(\omega_1,\omega_2)\Delta\omega(\omega_1,\omega_2)$along the curve of intersection in the $(\omega_{1},\omega_{2})$-plane for a statistically sufficient number of points.