© ESO, 2014

1. Introduction

Impact cratering is one of the most important processes that have shaped the planetary surfaces in the solar system. It plays a fundamental role in age determination thanks to the absolute calibration offered by radiometric dating of lunar surface strata. Thus, the rate of formation of craters is an essential parameter, which depends on the number, mass, and orbit distribution of the impactor population. Thanks to this dependence, studies of impact craters may also provide information on the dynamical evolution of the small body populations in the solar system and the evolution of the solar system itself.

A prime example from recent research concerns the Nice model and its explanation (Gomes et al. 2005) of the late heavy bombardment, an episodic spike in the lunar cratering rate that has been inferred but has also been disputed, based on the available evidence. For instance, Bottke et al. (2012) argued that the impactors should have mostly come from the E belt, a hypothesized continuation of the early asteroid main belt to heliocentric distances less than 2 AU. Gomes et al. (2005) had found earlier that the main belt itself would have provided roughly half the impactors; the other half coming from the trans-planetary, icy planetesimal disk. Morbidelli et al. (2010) applied a modified Nice Model to the asteroid belt and found its contribution to have been small, whereas the majority of the impactors could have been of the latter, cometary type.

Fundamentally, to judge the importance of different source populations, the issue is about their size, physical-dynamical lifetimes, and impact probabilities. The latter form the subject of the present paper. We consider the terrestrial planets as targets and small bodies on planet-crossing orbits as projectiles, and we are interested in the random probability that a projectile, during any of its orbits around the Sun, shall collide with a certain target. In Sect. 2, we describe exactly what these probabilities mean and how they are related to impact probabilities discussed in earlier literature. Our methods to compute them are described in Sect. 3. In Sect. 4, we investigate the behavior of the impact probability in the most important limiting cases, and in Sect. 5, we subject our methods to tests by comparing their results, verifying their general internal agreement, and identifying the few cases, where not all methods work equally well. In Sect. 6, we present results for probabilities of impacts with the Earth and Mars in the form of maps of parametric planes for the projectile orbits. We summarize our conclusions in Sect. 7.

2. Impact probabilities

Since impact probabilities have been under consideration in many previous papers and their actual meaning has not always been obvious, we first make sure that our concept is clearly defined and put it in a proper context by reviewing its background.

We consider projectiles of negligible dimensions and target planets that are massive enough to induce gravitational focusing at close encounters. Apart from this, both targets and projectiles are assumed to move on elliptic, Keplerian orbits around the Sun. Now, consider a planet with a given orbit. A projectile performing one orbital revolution may or may not collide with the planet, depending on the exact values of its orbital elements, and finding the outcome is then a deterministic problem. If the values are not known exactly, as is always the case in reality, we face the well-known problem of estimating impact risks for observed near-earth objects at upcoming encounters with our planet. However, our goal is different, since we are interested in a huge population of projectiles with a wide distribution of elements, and we want to find the random probability of an impact per orbital revolution, which corresponds to this distribution.

We take the semi-major axis of the projectile orbit, its eccentricity, and its inclination relative to the orbital plane of the target as given, and we denote them by a, e, and i, respectively. The three remaining elements (longitude of the ascending node Ω, argument of perihelion ω, and time of perihelion passage τ) are treated as stochastic variables with uniform probability distributions on the intervals [ 0,2π ] for Ω and ω, and [ 0,P_t ] for τ, where P_t is the orbital period of the target: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{t}}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathit{a}}_{\mathrm{t}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\left[}\mathrm{\right(}\mathit{G}\mathrm{(}{\mathit{M}}_{\mathrm{\odot }}\mathrm{+}{\mathit{M}}_{\mathrm{t}}\mathrm{)}{\mathrm{]}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(1)and M_⊙ and M_t are the masses of the Sun and target, respectively, with G being the gravitational constant. The average probability of collision with a specific target during one orbital revolution of the projectile for a uniform distribution of the angles is denoted P_coll(a,e,i).

2.1. Circular target orbit

The computation of P_coll is a long-standing problem. The first analytic results relevant to this case are due to Öpik (1951), who derived the following formula for a circular target orbit: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{coll}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }\sin \mathit{i}}\hspace{0.17em}{\left(\frac{{\mathit{R}}_{\mathrm{t}}}{{\mathit{a}}_{\mathrm{t}}}\right)}^{\mathrm{2}}\frac{\mathit{U}}{{\mathit{V}}_{\mathrm{t}}}\hspace{0.17em}{\left[\mathrm{2}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{t}}}{\mathit{a}}\mathrm{-}\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{t}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(2)The subscript “t” stands for the target, and R_t, a_t and V_t are its radius, orbital semi-major axis and orbital velocity, respectively. The variable U is the relative velocity between projectile and target upon close approach. In the formalism of the circular restricted three-body problem, we have $\begin{array}{ccc}\frac{\mathit{U}}{{\mathit{V}}_{\mathrm{t}}}\mathrm{=}{\mathrm{(}\mathrm{3}\mathrm{-}{\mathit{T}}_{\mathrm{t}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(3)where T_t is the Tisserand parameter of the projectile relative to the target: $\begin{array}{ccc}{\mathit{T}}_{\mathrm{t}}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{t}}}{\mathit{a}}\mathrm{+}\mathrm{2}{\left[\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{t}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\cos \mathit{i}\mathit{.}& & \end{array}$(4)Equation (2) has been widely used in different connections, but it has two potential singularities due to zero denominators. One is for coplanar orbits (i = 0), and the other occurs, if the projectile orbit is perihelion- or aphelion-tangent to the target orbit (the last factor tends to infinity). Moreover, it yields only a first-order approximation to the real case of the low-eccentricity elliptic orbits of the terrestrial planets. On the other hand, the formula is easily adaptable to gravitational focusing, since R_t can simply be replaced by the collisional radius: $\begin{array}{ccc}{\mathit{R}}_{\mathrm{coll}}\mathrm{=}{\mathit{R}}_{\mathrm{t}}{\left(\mathrm{1}\mathrm{+}\frac{{\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}}{{\mathit{U}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(5)where V_e is the escape velocity from the target surface, and U = U(a,e,i) is given by Eqs. (3) and (4).

2.2. Avoiding singularities

According to Kessler & Cour-Palais (1978) and Kessler (1981), the impact probability can be derived for a projectile with fixed (a,e,i) and uniformly changing angles Ω, ω, ℓ (the latter being the mean anomaly) from the corresponding spatial density function S(r,β) that expresses the probability to find the projectile at heliocentric distance r and ecliptic latitude β. However, this density becomes infinite, when r equals q or Q (the perihelion or aphelion distance) and when β = i = 0.

Finding P_coll implies an averaging of S over a torus surrounding the target orbit with radius R_coll, which we refer to as the “target ring”. For a circular target orbit with zero inclination, singularities then occur when using S(a_t,0) as an approximation, in the case of a_t = q, a_t = Q, or i = 0. Steel & Baggaley (1985) showed how this can be avoided by using relevant averages of S instead of the singular spatial density.

Milani et al. (1990) devised a method to regularize the integral involved, which works well even for projectile orbits that both have a very low inclination relative to the target orbit and are nearly tangent to it. Hence, the fundamental problem of the singularities in Öpik’s formula may clearly be dealt with analytically. However, from a pragmatic point of view, modern high-speed computers allow us to circumvent the problem by using the brute-force technique of massive numerical simulations, and this is the approach that we follow.

2.3. Elliptic target orbit

In a study of meteoroid collisions with asteroids, Wetherill (1967) developed a method for calculating the impact probability under the assumption of an elliptic target orbit. If put on a form similar to Eq. (2), the basic formula for the impact probability per orbital revolution of the projectile as a function of the target’s heliocentric distance ρ is $\begin{array}{ccc}{\mathit{P}}_{\mathrm{E}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }\sin \mathit{i}}\hspace{0.17em}{\left(\frac{{\mathit{R}}_{\mathrm{t}}}{{\mathit{a}}_{\mathrm{t}}}\right)}^{\mathrm{2}}\frac{\mathit{U}{\mathit{\rho }}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{\left[}\mathit{G}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\right(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{)}{\mathrm{]}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\hspace{0.17em}{\left[\mathrm{2}\mathrm{-}\frac{\mathit{\rho }}{\mathit{a}}\mathrm{-}\frac{\mathit{a}}{\mathit{\rho }}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}& & \end{array}$(6)where e_t is the eccentricity of the target orbit. We note that this formula assumes the nodal line to cross the target’s orbit at a distance ρ, which implies a condition on the longitude of the ascending node. It includes the probability that the argument of perihelion shall be within the range that allows a close enough approach to the target orbit at this particular distance. The factor four in the denominator can be interpreted as follows. The distance ρ may correspond to either the ascending or descending node, which corresponds to two different Ω values, and the orbital intersection may occur on the outward or inward branch of the projectile orbit, which corresponds to two different ω values. Equation (6) has to be evaluated for all four combinations, and the results have to be added to yield the total impact probability. While U has the same meaning as before, there is the further complication that it now varies along the target orbit and also varies depending on the outward or inward motion of the projectile (see below).

As shown by Wetherill, P_E can be used to derive P_coll(a,e,i) by integration, or averaging, over Ω using a uniform distribution of this angle. In accord with the previous paragraph for each value of Ω, the nodal line crosses the target orbit at two well-defined distances ρ, and the same values are obtained by replacing Ω by Ω + π. Each nodal distance corresponds to two points on the projectile orbit, which yields two arguments of perihelion (ω and − ω); and thus, the integration may extend over [ 0,π ] in Ω, adding the contributions of the four ω values, and finally dividing the result by π. In general, all four values of ω are associated with different values of P_E. Taking the target orbit to define the fundamental plane, the reference direction for measuring Ω may be taken as the target’s perihelion direction, so Ω can be expressed by the true anomaly ν of the target at the nodal line crossing: $\begin{array}{ccc}\mathit{\rho }\mathrm{=}\frac{{\mathit{a}}_{\mathrm{t}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{)}}{\mathrm{1}\mathrm{+}{\mathit{e}}_{\mathrm{t}}\cos \mathit{\nu }}\mathrm{·}& & \end{array}$(7)Integrating the Wetherill formula over ν has been practiced in previous works, but attention needs to be paid to the singularities appearing in Eq. (6). In Appendix A, we analyze in detail how they may be dealt with. However, they are artifacts of the underlying theory, which we prefer to avoid, and for this reason, we have chosen a different approach using Monte Carlo simulations (see Sect. 3). Moreover, as we account for gravitational focusing, there is the extra complication that, when replacing R_t by R_coll in Eq. (6), R_coll(U) is no longer a constant. Nevertheless, we have developed a computer code that performs the averaging by following the scheme described in Appendix A, by a standard Simpson rule integration, and in Sect. 5, we compare its results with those of our new methods.

Greenberg (1982) had already introduced a method that replaced the Monte Carlo type integration used for the averaging in Wetherill’s paper by a standard, numerical evaluation of the corresponding integral. However, he chose to integrate with respect to the target’s mean anomaly instead of the true anomaly, arguing that the integration has to reflect that encounters are more likely in the parts of the orbit, where the target spends more time. While this cannot be disputed, we note that the outcome is a fundamental change of the averaging that involves a non-uniform distribution of Ω.

Our goal is instead to find the average that characterizes a hypothetical population of projectiles with a uniform distribution. This is different from finding long-term temporal averages for individual objects undergoing uniform precession of Ω and ω. These impact probabilities, including the secular changes in e and i, were recently calculated by Pokorný & Vokrouhlický (2013). They should be used in temporal averaging, whenever secular effects are important, but we conjecture that, for instance, for Jupiter family comets that frequently undergo close encounters with Jupiter, the secular evolutions are quenched, and the relevant concern is instead to take account of the jumps in a, e, and i that are caused by the encounters. When doing this, a bias is likely introduced into the ω distribution, invalidating the use of a flat distribution for this angle.

Thus, the realism of the Wetherill-type averages that we consider may be questioned, but the aim of this paper is to present new methods to calculate statistical impact probabilities using arbitrary distributions of Ω and ω. For the purpose of illustrating the inherent features of how these probabilities vary with the orbital elements of the projectiles due to geometric and kinematic effects in the orbital motions, we choose to use flat distributions, but this choice does not follow from the methods themselves. In this paper, we introduce two new Monte Carlo simulation methods that are conceptually different from any previous ways to derive impact probabilities. Each of them exists in two independent variants. All of these avoid the above-mentioned singularities of the analytic formulae.

3. Simulation methods

As already mentioned, there are different combinations of Ω and ω, for which the target and projectile orbits actually intersect. In case of “perfect” timing, a head-on collision would result for those combinations. Considering the finite collisional radius of the target planet at each value of Ω, there are typically four finite ranges of ω, where collisions may occur. The impact probability depends on the width of these ranges and on the chance for the two objects to arrive at the meeting point simultaneously enough, implying a critical range for the perihelion time τ.

Our simulations aim to explore all these parameters and to integrate the relative measure of collision orbits as a function of (Ω,ω) over the whole domain, where this measure is non-zero. We have devised two methods to achieve this, and in the following subsections we will introduce and describe both these methods.

3.1. The Hill sphere method

This is the most straightforward method. For a given combination (a,e,i), we take a vast number of triplets (Ω,ω,τ), using the relevant uniform distributions for the three parameters as described above. Each of these triplets then defines a Keplerian orbit for the projectile, and we trace this motion along with that of the target over one projectile orbit, checking for close encounters. Here, we can choose between two different approaches. One uses the unperturbed Keplerian orbit of the projectile, while the other uses the elliptic restricted three-body problem (Sun-target-projectile). By comparing the two, we have a check on whether the long-range perturbations on the projectile due to the target have any important influence on the results.

It is well known that the impact probabilities with terrestrial planets are no larger than ~10^-8. If we were to estimate them by counting the number of hits onto the collisional sphere with a good statistical accuracy, we would need a sample size of ~10¹¹ or more. We therefore use a modified approach. The real target, which is the collisional sphere, is replaced by a much larger standard target, for which we choose the Hill sphere of radius, $\begin{array}{ccc}{\mathit{R}}_{\mathrm{H}}\mathrm{=}{\mathit{a}}_{\mathrm{t}}{\left(\frac{{\mathit{M}}_{\mathrm{t}}}{\mathrm{3}{\mathit{M}}_{\mathrm{\odot }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}& & \end{array}$(8)Using σ_coll and σ_H for the cross-sections of the collisional and Hill spheres, respectively, the number of hits onto the Hill sphere may be rescaled into the number of actual impacts using the factor, $\begin{array}{ccc}{\mathit{s}}_{\mathrm{coll}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{coll}}}{{\mathit{\sigma }}_{\mathrm{H}}}\mathrm{=}{\left(\frac{{\mathit{R}}_{\mathrm{t}}}{{\mathit{a}}_{\mathrm{t}}}\right)}^{\mathrm{2}}{\left(\frac{\mathit{M}}{\mathrm{3}{\mathit{M}}_{\mathrm{\odot }}}\right)}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}}{{\mathit{U}}^{\mathrm{2}}}\mathrm{·}\right)\mathrm{·}& & \end{array}$(9)Note that s_coll has to be evaluated separately for each hit, using the specific value of U. The reason is twofold, as mentioned above. On the one hand, the radial and transverse components of the target’s heliocentric velocity vary continuously along its orbit, and on the other hand, the radial velocity components of the target and the projectile may be directed either equally or oppositely, leading to two different values of U at each orbital position of the target.

The estimate obtained for the average impact probability is thus: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{coll}}\mathrm{\left(}\mathit{a,e,i}\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathrm{S}}^{-1}\sum _{\mathrm{\Omega }\mathit{,\omega ,\tau }}{\mathit{s}}_{\mathrm{coll}}\mathit{,}& & \end{array}$(10)where the sum is taken over all hits onto the Hill sphere (i.e., successful triplets), and N_S is the total sample size. We finally emphasize the basic assumption used in the above scheme, namely, that there is a flow of projectiles cut by the cross-sectional area from the vantage point of the target and that the density of points of intersection does not vary with distance from the center.

Figure 1 presents a graphic illustration of such a density variation as an example of tests we performed to verify the above assumption. As suggested by this plot, the result of these tests is generally in the affirmative. In particular, we verified that the surface density of intersection points does not vary significantly with distance from the center.

Fig. 1 Intersection of the Earth’s Hill sphere with the b-plane for projectiles having a = 3.5 AU, q = 0.5 AU, and i = 5°, showing the passage points of 3000 objects that entered within the Hill sphere. The center of the Earth is at the origin of the (ξ,ζ) coordinates. These coordinates are described in Sect. 3.2; see also Carusi et al. (1990).

We note a great deal of similarity between this method of calculating average impact probabilities and the technique of simulating actual impact probabilities from orbit integrations by “super-sizing” the terrestrial planets (Horner & Jones 2008; Horner et al. 2009). The idea is the same, and these authors noted that no unforeseen drawbacks had been found in detailed tests. We, however, describe some real drawbacks in Sect. 5.2, while acknowledging that these cases are of very rare occurrence.

3.2. The MOID method

In this method, we construct a sample of projectile orbits in a similar way as described above though with a different accuracy criterion for the sample size. In the first step, we compute the minimum orbit intersection distance (MOID) using the numerical method recently developed by Wiśniowski & Rickman (2013). Let us use ℳ to denote this quantity. Of all the objects in the sample, we retain only those with ℳ ≤ R_coll. These can be treated in two different manners, which we refer to as the MOID-chord method and the MOID-track method.

We first need to describe a preparatory step. Since the target orbits have low eccentricity, the critical ω values for which the orbits actually intersect are typically nearly insensitive to Ω, as illustrated in Fig. 2. Thus, rather than letting the sample of projectile orbits span the whole range of ω, we focus on the intervals that include all the orbits with ℳ ≤ R_coll. These are found empirically as follows.

Fig. 2 Map of the (ω,Ω) plane for projectiles with the same values of a, q, and i as in Fig. 1, showing all the combinations that yielded ℳ ≤ Rcoll with respect to the Earth. The four critical ranges of ω allowing collisions are easily recognized.

First, the critical values ω_c for which ℳ = 0 are found analytically (Greenberg 1982). Then, we choose a very small value of Δω and calculate ℳ for a maximum of 20 000 orbits with the given values of (a,e,i) and Ω at random, keeping ω = ω_c ± Δω. As soon as any of these proves to have ℳ ≤ R_coll, we increase the value of Δω and repeat the experiment. By the time none of the 20 000 orbits have small enough ℳ, we increase the value of Δω once more (the factor applied is 1.1) and keep the resulting range around ω_c, as the interesting one to search for collisions.

When the inclination is substantial, the interesting ω intervals are quite narrow and restricting the search to them means a drastic gain in computer time. On the other hand, if i is small, the four ranges may actually overlap, and we are forced to use their union rather than just adding them together. Another complication arises, when the perihelion distance of the projectile is intermediate between the perihelion and aphelion distances of the target. In this situation, the bands seen in Fig. 2 develop a significant curvature, meaning that Ω plays an important role along with ω to determine the value of ℳ. Then our procedure to select the interesting ω ranges remains the same, but the ranges become much broader. In any case, if A is the measure of the total range of ω considered, a sample size N_S within this range replaces a sample of size N_S·2π/A, which is uniformly distributed over the full range of ω.

Once the ω ranges have been selected, we proceed to creating a random sample of orbits within the proper ω limits, and in this investigation, we choose this sample size to obtain ≃10⁵ members with ℳ ≤ R_coll. Each of these is used to calculate an individual impact probability for a random timing of the encounter, as described below, and these probabilities are combined to get the total, statistical impact probability. With the chosen sample size, the statistical error of this result is about 1.5%.

As a final remark, it is shown below that there are projectile orbits that approach the target orbit to within R_coll at two different points under special circumstances. We refer to this as a double crossing of the target ring. In such cases, separate impact probabilities are calculated for both crossings and finally added together.

In the MOID-chord method, we use the b-plane (Greenberg et al. 1988) to describe the geometry of close encounters. This plane passes through the center of the target, and its normal is parallel to the unperturbed, relative velocity vector of the projectile at closest approach. During the course of an encounter, there is a continuous set of b-planes, which are all centered on the target. For each of these, there is a crossing point, where the projectile would pass through the plane for the relevant timing of the encounter. If the timing is “perfect” (τ = τ_c), this passage occurs at the point on the projectile’s orbit that is closest to the target, or the so-called MOID point. For τ ≠ τ_c, the b-plane is crossed before or after the passage of this point. By plotting all the b-planes on top of each other, we thus obtain a curve including the MOID point and extending on both sides of it. Since the MOID point is the closest possible one along the projectile orbit, this curve must be locally perpendicular to the line joining the MOID point and the origin.

If we let the timing of the encounter vary continuously across the perfect one, the actual b-plane crossing point thus moves along the curve, and this motion is a direct reflection of the target’s orbital motion projected onto the b-plane. The motion is usually considered rectilinear and uniform during the short timing range under consideration, and we follow the same approximation.

Fig. 3 Sketch of the b-plane, where the circular contour marks the intersection with the target’s collisional sphere. The MOID point of the projectile is marked by a small circle on the ξ axis, and the line through points E and E′ (also marked by small circles) is the locus of projectile crossing points for different timings of the encounter.

The geometry is illustrated in Fig. 3. This shows the b-plane with a circular contour marking the intersection with the collisional sphere (radius R_coll). The straight line, offset from the center by a distance equal to the MOID, is the above-described locus of projectile passage points for different timings of the encounter. The coordinates ζ and ξ, parallel and orthogonal, respectively, to this line are the same as those previously used by Greenberg et al. (1988) and Carusi et al. (1990). The interpretation in terms of the MOID was introduced by Valsecchi et al. (2003).

We let the delay of the projectile (or advance of the target) with respect to perfect timing increase in the direction of negative ζ. The timing range for impacts (Δτ) is therefore determined by the length of the chord between the entry point E and the exit point E′. Quantitatively, we have $\begin{array}{ccc}\mathrm{\Delta }\mathit{\tau }\mathrm{=}\frac{\mathrm{2}\sqrt{{\mathit{R}}_{\mathrm{coll}}^{\mathrm{2}}\mathrm{-}{\mathrm{ℳ}}^{\mathrm{2}}}}{{\mathit{V}}_{\mathrm{t}}\sin \mathit{\theta }}\mathit{,}& & \end{array}$(11)where V_t is the orbital velocity of the target and θ is the angle between the vectors V_t and U.

Each of the retained objects has its specific values of ℳ, V_t, and sinθ, which are uniquely defined by Ω and ω. The impact probability per orbital revolution is then $\begin{array}{ccc}{\mathit{P}}_{\mathrm{coll}}\mathrm{=}\mathrm{\Delta }\mathit{\tau }\mathit{/}{\mathit{P}}_{\mathrm{t}}& & \end{array}$(12)with Δτ given by Eq. (11). The average impact probability for the whole sample is: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{coll}}\mathrm{\left(}\mathit{a,e,i}\mathrm{\right)}\mathrm{=}\frac{\mathit{A}}{\mathrm{2}\mathit{\pi }{\mathit{N}}_{\mathrm{S}}}\sum _{\mathrm{\Omega }\mathit{,\omega }}{\mathit{P}}_{\mathrm{coll}}\mathit{,}& & \end{array}$(13)where the sum is taken over all the retained objects with their respective values of Ω and ω, and N_S is the total sample size contained within the restricted ω range of measure A.

Let us consider the MOID-chord method in the case where the target orbit is circular for a moment. There are then four equal ω ranges, where collisions may occur independent of Ω. Their common width is $\begin{array}{ccc}\mathrm{\Delta }\mathit{\omega }\mathrm{=}\frac{\mathrm{2}{\mathit{R}}_{\mathrm{coll}}}{{\mathit{a}}_{\mathrm{t}}}\mathrm{·}\frac{\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right){\cos }^{\mathrm{2}}\mathit{i}}}{\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right)}\sin \mathit{i}}\mathit{,}& & \end{array}$(14)where x = a/a_t. The background to this formula was presented by Valsecchi et al. (2005), and its explicit derivation is briefly explained in Appendix B. Note that the ranges are centered on four critical values ω_c that are situated one in each quadrant. In case, i is not too small; the intervals do not overlap, and the total range has width 4Δω. However, the formula is based on a linear approximation, which holds only in such situations. It should not be used when i → 0, where the total width approaches 2π and not infinity. Our MOID-chord method finds the total range empirically by Monte Carlo sampling and, thus, converges toward the proper value.

Within each ω range, ℳ takes all values from 0 at ω = ω_c to R_coll at the end points. Moreover, in these intervals ℳ scales linearly with ω − ω_c, so a uniform ω distribution yields a uniform distribution of ℳ in each collisional range.

The average impact probability over the whole (Ω,ω) plane is thus $\begin{array}{ccc}{\mathit{P}}_{\mathrm{coll}}\mathrm{\left(}\mathit{a,e,i}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Delta }\mathit{\omega }}{\mathit{\pi }\mathit{/}\mathrm{2}}\mathrm{·}\mathrm{⟨}{\mathit{P}}_{\mathrm{coll}}\mathrm{⟩}\mathit{,}& & \end{array}$(15)where P_coll is given by Eqs. (11) and (12), and the average is taken with a uniform distribution of ℳ ∈ [ 0,R_coll ] and constant values of V_t and sinθ. The result is $\begin{array}{ccc}\mathrm{⟨}{\mathit{P}}_{\mathrm{coll}}\mathrm{⟩}\mathrm{=}\frac{\mathit{\pi }\mathit{/}\mathrm{2}{\mathit{R}}_{\mathrm{coll}}}{{\mathit{P}}_{\mathrm{t}}{\mathit{V}}_{\mathrm{t}}\sin \mathit{\theta }}\mathrm{·}& & \end{array}$(16)Since we have P_tV_t = 2πa_t, we get from Eqs. (14)–(16): $\begin{array}{ccc}{\mathit{P}}_{\mathrm{coll}}\mathrm{\left(}\mathit{a,e,i}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\left(\frac{{\mathit{R}}_{\mathrm{coll}}}{{\mathit{a}}_{\mathrm{t}}}\right)}^{\mathrm{2}}\frac{\mathrm{1}}{\sin \mathit{i}}\frac{\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right){\cos }^{\mathrm{2}}\mathit{i}}}{\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right)}\sin \mathit{\theta }}\mathrm{·}& & \end{array}$(17)Note that the factor sinθ in the denominator may take the value zero in case the vectors V_t and U are parallel, which occurs when the projectile orbit has zero inclination and is tangent to the target orbit at perihelion or aphelion. This singularity of the MOID-chord method comes from the use of Eq. (11), which yields an infinite timing range due to the requirement for a chord of finite length to be traversed with zero velocity. This problem is caused by the assumption of rectilinear and uniform motion of the two objects and is thus an artifact that should be avoided.

Since it can easily be shown that $\begin{array}{ccc}\sin \mathit{\theta }\mathrm{=}\frac{\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right){\cos }^{\mathrm{2}}\mathit{i}}}{\mathit{U}\mathit{/}{\mathit{V}}_{\mathrm{t}}}\mathit{,}& & \end{array}$(18)we finally get $\begin{array}{ccc}{\mathit{P}}_{\mathrm{coll}}\mathrm{\left(}\mathit{a,e,i}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\left(\frac{{\mathit{R}}_{\mathrm{coll}}}{{\mathit{a}}_{\mathrm{t}}}\right)}^{\mathrm{2}}\frac{\mathrm{1}}{\sin \mathit{i}}\frac{\mathit{U}\mathit{/}{\mathit{V}}_{\mathrm{t}}}{\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right)}}\mathit{,}& & \end{array}$(19)thus recovering Eq. (2). This demonstrates the general consistency between the MOID-chord method and the Öpik formula in the case of a circular target orbit.

In the MOID-track method, we use the Keplerian orbits of both target and projectile close to the MOID configuration, thus, avoiding the linearization employed in the MOID-chord method. Just as in the latter, we consider only objects with ℳ ≤ R_coll. The principle is that we know that a central collision occurs for some τ = τ_c, and even without knowing the exact value of τ_c, we seek the interval of τ, enclosing τ_c, where collisions occur. We explore this interval and find its limits numerically by tracing the Keplerian motions of the two objects along their orbits and recording their minimum mutual distance D_min as a function of τ. The limits τ_min and τ_max are the values for which D_min = R_coll.

Note that the interval Δτ = τ_max − τ_min is always finite, even if the vectors V_t and U should be parallel. Thus, the MOID-track method allows to circumvent the above-mentioned problem of the MOID-chord method in cases, when the inclination goes to zero and the projectile orbit becomes tangent to the target orbit simultaneously.

Fig. 4 Minimum mutual distance of one of the objects represented in Fig. 2 from the Earth’s center is shown as a function of the difference between its time of perihelion passage and the critical one yielding the MOID. The values corresponding to the minimum and maximum times of perihelia, for which the mutual distance reaches below the Earth’s collisional radius, are marked.

In Fig. 4, we show a plot of D_min versus τ − τ_c for one of the orbits with ℳ ≤ R_coll in Fig. 2. In this case, since Δτ is found to be close to 400 s, the Earth impact probability per orbit for the orbital elements in question is about 1.3 × 10^-5, according to Eq. (12). In general, the average impact probability for the whole sample (for instance, the one shown in Fig. 2) is found using Eq. (13) as in the MOID-chord method.

4. Impact probability in limiting cases

4.1. Theoretical discussion

Let us now discuss the behavior of the impact probability close to the limits of i = 0 and q = 1 AU for the projectile elements, when considering the Earth as a target, moving on a circular orbit at r ≡ 1 AU with i = 0. We consider projectile orbits with a = 3.5 AU, whose perihelion distances are q ≤ 1 AU and inclinations i ≤ 90°. If the Earth’s collisional radius is R_coll, the three-dimensional target for impacts is a torus with radius R_coll centered on the Earth’s circular orbit. We call this the target ring. Impacts may occur only for projectile orbits cutting through the target ring.

Orbits crossing the central circle need to have one node at exactly r_N = 1 AU. Let q′ denote the value of q expressed in AU. The eccentricity is then $\begin{array}{ccc}\mathit{e}\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{2}{\mathit{q}}^{\mathrm{\prime }}}{\mathrm{7}}\mathit{,}& & \end{array}$(20)and the true anomaly f₁ at r = 1 AU is given by $\begin{array}{ccc}\mathrm{1}\mathrm{+}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{2}{\mathit{q}}^{\mathrm{\prime }}}{\mathrm{7}}\right)\cos {\mathit{f}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{2}{\mathit{q}}^{\mathrm{\prime }}}{\mathrm{7}}\left(\mathrm{7}\mathrm{-}{\mathit{q}}^{\mathrm{\prime }}\right)\mathit{,}& & \end{array}$(21)yielding two solutions ±f₁, where f₁ ∈ [ 0,π ]. These correspond to the node being situated on the outgoing branch (f₁) or the ingoing branch (− f₁). If the node is the ascending one, the critical argument of perihelion is ω_c = −f₁ or + f₁, respectively, and if it is the descending one, then ω_c = π − f₁ or π + f₁, respectively.

Thus, for a given value of q′, which yields a unique value of f₁, there are four values of ω_c, for which the orbit has a node at r_N = 1 AU. These four orbits have ℳ = 0. The ω_c values are symmetrically placed with respect to both the apsidal line (ω = 0 or π) and the perimetric line (ω = ± π/ 2). In the general case, they occupy one quadrant each. There are just a few special cases, where they coincide in pairs, so that the set of four values degenerates into two. These occur for f₁ = 0 (corresponding to q′ = 1) and f₁ = π (corresponding to q′ = 0), when ω_c is either 0 or π, and for f₁ = π/2 (corresponding to q′ = 0.54196), when ω_c equals ± π/2.

Fig. 5 The dependence of the first quadrant value of ωc on the perihelion distance, as expressed in AU, for an assumed semi-major axis of 3.5 AU and a targeted nodal distance of rN = 1 AU.

The dependence of the first quadrant ω_c value (${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{\prime }}$) on q′ is shown in the figure above. Close to the three above-mentioned q′ values, we note that the ω_c values are close to degeneracy: the four values form two tight pairs widely separated from each other. In addition, at q′ = 0.16735 and q′ = 0.87446, we have ω_c = ± π/4 and ±3π/4, representing the contrasting situation of equal spacings.

Each ω_c value is surrounded by an interval (of width Δω), where 0 ≤ ℳ ≤ R_coll, which means that the projectile orbit cuts through the target ring so that impacts may occur. At the limits of such an interval, the orbit is tangent to the ring, and ℳ = R_coll. If the inclination is large, Δω is very small, and the different intervals do not overlap, unless q′ is extremely close to one of the degeneracy values. The situation changes, if we consider the limit of very small inclination.

When the orbit cuts through the ring at a slant angle, it may do so even if the node is far away from the ring. Thus, with decreasing inclination, Δω increases at the same rate for any value of ω_c. As a consequence, the mean collision probability (P_coll) for the particular q′ in question increases, since it is averaged over the full range of 2π in ω. In general, this increase is monotonous but shaped by two phenomena, which occur at specific values of the inclination.

Considering the first quadrant and thus ${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{\prime }}$, the first is when the impacting ω interval extends to the nearest boundary (0 or π/2). As a consequence, the two neighboring ω intervals meet at this boundary and the opposite one. The boundary orbit is then tangent to the target ring at two different places. If the inclination is decreased beyond this critical value, there is a growing ω interval centered on the boundary, within which the orbits cross the ring twice. This causes a more rapid increase of P_coll.

To begin, there is still a range of ω′ for which the orbits do not cross the target ring at all. However, this range is bound to vanish, as the inclination reaches the next critical value, for which the impacting ω′ interval reaches the furthest boundary. Now, impacts may occur at all values of ω, and more and more orbits cross the ring twice, allowing for a sharp increase of the impact probability. The final increase of P_coll as the inclination drops to zero, and all orbits are double ring crossers, is again slower, and finally, all orbits reach two central crossings of the ring.

The two above-mentioned, critical inclinations may be computed using the sine rule of spherical trigonometry from $\begin{array}{ccc}{\mathit{i}}_{\mathrm{c}}\mathrm{=}\arcsin \left\{\frac{{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}}{\sin \mathit{\eta }}\right\}\mathit{,}& & \end{array}$(22)where ${\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$ is the value of R_coll expressed in AU and η is either ${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{\prime }}$ or $\mathit{\pi }\mathit{/}\mathrm{2}\mathrm{-}{\mathit{\omega }}_{\mathrm{c}}^{\mathrm{\prime }}$.

Let us also pay some attention to the special case of q′ = 1 and the neighboring q′ values. If the projectile orbits have a large inclination and if q′ is just slightly smaller than $\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$, P_coll increases sharply with q′, because the omega ranges that allow ring crossing widen rapidly. Their two limits are set by two cases of ring tangency: for the larger value of | ω |, the orbit touches the ring on its outer side, and for the smaller one, the touching occurs on the inner side of the ring. For ${\mathit{q}}^{\mathrm{\prime }}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$, the inner tangency occurs at perihelion, and thus, the range extends to ω = 0.

When q′ takes even larger values, the orbit crosses the target ring for a ω range including zero. This means that the impacting ω ranges have joined in pairs, but there can only be a single ring crossing at large inclinations. Hence, the further evolution of P_coll up to q′ = 1 is set by the further evolution of the combined ω range, which is slow.

Some new features appear, if we consider the case where q′ tends to unity with a very small inclination. These features are illustrated and analyzed in a later Subsection, but we first consider the above-described, general behavior for i tending to zero.

4.2. The i = 0 case

As mentioned above in the MOID method, we focus the search for small MOIDs on more or less narrow ranges of ω, and in Eq. (14) we gave an analytic formula for the width of an interval of ω where collisions may occur in the case of a circular target orbit. We also commented that this equation, like the Öpik formula in Eq. (2), diverges when i → 0 due to the use of a linear approximation. Let us now look closer at this problem.

We performed experiments, applying our MOID method, to numerically calculate the width of the interval Δω that yields ℳ ≤ R_coll for very small inclinations. We took a = 3.5 AU for the semi-major axis of the projectiles, and as discussed above, we used the Earth on a circular orbit with a = 1 AU as target. The Earth was treated as a sphere with a radius R_t = 6371 km, and we applied gravitational focusing, according to Eq. (5), when using both the MOID method and Eq. (14). For the projectile’s perihelion distance, we chose q′ = 0.16735 (eccentricity e = 0.9522), which yields four equidistant values of ω_c at ± π/4 and ± 3π/ 4, and, thus, Δω may reach all the way to π/2 without overlaps between the ranges. In this case, the two values of i_c given by Eq. (22) coincide at i_c = 0.0036 degrees.

We note that the case studied here is characterized by a projectile semilatus rectum of 1 AU equal to the radius of the circular target orbit. This allows collisions at both nodes in the case of ω = π/2 or ω = 3π/ 2, as studied by Hénon (1968), who regarded orbits that may bring an object from one collision to the next (for instance, a space probe launched from and returning to the Earth).

Fig. 6 Comparison of two different calculations of the width (Δω) of the ω range allowing collisions with the Earth by projectiles on low-inclination orbits. For details, see the text. The dashed curve shows results obtained from Eq. (14), and the solid curve shows numerically calculated values using the MOID-track method.

When comparing the results with the predictions by Eq. (14), we notice a perfect agreement for inclinations down to ~0.01 degrees. In Fig. 6, we show only data for the inclination range, where some difference can be seen. The numerical method yields somewhat larger Δω than the analytical formula, and for i = 0.0036 degrees, the limit of Δω = π/2 is reached. Thus, the MOID calculations are in perfect agreement with the prediction of Eq. (22). For even lower inclinations, the formal values of Δω>π/ 2 cannot be used to find the impact probability, so these are not plotted in the diagram. The analytic formula with its (sini)^-1 dependence would predict a further increase toward infinity.

For i<i_c, the orbit cuts through the target ring independent of ω, and as a consequence, any particle passing r = 1 AU spends a non-vanishing amount of time within the ring. This amount increases on the average with decreasing inclination. Moreover, for widening ω ranges, the orbits exhibit double ring crossings. Thus, in the above analysis, we concluded that the impact probability, averaged over all ω, should increase sharply as the inclination decreases. This is verified by our numerical results using the MOID-track method, as shown in Fig. 7 (upper panel). The maximum value in this case is P_coll(3.5 AU,0.9522,0) = 3.0 × 10^-5, while the Öpik formula in Eq. (2) diverges toward infinity.

Fig. 7 Mean impact probability with the Earth, as averaged over all arguments of perihelion, versus inclination as the latter goes to zero. Calculations have been made with the MOID-track method. The two panels correspond to two perihelion distances: q′ = 0.16735 and q′ = 0.3. Critical values of i, as defined in the text, are marked.

In the lower panel, we show the corresponding diagram for q′ = 0.3 as an example of a general case with two different i_c values. We have ${\mathit{\omega }}_{\mathrm{c}}^{\mathrm{\prime }}\mathrm{=}\mathrm{62.25}$ degrees, so the smallest of the two η values in Eq. (22) is 27.75 degrees. The same analysis as above then yields i_c = 0.00553 degrees, while we get i_c = 0.00291 degrees using the other η value (62.25 degrees). We see that the impact probability starts to rise rapidly as i decreases through this range, finally reaching a maximum of P_coll(3.5 AU,0.9143,0) = 2.9 × 10^-5. Generally, the features of both diagrams are the same. The increase of the impact probability with decreasing inclination is accelerated in the range, where i passes the i_c values, and slows down as the inclination drops to zero.

An interesting special case occurs for q′ = 0.54196, when the semilatus rectum of the orbit equals 1 AU and the node crossing occurs only for ω_c = ± π/2. In such a situation, an orbit crossing the target ring always does this twice for any inclination, and collisions may occur at both the ascending and descending nodes. However, the average impact probability is no larger than in the general case, because instead of the four ω ranges that normally allow impacts, there are now only two, within which the probability is twice as high.

4.3. The q′ = 1 case

There is a second possible singularity in Eq. (2), which corresponds to the factor $\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}$ in the denominator of Eq. (14). This tends to zero, when the perihelion distance of a projectile on an exterior orbit approaches the orbital radius of the target, irrespective of the inclination. We have studied this problem by numerical experiments using the MOID method and comparing the results with those predicted by Eq. (14).

To this end, we considered projectile orbits with a = 3.5 AU and q ≤ 1 AU. We fixed the inclination at different values and let q → 1 AU. Apart from this, everything was done as in the previous tests for i → 0. We applied our MOID method to numerically calculate the different ω intervals yielding ℳ ≤ R_coll for values of q′ equal to or slightly less than 1. As explained in Sect. 4.1, we expect to see the joining of two neighboring intervals for ${\mathit{q}}^{\mathrm{\prime }}\mathit{>}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$, and for this reason, we introduce a quantity, which we call Δω₍₂₎ with the following definition. For ${\mathit{q}}^{\mathrm{\prime }}\mathit{<}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$, the two intervals are separate, and they have the same width Δω. The variable Δω₍₂₎ is the width of the combined range, and thus, Δω₍₂₎ = 2·Δω. For ${\mathit{q}}^{\mathrm{\prime }}\mathit{>}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$, the combined range extends between the exterior limits of the joined intervals, and Δω₍₂₎< 2·Δω is the width of this range.

Fig. 8 Comparison of two different calculations of the width (Δω(2)) of two combined ω ranges that allow collisions with the Earth by projectiles on orbits with perihelia near the target orbit for i = 90°. For details, see the text. The dashed curve shows twice the Δω value that is obtained from Eq. (14), and the solid curve shows numerically calculated values using the MOID method.

In Fig. 8, we show Δω₍₂₎ vs. q′ for i = 90°. This case is extreme, yielding the narrowest ω intervals in general, but it serves well to illustrate the phenomena. It is analogous to Fig. 6, although we plotted only Δω in that case, because the joining of intervals meant double ring crossings. We have marked the critical value q′ = 0.999956, where the two ranges join for the current value of ${\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0.000044}$. We note that Δω₍₂₎ passes through a maximum close to 1 − q′ = 0.00003, which is about 0.5 degrees larger than the limit as q′ → 1. By contrast, if we used Eq. (14) to compute 2·Δω for all q′, the result would diverge toward infinity, as indicated by the dashed curve. This analytical prediction works well up to q′ = 0.9999, but it then falls victim to the singularity caused by the factor $\sqrt{\mathrm{2}\mathrm{-}{\mathit{x}}^{-1}\mathrm{-}\mathit{x}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}$ in the denominator.

In the following, we concentrate on a few cases, where i → 0 and q′ → 1 at the same time. We have found it convenient to plot the time spent by the projectile inside the target ring, as computed by the MOID-track method versus q′ and ω (the latter limited to the first quadrant) for constant values of the inclination. In Fig. 9, we show two such 3D plots for i = 1 and 0.1 degrees.

Fig. 9 Time of residence of the projectile (TR) within the target ring for different combinations of q′ and ω – the former is close to 1, the latter close to 0 – for i = 1 deg (upper panel), and i = 0.1 deg (lower panel). The MOID-track method has been used for the calculations. “Slices” showing TR vs. ω for different, constant q′ are shown by different colors, as explained in the diagrams.

If we estimate a maximum for the residence time, T_R, as the time it takes the projectile to move between − R_coll and + R_coll in terms of distance from the orbital plane of the target, we find about 30 000 s for i = 1 deg and about 300 000 s for i = 0.1 deg. In the upper panel, we see that T_R is constrained by this upper limit for the higher of the two inclinations. From the lower panel, we conclude that the lower inclination implies a different constraint, which is likely due to the different curvatures of the orbits.

We expect the integral of T_R with respect to ω to correspond very well to the mean impact probability, as averaged over ω, for the perihelion distance in question, and this expectation has been verified by our tests. This means that observing the areas of the colored “slices” yields an estimate of and an understanding of the reasons for the variation of P_coll with q′.

Fig. 10 Mean impact probability with the Earth, as averaged over all arguments of perihelion, versus perihelion distance q′ as the latter approaches 1. Calculations have been made with the MOID-track method. The two panels correspond to two inclinations: i = 1 deg and i = 0.1 deg. Critical values of q′, as defined in the text, are marked.

In Fig. 10, we show these variations as obtained with the MOID-track method. As seen in Fig. 9 for the smallest q′ values, the impacting ω range does not extend to zero. However, for any given inclination, there is a q′ value, where this range reaches exactly to zero with a vanishing residence time. These critical values are indicated in both panels of Fig. 10. However, as we explain in Appendix C, there are two different behaviors separated by a critical inclination amounting to 0.303 deg. Thus, the two panels in Figs. 9 and 10 refer to different cases.

In terms of geometry in the lowest inclination range (lower panels), there is an orbit with ω = 0 whose perihelion is situated in the target orbital plane slightly closer to the Sun than the inner edge of the ring (${\mathit{q}}^{\mathrm{\prime }}\mathit{<}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$). This orbit touches the ring on both sides of perihelion without penetrating into it. The specific value of q′ in the case shown is about 0.9997. By increasing this value just slightly and keeping ${\mathit{q}}^{\mathrm{\prime }}\mathit{<}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$, there are orbits with ω ≃ 0 that cross the ring twice before and after perihelion. This is observed in the lower panel of Fig. 9, where some of the colored slices (most notably, the blue one) show a bimodal structure with a doubling of T_R at ω ≃ 0.

In the higher inclination range, what was just described cannot happen, and orbits with ω = 0 and ${\mathit{q}}^{\mathrm{\prime }}\mathit{<}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathrm{coll}}^{\mathrm{\prime }}$ do not touch the ring at all. It is only when the perihelion point reaches the inner edge of the ring that a touch occurs. For i = 1°, the value of q′ where this happens is 0.999933, but this is just a single touch. When q′ is further increased, there is a single ring crossing for ω ≃ 0. As seen in the upper panel of Fig. 9, this allows T_R to reach considerable values when q′ is large enough, but there is no anomalous increase like the one seen in the lower panel.

The existence of the double crossings, occurring for the lowest inclinations, has some influence on the shape of the curves in Fig. 10, but this is not dramatic. The two curves are very similar, although they refer to the two different scenarios. They exhibit the same slight drop of P_coll in the highest range of q′ already seen for Δω₍₂₎ in Fig. 8, which shows that this narrowing of the ω range is the reason for the decrease of P_coll.

5. Tests and comparisons

Throughout this section, as in the whole paper, we use a projectile semi-major axis of 3.5 AU for the calculations.

5.1. Check of the two-body approximation

Impact probabilities are usually calculated on the basis of unperturbed, Keplerian motions of both projectile and target. Thus, in principle, the orbital elements used should be valid at the time of the encounter. Using elements with a different osculation epoch, one might possibly introduce a significant error. This could occur, in case the long-range perturbations on the projectile by the target planet would tend to modify the statistical distribution of the close encounter geometries.

Fig. 11 Collision probability with Mars (using its current orbital and physical properties) per orbital revolution as a function of the projectile perihelion distance (qp) divided by Mars’ perihelion distance (qt) according to the Hill sphere method. The probability peaks when this ratio equals unity, and the variation is the same for both dynamical models used: the two-body problem (2BP) and the restricted 3-body problem (R3BP).

Using the Hill sphere method, we decided to test this by comparing the results from the two dynamical models that we used. The two-body problem corresponds to the classical, Keplerian approach, while the restricted three-body problem (Sun-Target-Projectile) includes the long-range perturbations just mentioned. In Fig. 11, we show the results of a check, where impact probabilities were calculated with respect to Mars on its current elliptic orbit for projectiles on Mars-crossing orbits of cometary type. Clearly, the two methods yield equivalent results, thus verifying the validity of the Keplerian approximation in the present context.

5.2. Checks of consistency between the methods

We have made great efforts throughout this work to verify that our Monte Carlo methods agree internally when calculating the average impact probabilities and also that they agree with calculations using Wetherill’s formula in Eq. (6) by including gravitational focusing. The general result is illustrated by Fig. 12, which shows an example using the current orbit and physical properties of Mars for the target. A constant inclination of 1° has been used to calculate the probability as a function of the perihelion distance of the projectile.

Fig. 12 Collision probability with Mars (using its current orbital and physical properties) per orbital revolution as a function of the projectile perihelion distance (qp) divided by Mars’ perihelion distance (qt), according to different methods. The variations are the same to good accuracy for all methods, except in the immediate vicinity of qp = qt.

An excellent agreement between all methods is obvious under almost all the circumstances shown in the figure. The same holds for practically all inclinations including all values larger than 1°. However, for smaller inclinations, the Hill sphere method faces a severe problem. For a given projectile inclination i with respect to the target orbit and the target placed at heliocentric distance r as the projectile passes the same heliocentric distance, it can never be further than rsini from the target orbital plane. If the Hill sphere of the target has a radius exceeding this value, the crossing points cannot fill up the whole Hill sphere cross-section like illustrated in Fig. 1. As a consequence, by erroneously assuming that the cross-sectional coverage is uniform, the Hill sphere method is bound to underestimate the impact probability, since it underestimates the density of crossing points in the close vicinity of the target.

Fig. 13 Collision probability with the Earth (using a circular target orbit) per orbital revolution as a function of the projectile inclination (i), according to different methods. The upper panel shows the divergence between the MOID and Hill sphere methods for i between 0.1 and 1°, and the lower panel shows the divergence between the MOID and Wetherill averaging methods for i< 0.014 deg.

Fig. 14 Collision probability with Mars (using its current orbital and physical properties) per orbital revolution as a function of the projectile perihelion distance (qp) divided by Mars’ perihelion distance (qt), according to different methods. The variations are shown for the immediate vicinity of qp = qt, using an abscissa proportional to log (1 − qp/qt).

In the case under consideration, the target is Mars, and estimating the minimum inclination for the Hill sphere method to work as $\begin{array}{ccc}{\mathit{i}}_{\min }\mathrm{=}\arcsin {\left(\frac{{\mathit{m}}_{\mathrm{M}}}{\mathrm{3}{\mathit{M}}_{\mathrm{\odot }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}& & \end{array}$(23)the mass m_M of Mars yields i_min = 0.27 deg. Using the Hill sphere method for smaller inclinations for Mars as a target leads to erroneous results. For comparison, the minimum inclination for the Earth is i_min = 0.57 deg. In Fig. 13, we illustrate the divergence between Earth impact probabilities calculated by the Hill sphere and MOID methods, as the inclination drops from 1 to 0.1 deg. The problem is seen to start around the predicted value of i, and for inclinations of 0.01 deg or less, the Hill sphere method is off by orders of magnitude. To realize this, compare this to the lower panel, which illustrates the divergence between the results by Wetherill averaging and the MOID method. The MOID curve is the same, as was plotted in Fig. 7 (lower panel). This divergence is only of concern for inclinations below 0.002 degrees, but as expected, the Wetherill curve diverges toward infinity, while the MOID curve approaches a limiting value without any dramatic increase.

Though limited to very low inclinations, this is a serious shortcoming of the Hill sphere method. Thus, like the above-mentioned super-sizing technique (Horner & Jones 2008), it has to be used with care, even though it works well in almost all situations. Returning now to Fig. 12, there is another problematic situation that is not recognized in the plot but is seen in higher resolution for practically equal perihelion distances between target and projectile. An illustration is provided by a special set of calculations, as illustrated in Fig. 14. This plot shows how the average impact probability varies, when the projectile perihelion distance is close to, but less than, that of the target. An inclination of 0.1° has been used.

Here, some interesting features stand out. As we already saw in Sect. 4.3, the Wetherill formula cannot be used arbitrarily close to q′ = 1 in the circular case, and now, we show this also in the case of an eccentric target orbit, which verifies the result found in Appendix A. The impact probability predicted by this formula tends to infinity as q_p/q_t → 1, whereas the true value (given by the MOID methods) is about 1.6 × 10^-7. The reason for this discrepancy is thus understood, so we now discuss the behaviour of the Hill sphere results. The consequence of the low-inclination problem is seen to the left in the diagram, where there is a roughly constant difference between the Hill sphere curve (lower) and the others (higher). This holds true up to q_p/q_t ≃ 0.993, but the Hill sphere results fall increasingly short of those obtained by the MOID methods for larger values (up to q_p/q_t ≃ 0.999992). This phenomenon needs a different explanation.

The reason has to do with the special geometry of encounters occurring under the present circumstances, when the perihelion directions of the target and the projectile nearly coincide and the inclination is very low. Then, the two objects chase each other along nearly parallel trajectories, and a very wide timing range leads to collision, thus increasing the average impact probability. However, this particular phenomenon is lost, when considering the artificially enlarged cross-section of the Hill sphere method. The timing range in that case cannot increase in proportion to the increase of the cross-section because of the divergence of the orbits that follows from their different curvatures.

We illustrate this by Fig. 15. Here, we devised a numerical experiment, where the average impact probability with Mars was calculated by the MOID-track method for a very low inclination. Three different values of q_p/q_t were considered: 0.95, 0.99, and 1.00. We artificially changed the collisional radius as shown on the abscissa of the diagram (R denotes the real value). The ratios of the impact probabilities found with different perihelion distances are plotted vs. R_coll. With q_p = q_t, the impact probability is considerably higher than with a slightly lower q_p, as a consequence of the above-mentioned parallel chase of projectile and target within the collisional torus. However, the diagram shows that the amount of this increase depends strongly on the assumed collisional radius. If we overestimate this radius, we underestimate the q_p = q_t spike of the impact probability.

Fig. 15 Ratios of average impact probabilities with Mars for projectiles having different perihelion distances (in units of Mars’ perihelion distance) plotted versus the articially varied collisional radius of the planet.

Fig. 16 Collision probability with the Earth for very low inclination orbits, as the projectile perihelion distance approaches 1 AU from below. Computations have been made with the MOID-chord and MOID-track methods, and significant differences are seen as the ratio of perihelion distances exceeds 0.9999.

Fig. 17 Impact probability per orbital revolution by the projectile with respect to the Earth (top panel) and Mars (bottom panel), using color-coding, as explained to the right of the diagrams. All calculations were performed using the MOID-track method.

Our final comparison between methods refers to the two variants of the MOID method. In Fig. 14, one can notice slight differences between the MOID-chord and MOID-track results, when q_p/q_t is very close to unity. These are not caused by random noise in the plotted values, but they arise from the fundamental difference that the MOID-track method does not involve any approximation when finding the timing range that allows collisions, while the MOID-chord method assumes uniform, rectilinear motion of the two objects with respect to each other. In Fig. 16, we show the corresponding phenomenon for the case of the Earth as a target on a circular orbit, and we see that it stands out much more clearly. We interpret this to mean that the MOID-chord approximation (also used in the Öpik and Wetherill analyses) fails to some degree in the extreme situation at hand, where a parallel chase of the two objects occurs within the target ring. As a final result, we have shown that all methods with the likely exception of MOID-track have some shortcomings in the extreme situations considered here.

6. Earth and Mars impact probabilities

In view of what was found with regard to the performance of the methods, we selected the MOID-track to calculate the statistical impact probabilities of projectiles with a = 3.5 AU with respect to the Earth and Mars, using their current, eccentric orbits. In both cases, we surveyed the (q′,i) plane, where q′ = q_p/q_t by using a fine grid covering inclinations up to 30°. The two impact probability maps are shown in the two panels of Fig. 17.

With the use of a common color scale for both planets, Fig. 17 illustrates the difference in probabilities between them, which is mainly caused by the Earth’s larger size and mass. Both diagrams show a relatively flat impact probability over most of the parametric plane and the existence of ridges of increased probability associated with the critical cases examined above – that is, the very small inclinations (on top of the diagram) and the nearly coinciding perihelion distances (vertical ridge in the left part). The low-inclination ridge extends up to i ≃ 5°, while the vertical one is much narrower.

The highest probability is found where the two ridges intersect. The height of the ridge at i = 0 is not very sensitive to the perihelion distance, but the height of the q′ = 1 ridge drops with increasing inclination. We do not show here any diagrams with varied semi-major axis, since there is no particular structure associated with that variation, and the changes with a are very small. We finally note that the aim of these diagrams is only to show the general features of the variation of impact probability with projectile orbital elements. They are based on a particular assumption about the distributions of angular elements, in which these are uniform, but real populations of objects in the solar system will not behave according to this assumption. Therefore, the diagrams do not accurately represent the real impact probabilities characterizing small bodies like Jupiter family comets.

7. Conclusions

This investigation considered the statistical impact probability per orbital revolution by the projectile, P_coll(a,e,i) as a function of the projectile orbital elements a, e, i, by assuming uniform distributions of the angular elements. The target was considered to be one of the terrestrial planets, and the projectile orbits were taken as typical of Jupiter family comets for simplicity. Specifically, the Earth and Mars were selected for all the computations.

All our methods to calculate P_coll are numerical. Although the target orbit may be circular or elliptic, our methods deal with the general, elliptic case. The first is an average of the impact probability P_E (Wetherill 1967) over a uniform distribution of the target’s true anomaly, which corresponds to a uniform Ω distribution for the projectile. The second is a Monte Carlo simulation method, which counts the number of Hill sphere passages of a large population of test objects with random combinations of angular elements and scales with respect to the size of the target’s collisional cross-section. The third uses a similar search for random test objects, whose MOIDs with respect to the target are less than the target’s collisional radius, after which the target–object encounter timing range for a collision is computed and compared to the orbital period of the target.

Of these, only the MOID method (specifically, the MOID-track variant) involves no assumptions or approximations regarding the properties of the encounters. Using this method, we explored the vicinity of the singularities appearing in Wetherill’s formula, namely, for projectile elements, i_p → 0 or q_p → q_t. The behavior of P_coll in these situations was seen to be shaped by the appearance of double crossings of the target ring (a torus surrounding the target orbit with a cross-section equal to the

collisional one) for very small values of i_p. For nearly coinciding perihelion distances, the critical feature is instead the amount of overlap between the projectile orbit and the target ring, which governs the range of ω values that allow single or double crossings. The width of this range, and thus P_coll, has a maximum when q_p/q_t is slightly less than unity.

Detailed comparisons of the different methods showed that there is a nearly perfect agreement between all of them in most situations. However, in all the above-mentioned critical situations, the Hill sphere method (and other super-sizing methods based on the same principle) underestimates the statistical impact probability, while the Wetherill averaging yields results tending to infinity.

We thus used the MOID-track method to present the overall picture of statistical impact probability with respect to the Earth and Mars on their current orbits in the form of maps of the (q,i) parametric plane, using a constant semi-major axis of 3.5 AU for the projectiles. The maps are hence useful to provide insight into the general variation of planetary impact probabilities of Jupiter family comets. However, they should not be used for the real population of comets except for an order-of-magnitude estimate, since the angular elements are not uniformly distributed.

Online material

Appendix A: Appendix A

Equation (6) contains two factors raised to the power of − 1, which may take the value zero. The first is sini, which implies that an integration of the formula cannot be performed, if the inclination is zero. The other factor is more benign. It can be written as $\begin{array}{ccc}\sqrt{\mathrm{2}\mathrm{-}\frac{\mathit{\rho }}{\mathit{a}}\mathrm{-}\frac{\mathit{p}}{\mathit{\rho }}}\mathit{,}& & \end{array}$where a and p are the semi-major axis and semilatus rectum of the projectile, and the heliocentric distance ρ of the target is expressed by Eq. (7). We use ${\mathit{P}}_{\mathrm{t}}\mathrm{=}{\mathit{a}}_{\mathrm{t}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{t}}^{\mathrm{2}}\mathrm{)}$ for the semilatus rectum of the target. The reciprocal of the factor under consideration is $\begin{array}{ccc}\sqrt{\frac{\mathit{\rho }}{\mathrm{-}\frac{{\mathit{\rho }}^{\mathrm{2}}}{\mathit{a}}\mathrm{+}\mathrm{2}\mathit{\rho }\mathrm{-}\mathit{p}}}\mathrm{·}& & \end{array}$The non-singular nominator is easily dealt with by numerical integration, so we need to concentrate only on the integral, $\begin{array}{ccc}\mathit{I}\mathrm{=}{\mathrm{􏽉}}_{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}\mathit{\nu }}{\sqrt{\mathrm{-}\frac{{\mathit{\rho }}^{\mathrm{2}}}{\mathit{a}}\mathrm{+}\mathrm{2}\mathit{\rho }\mathrm{-}\mathit{p}}}\mathit{,}& & \end{array}$(A.1)where ν is the true anomaly of the target. The denominator can be written as $\sqrt{\mathrm{(}\mathit{\rho }\mathrm{-}\mathit{q}\mathrm{)}\mathrm{(}\mathit{Q}\mathrm{-}\mathit{\rho }\mathrm{)}\mathit{/}\mathit{a}}$, where q and Q are the perihelion and aphelion distances of the projectile, respectively. The integrand may thus become singular, only if there is a value of ν such that ρ = q or ρ = Q. In the present case, we are only interested in the occurrence of ρ = q, when the perihelion distance of the projectile is spanned by the target’s orbit. By ν_c we denote the positive true anomaly of the target, for which ρ = q. If | ν | <ν_c, collisions are excluded, and we may neglect this part of the integral, considering only $\begin{array}{ccc}{\mathit{I}}_{\mathrm{c}}\mathrm{=}\sqrt{\mathit{a}}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathrm{c}}}^{\mathrm{2}\mathit{\pi }\mathrm{-}{\mathit{\nu }}_{\mathrm{c}}}\frac{\mathrm{d}\mathit{\nu }}{\sqrt{\mathrm{(}\mathit{\rho }\mathrm{-}\mathit{q}\mathrm{)}\mathrm{(}\mathit{Q}\mathrm{-}\mathit{\rho }\mathrm{)}}}\mathit{.}& & \end{array}$(A.2)The non-singular factor $\mathrm{1}\mathit{/}\sqrt{\mathit{Q}\mathrm{-}\mathit{\rho }}$ can again be disregarded, so we concentrate on the expression $\begin{array}{ccc}& & \int \begin{array}{c}\mathrm{2}\mathit{\pi }\mathrm{-}{\mathit{\nu }}_{\mathrm{c}}\\ {\mathit{\nu }}_{\mathrm{c}}\end{array}\frac{\mathrm{d}\mathit{\nu }}{\sqrt{\mathrm{\left[}\mathit{\rho }\mathrm{\right(}\mathit{\nu }\mathrm{)}\mathrm{-}\mathit{\rho }\mathrm{(}{\mathit{\nu }}_{\mathrm{c}}\mathrm{\left)}\mathrm{\right]}}}\mathrm{=}\\ & & \frac{\mathrm{1}}{\sqrt{{\mathit{e}}_{\mathrm{t}}{\mathit{p}}_{\mathrm{t}}}}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathrm{c}}}^{\mathrm{2}\mathit{\pi }\mathrm{-}{\mathit{\nu }}_{\mathrm{c}}}\frac{\sqrt{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{e}}_{\mathrm{t}}\cos \mathit{\nu }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{e}}_{\mathrm{t}}\cos {\mathit{\nu }}_{\mathrm{c}}\mathrm{)}}}{\sqrt{\cos {\mathit{\nu }}_{\mathrm{c}}\mathrm{-}\cos \mathit{\nu }}}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathrm{·}\end{array}$

We note that the above manipulation fails in case e_t = 0. Hence, the integration cannot be performed, if the target orbit is circular with a_t = q. Excluding this case. the non-singular nominator can again be disregarded, so we concentrate on the integral, $\begin{array}{ccc}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathrm{c}}}^{\mathrm{2}\mathit{\pi }\mathrm{-}{\mathit{\nu }}_{\mathrm{c}}}\frac{\mathrm{d}\mathit{\nu }}{\sqrt{\cos {\mathit{\nu }}_{\mathrm{c}}\mathrm{-}\cos \mathit{\nu }}}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathrm{c}}}^{\mathrm{2}\mathit{\pi }\mathrm{-}{\mathit{\nu }}_{\mathrm{c}}}\frac{\mathrm{d}\mathit{\nu }}{\sqrt{{\cos }^{\mathrm{2}}\frac{{\mathit{\nu }}_{\mathrm{c}}}{\mathrm{2}}\mathrm{-}{\cos }^{\mathrm{2}}\frac{\mathit{\nu }}{\mathrm{2}}}}\mathrm{·}& & \end{array}$Substituting $\begin{array}{ccc}\cos \frac{\mathit{\nu }}{\mathrm{2}}\mathrm{=}\cos \frac{{\mathit{\nu }}_{\mathrm{c}}}{\mathrm{2}}\cos \frac{\mathit{\alpha }}{\mathrm{2}}\mathit{,}& & \end{array}$(A.3)the integral is transformed into $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}\mathit{\alpha }}{\sin \frac{\mathit{\nu }}{\mathrm{2}}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\frac{\mathrm{d}\mathit{\alpha }}{\sqrt{\mathrm{1}\mathrm{-}{\cos }^{\mathrm{2}}\frac{{\mathit{\nu }}_{\mathrm{c}}}{\mathrm{2}}{\cos }^{\mathrm{2}}\frac{\mathit{\alpha }}{\mathrm{2}}}}\mathit{,}& & \end{array}$which is proper, when cos²(ν_c/2) ≠ 1, or ν_c ≠ 0. Our method to integrate the Wetherill formula starts by flagging an error, if i = 0 or q_t = q (i.e., ν_c = 0). Otherwise, we either have q_t>q, in which case a Simpson rule integration of the Wetherill expression is performed without problems (as long as Q_t<Q), or q_t<q, in which case the Simpson rule integration is easily carried out using the variable substitution of Eq. (A.3). We finally note that a different analytic proof of the regularity of the integral was presented by Wetherill (1967).

In summary, we have found the following cases, when the Wetherill integral diverges due to singularities:

• co-planar orbits of the target and projectile (i = 0);

• circular target orbit at the perihelion distance of the projectile;

• elliptic target orbit with the same perihelion distance as the projectile.

Although we did not consider it explicitly here, the third case has an analog for a different type of projectile orbit, if the two aphelion distances are equal.

Appendix B: Appendix B

Valsecchi et al. (2003) present an analytic theory of close encounters by extending the work of Öpik (1976), in which the effects of a close encounter of a small body with a planet moving in a circular heliocentric orbit is described making use of two b-plane coordinates, namely, ξ and ζ. From these quantities, the latter is related to the time difference with which the small body and the planet arrive at their closest approach, while the modulus of the former is the MOID (Valsecchi et al. 2003, Sect. 3.1). It is clear that an encounter within a given distance d can only take place if the MOID is smaller than d; that is, if | ξ | ≤ d.

Building on that theory and on refined expressions given in Valsecchi (2006), it is possible to express small displacements of ξ and ζ from the origin in terms of the partial derivatives of ξ and ζ with respect to the orbital elements a,e,i,ω,Ω, and f of the small body at collision with the planet (neglecting the gravitational influence of the latter). Valsecchi et al. (2005) give a suitable procedure for this computation.

In terms of these derivatives, keeping a,e,i constant and d suitably small, we have: $\begin{array}{ccc}\mathrm{\left|}\mathit{\xi }\mathrm{\right|}& \mathrm{=}& \left|\frac{\mathit{\partial \xi }}{\mathit{\partial \omega }}\right|\mathit{\delta \omega }\mathrm{=}\sin \mathit{i}\mathrm{\left|}\sin \mathit{\phi }\mathrm{\right|}\hspace{0.17em}\mathit{\delta \omega }\\ & \mathrm{=}& \frac{\sqrt{\mathrm{2}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{t}}}{\mathit{a}}\mathrm{-}\frac{\mathit{a}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right)}{{\mathit{a}}_{\mathrm{t}}}}}{\sqrt{\mathrm{2}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{t}}}{\mathit{a}}\mathrm{-}\frac{\mathit{a}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right){\cos }^{\mathrm{2}}\mathit{i}}{{\mathit{a}}_{\mathrm{t}}}}}\sin \mathit{i}\hspace{0.17em}\mathit{\delta \omega }\\ \mathit{\delta \omega }& \mathrm{=}& \frac{\sqrt{\mathrm{2}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{t}}}{\mathit{a}}\mathrm{-}\frac{\mathit{a}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right){\cos }^{\mathrm{2}}\mathit{i}}{{\mathit{a}}_{\mathrm{t}}}}}{\sqrt{\mathrm{2}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{t}}}{\mathit{a}}\mathrm{-}\frac{\mathit{a}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right)}{{\mathit{a}}_{\mathrm{t}}}}\sin \mathit{i}}\hspace{0.17em}\mathit{d}\mathit{.}\end{array}$

Appendix C: Appendix C

The target has a circular orbit around the Sun with i = 0, which is situated at a distance of 1 AU. We denote here its collisional radius by Δ. Consider a projectile orbit with ω = 0. Thus, its perihelion point coincides with the ascending node. We also assume that $\begin{array}{ccc}{\mathit{q}}^{\mathrm{\prime }}\mathrm{=}\mathrm{1}\mathrm{-}\mathrm{\Delta }\mathit{,}& & \end{array}$so the perihelion point is situated on the surface of the target ring, where this cuts the reference plane on the side facing the Sun.

If the projectile inclination is i, we can calculate the distance of the projectile from the target orbit at any true anomaly f, since we have $\begin{array}{ccc}\mathit{r}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{=}\frac{\mathit{p}}{\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}}& & \end{array}$(C.1)for the heliocentric distance expressed in AU, using p = a(1 − e²) for the semilatus rectum, also expressed in AU. Using cylindrical coordinates (s,ϑ,z) (see Fig. C.1), we get $\begin{array}{ccc}\begin{array}{c}\\ & \\ & \end{array}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}& & \end{array}$(C.2)and we can write the square of the distance between the projectile and the target orbit as $\begin{array}{ccc}\mathit{t}\mathrm{=}{\mathit{z}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{s}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathit{.}& & \end{array}$(C.3)Concerning the function t(f), we know that t(0) = Δ² and the first derivative, t′(0) = 0, since the motion of the projectile is orthogonal to the direction to the target ring at perihelion. Therefore, whether or not the projectile orbit enters into the ring to first order is determined by the value of the second derivative, t′′(0). We now set out to calculate this.

Fig. C.1 Geometric illustration for orbits with perihelia close to the target ring. The z axis is perpendicular to the orbital plane of the target, and a cross-section of the ring is shown at unit distance from the Sun, including its radius Δ.

From the above equations, we get $\begin{array}{ccc}\mathit{t}& \mathrm{=}& {\mathit{r}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{i}{\sin }^{\mathrm{2}}\mathit{f}\mathrm{+}{\left[\mathit{r}\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{i}{\sin }^{\mathrm{2}}\mathit{f}}\mathrm{-}\mathrm{1}\right]}^{\mathrm{2}}\\ & \mathrm{=}& \mathrm{1}\mathrm{+}{\mathit{p}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}{\mathrm{)}}^{\mathrm{2}}\\ & & \mathrm{-}\mathrm{2}\mathit{p}\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{i}{\sin }^{\mathrm{2}}\mathit{f}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}\mathrm{)}\mathit{,}\end{array}$which upon differentiation yields

$\begin{array}{ccc}{\mathit{t}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathrm{2}\mathit{p}\sin \mathit{f}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}{\mathrm{)}}^{\mathrm{2}}}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathit{pe}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}\mathrm{)}}\mathrm{-}\frac{\mathit{e}{\cos }^{\mathrm{2}}\mathit{i}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{i}\cos \mathit{f}}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{i}{\sin }^{\mathrm{2}}\mathit{f}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathrm{·}& & \end{array}$(C.4)We easily verify that t′(0) = 0. This equation can be written as $\begin{array}{ccc}{\mathit{t}}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{pA}\mathrm{·}\mathit{B,}& & \end{array}$where $\begin{array}{ccc}\mathit{A}\mathrm{=}\frac{\sin \mathit{f}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}{\mathrm{)}}^{\mathrm{2}}}& & \end{array}$and B is the expression in curly brackets, so we have $\begin{array}{ccc}{\mathit{t}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{2}\mathit{p}\left(\mathit{A}\frac{\mathrm{d}\mathit{B}}{\mathrm{d}\mathit{f}}\mathrm{+}\mathit{B}\frac{\mathrm{d}\mathit{A}}{\mathrm{d}\mathit{f}}\right)\mathrm{·}& & \end{array}$Obviously, A(0) = 0, so we get $\begin{array}{ccc}{\mathit{t}}^{\mathrm{\prime \prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{pB}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathit{A}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{·}& & \end{array}$We easily find that $\begin{array}{ccc}{\mathit{A}}^{\mathrm{\prime }}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}{\mathrm{)}}^{-2}\mathit{,}& & \end{array}$which, like p, is positive definite, so the sign of t′′(0) is the same as the sign of $\begin{array}{ccc}\mathit{B}\mathrm{\left(}\mathrm{0}\mathrm{\right)}& \mathrm{=}& \mathit{pe}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\mathrm{)}\mathrm{-}\mathit{e}{\cos }^{\mathrm{2}}\mathit{i}\mathrm{+}{\sin }^{\mathrm{2}}\mathit{i}\\ & \mathrm{=}& {\cos }^{\mathrm{2}}\mathit{i}\{\mathit{pe}\mathrm{(}\mathrm{1}\mathrm{+}{\tan }^{\mathrm{2}}\mathit{i}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\mathrm{)}\mathrm{-}\mathit{e}\mathrm{+}{\tan }^{\mathrm{2}}\mathit{i}\}\mathit{.}\end{array}$There is, hence, a unique inclination i_c, for which B(0) = 0, given by $\begin{array}{ccc}{\tan }^{\mathrm{2}}{\mathit{i}}_{\mathrm{c}}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{e}\mathrm{-}\mathit{p}}{\mathrm{1}\mathrm{+}\mathit{p}\mathrm{+}{\mathit{e}}^{-1}}\mathrm{·}& & \end{array}$(C.5)Since we have p = q′(1 + e) = (1 − Δ)(1 + e), we obtain $\begin{array}{ccc}{\tan }^{\mathrm{2}}{\mathit{i}}_{\mathrm{c}}\mathrm{=}\frac{\mathit{e}\mathrm{\Delta }}{\mathrm{1}\mathrm{+}\mathit{e}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Delta }\mathrm{)}}\mathit{,}& & \end{array}$and since Δ ≪ 1, the final result is $\begin{array}{ccc}{\mathit{i}}_{\mathrm{c}}\mathrm{\approx }\arctan \sqrt{\frac{\mathit{e}\mathrm{\Delta }}{\mathrm{1}\mathrm{+}\mathit{e}}}\mathit{.}& & \end{array}$(C.6)We easily verify that B(0) < 0 at i<i_c, while B(0) > 0 at i>i_c. Thus, for the particular orbit that we have chosen, t has a maximum for f = 0 in the first case and a minimum in the second case. In conclusion, for the lowest inclinations (i<i_c), the orbit enters into the target ring on both sides of perihelion, while it just touches the ring at perihelion for higher inclinations. If we decrease q′ slightly in the first case, we get a double ring crossing, and in the second case, we get no crossing at all. If, on the other hand, we increase q′ slightly, we always get a single ring crossing.

Acknowledgments

We thank the referee, Andrea Milani, for useful comments on the manuscript. This work was supported by the Polish National Science Center under Grant No. 2011/01/B/ST9/05442. HR is also indebted to Grant No. 74/10:2 of the Swedish National Space Board.

References

All Figures

	Fig. 1 Intersection of the Earth’s Hill sphere with the b-plane for projectiles having a = 3.5 AU, q = 0.5 AU, and i = 5°, showing the passage points of 3000 objects that entered within the Hill sphere. The center of the Earth is at the origin of the (ξ,ζ) coordinates. These coordinates are described in Sect. 3.2; see also Carusi et al. (1990).
In the text

	Fig. 2 Map of the (ω,Ω) plane for projectiles with the same values of a, q, and i as in Fig. 1, showing all the combinations that yielded ℳ ≤ Rcoll with respect to the Earth. The four critical ranges of ω allowing collisions are easily recognized.
In the text

	Fig. 3 Sketch of the b-plane, where the circular contour marks the intersection with the target’s collisional sphere. The MOID point of the projectile is marked by a small circle on the ξ axis, and the line through points E and E′ (also marked by small circles) is the locus of projectile crossing points for different timings of the encounter.
In the text

	Fig. 4 Minimum mutual distance of one of the objects represented in Fig. 2 from the Earth’s center is shown as a function of the difference between its time of perihelion passage and the critical one yielding the MOID. The values corresponding to the minimum and maximum times of perihelia, for which the mutual distance reaches below the Earth’s collisional radius, are marked.
In the text

	Fig. 5 The dependence of the first quadrant value of ωc on the perihelion distance, as expressed in AU, for an assumed semi-major axis of 3.5 AU and a targeted nodal distance of rN = 1 AU.
In the text

	Fig. 6 Comparison of two different calculations of the width (Δω) of the ω range allowing collisions with the Earth by projectiles on low-inclination orbits. For details, see the text. The dashed curve shows results obtained from Eq. (14), and the solid curve shows numerically calculated values using the MOID-track method.
In the text

	Fig. 7 Mean impact probability with the Earth, as averaged over all arguments of perihelion, versus inclination as the latter goes to zero. Calculations have been made with the MOID-track method. The two panels correspond to two perihelion distances: q′ = 0.16735 and q′ = 0.3. Critical values of i, as defined in the text, are marked.
In the text

	Fig. 8 Comparison of two different calculations of the width (Δω(2)) of two combined ω ranges that allow collisions with the Earth by projectiles on orbits with perihelia near the target orbit for i = 90°. For details, see the text. The dashed curve shows twice the Δω value that is obtained from Eq. (14), and the solid curve shows numerically calculated values using the MOID method.
In the text

	Fig. 9 Time of residence of the projectile (TR) within the target ring for different combinations of q′ and ω – the former is close to 1, the latter close to 0 – for i = 1 deg (upper panel), and i = 0.1 deg (lower panel). The MOID-track method has been used for the calculations. “Slices” showing TR vs. ω for different, constant q′ are shown by different colors, as explained in the diagrams.
In the text

	Fig. 10 Mean impact probability with the Earth, as averaged over all arguments of perihelion, versus perihelion distance q′ as the latter approaches 1. Calculations have been made with the MOID-track method. The two panels correspond to two inclinations: i = 1 deg and i = 0.1 deg. Critical values of q′, as defined in the text, are marked.
In the text

Fig. 11 Collision probability with Mars (using its current orbital and physical properties) per orbital revolution as a function of the projectile perihelion distance (qp) divided by Mars’ perihelion distance (qt) according to the Hill sphere method. The probability peaks when this ratio equals unity, and the variation is the same for both dynamical models used: the two-body problem (2BP) and the restricted 3-body problem (R3BP).

In the text

	Fig. 12 Collision probability with Mars (using its current orbital and physical properties) per orbital revolution as a function of the projectile perihelion distance (qp) divided by Mars’ perihelion distance (qt), according to different methods. The variations are the same to good accuracy for all methods, except in the immediate vicinity of qp = qt.
In the text

	Fig. 13 Collision probability with the Earth (using a circular target orbit) per orbital revolution as a function of the projectile inclination (i), according to different methods. The upper panel shows the divergence between the MOID and Hill sphere methods for i between 0.1 and 1°, and the lower panel shows the divergence between the MOID and Wetherill averaging methods for i< 0.014 deg.
In the text

	Fig. 14 Collision probability with Mars (using its current orbital and physical properties) per orbital revolution as a function of the projectile perihelion distance (qp) divided by Mars’ perihelion distance (qt), according to different methods. The variations are shown for the immediate vicinity of qp = qt, using an abscissa proportional to log (1 − qp/qt).
In the text

	Fig. 15 Ratios of average impact probabilities with Mars for projectiles having different perihelion distances (in units of Mars’ perihelion distance) plotted versus the articially varied collisional radius of the planet.
In the text

	Fig. 16 Collision probability with the Earth for very low inclination orbits, as the projectile perihelion distance approaches 1 AU from below. Computations have been made with the MOID-chord and MOID-track methods, and significant differences are seen as the ratio of perihelion distances exceeds 0.9999.
In the text

	Fig. 17 Impact probability per orbital revolution by the projectile with respect to the Earth (top panel) and Mars (bottom panel), using color-coding, as explained to the right of the diagrams. All calculations were performed using the MOID-track method.
In the text

	Fig. C.1 Geometric illustration for orbits with perihelia close to the target ring. The z axis is perpendicular to the orbital plane of the target, and a cross-section of the ring is shown at unit distance from the Sun, including its radius Δ.
In the text