Monte Carlo simulations of the Oort comet cloud

ICARUS

10’4-121 (1990)

88,

Monte Carlo Simulations

of the Oort Comet Cloud

JULIA HEISLER’ Steward Observatory. University of Arizono. Tucson, Arizonrr 85721 Received

January

5. 1990; revised

May 31. 1990

The evolution of the Oort comet cloud between 10,000 and 40,000 AU is simulated using Monte Carlo techniques. In total, over 40 million comets are evolved for a period of 250 million years and 1.5 million comets for 4.5 billion years. Perturbations from both the Galactic tidal field and passing stars are included. The loss cone is filled for semimajor axes ~30,000 AU, so that comets at greater distances do not contribute to showers. Over the age of the Solar System the comet cloud beyond -20,000 AU loses about half its members. If comets are distributed in density as rm3.’ beyond r = 3000 AU, then the background flux of new comets within 2 AU of the Sun should be -0.2 comets/year for a cloud of 10” comets beyond 3000 AU, with 2% of the time spent in showers, events where the comet flux is more than five times the background rate. Furthermore, the energy distribution of new comets outside of showers should be strongly peaked at l/a = 35 x 1O-6 AU-‘. cc)1990Academic Press, Inc.

1. INTRODUCTION

Cometary dynamics has become a very active field in the last several years. Problems currently being investigated include the formation of the Oort cloud, the possible existence of a massive inner cloud, and the evolution of comets to short period orbits (e.g., Bailey and Stagg 1989, Fernandez 1980a, Duncan et al. 1987, 1988). Tidal torques induced by the smooth distribution of matter in the Galactic disk are now recognized to be more important than stellar perturbations (Heisler and Tremaine 1986, Morris and Muller 1986, Smoluchowski and Torbett 1984, and Torbett 1986). Furthermore, the proposal by Alvarez et ~1. (1980) that the impact of an asteroid caused the extinction of the dinosaurs 65 million years ago has sparked much interest in comet showers. Others, (e.g., Hsu 1980) have suggested the impact of a comet instead of an asteroid. It has been speculated that such ’ Currently at the Lunar and Planetary sity of Arizona. Tucson, AZ 85721.

Lab,

Univer-

events can be caused by a companion star (Davis et al. 1984, Whitmire and Jackson 1984), a tenth planet (e.g., Whitmire and Matese 1985), giant molecular clouds (e.g., Bailey 1983, Clube and Napier 1982, Hut and Tremaine 1985, Napier and Staniucha 1982, Rampino and Stothers 1984, van den Bergh 1982), and finally, by random passing stars (e.g., Fernandez and Ip 1987, Heisler et al. 1987, Hills 1981, Hut and Tremaine 1985). Weissman’s reviews (1985a, 1988, 1989a) sort out the many papers that have been written on all of these subjects. Much of this work would not have been possible without computer simulations. In particular, Monte Carlo methods for modeling the Oort cloud have proved to be most fruitful. Schreur (1979) was the first to model the Oort cloud with a Monte Carlo scheme. Later, Weissman (1979, 1982, 198_5b), Fernandez ( 1980b, 198 1, 1982) and Remy and Mignard (1985) used similar approaches to evolve a cloud of comets over the age of the Solar System. These early papers were mostly concerned with the formation of the Oort cloud from comets that originated be104

0019-1035/90

_

Cowrieht

0

All nghts

of

$3.00

1990 bv Academic

reproduction in

any

Press. form

Inc. reserved

SIMULATIONS OF COMETS tween Neptune and Uranus, in an extended solar accretion disk, or in situ in the Oort cloud itself. However, only small numbers of comets could be treated. About 20,000 comets were evolved in Weissman's simulations, less than I0,000 in Fern~indez' work, and only 200 in Remy and Mignard's models. Also, these early works did not include the Galactic tidal field [Duncan et al. (1987) have included the Galactic tide in their treatment of the formation of the comet cloud.] Heisler et al. (1987; hereafter HTA) were able to treat a comet cloud that was more than 100 times larger than any previously considered by employing a "cloning" technique that had previously been used to describe stellar orbits in a globular cluster with a central black hole (Shapiro and Marchant 1978). The two problems are similar because in both one is interested in the evolution of orbits in a relatively small region of phase space, namely near the event horizon of a black hole or near the planetary region of the Solar System. One is less interested in the evolution of orbits away from these regions; however, since the orbit distribution in phase space is roughly uniform, in a straightforward treatment one would be forced to spend the majority of time following uninteresting orbits. The solution is to follow " t o k e n s " instead of comets, where a token represents a variable number of comets. The distribution in phase space can be uniform for comets but not for tokens, the majority of which are assigned orbital perihelia near the planetary region. With this technique HTA simulated a total of 6 million comets while following only about 39,000 orbits. HTA integrated comets over a 300 myr time span to investigate the temporal variation in the flux of new comets. Most of their calculations incorporated stellar perturbations only. They found that the flux of comets into the inner Solar System was highly variable. About every 100 myr a 1-3 M o star passed through the Oort cloud and induced a shower of comets 10 times more intense than the average background level, for semi-

105

major axes below 30,000 AU. Even outside of showers the background flux could fluctuate by about a factor of two. This paper is a sequel to HTA; it includes the Galactic tide within essentially the same Monte Carlo scheme and presents simulations with greatly improved statistics. The tide was included by HTA in only a small number of runs and only for two values of the semimajor axis. In this paper a new way is presented to incorporate the tide more efficiently into the calculation, and through the use of an SCS40 near supercomputer a comet cloud is evolved that is seven times larger than that computed by HTA. Several long-term runs that last the age of the Solar System have also been completed. This paper is organized as follows. Section II describes the numerical model and the changes made to the code used by HTA. Section Ilia describes the results for a 250 myr evolution of 40 million comets between I0,000 and 40,000 AU. The long term, 4.5 Gyr, evolution of a smaller cloud with initially 1.5 million comets is described in Section IIIb. Section IIIc discusses the expected new comet flux from the entire Oort cloud and the expected energy distribution of new comets. Section IV has conclusions. II. NUMERICAL MODEL The model follows a large ensemble of randomly oriented comets, evolved for typically 270 myr subject to the perturbations by random passing stars and the Galactic tidal field. " S t a n d a r d " interstellar clouds are unlikely to perturb the Oort cloud significantly due to their low density (Morris 1985). Encounters with giant molecular clouds (GMCs) every billion years or more can significantly alter the orbits of comets with semimajor axes a >~ 25,000 AU (Weinberg et al. 1986). Weinberg et al. (1986) find that the half-life of a wide binary with initial semimajor axis a 0 is -i-it2 60(a0/pc)- 1.4myr; for a0 = 25,000 AU, r~12 = 1.1 Gyr. Others (e.g., Bailey 1986, Shoemaker and Wolfe 1986, Wasserman 1988) arrive at similar conclusions, namely that the half-life for =

106

JULIA HEISLER

a = 25,000 AU is of order l to 2 Gyr. H o w ever, Hut and T r e m a i n e (1985) find that the half-life is s o m e w h a t longer, about 3 Gyr. G i v e n the large uncertainties in the parameters of such clouds, it is impossible to say precisely what the half-life is for Oort cloud comets. Yet, the observational argument that a large fraction of stars exist in wide binaries of similar size argues in f a v o r that G M C ' s have probably not disrupted the Oort cloud o v e r the age of the Solar System (Hut and T r e m a i n e 1985). F u r t h e r m o r e , o v e r the shorter timescales we mostly consider, - 2 7 0 m y r , the G M C s should be ignorable, although the unlikely possibility of a close e n c o u n t e r with a G M C during this interval cannot be completely excluded. Certainly o v e r long timescales, more than a billion years, G M C s must be included in any c o m p l e t e treatment of the dynamical evolution of the c o m e t cloud. We will c o m e back to this issue in Section I I I b which describes simulations lasting 4.5 Gyr. A brief description of the details of our numerical model follows, but for further details see H T A . a. P l a n e t a o , P e r t u r b a t i o n s a n d the Loss Cone

Once a c o m e t ' s orbital perihelion enters the inner Solar System its orbit is strongly affected by the giant planets and in particular, by Jupiter. While planetary perturbations are not followed in detail, they are crudely represented by removing from the simulation all comets that pass the Sun with a perihelion q < q , -= i0 AU. Fernfindez (1981) has shown that within this critical distance c o m e t s on their first passage through the inner Solar system, averaged o v e r inclination, are likely to be r e m o v e d from the Oort cloud either to hyperbolic orbits or to orbits with much shorter period. While this crude p r o c e d u r e introduces some inaccuracies [e.g., it underestimates the duration of c o m e t showers by about a factor of two (Fern/mdez and lp 1987) because c o m e t s not immediately ejected are able to m a k e subsequent multiple returns], it greatly simplifies the calculation and is a valuable shortcut.

In later sections it will be convenient to describe these c o m e t s as entering the " l o s s c o n e . " The loss cone is that part o f phase space where c o m e t s have perihelia < q , or angular m o m e n t u m less than, J , = [/x(2q, - q2,/a)] ~ 2q,/x, where /x -= GM o and a is the semimajor axis. The larger is q . , the larger is the loss cone and the stronger are the planetary perturbations. C o m e t s that pass through perihelion while in the loss cone are a s s u m e d to be r e m o v e d from the Oort cloud and hence are r e m o v e d from the simulation. The loss cone is said to be " f u l l " if Galactic and stellar perturbations are large enough to send c o m e t s directly to all parts of the loss cone, giving a uniform perihelion distribution; this tends to be true during showers or for c o m e t s with large semimajor axes. The loss cone is " e m p t y " for c o m e t s with small semimajor axes that are only slightly perturbed by the Galactic tide and stars and thus penetrate only a small distance into the loss cone. Inclination does affect the critical perihelion, q , . C o m e t s on high inclination orbits, with respect to the ecliptic plane, suffer smaller changes in energy than c o m e t s on low inclination orbits. Duncan et al. (1987) find that to yield a root mean square energy change of < A(l/a) 2 > i n = 10 4 AU I requires a perihelion of 14 AU for a c o m e t with inclination i -< 30 ° and 6 AU for a c o m e t with 150 °-< i-< 180 ° . The value q , = 10AU chosen here represents a mean o v e r inclinations. Allowing q , to vary would affect most the inclination distribution of c o m e t s that encounter a relatively full loss cone. The relative n u m b e r of low inclination c o m e t s would decrease because they must traverse a larger loss cone to be seen from the Earth whereas high inclination c o m e t s would increase for the opposite reason. b. C o m e t s a n d Cloning

The orbits of c o m e t s are described with respect to a standard nonrotating Cartesian set of axes (2, 9, ~) with the Sun at the origin and ~ perpendicular to the Galatic disk. The S u n ' s peculiar motion with respect to the Local Standard of Rest. which is equal to

SIMULATIONS OF COMETS 16.5 k m / s e c (Mihalas and Binney 1981) is neglected. (While including the S u n ' s motion would increase the f r e q u e n c y o f stellar encounters, it would also w e a k e n by a similar factor the effect of individual encounters. Such quantities as the rate o f change of energy and the half-life due to strong encounters are to first order unaffected by Solar motion. H o w e v e r , the cumulative effects of w e a k encounters do decrease with relative velocity (c.f. Bahcall, et al. 1985).) All c o m e t s begin with semimajor axis a = a 0, with 10,000 -< a 0 -< 40,000 A U , where a 0 is in discrete 10,000 A U steps. The squared eccentricity e 2 = 1 - JZ/(lxa) is assigned randomly from a uniform distribution b e t w e e n 0 and e2, -= 1 - 2q,/a, appropriate for a uniform distribution of c o m e t s on the energy h y p e r s u r f a c e o f phase space. The u p p e r limit ensures that no c o m e t is initially inside the loss cone. The inclination angle i is chosen so that cos i is uniformly distributed between - 1 and 1, appropriate for a spherical Oort cloud. The remaining three orbital elements relative to the Galactic plane, to, the argument o f perihelion, ~ , the longitude of ascending node, and e, the mean a n o m a l y , are chosen randomly from a uniform distribution b e t w e e n 0 and 2~-. As in H T A this p r o g r a m follows tokens instead of c o m e t s , where a token represents a variable n u m b e r of actual comets. T o k e n s near the loss cone represent a single c o m e t , whereas tokens far from the loss cone represent m a n y comets. In this w a y one can concentrate on the interesting orbits that are near the loss cone, while still maintaining a uniform distribution in e 2 for actual comets. The total range in 1 - e z outside the loss cone is divided into seven bins. C o m e t s in bin I are closest to the loss cone and h a v e a t o k e n / c o m e t ratio of 1 : 1. In e v e r y successive bin a token represents three times m o r e c o m e t s than in the previous bin. Thus, comets in bin 7 are the least eccentric and have a t o k e n / c o m e t ratio of 1:729. If a token crosses into a m o r e eccentric bin, it is " c l o n e d " into the three tokens; if three tokens pass into a less eccentric bin they are " d e c l o n e d " into one token. Cloned tokens

107

are given the same values of to and cos i as their " p a r e n t s " since these two variables are directly affected by the Galactic tide. Decloned tokens are given the same values o f a, j2, to and cos i as one, randomly selected, parent. H o w e v e r , new values of and 12 are randomly assigned to both new cloned and decloned tokens. (See H T A for a fuller explanation.) This p r o c e d u r e was adapted f r o m one described by Shapiro and Marchant (1978) to study stellar d y n a m i c s in a globular cluster with a central black hole. One simulation, or set, consists of 1400 tokens followed for 270 m y r , representing initially a total of 218,600 actual comets. Typically 50 sets were c o m p u t e d for each semimajor axis, corresponding to a total of 10,930,000 comets. c. P e r t u r b a t i o n s - G a l a c t i c

Tide a n d Stars

In H T A a few simulations incorporated the Galactic tide by directly integrating the exact equations of motion for all comets. In Cartesian coordinates these equations are: )~ =

--/x

r 3

x

y = ~y

(I)

"~ = - ~ - z - 47rGpz, w h e r e r 2 = x 2 + y2 + z 2 a n d p = 0 . 1 8 5 M o pc -3 (Bahcall 1984) is the m e a n density of matter (stars, gas, and dark matter) in the Solar neighborhood. (These equations are a p p r o x i m a t e in that they neglect the solar motion, radial derivatives of the galactic potential, and n o n a x i s y m m e t r i c terms.) This p r o c e d u r e p r o v e d to be costly in c o m p u t e r time. This difficulty has been o v e r c o m e in the following manner. Heisler and T r e m a i n e (1986) showed that if the orbital evolution due to the tide is slow c o m p a r e d to the orbital period, one can orbit a v e r a g e the equations of motion and derive the following differential equations for the evolution of angular m o m e n t u m , J, and argument of perihelion, to:

108 &=

JULIA HEISLER 27rGpLZ[Jj /,2 j2

( j 2 + 5(L 2 _ j2) sin2co)

+ ( j 2 _ j z ) ( j _ 5J sin 2 to)J

(2)

affect the temporal variation of flux into the loss cone. Excluding low mass stars increases the code's speed by an additional 33%. d. M o n t e Carlo P r o c e d u r e

J

--

57rGp L 2 /./.2

j 2 ( J 2 _ j ~ ) ( L 2 _ j 2 ) sin 2oJ,

where L ~ (p.a) J/2. When a token passes through perihelion we update its orbit according to Eq. (2) integrated for one orbital period. This procedure works remarkably well even for comets on highly eccentric orbits, even though the equations are singular when J = 0. Nevertheless, to ensure that the flux into the loss cone is determined accurately, the exact equations (1) are integrated for the most eccentric orbits, namely those tokens in bin 1 nearest to the loss cone. The remaining tokens, in bins 2 to 7, are handled using Eq. (2). This new procedure for incorporating the effects of the Galactic tide executes about three times faster than the old procedure of integrating the exact equations of motion. The importance of the Galactic tide in determining the evolution of the Oort cloud depends on the value of p. The value used here, based on Bahcall (1984), is not without uncertainty. Kuijken and Gilmore (1989) argue that p = 0.1 M~jpc ~. Lowering p reduces the comet flux in direct proportion, for an empty loss cone, and in turn raises somewhat the value of the semimajor axis where comets begin to see a filled loss cone. H o w e v e r , even ifp is as low as Kuijken and Gilmore (1989) claim, the Galactic tide will still dominate o v e r discrete stellar encounters in its ability to perturb Oort cloud comets by a factor of three (see Section llIa). The kinematic properties of the stars were derived from the standard Galaxy model of Bahcall and Soneira (1980) and from velocity data given by Mihalas and Binney (1981). These data are summarized in Table I of HTA. In this work stellar masses below 0.2 M o are excluded because encounters by stars of such low mass do not noticeably

The comet cloud was evolved according to the following steps (a more detailed description is contained in HTA): (1) Compute the time, At, to the next stellar encounter according to the probability distribution p(>At) = exp(--~ZAt), where the mean encounter frequency is sc = 8.8 myr ~. (2) Determine the particular species of star involved in the encounter by choosing a mass at random from the Solar neighborhood mass function. (3) Randomly choose the star's entry point 7~ into a 1 pc sphere centered on the Sun. (4) Choose the star's velocity vector. Let v 3 be the component of the velocity antiparallel to 7~; then v3 > 0 must have a probability density proportional to v 3 e x p [ - v ~ / (2cry.)], where o-2 is the one-dimensional velocity variance of stars of species i. The components perpendicular to F~, namely v~ and v2, are distributed with a density proportional to e x p [ - ( v ~ + v~)/(2cr2)]. (4) Advance each comet's orbit in time by At. Comets in bin 1 are advanced by directly integrating Eq. (1) while the rest are advanced by replacing the mean anomaly e by + 2~rAt/Por b, solving Kepler's equation and calculating the comet's new (x, y, z) position. The comet is removed from the simulation and recorded in a separate file if it passes through perihelion with q < q , . Comets not in bin 1 that pass through perihelion with q > q , are updated for Galactic tidal effects using Eq. (2) integrated for one orbital period. (5) Compute the change to the comet velocity relative to the Sun due to the stellar encounter using the impulse approximation.

109

SIMULATIONS OF COMETS TABLE I a(AU)

q < 2 AU Rb

10,000 20,000 3,000 40,000

0.003 0.31 1.61 0.77

o-t

S

R jR o

0.008 181. 8 × 10-5 0.23 0.60 0.044 0.33 -0.63 0.17 -0.62

q < 10 AU Weak

Strong

Rb

3 . 2 % 12.8% 3.55 3.6% 2 . 4 % 8 . 2 9 --8.11 --3.85

o-t

S

Rb/R o Weak

Strong

1.03 1.57 1.46 0.70

7.7 ----

0.018 1 . 6 % 0.23 -0.63 -0.62 --

0.4% ----

0.001 0.04 0.20 0.20

N o t e . R b gives the background rate in comets/year normalized to a total cloud population of 10II at each semimajor axis, and R0 is the rate for a filled loss cylinder (see text); o-t is the temporal dispersion in the flux, from bin to bin, in the same units. S is the ratio of the mean shower intensity to the background, and r is the ratio of the inside 2 and 10 AU fluxes. "Weak" and "strong" refer to the percentage of time spent in weak and strong showers; when no showers occur (as defined in the text) the columns for S, weak and strong are blank.

(6) C o m p u t e n e w o r b i t a l e l e m e n t s for e a c h c o m e t . R e m o v e f r o m the s i m u l a t i o n c o m e t s that h a v e e s c a p e d the S o l a r S y s t e m or h a v e a > 0.5 pc. (Orbits with s u c h large a are u n c o m m o n b u t are u s u a l l y s h o r t - l i v e d , a n d e v o l v e q u i c k l y e n o u g h that the u s e o f o r b i t - a v e r a g e d e q u a t i o n s is i n a p p r o p r i a t e . ) T h e a b o v e steps are f o l l o w e d r e p e a t e d l y until the final t i m e , t y p i c a l l y 270 m y r , is r e a c h e d . S i n c e the c o m e t s are m a s s l e s s , ind e p e n d e n t particles, sets with different c o m e t s b u t the s a m e s e q u e n c e o f stellar enc o u n t e r s c a n be a d d e d t o g e t h e r . T h u s , it is p o s s i b l e to s i m u l a t e a c o m e t c l o u d o f arbit r a r y size. All s i m u l a t i o n s d e s c r i b e d in this p a p e r are b a s e d o n a single s e q u e n c e o f stellar e n c o u n t e r s , l a s t i n g for 170 m y r ( S e c t i o n I l i a ) or 4.5 G y r ( S e c t i o n IIIb). tli. RESULTS a. Indioidual

Semimajor

Axes

Fifty sets (each o f 1400 t o k e n s r e p r e s e n t ing 218,600 c o m e t s ) w e r e c o m p u t e d at each v a l u e o f the initial s e m i m a j o r axis at) = 10,000, 20,000, 30,000, a n d 40,000 A U , using d i f f e r e n t initial c o m e t o r b i t s b u t the s a m e s e q u e n c e o f stellar e n c o u n t e r s . I n total, 10,930,000 c o m e t s o r 70,000 t o k e n s were e v o l v e d at e a c h a0.

T h e r e s u l t s are s h o w n in Figs. 1-4 a n d T a b l e 1. T h e figures give the n u m b e r o f c o m ets p e r m i l l i o n y e a r s that pass w i t h i n 2 or 10 A U o f the S u n . T h e right h a n d vertical axis gives the c o m e t flux p e r y e a r n o r m a l i z e d to a total c l o u d p o p u l a t i o n o f 10 lj c o m e t s in c o m e t s / y e a r , at that s e m i m a j o r axis. A c t u ally the figures are h i s t o g r a m s , g i v i n g the total n u m b e r o f c o m e t s f r o m all 50 sets arriving in s u c c e s s i v e m i l l i o n y e a r time i n t e r v a l s with 250 time i n t e r v a l s s h o w n . (The first 20 m y r o f e a c h s i m u l a t i o n is t h r o w n out to r e d u c e the effect o f initial t r a n s i e n t s . ) T h e e r r o r b a r s , d r a w n a b o u t e v e r y 20 m y r , indicate the ltr statistical e r r o r in the result. ] T h e statistical e r r o r is 5 - 1 5 % o f the b a c k g r o u n d flux (defined b e l o w ) , for q < 10 A U at all s e m i m a j o r axes, a n d 7 - 2 0 % o f the b a c k g r o u n d for q < 2 A U a n d a --- 30,000 AU. T h e statistical error, tr, in the flux at a g i v e n t i m e is larger t h a n e x p e c t e d for p u r e P o i s s o n statistics b e c a u s e o f the effects o f c l o n i n g . I f the n u m b e r o f c o m e t s / m y r that If ns equals the number of comets entering in a specific I myr interval in set i and ~-= (1/50) Y~°I ni then o-2 = (50/49)(~0= I n] - 50 n-').

110

JULIA HEISLER

enter the planetary system in a given time interval follows a Poisson distribution and ni is the number of comets entering in a particular myr interval in set i, then the ratio of V -= (1/49) ~ i = l ( n i - n ) 2 t o ~ - = ( I / 5 0 ) ~50 ~=jn~ should be unity. 2 Instead this ratio is typically between 2 and 4 for q < 2 AU and between 3 and 9 for q < 10 AU. During showers the ratio can be as high as 18. The cloning technique enhances the statistical error both during showers and in background, for two related but distinct reasons. During showers the stellar encounter responsible for the shower is strong enough to kick tokens from higher angular momentum bins directly into the loss cone, and these tokens can be cloned by as much as a factor of 33 = 27. Outside of showers there is another problem. Near the loss cone tokens evolve rapidly in their angular m o m e n t u m due to the Galactic tide. Tokens that are drawn into lower angular m o m e n t u m bins near the loss cone are cloned, and the newly made clones fall together into the loss cone before stellar perturbations have a chance to randomize their orbits. Despite these problems, the statistical errors in the results in Figs. 1-4, as determined directly from the set-to-set variation, are small enough that most of the temporal structure is real within that sequence of stellar encounters. Table I gives some statistics for each semimajor axis. The background rate, Rb(q < 10), is the number of comets from a 10 ~ comet cloud that pass within q = I0 AU of the Sun every year; it is computed by

2 Strictly speaking, ev en without cloning the flux is not e x p e c t e d to follow pure Poissonian statistics because the n u m b e r of c o m e t s near the loss cone varies, especially for small s e m i m a j o r axes. H o w e v e r , most of lhe variation o ccurs when c o m p a r i n g one million-year interval with another, and less variation occurs when c o m p a r i n g one set with a n o t h e r but for the same million-year interval b e c a u s e the same s e q u e n c e of stellar e n c o u n t e r s is used for each set. Thus, the ratio of V to c o m p u t e d for a particular myr interval will still app r o x i m a t e l y follow Poisson statistics, with some deviations due to the fact that the n u m b e r of comets near the loss cone does still v a r y s o m e w h a t from set to set.

averaging the 225 lowest fluxes from the 250 measured 1 m y r intervals and then converting the flux to years i and normalizing to a 1011 cloud population at this semimajor axis. A shower is defined as an interval where the comet flux over 1 myr exceeds 5Rb; a strong shower exceeds 10R d, and a weak shower is between 5 R b and 10R d. Rb(q < 2) is similarly defined as the background flux for q < 2 AU. The ratio, r, of the inside 2 AU and 10 AU background fluxes measures the fullness of the loss cone; the loss cone is empty if r = 0 and full if r = 2 AU/10 AU = 1/5 (i.e., comets are uniformly distributed in phase space; dn/dq is constant if the energy hypersurface is uniformly filled). During showers a larger fraction of phase space is filled and r is larger. For example, for showers w i t h a = 10,000AU, 0.03 ~ r ~< 0.13, whereas in background r = 0.001. The variability factor, S, is the ratio of the mean intensity during showers to the background. The temporal dispersion in the background flux is given by o-t; it is the dispersion in the 225 lowest fluxes and is thus a measure of how the flux varies from one million years to the next. 0-1 and to a lesser extent S are strongly dependent on the arbitrary choice of 1 myr bins. Larger bins would decrease the variation described by these variables. Also given in Table I is the fraction of time spent in strong and weak showers. Figure I shows the comet flux for a 0 = 10,000 AU. Over the course of the 250 million years shown for q < 2, among the comets that enter the planetary system, almost all (99.4%) enter during showers. (There are no error bars drawn in the q < 2 AU plot because there were too few background comets from which to determine the error with the algorithm given above.) For q < l0 AU, most comets enter outside of showers; shower comets account for only 12% of all new comets. For both q < 2 AU and q < i0 A U the strongest showers o c c u r at t - 86 myr and t = 192 myr. Each of these events is caused by a close stellar encounter. The first at t = 86 myr is due to the passage of a 2.27 M~:~star at a speed of 31 km/sec with

SIMULATIONS OF COMETS a o = 10,0o0 AU; 5 0 s e t s ; q < 2 500 i

'

I

'

'

r-~

,

,

,

~

r

a# = 10,000 AU; 5(~ s e t s

AU

,

111 q < 10 AU

b

r

'

'

'

'

C-''l

I

'

~ 1

'

'

'--'

I

'

'

'

'

r'

'

'

~

38.6

3,7 3000

27.4

L

2.7

18.3

t

10oo

-

~9.1

0.9 /

,

50

i

1O0

i

i

i

i~

150

200

~ J

200

T i m e (Myr)

i

50

I O0

150

:

i

i

;

200

i

i

i

i

I

2,50

T i m e (Myr)

FIG. 1. (a) Flux of new comets with perihelia q < 2 AU and semimajoraxis a0 = 10,000AU. Fifty sets were computed, amounting to a total cloud population of 10,930,000 comets at this semimajor axis, or 70,000 tokens. The left band axis gives the flux per myr from a cloud with 10,930,000 comets at this semimajor axis. The right hand axis converts this flux to comets/year normalized to a 10It cloud. Nearly all comets enter the inner Solar System during "showers." (b) q < 10 AU; error bars drawn about every 20 myr indicate the lcr statistical error.

a closest approach distance of 15,800 AU from the Sun. The second shower is caused by a 1.19 M G star travelling at 8 km/sec and reaching a minimum distance o f 22,000 AU from the Sun. As in HTA, the strongest showers are created by the close passages of stars of 1-3 M o. The background flux in Fig. la is nearly zero; almost never does a 10,000 A U comet pass within 2 A U of the Sun except in showers. For q < 10 AU the background flux is larger, but only 10% of the intensity of the strongest showers. The temporal variation, as indicated by o-t, is about three times greater than the statistical uncertainty, as indicated by the error bars, and is due to the irregular nature o f the stellar encounters. For q < 2 A U 16% of the time is spent in showers with 3.2% in weak and 12.8% in strong showers. H o w e v e r , as a percentage of all new comets, weak showers account for 0.8%, strong showers 98.6%, and the background for 0.6%. Hence, the overwhelming majority o f q < 2 A U , a = 10,000 AU comets enter during strong showers. Figure 2 shows the flux of comets with a0 = 20,000 AU. These comets have larger lever arms and are more easily perturbed

into the inner Solar System by external perturbers; hence, the flux is higher. The peaks at t = 86 and 192 myr are still visible in both plots. H o w e v e r , for q < 10 AU the intensity at these times is no longer at least five times the background level, and thus, the peaks are not classified as showers. The temporal dispersion, o-t is about 20% of the background (for q < 10 AU) and is about twice the statistical error. For q < 2 A U the variability factor S decreases from 181 for a = 10,000 AU to 0.60 for a = 20,000 AU. (S is not defined for q < 10 AU since there are no true showers.) The percentage of time spent in weak and strong showers at this semimajor axis (a = 20,000 AU) is 3.6% and 2.4% respectively, for q < 2 AU. As a percentage of new comets weak showers account for 15%, strong showers 19%, and the background the remaining 66%. Hence, by 20,000 A U , most comets enter the inner planetary system outside of showers, but a large fraction are shower comets. In Figs. 3 and 4, which display comets with initial semimajor axes a 0 = 30,000 AU and 40,000 AU, respectively, the showers are clearly missing. Instead, there is a vari-

112

JULIA HEISLER a o = ~o,ooo " T ~ T ,

1

r~

~ ] - r

AU; 50 l

T

sets

q

[ T--

% : 20.000 AU: so sets; q < lo AU

2 AU T

I

[ rT~.T.

]

I

[

r

]

r

f

I

I

•

T

~ T

~

I

I

T ~ T ~ ' ~

T

! 500

t

t

18.3

t

7

400

Z

3.7 13.7

300

9.1 200

'

50

I O0

150 Tir, Le ( M y r )

~00

il.8

850

50

tO0

150 Time (Mvr)

~5o

200

FIG. 2. Same as for Fig. I but for a0 = 20,000 AU. Showers disappear for q < 10 AU.

a b l e b a c k g r o u n d flux, a n d t h e t e m p o r a l disp e r s i o n is a b o u t t w i c e t h e s t a t i s t i c a l e r r o r . T h e p e a k s in t h e flux in t h e s e d i a g r a m s d o n o t c o r r e s p o n d in t i m e to p e a k s in Figs. 1 a n d 2. F o r a < 20,000 A U t h e loss c o n e is e m p t y , a n d t h u s r ~ 0.2. S h o w e r s a r e p r o m i n e n t a g a i n s t t h e b a c k g r o u n d flux, a n d t h e t e m p o ral d i s p e r s i o n is high (ort ~ 0.2R b f o r q < 10 A U a n d o-t >~ 0.75R b f o r q < 2 A U ) . B y a = 30,000 A U t h e l o s s c o n e fills, r = 0.2, the flux is at its m a x i m u m a n d t h e s h o w e r s h a v e d i s a p p e a r e d . T h i s r e s u l t is r o u g h l y c o n s i s tent with that obtained by Heisler and Trema i n e (1986) w h o p r e d i c t e d t h a t t h e loss c o n e

ao

30.00(~ Atr 5{) ~et~ q ,

s h o u l d fill b y a = 31,000 A U . A n o t h e r m e a s u r e o f t h e f u l l n e s s o f t h e l o s s c o n e is t h e r a t i o o f R b to R 0, w h e r e R 0 is t h e flux f o r a filled l o s s c o n e a s s u m i n g a u n i f o r m d e n s i t y in 1 - e 2 f r o m 0 to 1 ( s e e E q . 8 in H T A ) . O p e r a t i o n a l l y R 0 is c a l c u l a t e d f r o m t h e d e n sity o f n e a r c i r c u l a r o r b i t s , in o t h e r w o r d s , f r o m c o m e t s t h a t a r e f a r a w a y f r o m t h e loss c o n e . R a t i o s o f R b to R 0 a r e g i v e n in T a b l e I. B y a = 30,000 A U t h e r a t i o s a s y m p t o t e to 0.6. T h e r a t i o d o e s not r e a c h 1 as w o u l d be expected because the density even by a = 40,000 A U is n o t u n i f o r m o v e r t h e e n i r e range 0 ~ 1 e 2 --< 1. T h e c o m e t d e n s i t y is still d e p l e t e d n e a r t h e l o s s c o n e r e l a t i v e to -

ao

:~ AIJ

3o.ooo ALl 5(} sets: q < IO Au

a

L

1

50

Ioo

15o , (~ y I

zoo

25o

2O00

18.3

tsoo

13.7

IOOO

9.1

5oo

4.6

o

,

~ [ ~ j S0

~

_ 100

~

_ Time

i

150 (Mvr)

~

i

i

]

200

L _ ~ _ ~

250

i

FIG. 3. Same as for Fig. 1 but ibr a0 = 30,000 AU. No showers at either perihelion limit; the loss cone is filled.

SIMULATIONS OF COMETS ao

40,000

=

AU;

50

sets:

q

<

2

ao = 40.000

AU

b

200 i

I

I

I

'

,

k

I

[

'

,

r

,

I

'

,

,

I I ' l [ I

,

113

I

I

I

I

I

1

AU; 50 I

I

sets:

I

q < 1o AU I

I

r

,

I

i

,

,

I

' /'.3

F -

z

0.9

1oo

50

0

1.4

-

~

i

[ 50

i

L

i

i

[ i 100

i

i Time

i

I i 150

i

i

L

i 200

i

i

i

I 250

z

6O0

5.,5

4o0

3.7

0.5

i

1.8

, , q , , 50

,

,I

,, 100

, Time

] , is0 (Myr)

,,

I 200

,

,

L

,

I 25o

L

(~,ly r )

FIG. 4. Same as for Fig. 1 but for a o = 40,000 AU.

far away. A p p a r e n t l y the perturbations are not yet strong enough to give a totally unif o r m c o m e t distribution o v e r the entire range of 1 - e 2. N o n e t h e l e s s , the density is constant o v e r regions that are c o m p a r a b l e in size to the loss cone. Therefore, the loss cone fills uniformly b y a = 30,000 A U , but the total flux is only about 60% of what it would be if the density was truly uniform o v e r all eccentricities. H T A found that for 10,000 -< a -< 20,000 A U , the b a c k g r o u n d flux due to the Galactic tide and stellar perturbations is four times greater than that due to stars alone. T h e s e simulations, which treat 25-50 times m o r e c o m e t s , confirm this result: b a c k g r o u n d rate is larger than the rate for stars alone determined by H T A by a factor of 5.5 at a = 20,000 A U and by 5.2 at a = 10,000 AU. H e n c e , the b a c k g r o u n d flux due to the c o m bination o f the Galactic tide and stellar perturbations is about what H T A had predicted. The tide definitely d o m i n a t e s the evolution o f c o m e t s with s e m i m a j o r axes under about 30,000 AU. B e y o n d 30,000 A U the loss cone is filled and the source of the perturbation is irrelevant to the determination of the flux. b. 4.5 Billion-Year Simulations

This section investigates the long-term evolution of the Oort cloud and the flux f r o m

it. All sources of perturbations, including giant molecular clouds, other than passing stars and the Galactic tide will continue to be neglected. O v e r short peridos of time, as considered in the previous section, the only important sources o f perturbation are stars and the tide. O v e r longer periods, o f o r d e r 1-3 G y r , e n c o u n t e r s with giant molecular clouds are important and m a y strip a w a y c o m e t s with a ~> 25,000 A U (e.g., H u t and T r e m a i n e 1985, Weinberg et al. 1986) but are not likely to disrupt the majority o f the cloud which is at distances ~< 10,000 AU. A p r o p e r long-term treatment of the cloud must include giant molecular clouds, since e v e n the largest estimate o f the half-life is less than the age of the Solar S y s t e m . H o w ever, this complication has not been added to the model b e c a u s e the p a r a m e t e r s describing the clouds are uncertain and because one can still learn m u c h a b o u t the evolution of the Oort cloud without them. N o n e t h e l e s s , it should be u n d e r s t o o d that the simulations described below are not intended to be a true reflection of the longt e r m orbital evolution of the cloud, since G M C s are a vital but missing c o m p o n e n t . Seven sets h a v e been c o m p l e t e d at a = 20,000 A U , equivalent to an initial population o f 1,530,200 c o m e t s or 9800 tokens. Figure 5 shows the c o m e t flux for q < l0 AU. The b a c k r o u n d flux d e c r e a s e s o v e r time.

! 14

JULIA HEISLER ,0

~)O~)

.

1

r

[

AL.

200oo T

r

-

:

T

5'

sel~ r

q .

10 .

.

.

A/ [

i

400

:26.1

19.6

13.1

6.5

FiG. 5. 4.5 Gyr simulation of a cloud with initially 1,530,200 comets (or 9800 tokens). New comets with q < 10 AU are shown.

During the first 500 myr, the average a b s o l u t e background flux is !15 _+ 41 comets/ myr, while during the last 500 myr is only 39 _+ 18 c o m e t s / m y r . (The frequency of showers at least five times more intense than the current background remains constant and is zero.) This drop by a factor of nearly three is mostly caused by attrition in the comet population, from an initial population of i.5 million to a final one of about 811,000 comets, In addition to a population decrease there is also significant diffusion in energy (Fig. 6) due to stellar perturbations. (The Galactic tide has no effect upon energy integrated over a complete orbit and if radial gradients in the tide are ignored.) Shown in Fig. 6 is the average energy distribution in the cloud every billion years. (The last panel in Fig. 6 is the average energy distribution during the last 500 myr). All of the comets initially have 1/a = 50 × 10 6 A U J. H o w e v e r , the energy distribution rapidly broadens and shifts to lower I / a . (The diffusion rate to energies I / a <~ 20 is underestimated for reasons discussed below.) During the second billion years the peak of the distribution has moved to l / a ~ 45 and 50% of the comets have 32 <~ l / a ~ 52. By the end o f three billion years the peak has shifted to under 40. After 4 billion years, the energy distribution has

flattened and is nearly uniform from 15 ~< I / a <~ 40. Thirty-seven percent of comets have a ~> 35,000 AU, 22% have 25,000 ~< a ~< 35,000 AU, 32% have 15,000 ~< a ~ 25,000 AU~ and 9% have a ~< 15,000 AU. For 811,000 comets with this energy distribution, one would expect an absolute background flux of - 5 0 c o m e t s / m y r (or Rb, ~ = 6.2 c o m e t s / y e a r for a 10 ~ comet cloud), based on the results from the previous section; this is somewhat higher than the final background measured here of 39 _+ 18 comets/myr (or R b ± o " I = 4.8 -+ 2.2 comets/ year). The remaining discrepancy between the measured R b = 4.8 and the expected Rb, c = 6.2 is small and is probably not significant in light of the fact that the background flux is highly variable. Both population attrition and energy diffusion play a role in determining the final background flux; however, population attrition is the dominant effect. This diffusion in orbital energy caused solely by stellar perturbations, not the Galactic tide, has been discussed by others. Weissman (1982, 1985b) and Fern~indez (1982) found 7 and 10%, respectively, of an initial comet population between Neptune and Uranus diffuses to distances beyond 0.5 pc after 4.5 Gyr. The percentages were about 20 and 10% tot a cloud formed in s i t u . Remy and Mignard (1985) also noted the effect, and tbund that 43-80% of comets diffuse to such large distances. This simulation finds that 21% of comets with a 0 = 20,000 AU diffuse to a > 0.6 pc or reach distances greater than 1 pc and 5% are directly ejected by stars. These percentages reflect a total loss of about 695,000 comets, of which 396,000 are due to direct stellar ejection or diffusion. (Comets that diffuse to a > 0.6 pc or to distances beyond I pc are removed from the simulation because their orbits b e c o m e too difficult to compute. If such low energy comets had been retained they would have likely been directly ejected within a short amount of time.) Of the comets that are not removed due to falling into the loss cone, about 33%

SIMULATIONS OF COMETS -

'

'

*

I

'

~

'

I

'

'

,'"''"]'"'""r

,

,

I

'

115 '

'

I

'

'

2_

t

-

1.5xI0B ixlOs 5xlO4 0 20

40

80

60

120

I00

IxlO 5

1

_

8x104 6x10 a 4x104 2x104 0 20

40

60

BO

40

60

O0

I00

80

I00

8o

loo

100

IZO

lx105

5xlO 4

I

0

20

0 Lzxl°S

'

'

'

T

l '

'

~

I

~

'

1 Ill ~

~

'

12.0

I

BxlO4~-

4x104

2×1o4~_ 0

L_L 20

40

60

40000

20000 0

~ -~--~J

zo

~

40

1/~

~o~ ~ ( 1 o -o AU" )

leo

FIG. 6. Evolution in orbital energy in the cloud. E a c h panel s h o w s the average distribution in 1/a in the cloud during the first, s e c o n d , third, a n d fourth billion years, a n d the fifth panel is the average distribution during the final 500 myr, |nitially all c o m e t s had l/a = 40 × 10 -6 A U - L With time t h e distribution broadens, and by the final 500 m y r , 37% o f all surviving c o m e t s have lla < 29 x l0 ~ A U -~ (a = 35,000 AU).

116

JULIA HEISLER

are lost due to stellar ejection or diffusion, implying a half life for comets for these two processes alone of 7.8 Gyr. This loss in population can be compared with estimates made by Weinberg et al. (1986) of the survival probability of wide binaries subject to stellar perturbations. T h e y found that the half-life of binary stars separated by about 20,000 AU is about 3.1 Gyr, more than a factor of two less than that found here. H o w e v e r , it is not possible to compare the two results in detail because Weinberg et al. models a somewhat different situation, wide binaries without a loss cone, and with an entirely different method utilizing diffusion coefficients from a Fokk e r - P l a n c k treatment where results may be sensitive to certain parameter choices. Another comparison can be made to the semianalytic estimates of the half life made by Bahcall et al. (1985), also for the wide binary problem. They found a half-life due to disruption from visible stars alone of tln=/32y3.81 x I 0 3 G y r ( 0 - ~ ) + m2] (0.03 MQ2 pc-3 ) mlMO / \ P . . . . Mseen /

(3)

V. (50 k ~ s e c ) ' where P~eenMseen =

f0~085Mo n ( M ) M d M ,

n ( M ) is the n u m b e r of visible stars per cubic pc with mass M, roland m2 are the masses of the binary components, V, is the rms relative velocity of stars with respect to the binary, and finally/3 = 0.023 -+ 0.005 and 3' = 2.2 _+ 0.5 are constants that were numerically determined. In these simulations V, -~ 40 km/sec and p M 0.029 Mo 2 pc -3 due to omitting stars with masses less than 0.2 M o. Thus, tl/2 = 3.5 Gyr, but given the uncertainties in/3 and 3", tv2 could be as low as 2.1 G y r or as large as 6.5 Gyr. The half-life implied by these simulations, tl/2 = 7.8 Gyr, is higher than Eq. (3) predicts, even given the uncertainties in/3 and 3". Most ....

....

=

likely the reason why the diffusion and ejection rates here are low is that only stellar encounters within 1 pc are allowed. The 1 pc limit is adequate for comets with a ~< 40,000 AU, or for short time scale simulations where diffusion is unimportant. However, once comets reach semimajor axes ~> 50,000 AU the combined effect of weak encounters up to 2 or 3 pc from the Sun should probably be included to obtain a more accurate estimate of the diffusion and ejection rate for large a comets. By 4.5 Gyr, about 20% of the surviving comets have a --- 50,000 AU (Fig. 6). Hence, the loss rate due to ejection and diffusion found here is likely just a lower limit. H o w e v e r , the rate that comets enter the loss cone, as shown in Fig. 5, is unaffected by whether extremely weak encounters are included or not, since this rate is almost entirely determined by the number of small a comets (a ~< 30,000 AU). Nonetheless, at lower semimajor axes (i.e., - 1 0 , 0 0 0 AU) energy diffusion is much less important over the age of the Solar System. Since the distribution of comets in the Oort cloud is likely a strongly decreasing function of distance with - 7 5 % of the cloud inside 15,000 AU, diffusion is probably not a strong effect for the cloud as a whole. H o w e v e r , diffusion from lower semimajor axes may allow the outer region (~>20,000 AU) of the cloud to replenish itself following destructive GMC encounters if and when such encounters occur. c. Combining Semimajor Axes An interesting use of the simulated data is to combine semimajor axes to model the overall flux from the Oort cloud. This exercise yields the best theoretical description of the actual flux into the inner Solar System. The results of Duncan et al. (1987; hereafter DQT)are used to combine semimajor axes. They simulated the formation and subsequent evolution of the comet cloud from a source in the outer planetary system and found that from about 3000 to 50,000 AU the steady-state density in the cloud is a power law proportional to r - 3 5 . A representation

SIMULATIONS OF COMETS 10.000

z

< a < ,10.ooo AU; q < 2

AU

3.5

400

2.8

3OO

2.1

~

1.4

1oo

0.7

o 50

100

150 T i m e (Myr)

2O0

250

FIG. 7. Total n e w c o m e t flux inside q = 2 A U . C o m ets from different s e m i m a j o r axes were c o m b i n e d according to the distribution predicted by DQT.

o f this distribution in r is achieved by weighting these results at different semimaj o r axes as follows (taken f r o m Fig. 8 of DQT): f o r a = 10,000AU,

2, a > 20,000) -+ o-t 0.91 -+ 0.22 c o m e t s / years f r o m a l0 II (a > 20,000 A U ) cloud. C o m p a r e d to the o b s e r v e d new c o m e t flux of 2.1 c o m e t s / y e a r s for a < 1 A U and absolute magnitude H10 brighter than 11 (see W e i s s m a n 1989b, Bailey and Stagg 1988), this simulated flux implies an actual Oort cloud population of 5 × l0 H c o m e t s with a > 20,000 AU. Figure 7 is the best representation available of the new c o m e t a r y flux in the vicinity o f the Earth and its temporal variation. N e x t the energy distribution of new c o m ets in the b a c k g r o u n d flux was c o m p u t e d . C o m e t s with q < 2 A U that entered the Solar S y s t e m during the last 50 million years shown in Fig. 7 were examined. The c o m e t s were divided into two groups: those that entered during periods of relatively higher background flux (>0.94R b or N > 30 in Fig. 7) and those that entered during periods of lower flux (--<0.94Rb or N < 30), where 0.94R b is the median flux during the last 50 myr. The resulting distributions in 1/a for both groups of c o m e t s are s h o w n in Figs. 8a and 8b. During times of low c o m e t flux (i.e., 50% of the time), new c o m e t energies are very strongly p e a k e d at l / a ~ 35 (in units of 10 - 6 A U - 1), and virtually none have 1/a as large as 60. [A different result though with far f e w e r c o m e t s , was d e m o n s t r a t e d by R e m y and Mignard (1985). T h e y found that for c o m e t s with perihelia initially b e t w e e n N e p t u n e and Uranus, the i / a distribution after 4.5 G y r was p e a k e d at l / a ~ 20 × l0 6 A U - t . ] During periods of higher background flux, the distribution in l / a is very similar, except there are s o m e w h a t m o r e c o m e t s with 1/a > 40, and some, though relatively few, with l / a - 100. In contrast, during the showers depicted in Fig. 7 the majority (>70%) have l / a = 100. (Note that the discrete sampling in l / a m a k e s the distributions in Figs. 8a and 8b a p p e a r clumpy w h e n they would not otherwise be had l / a been sampled in a continuous fashion.) One-half of the c o m e t s that enter during periods o f high background flux have 29 ~< 1/a <<, 39; during periods o f low flux the :

500

200

117

Wlo = 1

a = 20,000AU,

W2o = 0.12

a = 30,000AU,

W3o = 0 . 1 1

a = 40,000 A U ,

w40 = 0.09.

(4)

The total flux of c o m e t s / m y r into the inner Solar S y s t e m at time t is then given by Ntot(t) = Njo(t)Wlo + N2o(t)W2o + Nao(t)W3o + N4o(l)W4o , where Ni0, N20, N30, and N4o are

the c o m e t fluxes at a = 10,000, 20,000, 30,000, and 40,000 A U , respectively. Unfortunately, this realization o f the D Q T distribution is crude since the simulated data begins at 10,000 AU and extends only to 40,000 AU. The total c o m e t flux within 2 A U o f the Sun was c o m p u t e d by applying the a b o v e weighting s c h e m e to the results f r o m Section I l i a and is shown in Fig. 7. The background rate Rb( q < 2) ---+o"t = 0.22 ---+0.05 c o m e t s / y e a r f r o m a 10 jl c o m e t cloud with a -> 3000 A U ; 1.2% of the time is spent in w e a k showers, and 0.8% in strong showers. I f one considers only the flux f r o m a -> 20,000 A U , then the b a c k g r o u n d is Rb(q <

118

JULIA HEISLER DQT high f l u x

last 50 m y r

DQ1 Low flux, t~lsl 5O r x y r

h tSO

~

'

'

~

l

"

~

'

,

]

r

T

F

T

T

"

,

r

r

r

T

T

.

.

!: I oo I

t

m)

i, i

!

i

60 ff 60 ~

t 40i

t 40 i

I

t

l,!/!

20

:° ao

4o

t/,

6 06 (10 AU 1)

~(J

I ~)o

l :I()

ill

I

[

o

20

" ~f}

60 I :, ( t Q 6 M r

I

8O

10Q

FIG. 8. Distribution in I / a of new c o m e t s arriving during the last 50 m y r simulated (220 < t < 270 myr). (a) intervals of high b a c k g r o u n d flux and (b) intervals of low b a c k g r o u n d flux. E s s e n t i a l l y no ne w c o m e t s should have l / a > 60 x 10 6 AU ~ outside of showers.

interval is only slightly narrower. This strong peak at 1/a ~ 35 in the simulations is m u c h n a r r o w e r than the peak in the data. H a l f of 39 class 1 and 2 c o m e t s with 0 < I/aorig < 120 and q < 2 A U have 20 < I/aorig < 50 (Marsden 1986). H o w e v e r , this disc r e p a n c y could well be due to observational errors. Errors as large as 30 × 10 6 AU would m a s k differences in the l / a distribution as a function of flux. F u r t h e r m o r e , the existence of a small n u m b e r of mildly hyperbolic orbits, most with - 1 5 0 ~< l /a % 0 but one with l / a = - 4 2 0 0 x 10 6 A U (Yabushita and H a s e g a w a 1989), probably implies errors due to neglecting nongravitational forces, although it has been suggested (e.g., Clube and N a p i e r 1982, Bailey 1983) that these c o m e t s could be evidence of an e n c o u n t e r with a G M C and capture from interstellar space. It is also interesting that a large fraction, 20%, of the data have large l/a: Out of the same 39 class 1 and class 2 comets, 7 have 65 < l / a < 89; in the simulated distributions almost none have 1/a this large. This d i s c r e p a n c y can be entirely explained as due to positive or negative error in the determination of 1/a for o b s e r v e d c o m e t s and contamination by c o m e t s with smaller original

1/a that are not on their first passage through the inner Solar System. Failure to include nongravitational forces is another source of error, making a c o m e t a p p e a r less bound than it truly is. A further complication c o m e s f r o m c o m e t a r y fading of older c o m e t s that are not truly new ( K r e s a k 1977, Everhart 1982, Bailey 1984), making the comparison with the data e v e n m o r e problematic. Nevertheless, the clustering at l / a ~ 75, taken at face value, is interesting because assuming no errors in I/a or contamination by old c o m e t s , it implies that we are in a weak c o m e t shower, an idea discussed also in the past by Oort (1950) and others. (The s h o w e r would have to be w e a k because there are no c o m e t s with 89 < l/a < 120, which would be present in a strong shower.) Thus, to improve o n e ' s understanding of what should be e x p e c t e d at these values of I/a, an additional 50 sets were c o m p u t e d at l / a = 75, the median 1/a of the seven clustered comets. The background flux within 2 AU of the Sun was 1.4 c o m e t s / myr, and if about 9% of the c o m e t cloud is between about 11,500 and 15,000 AU, then only 0.5% of the background c o m e t flux should be at l / a ~ 75. The o b s e r v e d percentage is m o r e than 33 times larger than

SIMULATIONS OF COMETS this. H e n c e , since we are p r o b a b l y not in a w e a k s h o w e r the d i s c r e p a n c y b e t w e e n model and data only u n d e r s c o r e s the observational uncertainties, namely that they are large and cannot be neglected. IV. C O N C L U S I O N S

Monte Carlo simulations of the Oort c o m e t cloud have been c o m p l e t e d at semimajor axes o f 10,000, 20,000, 30,000, and 40,000 AU. At each s e m i m a j o r axis a cloud of o v e r 10 million c o m e t s was simulated and their orbits were evolved o v e r a period of 270 m y r subject to Galactic tide and stellar perturbations. A m u c h longer simulation of 4.5 billion years treated a cloud with initially 1.5 million c o m e t s at a 0 = 20,000 AU. The principal results can be s u m m a r i z e d as follows: (1) The b a c k g r o u n d flux rate within 2 A U of the Sun was (assuming 10 jt c o m e t s in the cloud): - 0 . 0 0 3 c o m e t s / y e a r at 10,000 AU 0.31 c o m e t s / y e a r at 20,000AU

(5)

1.61 c o m e t s / y e a r at 30,000 AU 0.77 c o m e t s / y e a r at 40,000 AU. The temporal variation of this flux is shown in Figs. 1-4. (2) At 10,000 AU the flux varies tremendously o v e r time; the variability factor, S, which is the ratio of the mean s h o w e r flux to the background equals 181 for q < 2 AU. As a increases, the temporal variation decreases quickly. By a = 20,000 AU, S = 0.6 for q < 2 AU. (3) By 30,000 A U the loss cone is already full. B e y o n d this point, no showers are present and we expect the fractional background flux to vary as a-5/2 (HTA). (4) O v e r the course o f the age of the Solar System, the Oort cloud b e y o n d a = 20,000 A U likely reduces in population by m o r e than a factor of two. (About 20% o f c o m e t s with initial a = 20,000 A U are lost by entering the loss cone.) C o m e t s that initially have l / a = 50 (in units o f 10 -6 A U - I ) diffuse in

119

1/a and after 4.5 G y r are roughly uniformly distributed b e t w e e n 10 and 40. (5) Combining these results f r o m different semimajor a x e s according to the distribution in a predicted by D Q T , a representation was m a d e o f the actual new c o m e t a r y flux near the E a r t h (Fig. 7). Figure 7 should be the best such description of this flux produced so far. The fractional b a c k g r o u n d flux is 0.2 c o m e t s / y e a r inside q = 2 AU for a 10 ~1 c o m e t cloud and 2% o f the time is spent in showers (rate > 5 times background). Furthermore, the distribution in l / a for new c o m e t s in b a c k g r o u n d is strongly p e a k e d at l / a ~ 35 (Fig. 8), and 50% of all new c o m e t s have 25,000 ~< a ~< 33,000 A U (30 <~ 1/a <<40). Essentially no new c o m e t s have I / a > 70, unless in a shower. The fact that a significant n u m b e r of c o m e t s have been observed with 65 < 1/a < 90 e m p h a s i z e s either that there must be a large error in the determination of 1/a or a large contamination by older c o m e t s b e c a u s e if not, the d i s c r e p a n c y would imply that we are in a weak c o m e t shower. The possibility that we are in a s h o w e r is an exciting one. Stricter limits should be possible by examining other orbital elements, such as the argument of perihelion and Galactic latitude, and this will be the subject of a future p a p e r (Heisler et al. 1990). ACKNOWLEDGMENTS 1 thank Brian Marsden, Jan Oort, and Eugene Shoemaker for the suggestions and information they provided through various correspondences, and Simon White and Paul Weissman for their helpful comments, suggestions, and discussions. I am particularly grateful to Scott Tremaine for his many recommendations and critical review of this manuscript. This work was supported by NSF Grant AST 8352062. REFERENCES ALVAREZ, L. W., W. ALVAREZ, F. ASARO, AND H. MICHEL 1980. Extraterrestrial cause for the Cretaceous-Tertiary extinction: Experimental results and theoretical interpretation. Science 208, 1095-1108. BAHCALL, J. N. 1984. Self-consistent determinations of the total amount of matter near the Sun. Astrophys. J. 276, 169-181. BAHCALL, J. N., P. HUT, AND S. TREMAINE 1985.

120

JULIA

M a x i m u m m a s s o f objects that constitute u n s e e n disk material. Astrophys. J. 290, 15-20. BAHCALL, J. N., AND R. M. SONEIRA 1980. T h e universe at faint magnitudes. I. Models for the Galaxy and the predicted star counts. Astrophys. J. Suppl. 44, 73-110. BAILEY, M. E. 1977. Some c o m m e n t s on the Oort cloud. Astrophys. Space Sci. 50, 3-22. BAILEY, M. E. 1983. T h e structure and evolution of the Solar S y s t e m c o m e t cloud. Mort. Not. R. Astron. Soc. 204, 603-633. BAILEY, M. E. 1984. The steady-state l/a distribution and the problem of c o m e t a r y fading. Mon. Not. R. Astron. Soc. 211, 347-368. BAILEY, M. E. 1986. T h e m e a n energy transfer rate to c o m e t s in the Oort cloud and implications for c o m e t a r y origins. Mon. Not. R. Astron. Soc. 218, 1-30. BAILEY, M. E., AND C. R. STAGG 1988. Cratering constraints on the inner Oort cloud: Steady-state models. Mort. Not. R. Astron. Soc. 235, 1. CLUSE, S. V. M., AND W. M. NAPIER 1982. Spiral a r m s , c o m e t s , and terrestrial catastrophism. Quart. J. R. Astron. Soc. 23, 45-66. DAVIS, M., HUT, P., AND R. A. MULLER 1984. Extinction o f species by periodic c o m e t showers. Nature (London) 308, 715-717. DUNCAN, M., QUINN, T., AND S. TREMAINE 1987. T h e formation and extent o f the Solar S y s t e m c o m e t cloud. Astron. J. 94, 1330-1338. DUNCAN, M., QUINN, T., AND S. TREMAINE 1988. The origin o f short period c o m e t s . Astrophys. J. 328, L69-L73. EVERHART, E. 1982. Evolution o f long and short period orbits. In Comets (L. L. Wilkening, Ed.), pp. 659-664. Univ. of Arizona Press, T u c s o n . FERNANDEZ, J. A. 1980a. On the existence of a c o m e t belt b e y o n d N e p t u n e . Mon Not. R. Astron. Soc. 192, 481-491. FERNANDEZ, J. A. 1980b. Evolution of c o m e t orbits u n d e r the perturbing influence o f the giant planets and nearby stars. Icarus 42, 406-421. FERNANDEZ, J. A. 1981. N e w and evolved c o m e t s in the Solar S y s t e m . Astron. Astrophys. 96, 26-35. FERNANDEZ, J. A. 1982. Dynamical a s p e c t s of the origin of c o m e t s . Astron. J. 87, 1318-1332. FERNANDEZ, J. A. AND W . - H . IP 1987. Time d e p e n d e n t injection o f Oort cloud c o m e t s into Earth-crossing orbits. Icarus 71, 46-56. HEISLER, J., AND S. TREMAINE 1986. The influence of the Galactic tidal field on the Oort c o m e t cloud. Icarus 65, 13-26. HEISLER, J., TREMAINE, S., AND C. ALCOCK 1987. The f r e q u e n c y and intensity o f c o m e t s h o w e r s from the Oort cloud. Icarus 70, 269-288. HEISLER, J., WEISSMAN, P., AND S. TREMAINE 1990. Are we in a c o m e t s h o w e r ? In progress. HILLS, J. G. 1981. C o m e t s h o w e r s and the steady-state

HEISLER infall of c o m e t s from the Oort cloud. Astron. J. 86, 1730-1740. H s u , K. J. 1980. Terrestrial c a t a s t r o p h e c a u s e d by c o m e t a r y impact at the end of the Cretaceous. Nature (London) 285, 201-203. HUT, P., AND S. TREMAINE 1985. H a v e interstellar clouds disrupted the Oort comet cloud? Astron. J. 90, 1548-1557. KRESAK, L. 1977. In Comets, Asteroids, Meteorites. 1AU Coll. 39 (A. H. D e l s e m m e , Ed.) University of Toledo Press, Toledo, p. 97. KUIJKEN, K., AND G. GILMORE 1989. The m a s s distribution in the galactic disk- I11. The local volume m a s s density. Mon. Not. R. Astron. Soc. 239, 651-664. MIHALAS, D., AND J. BINNEY 1981. Galactic Astronomy: Structure and Kinematics, p. 423. F r e e m a n , San Francisco. MARSDEN, B. 1986. Catalogue o f Cometary Orbits, 5th ed. Cambridge, MA. MORRIS, D. E. 1985. Comet Showers Are Not Induced by Interstellar Clouds. L a w r e n c e Berkeley Lab Report LBL-20650. MORRIS, D. E. AND R. A. MULLER 1986. Tidal gravitational forces: The infall of " n e w " c o m e t s and c o m e t showers. Icarus 65, 1-12. NAPIER, AND STANIUCHA 1982. Interstellar planetesimals. 1. dissipation o f a primordial cloud of c o m e t s by tidal e n c o u n t e r s with m a s s i v e nebulae. Mon. Not. R. Astron. Soc. 198, 723-735. OORT, J. H. 1950. The structure o f the cloud of c o m e t s surrounding the solar s y s t e m , and a hypothesis concerning its structure. Bull. Astron. Inst. Neth. 11, 91-110. RAMPINO, M. R., AND R. B. STOTHERS 1984. Terrestrial m a s s extinctions, c o m e t a r y impacts and the S u n ' s motion perpendicular to the galactic plane. Nature (London) 308, 709-712. REMY, F., AND F. MIGNARD 1985. Dynamical evolution of the Oort cloud. I. A. Monte Carlo simulation. Icarus 63, 1-19. SCHREUR, B. 1979. Doctoral dissertation, Florida State Univ. SHAPIRO, S. L., AND A. B. MARCHANT 1978. Star clusters containing m a s s i v e , central black holes: Monte Carlo simulations in two-dimensional phase space. Astrophys. J. 225, 603-624. SHOEMAKER, E. M., AND R. F. WOLFE 1984. Evolution of the U r a n u s - N e p t u n e planetisimal swarm. Lunar Planet. Sci. XV, 780-781. SHOEMAKER, E. M., AND R. F. WOLFE 1986. M a s s extinctions, crater ages and c o m e t showers. In The Galaxy and the Solar System (M. S. M a t t h e w s , J. Bahcall, and R. S m o l u c h o w s k i , Eds.), pp. 338-386. Univ. of Arizona Press, T u c s o n . SMOLUCHOWSKI, R., AND M. TORBETT 1984. The b o u n d a r y of the Solar S y s t e m . Nature (London) 311, 38-39. STAGG, C. R. AND M. E. BAILEY 1989. Stochastic

SIMULATIONS OF COMETS capture of short-period comets. Mon. Not. R. Astron. Soc. 241, 507. ToaBE-rr, M. 1986. Dynamical influence of Galactic tides and molecular clouds on the Oort cloud of comets. In The Galaxy and the Solar System (M. S. Matthews, J. Bahcall, and R. Smoluchowski, Eds.), pp. 147-172. Univ. of Arizona Press, Tucson. VAN DEN BERGH, S. 1982, Giant molecular clouds and the solar system comets. J. R. Astron. Soc. 76, 303-308. WASSERMAN, I. 1988. Theoretical models for wide binary evolution. Astrophys. Space Sci. 142, 267276. WEINBERG, M. D., SHAPIRO, S. L., AND I. WASSERMAN 1986. The dynamical fate of wide binaries in the solar neighborhood: Encounters with stars and molecular clouds. Icarus 65, 27-36. WEISSMAN,P. 1979. Physical and dynamical evolution of long-period comets. In Dynamics o f the Solar System, pp. 277-282. (R. L. Duncombe, Ed.) Reidel, Dordrecht. WEISSMAN, P. 1982. Dynamical history of the Oort

121

cloud. In Comets, (L. L. Wilkening, Ed.), pp. 637-658. Univ. of Arizona Press, Tucson. WElSSMAN, P. 1985a. Cometary Dynamics. Space Sci. Rev. 41, 299-349. WEISSMAN, P. 1985b. In Dynamics o f Comets (A. Carusi and G. B. Valsecchi, Eds.) Reidel, Dordrecht. WEISSMAN, P. 1988. The impact history of the solar system: Implications for the origin of atmospheres. In Origin and Evolution o f Planetary and Satellite Atmospheres. WE|SSMAN, P. 1989a. The Oort cloud. Preprint. WEISSMAN, P. 1989b. The cometary impactor flux at the Earth. Preprint. WHITMIRE, D. P., AND A. A. JACKSON, JR. 1984. Are periodic mass extinctions driven by a distant solar companion? Nature (London) 308, 713-715. WHITMIRE, D. P., AND J. J. MATESE 1985. Periodic comet showers and planet X. Nature (London) 313, 36-38. YABUSHITA, S., AND I. HASEGAWA 1989. Original and future orbits of ten hyperbolic comets. Observatory 109, 189.