1 Introduction

The problem of working out collision rates between solar-system bodies has many applications in theories concerning the long-term physical evolution of objects which make up the solar system. Collision rates are of particular importance in the main asteroid belt (Farinella & Davis 1992), where it has long been recognized that collisions play a key role; also in the Edgeworth-Kuiper belt (Farinella & Davis 1996), in meteoritic astronomy, and in all studies concerning the cratering history of the terrestrial planets. Collision rates are also a key factor in understanding the evolution of the zodiacal cloud and in modelling the development of asteroidal surface regoliths.

Algorithms for assessing collision probabilities are also demanded in fields such as NEO studies, and where there are suggestions that previous values of collision probabilities in the main belt (Farinella & Davis 1992; Vedder 1996, 1998), based on the observed population of bodies larger than 50 km, may not be representative of the collision probability for smaller bodies (Cellino et al. 1996). Resolving the latter question requires a new computational effort based on a much larger sample extending to include the small bodies.

One of the first to provide an analytic expression for the collisional probability between a pair of orbits was Öpik (1951). His approach has serious limitations (e.g. one orbit must be circular), but it has been widely used in the past and is still ideal as a quick estimator of collision rates against the planets, whose orbits are roughly circular. Unfortunately Öpik's method, although elegant and purely analytical, has singularities at the apsides (i.e. when the perihelia or aphelia just touch the circular orbit) and whenever the orbits are coplanar.

An improved technique by Wetherill (1967) allows both orbits to be elliptical, and allows for an assumed constant perihelion precession by assuming that the distribution of the value of the true anomaly of one orbit at the point of intersection with the other is uniform. This assumption allows the probability function to be integrated by Monte Carlo sampling, but for such low-dimensional problems Monte Carlo integrations are relatively inefficient.

Greenberg's (1982) method (cf. Bottke & Greenberg 1993), which can be viewed as a variation of the Öpik-Wetherill schemes, is more comprehensive that either of the above, but the geometry unfortunately makes it difficult to implement. The algorithm which we describe below is based on these methods, but uses different reference geometries in order to reduce programming complexity and improve performance. Following Bottke et al. (1994), our algorithm correctly weights the collision velocities between any pair of particles, and hence provides a more accurate determination of the frequency distribution of collision velocities in any particular case.

Another method, implemented by Kessler (1981) and Steel & Baggeley (1985), evaluates the probability density of the particles around their respective orbits using a simple kinetic theory to determine the collision probability. This method generally produces excellent results, but is comparatively slow. The more basic "particle in a box" methods make no allowance for the specific orbital geometries of individual colliding particles, and while useful for simulating large systems are comparatively inaccurate.

More recent innovations include those by Dell'Oro & Paolicchi (1997) and Vedder (1996). The former authors introduce the novel concept of estimating the collision probability between a target and a large population of field particles, by selecting only those orbits able to collide with the target, and then adjusting the distribution of such orbits to match the statistical properties of the underlying ensemble. This has the advantage of focusing attention solely on orbits which might collide with a given target, allowing the investigation of different field populations with little additional computational effort. By contrast, Vedder (1996) develops a probabilistic method based on the frequency distribution of close approaches between elliptical orbits. Whereas this provides an interesting alternative to other methods, avoiding singularities associated with some methods and orbital geometries, it is relatively computer intensive and not ideally suited as a quick collision probability estimator.

All these methods assume that the collision probability per unit time is sufficiently small that the interacting particles fully sample the orbital phase space before any collision can take place. However, some pairs of particles may be restricted (e.g. by resonances) to certain regions of phase space, whereas other pairs may undergo significant evolution in elements (e.g. semi-major axis) which are assumed in most methods to be constant or slowly varying. The alternative approach, namely direct numerical integration of the relevant particles, in principle allows a very good determination of the collision probability between two objects, but is very computationally expensive. Such an approach has been used, for example, by Yoshikawa & Nakamura (1994), and the numerical integration approach has been used in different regimes by other authors (e.g. Michel et al. 1996). There are some problems for which this is the only valid approach, a good example being the case of the Trojan swarms, where all the objects are librating about the Jovian L₄ and L₅ Lagrange points (Marzari et al. 1996).

Here we present a numerical-based method for assessing the collision probability between pairs of particles which is applicable to various types of problem. It is fast, depends on simple assumptions, and in general allows an accurate estimation of the collision probability between arbitrary pairs of elliptical orbits.