Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://star.arm.ac.uk/preprints/426.pdf Äàòà èçìåíåíèÿ: Fri Oct 29 17:38:58 2004 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 21:52:35 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: http astrokuban.info astrokuban

Mon. Not. R. Astron. Soc. 354, 798­810 (2004)

doi:10.1111/j.1365-2966.2004.08240.x

Simulations of the population of Centaurs ­ I. The bulk statistics
J. Horner,1
1 2 3 4

,2

N. W. Evans

2,3

and M. E. Bailey

4

¨ Physikalisches Institut, Universitat Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland Theoretical Physics, Department of Physics, 1 Keble Rd, Oxford OX1 3NP Institute of Astronomy, Madingley Rd, Cambridge CB3 0HA Armagh Observatory, College Hill, Armagh BT61 9DG

Accepted 2004 July 19. Received 2004 July 19; in original form 2004 April 28

ABSTRACT

Large-scale simulations of the Centaur population are carried out. The evolution of 23 328 particles based on the orbits of 32 well-known Centaurs is followed for up to 3 Myr in the forward and backward direction under the influence of the four massive planets. The objects exhibit a rich variety of dynamical behaviour with half-lives ranging from 540 kyr (1996 AR20) to 32 Myr (2000 FZ53). The mean half-life of the entire sample of Centaurs is 2.7 Myr. The data are analysed using a classification scheme based on the controlling planets at perihelion and aphelion, previously given in Horner et al. Transfer probabilities are computed and show the main dynamical pathways of the Centaur population. The total number of Centaurs with diameters larger than 1 km is estimated as 44 300, assuming an inward flux of one new short-period comet every 200 yr. The flux into the Centaur region from the Edgeworth­Kuiper Belt is estimated to be one new object every 125 yr. Finally, the flux from the Centaur region to Earth-crossing orbits is one new Earth-crosser every 880 yr. Key words: stellar dynamics ­ celestial mechanics ­ Kuiper belt ­ minor planets, asteroids ­ planets and satellites: general.

1 INTR ODUCTION The dynamical behaviour of Centaurs is still poorly understood. It is possible for a Centaur to work its way slowly inwards through the outer Solar system, leading to eventual capture by Jupiter and designation as a short-period comet. It is also possible for Centaurs to drift outwards to join the Edgeworth­Kuiper Belt, to be ejected from the Solar system in an encounter with one of the massive, outer planets, or even to be captured by these planets into temporary satellite orbits. A small number may even impact upon the planets. Therefore, the Centaurs potentially hold the key to understanding the mechanisms by which the short-period comet population is maintained, to explaining the distant, retrograde satellites of the massive planets, and to allowing us a glimpse of objects newly introduced to the Solar system from the Edgeworth­Kuiper Belt. As early as 1990, when the only known Centaur was Chiron, it was realized that such objects may lie on very unstable orbits. Numerical integrations by Hahn & Bailey (1990) found that Chiron had a halflife for ejection of around 1 Myr, but that the half-life to become a short-period comet for the object was around only 200 000 yr, implying that Chiron could well have been a short-period comet in the past and could possibly become one in the future. This is of particular interest given the size of Chiron (d 140 ­180 km; Groussin, Lamy
E-mail: jonathan.horner@phim.unibe.ch (JH); nwe@ast.cam.ac.uk (NWE); meb@arm.ac.uk (MEB)

& Jorda 2004 and references therein) and other Centaurs, because objects that large entering the inner Solar system would be both spectacular and dangerous. In fact, the idea has been mooted that objects of such size arrive in the inner Solar system with some frequency, and then fragment, leading to swarms of debris that have the potential to encounter the Earth. The Kreutz sun-grazer family may represent one example of such hierarchical fragmentation, whereas other cases in which comets of more ordinary size have undergone catastrophic fragmentation include 3D/Biela, D/1994 (Shoemaker­ Levy 9) and C/1994 S4 (LINEAR). Whether such a decay mode represents a generic process in determining the number of short-period comets can, in principle, be tested by examining the differences in the size distribution of Jupiter-family comets from those of their probable source objects, namely Centaurs, Edgeworth­Kuiper Belt objects, long-period comets and so on (cf. Lowry, Fitzsimmons & Collander-Brown 2003; Lamy et al. 2004). Clube & Napier (1984) have suggested that the Taurid meteoroid swarm may be the relic of the last large object to undergo such a decay. Despite the importance of the Centaurs, there has been little systematic study of the population using numerical simulations. Early calculations on Chiron (Oikawa & Everhart 1979; Hahn & Bailey 1990) and Pholus (Asher & Steel 1993) identified the chaotic nature of these two objects, though only small numbers of clones and modest integration times (<1 Myr) were used. Work by Dones, Levison & Duncan (1996) looked at the behaviour of four Centaurs (Chiron, Pholus, Nessus and 1994 TA) and two Jupiter-family
C

2004 RAS


Simulations of Centaurs ­ I. Statistics
comets (29P/Schwassmann­Wachmann 1 and 39P/Oterma). They found that the number of surviving objects decays exponentially during the early part of the integrations, whilst the decay becomes flatter after a number of half-lives. They also noted that Centaur halflives inferred from numerical integrations are smaller than those ¨ deduced from approximations such as Opik's (1976) theory and diffusion equations (e.g. van Woerkom 1948). Levison & Duncan (1997) ran orbital integrations of 2200 test particles evolving from the Edgeworth­Kuiper Belt to short-period comets, passing through the Centaur region in the process. The study of the integrations was mainly focused on the behaviour of the objects both in the Edgeworth­Kuiper Belt and the cometary region, rather than in the Centaur region itself. Hence, the dynamics of the Centaur population remains largely unexplored. In an earlier paper, Horner et al. (2003) introduced a new method of classifying objects in the Solar system. This was based on the idea that the dynamical evolution is largely determined by the planets that control the perihelion and the aphelion of the object. This classification is particularly useful for the Centaurs, as it breaks down the trans-Jovian region into 20 categories given in the upper panel of Table 1. Objects are labelled according to the controlling giant planet, so, for example, a JS object has perihelion under the control of Jupiter and aphelion under the control of Saturn. Objects with perihelion distance q 33.5 au are designated as either members of the Edgeworth­Kuiper Belt (EK) or trans-Neptunian disc (T). In addition, objects with q 4 au are designated as comets and are subdivided into Encke types (E), short-period (S), intermediate (I) and long-period (L), as summarized in the lower panel of Table 1.
Table 1. The classification scheme introduced by Horner et al. (2003). In the upper table, the first letter designates the planet controlling the perihelion, the second letter the planet controlling the aphelion or the region in which the aphelion lies, with the final two classes EK and T being beyond all the giant planets. (J, Jupiter; S, Saturn; U, Uranus; N, Neptune; EK, Edgeworth­ Kuiper Belt; T, trans-Neptunian belt). The lower table refers to cometary bodies (E, Encke; SP, short-period; I, intermediate and L, long-period). Object J JS JU JN JE JT S SU SN SE ST U UN UE UT N NE NT EK T E SP I L
C

799

The aim of this paper is to provide results from a detailed set of simulations, exploiting the new classification scheme. The orbits of 32 of the best known Centaurs were used to create an ensemble of 23 328 clones. The clones were integrated in the presence of the four massive outer planets in both forward and backward directions for a period of 3 â 106 yr, giving a vast data set with which to examine the dynamics of their orbits. Section 2 describes the details of the numerical simulations, whilst Section 3 provides half-lives for individual Centaurs. Section 4 gives transition probabilities, which allow the main dynamical pathways through this region of the Solar system to be identified. The simulations are used to estimate the total population of Centaurs in Section 5, together with the typical fluxes inwards from the Edgeworth­Kuiper Belt. Finally, Section 6 considers possible correlations between the dynamics and the colours of the Centaurs. 2 INTEGRA TIONS In order to study the bulk statistics of Centaurs, 32 objects were selected from the list of Centaurs given on the website of the Minor Planet Center. 1 The objects were restricted to those with an observed arc of at least 30 d and an aphelion distance of less than 40 au. This ensures that only Centaurs with moderately well-determined orbits were included in our sample. The list of objects is given in Table 2. Over time, as the Centaurs are observed over longer arcs, the accuracy with which the orbits are known increases, and the orbits given on the Minor Planet website change accordingly. The orbits used in these integrations therefore represent the best available information as of 2002 June. Table 2 also gives the absolute magnitude H of each Centaur, which is defined as the apparent magnitude the object would have, if it were placed at both 1 au from the Earth and 1 au from the Sun and was observed at zero phase angle. This is calculated in ignorance of any out-gassing that might occur. We can estimate the maximum and minimum diameter, assuming values of the albedo between 0.15 and 0.02. This gives a crude calculation of the size, though photometric work is required to obtain any more detail. Of the objects studied in these integrations, the one with the brightest absolute magnitude is 1995 SN55, with a value of H = 6.0, which corresponds to a diameter of between 220 and 590 km. The object with the faintest absolute magnitude is 2000 GM137, with H = 14.3 giving a diameter of between 5 and 13 km. This is similar to the size determined for some cometary nuclei (e.g. Sanzovo et al. 2001; Lamy et al. 2004). Hence, the Centaurs come in a wide range of sizes, from very large (1995 SN55 and Chiron) to those comparable in size with normal comets (2000 GM137). The orbital elements of each object were used to create a swarm of 729 clones, distributed through a small cube of a­e­i (semimajor axis, eccentricity and inclination) space, centred on the original orbit. The clones of the objects were created by incrementally increasing (and decreasing) the semimajor axis of the object by 0.005 au, the eccentricity by 0.005, and the inclination by 0.01 . These increments are sufficiently small that the clones can be considered as initially essentially identical to one another, yet they are large enough to ensure that the subsequent chaotic dynamical evolution following close planetary encounters rapidly disperse their orbits through phase space. So, nine values were used for each of these elements, with the central (fifth) value of the nine having the original orbital elements. This gives 729 clones of each of the 32 Centaurs, giving a grand total of 23 328 objects.
1

Perihelion (in au) 4 4 4 4 4 4 6.6 6.6 6.6 6.6 6.6 12.0 12.0 12.0 12.0 22.5 22.5 22.5 33.5 33.5 q q q q q q q q q q q q q q q q q q q q q q q q 4 4 4 4 6.6 6.6 6.6 6.6 6.6 6.6 12.0 12.0 12.0 12.0 12.0 22.5 22.5 22.5 22.5 33.5 33.5 33.5 60.0 60.0

Aphelion (in au) Q 6.6 6.6 Q 12.0 12.0 Q 22.5 22.5 Q 33.5 33.5 Q 60.0 Q 60.0 Q 12.0 12.0 Q 22.5 22.5 Q 33.5 33.5 Q 60.0 Q 60.0 Q 22.5 22.5 Q 33.5 33.5 Q 60.0 Q 60.0 Q 33.5 33.5 Q 60.0 Q 60.0 Q 60.0 Q 60.0 Q4 4 Q 35 35 Q 1000 Q 1000

http://cfa-www.harvard.edu/iau/lists/Centaurs.html

2004 RAS, MNRAS 354, 798­810


800

J. Horner, N. W. Evans and M. E. Bailey
Table 2. The names of the objects simulated, arranged in order of increasing semimajor axis, together with their orbital elements as of 2002 June. Here, a is the semimajor axis measured in astronomical units (au), e is the eccentricity, i is the inclination (in degrees), is the argument of perihelion (in degrees), is the longitude of the ascending node of the orbit (in degrees) and H is the absolute visual magnitude. D min and D max are the values for the diameter (in km) assuming albedos of 0.15 and 0.02, respectively, which are typical upper and lower limits from albedo measurements of comets and Centaurs performed to date. `Class' is the classification of the object, using the scheme given in Horner et al. (2003). (The data are compiled from the Minor Planet Center.) Object 2000 GM137 1998 SG35 2001 BL41 2001 PT13 2000 EC98 1999 UG5 Chiron 1996 AR20 Chariklo 2001 XZ255 2000 QC243 1994 TA 2001 SQ73 2000 CO104 1999 XX143 Asbolus 2002 GO9 1998 QM107 Pholus 2002 CA249 1999 HD12 2002 DH5 2002 GZ32 1995 SN55 2000 FZ53 Nessus Hylonome 2002 GB10 2001 KF77 1998 TF35 2002 FY36 2000 QB243 a 7.853 8.420 10.071 10.624 10.651 12.778 13.601 15.197 15.775 16.039 16.560 16.849 17.485 17.497 17.886 17.938 19.418 20.042 20.265 20.713 21.322 22.433 23.081 23.564 23.765 24.404 24.909 25.139 25.992 26.429 28.969 28.953 e 0.118 0.307 0.267 0.197 0.471 0.415 0.379 0.627 0.171 0.043 0.203 0.301 0.177 0.256 0.458 0.619 0.277 0.136 0.573 0.385 0.583 0.384 0.216 0.663 0.479 0.517 0.243 0.396 0.240 0.383 0.114 0.381 i 15.9 15.6 11.5 20.4 4.4 5.6 6.9 6.2 23.4 2.6 20.7 5.4 17.4 4.0 6.8 17.6 12.8 9.4 24.7 6.4 10.1 22.5 15.0 5.0 34.9 15.7 4.2 13.3 4.4 12.6 5.4 6.5 123.5 337.5 319.6 86.6 163.5 289.4 339.1 107.9 241.4 294.2 150.0 154.9 304.2 339.2 214.9 290.3 92.0 154.9 354.6 182.4 288.8 323.7 154.4 49.3 290.8 170.1 5.5 238.9 266.4 301.8 194.1 339.4 89.7 173.2 280.1 205.3 173.2 87.4 209.4 330.1 300.4 77.8 337.9 137.7 16.3 346.8 103.8 6.1 117.4 127.2 119.3 313.6 177.7 157.0 107.2 144.6 202.4 31.4 178.2 315.5 14.6 52.0 332.8 331.1 H 14.3 11.3 11.7 9.0 9.5 10.1 6.5 14.0 6.4 11.1 7.6 11.5 9.6 10.0 8.5 9.0 8.5 10.4 7.0 12.0 12.8 10.4 6.9 6.0 11.4 9.6 8.0 7.8 9.4 9.3 8.4 8.2 D
min

D

max

Class S JS SU SU JU SU SU JN U U U SU U U SN SN UN UN SN UN SE UN UN SE UE SE UN UE UN UE N UE

4.7 19 16 54 43 33 170 5 180 21 100 17 41 34 68 54 68 28 140 14 9.4 28 140 220 18 41 86 95 45 47 72 79

13 52 43 150 120 90 470 15 490 57 280 47 110 94 190 150 190 78 370 27 26 78 390 590 49 110 240 260 120 130 200 220

The use of multiple clones of an individual object in the study of its behaviour over time is desirable for a number of reasons. First, the observations contain some uncertainty, which means that the orbit itself is not known beyond a certain degree of precision. This alone would be enough to promote the use of a cluster of orbits with slightly different parameters. In addition, the chaotic nature of the orbits implies that an infinitesimally small change in the initial parameters may lead to a major difference in the final outcome of the simulation. This means that, beyond a certain time in the future, an object could be anywhere within the Solar system, or even beyond, as a result of a tiny change in the initial elements. These two facts taken together suggest that the best means of examining the future or past behaviour of an object is to integrate a large number of clones, and to examine the statistical properties of the data set (Hahn & Bailey 1990; Dones et al. 1996). The number of clones used in such a simulation is chosen to maximize the size of the data set available for analysis, without requiring an excessive amount of time for the simulations to run. The simulations described here took approximately three months to run on a desktop workstation. The clusters of Centaurs were then integrated for 3 Myr in both the forward and backward directions. The gravitational influence of the four Jovian planets (Jupiter, Saturn, Uranus and Neptune)

was included in the integrations, which were all carried out using the hybrid integrator within the MERCURY (Chambers 1999) software package. This is a symplectic integrator, which makes use of a turnover function to switch to an accurate Bulirsch­Stoer algorithm for close encounters. The terrestrial planets (Mercury, Venus, Earth and Mars) were all omitted from the integrations, and their masses added to that of the Sun. The only slight detriment that this causes is loss of accuracy when objects are captured into orbits crossing those of the terrestrial planets. Even in this case, however, the effects of Jupiter (and the other giant planets, if the aphelion of the object lies sufficiently far from the Sun) are generally much greater than those of the terrestrial planets. After some trials, a time-step of 120 d was used for the integrations, as this was found to give a good compromise between speed and accuracy. In order to determine the most efficient time-step, an object was placed on a typical short period cometary orbit with perihelion near the Earth and aphelion near Jupiter. A number of clones were created, and the ensemble was integrated for 105 yr with time-steps of 30, 60, 120, 240 and 360 d. The resulting orbital elements were then compared, and it was found that the results for time-steps of 30, 60, 120 and 240 d gave consistent results, while 360 d was too long a time-step. After a number of such trials, a
C

2004 RAS, MNRAS 354, 798­810


Simulations of Centaurs ­ I. Statistics
7 600 6 400

801

5

200

4

0 0 Time (years)

3 0

1996 AR20 Pholus Nessus 2000 FZ53 Time (years)

Figure 1. The number of surviving clones (left-hand panel) and its natural logarithm (right-hand panel) as a function of time for 1996 AR20, Pholus, Nessus and 2000 FZ53. In all cases, the graphs are deduced from the simulations in the forward direction.

time-step of 120 d was used for the integrations, as this was found to give a good compromise between speed and accuracy, the midrange value of 120 d being chosen so as to err on the side of accuracy where possible. An ejection distance of 1000 au from the Sun was used, following Levison & Duncan (1994). Any object that reached this distance was removed from the integration.2 Also removed were those objects that impacted upon the surface of the Sun (q < 0.005 au), or on any one of the giant planets (the separation is less than the physical radius of the planet). On completion of the integration, data files were created for each clone that gave the values of the orbital elements at 100-yr intervals for the entire lifetime within the integration of the clone. It wason these files that the analysis was carried out. Each one of these files was approximately 5 Mb in size (the exact size varied, as the file terminated with the ejection of the clone from the simulation). Hence, the 23 328 clones in total occupied 120 Gb of disc space for their orbital elemental evolution alone, prior to any analysis ­ a daunting data set by any standards!

The number of clones, N , that remain after time, t, is therefore given by N = N0 e-t , = 0.693 , T1/2 (1)

where N 0 is the initial number of clones, is the decay constant, and T 1/2 is the half-life for the object. To calculate the half-life, the simulation data were analysed with the help of least-square fitting routines from Press et al. (1992). The software provides the value of the 2 function
N

2 () =
i =1

ln Ni - ln N0 + ti i

2

,

(2)

3 HALF-LIVES It is straightforward to calculate the value of the half-life for each Centaur. As illustrated in Fig. 1, the number of clones remaining within the simulation decays in a roughly exponential manner as a function of time. The four objects where the decay is shown are 1996 AR20 (the object with the shortest lifetime of those studied), Pholus (an object with a moderately short lifetime), Nessus (a relatively long-lived object) and 2000 FZ53 (the object with the longest lifetime). It has often been noted in long Solar system integrations that the number of clones remaining within a simulation decays exponentially with time (e.g. Dones et al. 1996; Holman 1997; Evans & Tabachnik 1999). The trait is less obvious with the longer lived objects. For example, Dones et al. (1996) found that the number of surviving clones in their integrations decayed exponentially at early times, while at later times (generally greater than twice the half-life) the decay was slower. This is because those objects surviving the longest are those transferred to the most stable areas of the Solar system.

2

The long-period comets (the L class in Table 1) are defined to have aphelion Q in excess of 1000 au. Our choice of ejection distance means that the statistical properties of the L class cannot be reliably computed from our simulations.
C

where N i is the number of clones remaining at time t i and the i are the individual standard deviations on the data points. As the i are unknown, we proceed by first assigning uniform errors, fitting for the model parameters by minimizing the 2 and then rescaling the errors using equation (15.1.6) of Press et al. (1992). Of course, this well-known procedure precludes an independent estimate of the goodness of fit. The overall data set of 23 328 objects has an ensemble half-life of 2.76 Myr in the forward direction and 2.73 Myr in the backward direction. This gives us an estimate of the mean lifetime of a typical Centaur. Note that this lifetime adds weight to the argument that the population of the Centaur region is in a steady state (a reservoir of objects constantly being drained by Jupiter, and refilled from a long-lived source such as the Edgeworth­Kuiper Belt). Dones et al. (1996) calculated the half-lives of Chiron, Pholus, Nessus and 1994 TA using 100 clones. Their results for Chiron and Nessus are in excellent agreement with ours, but they found T 1/2 = 2.1 and 2.4 Myr for Pholus and 1994 TA ­ somewhat larger than our results. The most likely cause of the discrepancy is in the different algorithms used to populate the clones. Dones et al. carried out a ¨ comparison with approximations such as Opik's (1976) theory and diffusion equations (e.g. van Woerkom 1948) and concluded that both methods significantly overestimate the lifetimes by factors of between 2 and 5. It seems that numerical simulations with large numbers of clones are the only reasonably reliable method for halflife estimation. The half-lives of individual Centaurs are given in Table 3. The value of the Poisson uncertainty is calculated as = T1/2 / N0 , where N 0 is the initial number of clones. This uncertainty is added in quadrature to the uncertainty in the fitted half-life, as judged from the 2 surface in the space of model parameters and as returned

2004 RAS, MNRAS 354, 798­810


802

J. Horner, N. W. Evans and M. E. Bailey
Table 3. Half-lives T 1/2 (in Myr) of the simulated objects, together with an error estimate (in Myr). The column labelled C gives the class of the Centaur, while the column labelled D gives the direction of integration (F for forward, B for backward). Object 1996 AR20 1996 AR20 2000 EC98 2000 EC98 1998 SG35 1998 SG35 2000 GM137 2000 GM137 1995 SN55 1995 SN55 1999 UG5 1999 UG5 Asbolus Asbolus 2001 PT13 2001 PT13 2001 BL41 2001 BL41 Chiron Chiron 1999 XX143 1999 XX143 1999 HD12 1999 HD12 Pholus Pholus 1994 TA 1994 TA 2000 CO104 2000 CO104 2001 XZ255 2001 XZ255 C JN JN JU JU JS JS S S SE SE SU SU SN SN SU SU SU SU SU SU SN SN SE SE SN SN SU SU U U U U D F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B T
1/2

0.02 0.02 0.02 0.02 0.03 0.02 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.04 0.03 0.04 0.04 0.04 0.04 0.04 0.05 0.05 0.04 0.05 0.05 0.07 0.06 0.07 0.08 0.11 0.09

Object 2001 SQ73 2001 SQ73 2002 GO9 2002 GO9 2000 QC243 2000 QC243 2002 CA249 2002 CA249 1998 QM107 1998 QM107 Nessus Nessus Hylonome Hylonome 2001 KF77 2001 KF77 2002 DH5 2002 DH5 2002 GZ32 2002 GZ32 Chariklo Chariklo 1998 TF35 1998 TF35 2002 GB10 2002 GB10 2002 FY36 2002 FY36 2000 QB243 2000 QB243 2000 FZ53 2000 FZ53

C U U UN UN U U UN UN UN UN SE SE UN UN UN UN UN UN UN UN U U UE UE UE UE N N UE UE UE UE

D F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B

T

1/2

0.11 0.10 0.11 0.14 0.12 0.13 0.15 0.09 0.18 0.21 0.18 0.24 0.24 0.27 0.33 0.4 0.34 0.5 0.4 0.28 0.4 0.35 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.7 1.0 1.2

0.54 0.59 0.61 0.63 0.67 0.65 0.72 0.68 0.70 0.80 0.74 0.85 0.86 0.75 0.94 0.87 0.95 0.95 1.03 1.07 1.06 1.38 1.22 1.13 1.28 1.39 1.78 1.52 1.89 2.24 2.94 2.43

2.86 2.73 2.93 3.67 3.18 3.44 4.06 2.54 4.87 5.65 4.91 6.40 6.37 7.30 8.89 10.1 9.08 12.8 11.3 7.78 10.3 9.38 11.5 10.8 11.1 13.1 12.5 13.5 13.0 17.8 26.8 32.3

by our fitting software. Note that our ignorance of individual error bars on our simulation data points may lead to an underestimate of the latter quantity ­ it is generally smaller than by an order of magnitude ­ as our algorithm is tantamount to assuming a good fit to the exponential decay law. The object with the shortest half-life is 1996 AR20, a JN object, which has a half-life of approximately 540 kyr in the forward integration and 590 kyr in the backward integration. The object with the longest half-life is 2000 FZ53, a UE object, with half-lives of 26.8 Myr (forwards) and 32.3 Myr (backwards). On comparison with the orbital elements (Table 2), a correlation can be seen between the position of the orbit of a Centaur within the Solar system and its half-life. The further from the Edgeworth­ Kuiper Belt, the shorter the half-life. This is not unexpected ­ Jupiter is significantly more massive than Saturn, and Saturn in turn is more massive than either Uranus or Neptune. So, the further from the Edgeworth­Kuiper Belt, the more massive are the planets with which the Centaur interacts and the more frequently do such encounters occur. This effect is evident when the orbital elements of the objects are plotted against the logarithm of the half-life, as in Fig. 2. There is only a rough correlation with semimajor axis, but the data indicate a lower bound to the half-life as a function of perihelion q and an upper bound as a function of aphelion Q. Specifically,

we find that 0.392 exp(0.135q ) T1
/2

0.064 exp(0.275 Q ),

(3)

where q and Q are in au and T 1/2 is in Myr. This holds for all the Centaurs in our sample, but it is conceivably possible that loweccentricity orbits between the planets are extremely long-lived (e.g. Holman 1997; Evans & Tabachnik 1999). Additionally, in the plot of perihelion distance versus half-life, there seem to be three rough groupings of objects. The first are those in a band along the line from q = 6, log T 1/2 = 5.75 to q = 26, log T 1/2 = 7.2, which accounts for the bulk of the objects. A second group comprises six objects, which lie roughly on a parallel track at values of log T 1/2 greater by 0.7. Finally, 2000 FZ53 sits alone far above either of these groups. The objects that are in the first group currently lie away from the positions of any major resonances and so tend to have short lifetimes. The objects in the second group tend to lie nearer to stable mean motion or secular resonances. In fact, the orbit of 2002 GB10 lies within 0.009 au of the 3:4 resonance of Uranus, while the orbit of Chariklo lies within 0.09 au of the 4:3 resonance of Uranus. Clones of 2000 FZ53 quite frequently display resonant behaviour during the course of the simulations, although it does not currently lie near any major mean motion resonances. A possible cause of the exceptionally long half-life of 2000 FZ53 is its
C

2004 RAS, MNRAS 354, 798­810


Simulations of Centaurs ­ I. Statistics
8 7.5 7 6.5 6 5.5 0 10 20 30 Semi-major axis (au) 8 7.5 7 6.5 6 5.5 0 10 20 30 Perihelion distance (au) 8 7.5 7 6.5 6 5.5 0 10 20 30 40 Inclination (deg) 8 7.5 7 6.5 6 5.5 0 10 20 30 40 Aphelion distance (au)

803

Figure 2. The relationship of the logarithm of the half-lives of the Centaurs with their semimajor axis, perihelion and aphelion distances, and their inclination. The filled triangles show the results of the forward direction, and the open squares show the results of the backward direction integrations.

abnormally large inclination ­ higher than that of any other Centaur studied by over 10 . Any correlations of half-life with inclination and eccentricity are less clear-cut than those with position. However, there is a lack of long-lived objects at large e. This is a consequence of the fact that Centaurs with large e must cross the orbits of several of the outer planets, and so inevitably are more unstable than bodies where the close approaches are restricted to just one or two planets. As the equations of motion are time-reversible, it might naively be expected that the forward and backwards half-lives should be the same. In fact, it is often found that inward evolution of the orbits of minor bodies is more likely than outward evolution (e.g. Oikawa & Everhart 1979). This is because the Edgeworth­Kuiper Belt provides a source, while Jupiter provides a sink, so that clones ejected by Jupiter do not have an opportunity to return. Of all the objects, 19 have shorter half-lives in the forward direction, and 13 have shorter half-lives in the backward direction. The least discrepancy occurs for 2001 BL41, with the two half-lives agreeing to within 0.25 per cent. The greatest discrepancy occurs for 2002 GZ32, where the backward half-life is 46 per cent longer than the forward halflife. These differences can be visualized as the effect of the first encounters the clones have with the major planets. If an encounter with one of the planets occurs early enough in the simulation, the entire ensemble of clones will be perturbed, moving the objects on to slightly different orbits, with slightly different half-lives. For the entire data set, the forward and backward half-lives only diverge by a matter of 30 000 yr out of 2.76 Myr ­ a discrepancy of only just over 1 per cent. From the data set, it is also possible to calculate the half-lives of the starting class of the objects. The results of this calculation are given in Table 4. The number of clones in any particular class is not necessarily an exact multiple of 729. This is because a number
C

of the objects have outlying clones that actually fall into a different class at the start of the integration as compared with the seed. The results of this analysis again show the dependence of half-life on perihelion position ­ objects in the Jupiter classes have half-lives shorter than those in the Saturn classes, and these in turn are more short-lived than the objects in the Uranus classes. There are also hints in the table that more eccentric objects under the control of any particular planet may be more long-lived than their less eccentric counterparts (compare, for example, the half-lives of U, UN and UE objects). Orbital periods of the more eccentric objects are greater than those on near-circular orbits with similar perihelion distances, and hence encounters with the giant planets happen less frequently. Over the course of the integrations, the clones of each object are repeatedly transferred between classification bins. This allows us to evaluate the amount of time that is spent in each of these classes over the simulation, together with the number of times the object is transferred into that class. From this, we can calculate the mean time that an object spends in any particular class before being transferred into another. The results of such calculations are presented in Table 5. In this table, the value of the mean lifetime for the L or long-period comet class has been ignored, as objects that enter this classification are then removed from the simulation as they pass 1000 au from the Sun. This means that the value of mean lifetime for objects of class L is unrealistically small. It is also noteworthy that the EK and T classes have particularly short mean lifetimes. A stable orbit in these regions requires decoupling from Neptune, and there are no non-gravitational forces within the integrations that could allow this to happen. Hence, the very small number of objects that attain these two classes only do so at the extremes of a series of perturbations and are immediately perturbed back into classes under the control of Neptune.

2004 RAS, MNRAS 354, 798­810


804

J. Horner, N. W. Evans and M. E. Bailey
half-lives of the individual classification bins across which at the start of the integrations. N is the number of clones in class at the start of the integrations, T 1/2 is the half-life (in again the uncertainty on the half-life (in Myr). Direction Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward N 729 729 729 729 891 891 729 729 3519 3519 2205 2205 1917 1917 3780 3780 5094 5094 3006 3006 729 729 T
1/2

Table 4. The the objects fall that particular Myr), and is C JS JU JN S SU SN SE U UN UE N

0.03 0.02 0.02 0.02 0.02 0.02 0.03 0.03 0.02 0.02 0.02 0.03 0.04 0.04 0.06 0.06 0.08 0.09 0.23 0.26 0.46 0.50

Table 5. The mean time (in yr) spent in each classification bin, before a Centaur clone is transferred to another bin. The shortness of this mean time is understandable, as objects close to the boundary are often transferred to and fro. Note that EC stands for Earth-crossing objects and MC for Mars-crossing. Class EC MC E SP I J JS JU JN JE JT S SU SN SE ST U UN UE UT N NE NT EK T Forward lifetime 1960 950 4020 1680 570 290 890 1320 1110 1470 830 1400 3670 2710 3060 2960 6650 4860 4710 1560 2600 4640 4100 330 100 Backward lifetime 1990 890 3480 1630 630 300 890 1310 1120 1470 860 1350 3630 2710 3010 3090 6150 5040 4870 1480 2500 4350 4440 320 100

0.67 0.65 0.61 0.64 0.59 0.63 0.72 0.68 1.03 1.02 1.11 1.22 1.75 1.77 3.72 3.36 5.84 6.19 12.5 14.2 12.5 13.5

4 TRANSFER PR OB ABILITIES It is also straightforward to calculate the probability of an object being transferred from one class to another. This can be visualized by constructing a 24 â 24 grid with the initial class on the vertical axis and the final class on the horizontal (these classes are just initial and final with respect to a single transfer, not for the entire integration). This is performed by recording every transfer that occurs within the integrations, and hence calculating the fraction of objects that, for example, are transferred from class J to class JS. The results are shown in Table 6. The numbers have been normalized so that the sum along any row is unity. For any class, the probabilities give the relative likelihood of leaving from that class to the target classes given on the horizontal axis. As an example, let us take a typical result from one of these tables, namely that the value of the probability of transfer from class J to class JS is 0.49 (Table 6). This means that, for an object in the J class, there is a 49 per cent chance of the object being transferred directly into the JS class the next time the classification changes. We can see that for such an object, the two most likely transfers are to the JS class or to a short-period (SP) comet, and between them, these two possibilities make up the great bulk of transfers for all J objects. Table 6 shows a number of interesting features. Whenever an object is controlled by two planets (one at perihelion and one at aphelion), the classes to which the object is most likely to move correspond to transitions at the perihelion and aphelion of the planet. For instance, an SU object is controlled at perihelion by Saturn and at aphelion by Uranus. It is most likely to be transferred to one of the classes JU or U by an encounter at aphelion. These cases corresponds to an encounter with Uranus either increasing the eccentricity of the orbit, and hence pushing the perihelion down to the control of Jupiter, or decreasing the eccentricity, pulling the perihelion away from the control of Saturn. For encounters at perihelion, the most likely classes are S or SN corresponding either to a circularization of the orbit at Saturn, or to a pumping of the eccentricity of the orbit, as the aphelion moves control from Uranus to Neptune. These most popular transfers can be traced diagonally down the tables, around the empty diagonal corresponding to the same initial as final class. These four parallel lines of high probabilities give the appearance of two sets of `tram lines' running through the tables. After these possibilities, other transfers are also viable, albeit with lower probabilities ­ for example, an SU object can suffer a perihelion­aphelion interchange at Saturn, moving to the JS class. However, the fact that the four classes most likely to be reached in a transfer lie along the `tram lines' vindicates the classification scheme, which is based on the idea of transfers by interaction primarily at perihelion and aphelion. Also of interest is an effect that can be seen on comparing Table 6 with the equivalent results for each individual Centaur (given in Appendix A of Horner 2003). From any class, the probability of transfer to another class is roughly constant, regardless of the direction of integration or the object in question. The main discrepancies lie in very low probability transfers, where the uncertainty is large because of the small numbers involved. This means that for a newly discovered object, it is possible to give the probabilities of its transfer to any new class, as long as the initial class can be computed. It also permits insight into the main dynamical pathways followed by a Centaur. For example, using the values in Table 6 and assuming an initial population of 1000 short-period comets or SP objects, it can be seen that 27 of these objects become E types, all of which would
2004 RAS, MNRAS 354, 798­810

C


C

Table 6. Transfer probabilities for the entire simulation in the forward (upper table) and backward (lower table) direction. Note that the probabilities are nearly the same regardless of the direction of integration. An empty entry means that the transfer probability is <10-3 . The leading diagonal is empty by definition.

Pre J 0.002 0.001 0.304 0.004 0.498 0.144 0.004 0.381 0.003 0.001 0.415 0.007 0.340 0.001 0.051 0.628 0.001 0.045 0.001 0.269 0.001 0.289 0.109 0.336 0.001 0.001 0.303 0.001 0.001 0.673 0.658 0.222 0.481 0.487 0.005 0.471 0.059 0.102 0.001 0.006 0.201 0.004 0.151 0.001 0.077 0.001 0.001 0.606 0.136 0.001 0.320 0.321 0.003 0.002 0.353 0.001 0.006 0.885 0.006 0.098 0.002 0.006 0.366 0.001 0.363 0.006 0.001 0.004 0.251 0.016 0.468 0.174 0.002 0.010 0.198 0.003 0.001 0.001 0.007 0.357 0.018 0.001 0.004 0.343 0.043 0.037 0.022 0.004 0.521 0.038 0.014 0.011 0.009 0.422 0.146 0.001 0.046 0.009 0.083 0.002 0.835 0.271 0.001 0.001 0.001 0.001 JS JU JN JE JT S SU SN SE ST U UN UE UT N NE

Post E

SP

I

L

NT

0.993

0.027

0.001 0.041

0.071

0.005

2004 RAS, MNRAS 354, 798­810 0.006 0.189 0.004 0.001 0.001 0.435 0.006 0.477 0.001 0.001 0.274 0.002 0.151 0.003 0.508 0.167 0.005 0.496 0.002 0.001 0.001 0.003 0.567 0.358 0.430 0.156 0.136 0.308 0.001 0.149 0.004 0.384 0.003 0.001 0.419 0.007 0.345 0.001 0.055 0.641 0.001 0.041 0.001 0.001 0.076 0.001 0.136 0.269 0.001 0.001 0.563 0.330 0.001 0.006 0.094 0.002 0.349 0.006 0.001 0.007 0.362 0.001 0.002 0.004 0.253 0.163 0.002 0.001 0.499 0.017 0.016 0.471 0.005 0.001 0.003 0.346 0.016 0.026 0.010 0.192 0.003 0.421 0.001 0.147 0.044 0.009 0.068 0.004 0.002 0.861 0.256 0.002 0.001 0.001 0.001 0.001 0.038 0.026 0.004 0.540 0.014 0.021 0.013 0.008 0.001 0.001 0.006 0.360 0.002 0.001 0.312 0.292 0.111 0.664 0.218 0.001 0.701 0.470 0.501 0.005 0.449 0.061 0.001 0.304 0.001 0.117 0.001 0.005 0.212 0.003 0.152 0.006 0.868 0.311 0.003 0.001 0.366 0.001 0.005 0.181 0.004 0.001 0.001 0.449 0.005 0.475 0.001 0.001 0.314 0.002 0.159 0.002 0.494 0.159 0.005 0.492 0.001 0.002 0.566 0.374 0.432 0.137 0.125

0.012

0.484 0.198 0.056 0.022 0.004

0.022 0.050

0.009

0.001

E SP I L J JS JU JN JE JT S SU SN SE ST U UN UE UT N NE NT

0.009

0.992

0.027

0.002 0.039

0.005

0.060 0.001 0.483 0.201 0.057 0.022 0.004

0.011

0.020 0.063

0.009

0.001

Simulations of Centaurs ­ I. Statistics

E SP I L J JS JU JN JE JT S SU SN SE ST U UN UE UT N NE NT

0.010

805


806

J. Horner, N. W. Evans and M. E. Bailey
2004). This calculation also only takes into account objects becoming short-period comets for the first time. Given that the mean fade time is 12 kyr, it is reasonable to assume that objects captured for the first time are significantly brighter as short-period comets than those that have experienced a number of prolonged stays in the region. Therefore, the objects that in the simulations display a number of prolonged periods as short-period comets, would actually exhaust all their volatiles early on and should not contribute to the new short-period comet flux in later passages through the inner Solar system. Note that the usage of the flux of short-period comets to normalize the source populations has also been exploited recently by Emel'yanenko, Asher & Bailey (2004) to estimate the total population of trans-Neptunian objects on highly eccentric orbits, which reside at still greater heliocentric distances than the Centaurs. Given our estimate of the total population of Centaurs and the knowledge that the half-life is 2.75 Myr, we can estimate the influx of new Centaurs from the Edgeworth­Kuiper Belt. Neglecting those few objects that could be captured on to Centaur-like orbits from high-eccentricity orbits from the Oort cloud, we can see that 22 150 objects must be replaced every 2.75 Myr. This is equivalent to one object transferred to the Centaur region from the Edgeworth­ Kuiper Belt every 125 yr. This calculation also ignores the small flux of objects from Centaur-type orbits to the Edgeworth­Kuiper Belt. Our simulations allow neither for non-gravitational effects, such as collisions, nor for the gravitational perturbations between Edgeworth­Kuiper Belt objects. So, it is impossible to determine the flow of objects from orbits that encounter Neptune to those that are stable beyond. To drive the inward flux of Centaurs, the effects of collisions and of perturbations between the Edgeworth­Kuiper Belt objects must be considered. Durda & Stern (2000) suggest that collisions of objects greater than 4 m in diameter on to comet-sized bodies within the Edgeworth­Kuiper Belt occur every few days. In addition, they calculate that the time-scale for the disruption of 1 km objects is 1 Gyr. These two facts in concert imply a high rate of collision within the Edgeworth­Kuiper Belt, sufficiently high that an inward flux of one new Centaur every 125 yr seems reasonable. From the forward integrations, we can deduce that 1799 clones (7.7 per cent of the sample) became Earth-crossing and 3799 clones (16 per cent) became Mars-crossing. Very similar numbers are yielded by the backward integrations. Therefore, we expect typically one Centaur to become Earth-crossing for the first time approximately every 880 yr, and one new Mars-crosser every 420 yr. Most of the known population of near-Earth objects (NEOs) is asteroidal in nature, as the Main Belt provides the great majority of NEOs. Morbidelli et al. (2002) state that only 6 per cent of the NEOs are ultimately of Edgeworth­Kuiper Belt origin. NEOs originating in the Main Belt are expected to survive for far longer times within the inner Solar system than the cometary bodies, due to the fact that they lie on orbits significantly decoupled from direct perturbations by Jupiter. An approximate flux of one new Earth-crosser of Centaur (and hence, originally Edgeworth­Kuiper Belt) origin per millennium seems to be in reasonable agreement with such work. Though the number of NEOs of Edgeworth­Kuiper Belt origin is small, they are exceptionally important in judging potential hazards. The size distribution of NEOs of asteroidal origin is heavily weighted towards small particles, consistent with the idea that they are collision fragments. The largest NEO is (1036) Ganymede with an absolute magnitude H = 9.45 and a diameter d 50 km. Very few NEOs are larger than 10 km across. However, for the Centaurs, the upper end of the size distribution is well-populated, with 16 objects of the 32 listed in Table 2 intrinsically brighter than Ganymede. The passage of a large Centaur such as Chiron or Pholus into the
C

return to the SP bin, 41 enter the I class, three of which would on average then return, 304 travel to the J class, 147 making the return trip. For the other bins, we find 84 returning objects from JS, eight from the JU class and one from JN. Therefore, 270 of 1000 objects that leave the short-period class return immediately the next time that their classification changes. Given that the total time spent in each class and the number of times that class becomes occupied are calculable from the simulation data, it is also possible to compute the probability per unit time of a transfer. We already know the probabilities that an object in a particular class will be transferred to any other. Dividing the mean time spent in any class by this probability, we obtain the mean transfer time from one class to any target class. The inverse of this is the probability per unit time of the transfer. The values obtained in this type of analysis are given in Horner (2003) for the individual Centaurs, while the probabilities for the entire data set are given in Table 7. The results are given as probabilities per Myr, so that a value of 0.1 in a particular box means that an object making the relevant transfer would have a mean transfer time of approximately 10 Myr. This means that the population in any class N evolves according to dN = Pi Ni - Pj N j , (4) dt
i j

where P j and N j are the probabilities per unit time and the populations in the bins along the row, while the P i and N i are the corresponding quantities along the column. In other words, the ingress to a particular class is governed by the numbers in the column, and the egress by the row. Mathematically speaking, this gives us coupled sets of exactly solvable linear first-order differential equations that govern the evolution of Centaur clones. We will return to this in a later publication, but these ideas are already prefigured in Bailey et al. (1992). 5 CENT A U R FLUXES AND POPULA TION During the simulations, we also record the numbers of clones that become Earth-crossing objects, Mars-crossing objects and shortperiod comets. This gives us a means to calculate the total population of Centaurs. Fernandez (1985) suggested a flux of 10-2 new ´ short-period comets per year, with a mean lifetime of 6kyr. More recently, Levison & Duncan (1994) find that the mean fade time for a short-period comet is 12 kyr. Under the assumption that the current population of short-period comets is in steady state, the work of Levison & Duncan implies a flux of 0.5 â 10-2 new short-period comets per year. This is equivalent to one new short-period comet being captured, on average, every 200 yr. If we assume that the entirety of this flux comes from the Centaur region, then this allows us to estimate the total population of the Centaur region using the simulations. From the simulation data, a total of 7900 out of 23 328 clones (34 per cent of the initial population) become short-period comets at some point during the forward integrations, and 8068 (again 34 per cent) become short-period comets during the backward integrations. This is a flux of one new short-period comet every 380 yr. If we assume all short-period comets are captured from the Centaur region, an estimate of the total number of Centaurs (with perihelion distance q 4 au and aphelion distance Q 60 au) is 44 300. This represents the population of objects bright enough to be visible as short-period comets, were they to be captured into such an orbit. An effective nuclear diameter d greater than 1 km seems a reasonable limit to place for this value, though there are an increasing number of comets with d 0.5 km (e.g. Lamy et al.

2004 RAS, MNRAS 354, 798­810


C

Table 7. Transfer probabilities per Myr for the entire simulation in the forward (upper table) and backward (lower table) direction. An empty entry means that the transfer probability per Myr is <10 diagonal is empty by definition. . The leading

-2

Pre J 181.12 10.45 86.67 0.29 JS JU JN JE JT S SU SN SE ST U UN UE UT N NE

Post E

SP

I

L

N

247.24 5.24 144.25 1.88 1.35 5.48 323.20 61.89 0.01 36.99 0.08 0.01 0.03 0.10 0.14 3.88 1063.87 0.04 0.07 0.15 49.43 0.01 2.22 156.10 0.28 0.02 96.93 45.80 0.01 0.01 0.14 430.55 102.28 0.01 187.82 1.08 0.01 106.75 1.94 107.07 0.38 0.03 0.03 0.02 0.03 0.10 0.50 193.99 1.84

16.19 0.06 55.36 521.23

0.15 24.33 251.13 0.08 4.89 1711.44 0.02 6.65 627.59 0.07 0.11 13.56

122.96

0.10 9.54

2004 RAS, MNRAS 354, 798­810 0.31 27.69 0.74 21.79 3.95 3.96 260.57 1.29 1454.53 326.64 0.50 0.44 1.10 4.00 170.39 0.01 0.11 160.84 3.02 0.21 0.05 0.01 0.16 0.23 0.57 0.63 51.42 0.23 11.45 150.93 2.98 0.10 0.50 450.08 0.23 0.32 0.74 45.14 0.32 0.34 4.40 170.75 2.76 0.63 0.79 118.60 0.08 0.20 0.74 26.68 0.18 0.19 0.28 5.76 248.55 1.77 0.20 0.35 118.33 0.17 92.48 290.06 2.47 0.50 0.10 260.30 1.62 0.06 0.03 0.01 375.78 4.75 0.24 4.39 26.65 0.65 0.03 0.01 231.08 0.92 0.01 0.34 16.48 0.42 0.01 0.33 25.23 0.28 0.37 204.61 0.39 0.22 0.85 46.19 0.05 0.04 0.05 0.06 0.09 0.38 0.01 87.64 1.23 0.01 0.03 0.31 0.17 0.25 16.45 0.05 0.02 0.02 0.02 0.06 0.18 0.01 0.46 130.33 0.48 0.04 98.95 12.46 21.08 0.22 0.01 0.01 0.04 0.03 0.05 0.11 0.01 99.38 0.24 0.09 50.49 0.15 0.04 0.12 0.02 0.03 0.05 94.71 0.11 0.85 41.35 0.75 0.05 0.01 0.01 0.04 36.76 0.01 0.19 64.34 0.56 0.02 0.04 0.07 0.32 32.16 1.63 195.65 0.64 138.48 77.16 0.01 104.97 0.04 99.46 0.01 29.37 259.35 9.11 1697.47 189.36 0.02 4.93 0.18 90.64 0.55 140.07 55.66 526.67 27.20 0.21 39.48 2.78 3.90 264.54 1.33 1372.92 2100.71 1.21 0.56 1.21 3.86 172.20 0.01 0.03 6.08 4424.25 0.09 0.11 0.30 107.05 0.01 0.02 61.50 0.01 0.02 7.97 6.83 0.76 0.09 5.50 109.06 36.06 1.12 1.33 5.65 331.59 64.19 166.55 2.72 0.19 0.05 0.01 0.21 0.01 293.69 2.50 0.49 0.08 258.23 1.58 0.06 0.02 235.46 0.76 0.02 0.28 15.18 0.41 0.29 25.10 0.29 0.36 182.06 0.09 0.26 0.23 134.76 0.81 182.86 1.89 0.03 0.04 0.21 0.54 0.38 0.46 213.72 0.18 11.37 146.54 2.56 0.06 0.35 474.09 0.18 0.39 0.80 311.65 0.18 0.25 4.09 165.74 2.61 0.65 0.84 123.64 0.01 0.27 0.22 0.48 214.09 0.13 0.14 0.23 5.66 245.12 1.82 0.22 0.36 115.01 37.42 0.08 0.02 0.02 35.84 0.04 0.10 0.14 3.96 1009.83 0.06 0.07 0.14 50.43 0.01 0.09 0.04 0.16 113.88 0.04 0.02 0.02 0.02 0.04 0.18 0.01 85.52 1.08 0.01 0.04 43.21 0.01 0.19 474.46 0.27 0.14 0.27 171.58 0.04 0.04 0.03 0.06 0.06 0.24 0.01 0.40 135.05 0.49 0.03 108.00 0.04 0.03 0.03 0.04 0.14 0.02 0.03 1.96 157.74 0.34 0.05 89.22 99.24 0.23 0.06 53.71 0.12 0.01 0.03 0.17 62.31 0.75 0.23 101.61 12.54 23.31 0.23 0.32 32.71 1.11 197.29 0.01 0.01 0.02 0.01 0.01 96.50 0.08 200.31 1.15 0.03 107.39 1.84 113.18 0.24 0.45 127.54 86.08 0.01 97.38 0.03 92.84 0.01 28.74 0.01 384.13 4.84 0.50 4.56 25.99 0.61 0.04 0.02 0.01 0.01 0.02 0.07 0.01 0.02 0.04 0.01 0.01 0.04 0.01 0.03 0.07 96.90 0.13 0.85 42.19 0.70 0.13 0.02 0.05 0.09

13.90

1661.87 220.88 42.57 19.51 2.62 0.08 0.87 0.11 0.02 0.01

0.02 14.83 60.54

0.09 0.02 0.01 0.03 0.10 10.71

0.05

3.07

E SP I L J JS JU JN JE JT S SU SN SE ST U UN UE UT N NE NT

0.18

284.75

16.43 0.05

0.54 23.88

0.11 8.36

93.76

Simulations of Centaurs ­ I. Statistics

94.91 4.56 1641.73 224.80 43.74 19.99 2.65 0.14 0.84 0.09 0.03 0.01 0.01

0.01 0.02 13.92 72.87

0.01 0.04 0.06 0.11 10.80

0.05

3.11

E SP I L J JS JU JN JE JT S SU SN SE ST U UN UE UT N NE NT

0.15

0.01 0.02

807


808

J. Horner, N. W. Evans and M. E. Bailey
be underestimates of the impact rate on the planets, given that the errors in integration are at their largest when the clone is closest to a massive body. Of course, all the numbers derived in this section are dependent on the assumed flux of one new short-period comet every 200 yr. If the true flux is higher or lower, then the total populations of Centaurs, Earth-crossers and Mars-crossers would need to be correspondingly adjusted upwards or downwards. Although the flux of new shortperiod comets is perhaps uncertain by a factor of 2, it is not uncertain by a factor of 10, and so our population estimates are surely correct to an order of magnitude. 6 D YN AMICS AND PHO T OMETRIC COLOURS As the Centaurs range over a large area of the outer Solar system, it is possible that there could be variations in their observable characteristics (such as colour and light curves) as a function of their position. Observations have been carried out by a number of groups (e.g. Weintraub, Tegler & Romanishin 1997; Peixinho et al. 2001; Hainaut & Delsanti 2002; Bauer et al. 2003). Although the number of Centaurs with well-determined colours is restricted, the situation is rapidly improving. Colours for 19 Centaurs are given in Table 9. They range from some of the reddest objects in the Solar system (e.g. Pholus) to much bluer objects (e.g. Chiron). There appears to be no correlation between colour and half-life within the data set, or indeed between colour and classification (for example the reddest and the bluest objects on the list, Pholus and Chiron, are both controlled by the same planet at perihelion). We expect objects closest to the Sun to have undergone some outgassing, and this would lead to a resurfacing of the object, covering the older, darker and redder material with fresh material from the interior of the object. Conversely, those objects that have not displayed cometary behaviour since the formation of the Solar system are expected to be darker and redder. Other factors that could alter the colours of the Centaurs include impacts between objects. Collisions might activate areas of the surface of the objects, exposing fresh material from the interior, manifesting itself as bluer colours. It follows that photometric observations of the Centaurs do have the potential to provide important clues as to the dynamical history of the object, which cannot be determined by integration alone. A good example of this is the case of Chiron and Pholus. These two objects have similar half-lives, and are controlled by the same planet. We can easily calculate that 60 per cent of the clones of Chiron become short-period comets in the forward integrations, against 30 per cent of the clones of Pholus. Yet, it is impossible to determine whether either object has been a short-period comet in the past, purely by means of the integrations. All that we can calculate are percentage probabilities. However, it is possible that the fact that Chiron is the bluest of the Centaurs listed in Table 9 is caused by its well-known cometary activity, which may indicate that at some time in the past Chiron was a short-period comet. The red colour of Pholus may similarly indicate it has not yet been transferred into an orbit that would activate the surface sufficiently for resurfacing to take place. Hence, we may be observing the surface of an object that has not been active since the birth of the Solar system. Bauer et al. (2003) find that the Centaurs display a wider colour distribution than both Jupiter-family short-period comets and Edgeworth­Kuiper Belt objects. The members of the former population are all active, while the members of the latter are not. An intermediate population such as the Centaurs should at least show
C

inner Solar system would provide a very significant environmental disturbance (Hahn & Bailey 1990), as its fragmentation and possible decay could overwhelm the local space environment with debris and dust. Another interesting set of statistics is the total number of times a class becomes occupied during the integrations. In the forward direction, the Earth-crossers become occupied 23 300 times over 3 â 106 yr. However, we already know that only 1799 of the clones enter this class, so it is obvious that the bulk of objects that become Earth crossing at least once will in fact enter, leave and re-enter the class a number of times (in this case, 13 times). Table 8 summarizes the number of visits to the short-period (SP), Earth-crossing (EC), Mars-crossing (MC) and Encke-type (E) classes that a typical object entering these classes would have. The flux of objects into the Encke class of comets can also be calculated. This is subject to greater errors ­ because of the much lower numbers of objects that become Encke-type comets than become Earth-crossers and because a time-step of 120 d is insufficient to track the orbits of such objects with reasonable accuracy and because the gravitational influence of the neglected terrestrial planets is now significant. For objects in the forward integrations, 303 clones became Encke types, implying a capture rate of one every 5200 yr. The backward integration yields the same capture rate. Typically, an object will become an Earth- or Mars-crosser, while classified as a short-period comet, having been transferred to the inner Solar system by an encounter with Jupiter. Although some objects can be perturbed on to Earth and Mars crossing orbits by the other outer planets (for example, the comets P/1997 T3 Lagerkvist­ Carsenty and P/1998 U3 Jager, discussed in Horner et al. 2003), the ¨ vast majority of such objects are under the control of Jupiter prior to a period of cometary behaviour. This is evident in Table 6. In both the forward and backward integrations, the only classes from which the likelihood of transfer to short-period orbits is greater than 0.001 are the other cometary classes (E, I, L), the classes in which the object is controlled by Jupiter (J, JS, etc.) and the S class. Finally, we can use the simulation data to estimate the impact rate on the giant planets from Centaurs. In the forward integration, we find that 144 objects hit Jupiter, 53 hit Saturn, five hit Uranus and a further five hit Neptune. In the backward integration, these numbers are 135, 48, 5 and 1, respectively. Given that the estimated population of the Centaur region is 44 300, then we expect one impact per 10 kyr on Jupiter, one per 28 kyr on Saturn, and one per 300 kyr on Uranus and Neptune. These numbers are likely to
Table 8. Table showing the mean number of times an object that enters a cometary class at least once will go on to enter and re-enter that class through its lifetime. D gives the direction of the integration (F, forward; B, backward), while N gives the mean number of visits to that class. Here, again EC is Earth-crossing, MC is Marscrossing, SP is short-period and E is Encke-type. Class EC MC SP E D F B F B F B F B N 12.9 13.7 13.7 14.6 28.0 28.4 16.5 18.3

2004 RAS, MNRAS 354, 798­810


Simulations of Centaurs ­ I. Statistics
Table 9. List of the Centaurs with colour information. D is the direction of integration and T are taken from Hainaut & Delsanti (2002), Boehnhardt et al. (2002) and Bauer et al. (2003). Object 2000 EC98 1998 SG35 1999 UG5 Class JU JS SU D F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B F B T
1/2 1/2

809

is the half-life in Myr. The colours

B-V 0.854 ± 0.081 0.725 ± 0.089 0.88 ± 0.18 0.964 ± 0.085 0.75 ± 0.04 0.750 ± 0.040

Colour V -R 0.47 ± 0.06 0.466 ± 0.05 0.42 ± 0.08 0.456 ± 0.050 0.60 ± 0.08 0.607 ± 0.060 0.68 ± 0.02 0.47 ± 0.04 0.513 ± 0.068 0.53 ± 0.02 0.44 ± 0.03 0.53 0.34 0.359 0.36 ± ± ± ± 0.05 0.03 0.027 0.10

R-I 0.50 0.439 0.46 0.546 0.72 0.625 0.58 ± ± ± ± ± ± ± 0.06 0.076 0.04 0.063 0.13 0.042 0.02

0.61 0.63 0.67 0.65 0.74 0.85 0.86 0.75 0.94 0.87 0.95 0.95 1.03 1.07 1.06 1.38 1.28 1.39 1.78 1.52 2.93 3.67 3.18 3.44 4.87 5.65 4.91 6.40 6.37 7.30 11.1 13.1 10.3 9.38 11.5 10.8 11.1 13.1

Asbolus

SN

0.523 ± 0.045 0.47 ± 0.02 0.48 ± 0.02 0.56 ± 0.04 0.356 ± 0.037 0.32 ± 0.10 0.70 ± 0.06 0.81 0.814 ± 0.056 0.740 ± 0.210 0.66 0.41 0.397 0.368 ± ± ± ± ­ 0.64 ± ± ± ± ± ± 0.05 0.06 0.069 0.102 0.10 0.066 0.07 0.18 0.122 0.06

2001 PT13 2001 BL41 Chiron

SU SU SU

0.67 ± 0.06 0.679 ± 0.039

1999 XX143 Pholus

SN SN

1.35 1.19 ± 0.10 1.299 ± 0.099 1.261 ± 0.139

0.67 ± 0.07 0.71 0.78 ± 0.04 0.794 ± 0.032 0.672 ± 0.080 0.74 0.38 0.448 0.474 0.520 0.63 0.77 0.793 0.74 0.41 0.464 0.50 ± ± ± ± ± ± ± ± ± ± ± ± 0.06 0.06 0.044 0.095 0.030 0.12 0.05 0.041 0.08 0.10 0.059 0.07

1994 TA 2002 GO9 2000 QC243 1998 QM107

SU UN U UN

0.724 ± 0.062 0.771 ± 0.100 0.730 ± 0.060 0.88 ± 0.07 1.090 ± 0.040

Nessus

SE

Hylonome

UN

0.643 ± 0.082

0.695 0.64 0.46 0.490 0.52

2002 GB10 Chariklo

UE U

0.802 ± 0.049

1998 TF35 2002 GB10

UE UE

1.085 ± 0.111

0.58 ± 0.07 0.47 0.479 ± 0.029 0.49 ± 0.02 0.65 ± 0.08 0.697 ± 0.064 0.71 ± 0.02

0.62 ± 0.07 0.55 0.542 ± 0.030 0.51 ± 0.02 0.66 ± 0.07 0.651 ± 0.119 0.65 ± 0.02

properties varying between one extreme and the other, and this seems to be borne out by the data in Table 9. 7 CONCLUSIONS Detailed simulations of the evolution of the Centaur population under the gravitational influence of the Sun and the four giant planets (Jupiter, Saturn, Uranus and Neptune) have been carried out. 23 328 test particles were created by taking the orbits of 32 well-known
C

Centaurs and producing nine clones in each of semimajor axis (a), eccentricity (e) and inclination (i), giving a total of 729 clones of each object. The clones were then integrated in both the forward and backward directions for a period of 3 Myr. Any clone that reached a heliocentric distance of 1000 au was deemed ejected from the Solar system and removed from the integration. The half-lives for ejection of the Centaurs were calculated from the simulation data. These ranged from 540 kyr for 1996 AR20 (in the forward direction) to 32 Myr 2000 FZ53 (in the backward

2004 RAS, MNRAS 354, 798­810


810

J. Horner, N. W. Evans and M. E. Bailey
Bauer J. M., Meech K. J., Fernandez Y. R., Pittichova J., Hainaut O. R., ´ Boenhnhard H., 2003, Icarus, 166, 195 Boehnhardt H. et al., 2002, A&A, 395, 297 Chambers J. E., 1999, MNRAS, 304, 793 Clube S. V. M., Napier W. M., 1984, MNRAS, 211, 953 Dones L., Levison H. F., Duncan M., 1996, ASP Conf. Ser. 107: Completing the Inventory of the Solar System. Astron. Soc. Pac., San Francisco, p. 233 Durda D. D., Stern S. A., 2000, Icarus, 145, 220 Emel'yanenko V. V., Asher D. J., Bailey M. E., 2004, MNRAS, 350, 161 Evans N. W., Tabachnik S., 1999, Nat, 399, 41 Fernandez J. A., 1985, Icarus, 64, 308 ´ Groussin O., Lamy P., Jorda L., 2004, A&A, 413, 1163 Hahn G., Bailey M. E., 1990, Nat, 348, 132 Hainaut O. R., Delsanti A. C., 2002, A&A, 389, 641 Holman M. J., 1997, Nat, 387, 785 Horner J., 2003, D. Phil. thesis, Univ. Oxford Horner J., Evans N. W., Bailey M. E., Asher D. J., 2003, MNRAS, 343, 1057 Horner J., Evans N. W., Bailey M. E., 2004, MNRAS, in press Lamy P. L., Toth I., Fernandez Y. R., Weaver H. A., 2004, Comets II, M. C. ´ Festou U. Keller H. A. Weaver (eds), Univ. Arizona Press, Tucson, in press Levison H. F., Duncan M. J., 1994, Icarus, 108, 18 Levison H. F., Duncan M. J., 1997, Icarus, 127, 13 Lowry S. C., Fitzsimmons A., Collander-Brown S., 2003, A&A, 397, 329 Morbidelli A., Bottke W. F., Froeschle C., Michel P., 2002, in Bottke W. F., ´ Jr, Cellino A., Paolicchi P., Binzel R. P., eds, Asteroids III, Univ. Arizona Press, Tucson, p. 409 Oikawa S., Everhart E., 1979, AJ, 84, 134 ¨ Opik E. J., 1976, Interplanetary Encounters. Elsevier, New York Peixinho N., Lacerda P., Ortiz J. L., Doressoundiram A., Roos-Serote M., Gutierrez P. J., 2001, A&A, 371, 753 ´ Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical Recipes. Cambridge Univ. Press, Cambridge, ch. 15 Sanzovo G. C., de Almeida A. A., Misra A., Torres R. M., Boice D. C., Huebner W. F., 2001, MNRAS, 326, 852 van Woerkom A. J., 1948, Bull. Astron. Inst., 10, 445 Weintraub D. A., Tegler S. C., Romanishin W., 1997, Icarus, 128, 456

direction). The half-life of the ensemble of Centaurs was 2.7 Myr, irrespective of the direction of integration. To analyse the simulation data we exploited a classification scheme introduced by Horner et al. (2003), which breaks down the orbits of objects in the Solar system into different classes according to the planets controlling the perihelion and aphelion. The use of the new classification scheme allowed us to determine the dynamical pathways through which the Centaur population evolves. Transfer probabilities were calculated between the different classes. These were found to be remarkably independent of the object under study, offering the prospect of determining possible future histories for objects from the knowledge of their current classification. Both observations and simulations suggest that one new shortperiod comet is produced every 200 yr. This is used to normalize our simulation results. In this case, the total population of the Centaur region (from q 4 au to Q 60 au) with nuclei large enough to provide visible comets (d 1 km) is estimated to be 44 300. A flux of one new object into the Centaur region from the Edgeworth­ Kuiper Belt every 125 yr is required to maintain the population in a steady-state. Additionally, one fresh Earth-crossing object is expected to arise from the Centaur region every 880 yr. This is both of interest and concern, as large Centaurs entering the inner Solar system are likely to fragment with the production of much dangerous dust and debris. In a companion paper (Horner, Evans & Bailey 2004), individual Centaur clones are discussed, showing examples of stable resonant behaviour, Kozai instabilities, capture into satellite orbits, and evolution into objects spending long time periods in the inner Solar system, amongst other things. This gives a feel for the rich variety of behaviour of the diverse objects known as `Centaurs'. A CKNO WLEDGMENTS JH acknowledges financial support from the Particle Physics and Astronomy Research Council. We thank the referee S.C. Lowry for a careful reading of the manuscript. REFERENCES
Asher D. J., Steel D. I., 1993, MNRAS, 263, 179 Bailey M. E., Chambers J. E., Hahn G., Scotti J., Tancredi G., 1992, Brahic A., Gerard J.C., Surdej J., eds, Proc. 30th Liege Int. Astrophysical Coll., Univ. Liege, Liege, p. 285 ` `

A This paper has been typeset from a TEX/L TEX file prepared by the author.

C

2004 RAS, MNRAS 354, 798­810