Issue 
A&A
Volume 569, September 2014



Article Number  A47  
Number of page(s)  15  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201423966  
Published online  18 September 2014 
Monte Carlo methods to calculate impact probabilities^{⋆}
^{1}
P.A.S. Space Research Center, Bartycka 18A, 00716
Warszawa, Poland
email:
wajer@cbk.waw.pl
^{2}
Dept. of Physics and Astronomy, Uppsala University,
Box 516, 75120
Uppsala,
Sweden
^{3}
IAPSINAF, via
Fosso del Cavaliere 100, 00133
Roma,
Italy
^{4}
IFACCNR, via
Madonna del Piano 10, 50019
Sesto Fiorentino (FI),
Italy
Received:
9
April
2014
Accepted:
28
June
2014
Context. Unraveling the events that took place in the solar system during the period known as the late heavy bombardment requires the interpretation of the cratered surfaces of the Moon and terrestrial planets. This, in turn, requires good estimates of the statistical impact probabilities for different source populations of projectiles, a subject that has received relatively little attention, since the works of Öpik (1951, Proc. R. Irish Acad. Sect. A, 54, 165) and Wetherill (1967, J. Geophys. Res., 72, 2429).
Aims. We aim to work around the limitations of the Öpik and Wetherill formulae, which are caused by singularities due to zero denominators under special circumstances. Using modern computers, it is possible to make good estimates of impact probabilities by means of Monte Carlo simulations, and in this work, we explore the available options.
Methods. We describe three basic methods to derive the average impact probability for a projectile with a given semimajor axis, eccentricity, and inclination with respect to a target planet on an elliptic orbit. One is a numerical averaging of the Wetherill formula; the next is a Monte Carlo supersizing method using the target’s Hill sphere. The third uses extensive minimum orbit intersection distance (MOID) calculations for a Monte Carlo sampling of potentially impacting orbits, along with calculations of the relevant interval for the timing of the encounter allowing collision. Numerical experiments are carried out for an intercomparison of the methods and to scrutinize their behavior near the singularities (zero relative inclination and equal perihelion distances).
Results. We find an excellent agreement between all methods in the general case, while there appear large differences in the immediate vicinity of the singularities. With respect to the MOID method, which is the only one that does not involve simplifying assumptions and approximations, the Wetherill averaging impact probability departs by diverging toward infinity, while the Hill sphere method results in a severely underestimated probability. We provide a discussion of the reasons for these differences, and we finally present the results of the MOID method in the form of probability maps for the Earth and Mars on their current orbits. These maps show a relatively flat probability distribution, except for the occurrence of two ridges found at small inclinations and for coinciding projectile/target perihelion distances.
Conclusions. Our results verify the standard formulae in the general case, away from the singularities. In fact, severe shortcomings are limited to the immediate vicinity of those extreme orbits. On the other hand, the new Monte Carlo methods can be used without excessive consumption of computer time, and the MOID method avoids the problems associated with the other methods.
Key words: celestial mechanics / comets: general / planets and satellites: terrestrial planets / minor planets, asteroids: general
Appendices are available in electronic form at http://www.aanda.org
© 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 transplanetary, 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, physicaldynamical lifetimes, and impact probabilities. The latter form the subject of the present paper. We consider the terrestrial planets as targets and small bodies on planetcrossing 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 wellknown problem of estimating impact risks for observed nearearth 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 semimajor 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: (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 longstanding problem. The first analytic results relevant to this case are due to Öpik (1951), who derived the following formula for a circular target orbit: (2)The subscript “t” stands for the target, and R_{t}, a_{t} and V_{t} are its radius, orbital semimajor 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 threebody problem, we have (3)where T_{t} is the Tisserand parameter of the projectile relative to the target: (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 apheliontangent to the target orbit (the last factor tends to infinity). Moreover, it yields only a firstorder approximation to the real case of the loweccentricity 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: (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 & CourPalais (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 highspeed computers allow us to circumvent the problem by using the bruteforce 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 (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 welldefined 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: (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 nonuniform 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 longterm 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 Wetherilltype 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 abovementioned 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 headon 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 nonzero. 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 threebody problem (Suntargetprojectile). By comparing the two, we have a check on whether the longrange 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^{11} 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, (8)Using σ_{coll} and σ_{H} for the crosssections 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, (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: (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 crosssectional 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 bplane 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). 

Open with DEXTER 
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 “supersizing” 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 MOIDchord method and the MOIDtrack 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 ℳ ≤ R_{coll} with respect to the Earth. The four critical ranges of ω allowing collisions are easily recognized. 

Open with DEXTER 
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^{5} 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 MOIDchord method, we use the bplane (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 bplanes, 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 socalled MOID point. For τ ≠ τ_{c}, the bplane is crossed before or after the passage of this point. By plotting all the bplanes 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 bplane crossing point thus moves along the curve, and this motion is a direct reflection of the target’s orbital motion projected onto the bplane. 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 bplane, 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. 

Open with DEXTER 
The geometry is illustrated in Fig. 3. This shows the bplane 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 abovedescribed 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 (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 (12)with Δτ given by Eq. (11). The average impact probability for the whole sample is: (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 MOIDchord 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 (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 MOIDchord 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 (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 (16)Since we have P_{t}V_{t} = 2πa_{t}, we get from Eqs. (14)–(16): (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 MOIDchord 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 (18)we finally get (19)thus recovering Eq. (2). This demonstrates the general consistency between the MOIDchord method and the Öpik formula in the case of a circular target orbit.
In the MOIDtrack method, we use the Keplerian orbits of both target and projectile close to the MOID configuration, thus, avoiding the linearization employed in the MOIDchord 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 MOIDtrack method allows to circumvent the abovementioned problem of the MOIDchord 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. 

Open with DEXTER 
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 MOIDchord 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 threedimensional 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 (20)and the true anomaly f_{1} at r = 1 AU is given by (21)yielding two solutions ±f_{1}, where f_{1} ∈ [ 0,π ]. These correspond to the node being situated on the outgoing branch (f_{1}) or the ingoing branch (− f_{1}). If the node is the ascending one, the critical argument of perihelion is ω_{c} = −f_{1} or + f_{1}, respectively, and if it is the descending one, then ω_{c} = π − f_{1} or π + f_{1}, respectively.
Thus, for a given value of q′, which yields a unique value of f_{1}, 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_{1} = 0 (corresponding to q′ = 1) and f_{1} = π (corresponding to q′ = 0), when ω_{c} is either 0 or π, and for f_{1} = π/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 semimajor axis of 3.5 AU and a targeted nodal distance of r_{N} = 1 AU. 

Open with DEXTER 
The dependence of the first quadrant ω_{c} value () on q′ is shown in the figure above. Close to the three abovementioned 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 , 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 abovementioned, critical inclinations may be computed using the sine rule of spherical trigonometry from (22)where is the value of R_{coll} expressed in AU and η is either or .
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 , 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 , 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 abovedescribed, 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 semimajor 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 lowinclination 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 MOIDtrack method. 

Open with DEXTER 
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 nonvanishing 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 MOIDtrack 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 MOIDtrack 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. 

Open with DEXTER 
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 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 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 , and for this reason, we introduce a quantity, which we call Δω_{(2)} with the following definition. For , the two intervals are separate, and they have the same width Δω. The variable Δω_{(2)} is the width of the combined range, and thus, Δω_{(2)} = 2·Δω. For , the combined range extends between the exterior limits of the joined intervals, and Δω_{(2)}< 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. 

Open with DEXTER 
In Fig. 8, we show Δω_{(2)} 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 . We note that Δω_{(2)} 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 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 MOIDtrack 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 (T_{R}) 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 MOIDtrack method has been used for the calculations. “Slices” showing T_{R} vs. ω for different, constant q′ are shown by different colors, as explained in the diagrams. 

Open with DEXTER 
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 MOIDtrack 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. 

Open with DEXTER 
In Fig. 10, we show these variations as obtained with the MOIDtrack 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 (). 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 , 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 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 Δω_{(2)} 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 semimajor axis of 3.5 AU for the calculations.
5.1. Check of the twobody 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 longrange 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 (q_{p}) divided by Mars’ perihelion distance (q_{t}) 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 twobody problem (2BP) and the restricted 3body problem (R3BP). 

Open with DEXTER 
Using the Hill sphere method, we decided to test this by comparing the results from the two dynamical models that we used. The twobody problem corresponds to the classical, Keplerian approach, while the restricted threebody problem (SunTargetProjectile) includes the longrange 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 Marscrossing 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 (q_{p}) divided by Mars’ perihelion distance (q_{t}), according to different methods. The variations are the same to good accuracy for all methods, except in the immediate vicinity of q_{p} = q_{t}. 

Open with DEXTER 
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 crosssection like illustrated in Fig. 1. As a consequence, by erroneously assuming that the crosssectional 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. 

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

Open with DEXTER 
In the case under consideration, the target is Mars, and estimating the minimum inclination for the Hill sphere method to work as (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 abovementioned supersizing 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 lowinclination 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 crosssection of the Hill sphere method. The timing range in that case cannot increase in proportion to the increase of the crosssection 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 MOIDtrack 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 abovementioned 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. 

Open with DEXTER 
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 MOIDchord and MOIDtrack methods, and significant differences are seen as the ratio of perihelion distances exceeds 0.9999. 

Open with DEXTER 
Fig. 17 Impact probability per orbital revolution by the projectile with respect to the Earth (top panel) and Mars (bottom panel), using colorcoding, as explained to the right of the diagrams. All calculations were performed using the MOIDtrack method. 

Open with DEXTER 
Our final comparison between methods refers to the two variants of the MOID method. In Fig. 14, one can notice slight differences between the MOIDchord and MOIDtrack 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 MOIDtrack method does not involve any approximation when finding the timing range that allows collisions, while the MOIDchord 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 MOIDchord 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 MOIDtrack 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 MOIDtrack 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 lowinclination 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 semimajor 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 crosssection. 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 MOIDtrack 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 crosssection 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 abovementioned critical situations, the Hill sphere method (and other supersizing methods based on the same principle) underestimates the statistical impact probability, while the Wetherill averaging yields results tending to infinity.
We thus used the MOIDtrack 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 semimajor 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 orderofmagnitude estimate, since the angular elements are not uniformly distributed.
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
 Bottke, W. F., Vokrouhlický, D., Minton, D., et al. 2012, Nature, 485, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Carusi, A., Valsecchi, G. B., & Greenberg, R. 1990, Celest. Mech. Dyn. Astron., 49, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Greenberg, R. 1982, Astron. J., 87, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Greenberg, R., Carusi, A., & Valsecchi, G. B. 1988, Icarus, 75, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hénon, M. 1968, Bull. Astron., 3, 3, 377 [Google Scholar]
 Horner, J., & Jones, B. W. 2008, Int. J. Astrobiology, 7, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Horner, J., Mousis, O., Petit, J.M., & Jones, B. W. 2009, Planet. Space Sci., 57, 1338 [NASA ADS] [CrossRef] [Google Scholar]
 Kessler, D. J. 1981, Icarus, 48, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Kessler, D. J., & CourPalais, B. G. 1978, J. Geophys. Res., 83, 2637 [NASA ADS] [CrossRef] [Google Scholar]
 Milani, A., Carpino, M., & Marzari, F. 1990, Icarus, 88, 292 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Brasser, R., Gomes, R., Levison, H. F., & Tsiganis, K. 2010, AJ, 140, 1391 [NASA ADS] [CrossRef] [Google Scholar]
 Öpik, E. J. 1951, Proc. R. Irish Acad. Sect. A, 54, 165 [Google Scholar]
 Öpik, E. J. 1976, Interplanetary encounters: closerange gravitational interactions (New York: Elsevier) [Google Scholar]
 Pokorný, P., & Vokrouhlický, D. 2013, Icarus, 226, 682 [NASA ADS] [CrossRef] [Google Scholar]
 Steel, D. I., & Baggaley, W. J. 1985, MNRAS, 212, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Valsecchi, G. B. 2006, in Dynamics of Extended Celestial Bodies and Rings, ed. J. Souchay, Lect. Not. Phys. (Berlin: Springer Verlag), 682, 145 [Google Scholar]
 Valsecchi, G. B., Milani, A., Gronchi, G. F., & Chesley, S. R. 2003, A&A, 408, 1179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valsecchi, G. B., Rossi, A., Milani, A., & Chesley, S. R. 2005, in Dynamics of Populations of Planetary Systems, eds. Z. Knežević & A. Milani, IAU Colloq. 197, 249 [Google Scholar]
 Wetherill, G. W. 1967, J. Geophys. Res., 72, 2429 [NASA ADS] [CrossRef] [Google Scholar]
 Wiśniowski, T., & Rickman, H. 2013, Acta Astron., 63, 293 [NASA ADS] [Google Scholar]
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 where a and p are the semimajor axis and semilatus rectum of the projectile, and the heliocentric distance ρ of the target is expressed by Eq. (7). We use for the semilatus rectum of the target. The reciprocal of the factor under consideration is The nonsingular nominator is easily dealt with by numerical integration, so we need to concentrate only on the integral, (A.1)where ν is the true anomaly of the target. The denominator can be written as , 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 (A.2)The nonsingular factor can again be disregarded, so we concentrate on the expression
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 nonsingular nominator can again be disregarded, so we concentrate on the integral, Substituting (A.3)the integral is transformed into which is proper, when cos^{2}(ν_{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:

coplanar 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 bplane 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:
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 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 (C.1)for the heliocentric distance expressed in AU, using p = a(1 − e^{2}) for the semilatus rectum, also expressed in AU. Using cylindrical coordinates (s,ϑ,z) (see Fig. C.1), we get (C.2)and we can write the square of the distance between the projectile and the target orbit as (C.3)Concerning the function t(f), we know that t(0) = Δ^{2} 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 crosssection of the ring is shown at unit distance from the Sun, including its radius Δ. 

Open with DEXTER 
From the above equations, we get which upon differentiation yields
(C.4)We easily verify that t′(0) = 0. This equation can be written as where and B is the expression in curly brackets, so we have Obviously, A(0) = 0, so we get We easily find that which, like p, is positive definite, so the sign of t′′(0) is the same as the sign of There is, hence, a unique inclination i_{c}, for which B(0) = 0, given by (C.5)Since we have p = q′(1 + e) = (1 − Δ)(1 + e), we obtain and since Δ ≪ 1, the final result is (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.
All Figures
Fig. 1 Intersection of the Earth’s Hill sphere with the bplane 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). 

Open with DEXTER  
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 ℳ ≤ R_{coll} with respect to the Earth. The four critical ranges of ω allowing collisions are easily recognized. 

Open with DEXTER  
In the text 
Fig. 3 Sketch of the bplane, 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. 

Open with DEXTER  
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. 

Open with DEXTER  
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 semimajor axis of 3.5 AU and a targeted nodal distance of r_{N} = 1 AU. 

Open with DEXTER  
In the text 
Fig. 6 Comparison of two different calculations of the width (Δω) of the ω range allowing collisions with the Earth by projectiles on lowinclination 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 MOIDtrack method. 

Open with DEXTER  
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 MOIDtrack 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. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 9 Time of residence of the projectile (T_{R}) 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 MOIDtrack method has been used for the calculations. “Slices” showing T_{R} vs. ω for different, constant q′ are shown by different colors, as explained in the diagrams. 

Open with DEXTER  
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 MOIDtrack 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. 

Open with DEXTER  
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 (q_{p}) divided by Mars’ perihelion distance (q_{t}) 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 twobody problem (2BP) and the restricted 3body problem (R3BP). 

Open with DEXTER  
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 (q_{p}) divided by Mars’ perihelion distance (q_{t}), according to different methods. The variations are the same to good accuracy for all methods, except in the immediate vicinity of q_{p} = q_{t}. 

Open with DEXTER  
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. 

Open with DEXTER  
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 (q_{p}) divided by Mars’ perihelion distance (q_{t}), according to different methods. The variations are shown for the immediate vicinity of q_{p} = q_{t}, using an abscissa proportional to log (1 − q_{p}/q_{t}). 

Open with DEXTER  
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. 

Open with DEXTER  
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 MOIDchord and MOIDtrack methods, and significant differences are seen as the ratio of perihelion distances exceeds 0.9999. 

Open with DEXTER  
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 colorcoding, as explained to the right of the diagrams. All calculations were performed using the MOIDtrack method. 

Open with DEXTER  
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 crosssection of the ring is shown at unit distance from the Sun, including its radius Δ. 

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