Free Access
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/0004-6361/201423966
Published online 18 September 2014

© 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,Pt ] for τ, where Pt is the orbital period of the target: (1)and M and Mt 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 Pcoll(a,e,i).

2.1. Circular target orbit

The computation of Pcoll 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: (2)The subscript “t” stands for the target, and Rt, at and Vt 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 (3)where Tt 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 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 Rt can simply be replaced by the collisional radius: (5)where Ve 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 Pcoll implies an averaging of S over a torus surrounding the target orbit with radius Rcoll, which we refer to as the “target ring”. For a circular target orbit with zero inclination, singularities then occur when using S(at,0) as an approximation, in the case of at = q, at = 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 (6)where et 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, PE can be used to derive Pcoll(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 PE. 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 Rt by Rcoll in Eq. (6), Rcoll(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 ~1011 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 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, (9)Note that scoll 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 NS 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.

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

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 “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 ℳ ≤ Rcoll. 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 ℳ ≤ Rcoll. These are found empirically as follows.

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

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 ℳ ≤ Rcoll, 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 NS within this range replaces a sample of size NS·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 105 members with ℳ ≤ Rcoll. 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 Rcoll 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.

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

Open with DEXTER

The geometry is illustrated in Fig. 3. This shows the b-plane with a circular contour marking the intersection with the collisional sphere (radius Rcoll). 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 (11)where Vt is the orbital velocity of the target and θ is the angle between the vectors Vt and U.

Each of the retained objects has its specific values of , Vt, 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 NS 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 (14)where x = a/at. 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 ω. 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 Rcoll 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 Pcoll is given by Eqs. (11) and (12), and the average is taken with a uniform distribution of ℳ ∈ [ 0,Rcoll ] and constant values of Vt and sinθ. The result is (16)Since we have PtVt = 2πat, we get from Eqs. (14)–(16): (17)Note that the factor sinθ in the denominator may take the value zero in case the vectors Vt 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 (18)we finally get (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 ℳ ≤ Rcoll. 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 Dmin as a function of τ. The limits τmin and τmax are the values for which Dmin = Rcoll.

Note that the interval Δτ = τmaxτmin is always finite, even if the vectors Vt 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.

thumbnail 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 Dmin versus ττc for one of the orbits with ℳ ≤ Rcoll 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 Rcoll, the three-dimensional target for impacts is a torus with radius Rcoll 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 rN = 1 AU. Let q denote the value of q expressed in AU. The eccentricity is then (20)and the true anomaly f1 at r = 1 AU is given by (21)yielding two solutions ±f1, where f1 ∈ [ 0 ]. These correspond to the node being situated on the outgoing branch (f1) or the ingoing branch (f1). If the node is the ascending one, the critical argument of perihelion is ωc = −f1 or + f1, respectively, and if it is the descending one, then ωc = πf1 or π + f1, respectively.

Thus, for a given value of q, which yields a unique value of f1, there are four values of ωc, for which the orbit has a node at rN = 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 f1 = 0 (corresponding to q′ = 1) and f1 = π (corresponding to q′ = 0), when ωc is either 0 or π, and for f1 = π/2 (corresponding to q′ = 0.54196), when ωc equals ± π/2.

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

Open with DEXTER

The dependence of the first quadrant ωc value () 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 ≤ ℳ ≤ Rcoll, 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 ℳ = Rcoll. 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 (Pcoll) 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 Pcoll.

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 Pcoll 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 (22)where is the value of Rcoll 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 , Pcoll 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 Pcoll 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 ℳ ≤ Rcoll 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 Rt = 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 ic given by Eq. (22) coincide at ic = 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).

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

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<ic, 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 Pcoll(3.5 AU,0.9522,0) = 3.0 × 10-5, while the Öpik formula in Eq. (2) diverges toward infinity.

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

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 ic 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 ic = 0.00553 degrees, while we get ic = 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 Pcoll(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 ic 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 ℳ ≤ Rcoll 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.

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

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

Open with DEXTER

If we estimate a maximum for the residence time, TR, as the time it takes the projectile to move between Rcoll and + Rcoll 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 TR 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 TR 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 Pcoll with q.

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

Open with DEXTER

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 (). 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 TR 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 TR 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 Pcoll 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 Pcoll.

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.

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

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 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 has been used to calculate the probability as a function of the perihelion distance of the projectile.

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

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

thumbnail 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 , and the lower panel shows the divergence between the MOID and Wetherill averaging methods for i< 0.014 deg.

Open with DEXTER

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

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 mM of Mars yields imin = 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 imin = 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 qp/qt → 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 qp/qt ≃ 0.993, but the Hill sphere results fall increasingly short of those obtained by the MOID methods for larger values (up to qp/qt ≃ 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 qp/qt 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. Rcoll. With qp = qt, the impact probability is considerably higher than with a slightly lower qp, 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 qp = qt spike of the impact probability.

thumbnail 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

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

Open with DEXTER

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

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 MOID-chord and MOID-track results, when qp/qt 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′ = qp/qt 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, Pcoll(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 Pcoll 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 PE (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, ip → 0 or qpqt. The behavior of Pcoll 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 ip. 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 Pcoll, has a maximum when qp/qt 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.

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

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 semi-major 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 non-singular 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 non-singular factor can again be disregarded, so we concentrate on the expression

We note that the above manipulation fails in case et = 0. Hence, the integration cannot be performed, if the target orbit is circular with at = q. Excluding this case. the non-singular nominator can again be disregarded, so we concentrate on the integral, Substituting (A.3)the integral is transformed into which is proper, when cos2(νc/2) ≠ 1, or νc ≠ 0. Our method to integrate the Wetherill formula starts by flagging an error, if i = 0 or qt = q (i.e., νc = 0). Otherwise, we either have qt>q, in which case a Simpson rule integration of the Wetherill expression is performed without problems (as long as Qt<Q), or qt<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:

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 − e2) 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.

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

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 ic, 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<ic, while B(0) > 0 at i>ic. 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<ic), 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

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

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

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

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

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

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

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

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

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

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

Open with DEXTER
In the text
thumbnail 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 , 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
thumbnail 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).

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

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

Open with DEXTER
In the text
thumbnail 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 Δ.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.