Embedded star clusters and the formation of the Oort cloud

Icarus 191 (2007) 413–433 www.elsevier.com/locate/icarus

Embedded star clusters and the formation of the Oort cloud II. The effect of the primordial solar nebula R. Brasser a,b,∗ , M.J. Duncan a , H.F. Levison c a Department Physics and Astronomy, Stirling Hall, Queen’s University, Kingston, K7L 3N6, Ontario, Canada b Canadian Institute for Theoretical Astrophysics, 60 St. George Street, Toronto, M5S 3H8, Ontario, Canada c Southwest Research Institute, 1050 Walnut Street, Boulder, 80302 CO, USA

Received 8 March 2007; revised 2 May 2007 Available online 29 May 2007

Abstract This paper deals with Oort cloud formation while the Sun was in an embedded cluster and surrounded by its primordial nebula. This work is a continuation of Brasser et al. [Brasser, R., Duncan, M., Levison, H., 2006. Icarus 184, 59–82], building on the model presented therein, and adding the aerodynamic drag and gravitational potential of the primordial solar nebula. Results are presented of numerical simulations of comets subject to the gravitational influence of the Sun, Jupiter, Saturn, star cluster and primordial solar nebula; some of the simulations included the gravitational influence of Uranus and Neptune as well. The primordial solar nebula was approximated by the minimum-mass Hayashi model [Hayashi, C., Nakozawa, K., Nakagawa, Y., 1985. In: Black, D.C., Matthews, M.S. (Eds.). Protostars and Planets II. Univ. of Arizona Press, Tucson, AZ] whose inner and outer radii have been truncated at various distances from the Sun. A comet size of 1.7 km was used for most of our simulations. In all of our simulations, the density of the primordial solar nebula decayed exponentially with an e-folding time of 2 Myr. It turns out that when the primordial solar nebula extends much beyond Saturn or Neptune, virtually no material will end up in the Oort cloud (OC) during this phase. Instead, the majority of the material will be on circular orbits inside of Jupiter if the inner edge of the disk is well inside Jupiter’s orbit. If the disk’s inner edge is beyond Jupiter’s orbit, most comets end up on orbits in exterior mean-motion resonances with Saturn when Uranus and Neptune are not present. In those cases where the outer edge of the disk is close to Saturn or Neptune, the fraction of material that ends up in the subsequently formed OC is much less than that found in Brasser et al. [Brasser, R., Duncan, M., Levison, H., 2006. Icarus 184, 59–82] for the same cluster densities. This implies that for comets of roughly 2 km in size, the presence of the primordial solar nebula hinders OC formation. A byproduct of some of our simulations are endresults with a substantial fraction of the comets in the Uranus–Neptune scattered disk. A subsequent followup of this material is planned for the near future. In order to determine the effect of the size of the comets on OC formation efficiency, a set of runs with the same initial conditions but different cometary radii have been performed as well, from which it is determined that the threshold comet size to begin producing significant Oort clouds is roughly 20 km. This implies that the presence of the primordial solar nebula acts as a size-sorting mechanism, with large bodies unaffected by the gas drag and ending up in the OC while small bodies remain trapped in the planetary region, in the models studied. © 2007 Elsevier Inc. All rights reserved. Keywords: Planetesimals; Comets, dynamics; Comets, origin; Origin, Solar System

1. Introduction The reader is referred to Brasser et al. (2006), for a summary of recent studies in the formation and dynamics of the * Corresponding author.

E-mail addresses: [email protected] (R. Brasser), [email protected] (M.J. Duncan), [email protected] (H.F. Levison). 0019-1035/$ – see front matter © 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2007.05.003

Oort cloud (Oort, 1950). In Brasser et al. (2006), simulations were performed of the scattering of comets by Jupiter and Saturn while the Sun was still in its putative birth cluster. The Oort cloud (OC) was formed in such simulations because, at sufficiently large distances, the stars and cluster gas from the Sun’s birth cluster torqued the pericenter of the comet’s orbit away from these planets, thereby saving the comets from ejection. They examined cluster densities ranging from 100 to 106 M pc−3 and found that typical efficiencies at populating

414

R. Brasser et al. / Icarus 191 (2007) 413–433

the OC ranged from 1–15%, with a mean efficiency rate of 8% after 3 million years of simulations. By that time the Sun was assumed to have left the cluster and the computations were terminated. The simulations in Brasser et al. (2006) were able to reproduce the orbits of the trans-neptunian bodies (90377) Sedna and 2000 CR105 if the mean density of stars encountered by the Sun in the cluster was of order 104 M pc−3 or larger. The model in Brasser et al. (2006) used a number of simplifying assumptions. First, the comets were started on cold orbits in the region of fully formed Jupiter and Saturn. It is well known that such comets are short-lived and that their eccentricities increase within a few tens of thousands of years (Gladman and Duncan, 1990). Second, there was no gas from the primordial solar nebula present in the computations. Third, the computations were terminated after 3 Myr since it was decided that by that time or earlier the Sun had left the cluster, given that embedded star clusters do not last much longer than 5 Myr (Lada and Lada, 2003). However, the computations could have run for longer than that, and it is unknown whether or not the Sun remained as part of a subsequently formed open cluster or what precisely happened after its escape. In this paper we investigate the role of the gas from the primordial solar nebula on the formation of the OC. The other assumptions listed above shall be tackled in forthcoming publications. The first author to study the effect of a rotating gas nebula on the elliptic orbits of planets was performed by Kiang (1962), who investigated the two-dimensional case for both the uniform and Keplerian rotation of the gas disk assuming a force law of the form F = kvrl s m , where vr is the relative velocity between the gas and the planet and s 2 = x 2 + y 2 . Kiang (1962) obtained expressions for a˙ and e˙ for the cases (l, m) = (1, 2), (3, 0), (1, −2) and (3, −2). Unfortunately, the case l = 2, which is the most important one in planetary formation, was not studied. This led Adachi et al. (1976) to investigate the three-dimensional problem with a disk in circular motion with a rotational velocity different from the Kepler velocity due to pressure in the gas, using a law of resistance with l = 2. They obtained formulae for the orbit-averaged variations of a, e and i for small e and i incorporating a parameter η, which is a measure of the difference in angular velocity between the gas and the Kepler velocity. Their Eq. (4.21) allow the reader to compute characteristic timescales for the damping in each of these quantities in the quasi-circular case. They do not investigate the case where the eccentricity is moderate or the orbit is nearly parabolic, which is of interest to us. Around the same time, Weidenschilling (1977) used a similar force law to compute the decay times of orbits and lifetime vs spiraling into the Sun. It should be noted that none of the above studies included the effect of a planet stirring up any planetesimals. Ida and Makino (1993) performed a semi-analytical study on the stirring of planetesimals in the vicinity of a protoplanet. Their analytical formulae were checked with numerical experiments. They calculated that the protoplanet clears out a region of 4 Hill radii on either side of it. The same results are found by Tanaka and Ida (1996) who, like Ida and Makino (1993) apply the re-

sults to protoplanetary growth and compute the characteristic damping timescales of the eccentricities. A first attempt to investigate the possibility of forming a comet cloud around the Sun by scattering comets off a protoplanet in the presence of gas was performed by Higuchi et al. (2002). However, the comets are not followed for many orbits. Instead, they investigate how efficient a comet is scattered by a protoplanet as a function of the orbital spacing in units of Hill radii, and argue that repeated close encounters will ultimately send the comet to the Oort cloud. However, at least in the case without gas drag Duncan et al. (1987) have shown that it is not the close approaches but the repeated small perturbations by the planets that diffuse most comets to the OC. This paper extends the work of Brasser et al. (2006) in investigating the evolution of comets in the Jupiter–Saturn zone. Here we include the effects of gas drag and the gravitational potential of the solar nebula. Section 2 deals with the analytics of aerodynamic drag acting on the comets. Section 3 discusses the modeling of the gravitational potential of said disk. In Section 4 some timescale arguments are presented from which comet evolution can be predicted. In Section 5 the numerical methods employed to simulate the OC formation are discussed. Section 6 contains the results; and the summary and conclusions are presented in Section 7. 2. Aerodynamic drag The concept of aerodynamic drag has first been studied in detail by Adachi et al. (1976). Therefore, our results are based on their theory, which is quickly summarized below. It is assumed that the gas disk is rotating in circular motion around the Sun and that the disk is not globally turbulent. There generally is a radial pressure gradient in the nebula and accordingly the angular velocity of the gas is slower than the Keplerian velocity. The balance of forces requires that 2 sνg2 = sνK +

1 dPg , ρg ds

(1)

where νg is the angular velocity of the gas, νK is the Keplerian angular velocity, ρg and Pg are the density and pressure of the gas in the plane and s 2 = x 2 + y 2 . If the z-dependence of the gas is temporarily ignored, then one can parametrize the radial dependence of the density and the temperature as power laws −α  s , ρg (s) = ρ0g 1 AU −β  s T (s) = T0 (2) , 1 AU where T0 is the temperature and ρg0 is the density at 1 AU. Then one can write for Eq. (1)  νg (s) = νK (s) 1 − 2η(s), (3) where η(s) =

 2 1 cs , (α + β) 2γ vK

(4)

Star clusters and comet cloud II

where vK is the Kepler velocity, cs is the sound speed and γ is the ratio of the specific heat capacities of the gas. Throughout this paper the minimum-mass Hayashi disk (Hayashi et al., 1985) will be used as a reference model. For this disk α = 11/4 and β = 1/2 so that −11/4  s 2 2 e−z /zs , ρg (s, z) = ρ0g (5) 1 AU where s2 = x2 + y2,  5/4 s AU zs = 0.047 1 AU

(6) (7)

and ρ0g is the gas density at 1 AU. For the minimum-mass model, this is ρ0g = 1.4 × 10−6 kg m−3 (Hayashi et al., 1985). In addition −1/2  s K T = 280 (8) 1 AU and η(s) = 2.13 × 10

−3



s 1 AU

1/2 (9)

.

In general the magnitude of the drag force on a spherical object is expressed as 1 FD = CD (R)πrc2 ρg vr2 ≡ Ag ρg vr2 , (10) 2 where rc is the comet radius, vr is the relative velocity between the comet and the gas and CD (R) is the drag coefficient, which for a spherical object is a function of the Reynolds number, R, only. Due to the dependence of CD on R, the drag force can sometimes take the form FD ∝ vr instead of FD ∝ vr2 (Landau and Lifschitz, 1987). We shall look at this in more detail in the next subsection, after stating the drag force per unit mass, as being FD =

3 CD ρ g 2 v . 8 rc ρc r

(11)

415

First we shall estimate the Mach number. The sound speed in an ideal gas of temperature T is given by  −β/2  γ RT s 3 ≈ 1.43 × 10 m/s, cs = (12) mH 1 AU where R is the gas constant, γ is the ratio of specific heats and mH is the mass of a mole of hydrogen molecules, and Eq. (2) was used and a temperature of 280 K at 1 AU (Hayashi et al., 1985). For small eccentricity the relative velocity vr ∼ evK and for planar comets −1/2  s vK ≈ 3.0 × 104 (13) m/s, 1 AU so that   β−1 2 vr s M= ≈ 21e . cs 1 AU

(14)

Note that M > 1 provided e  0.1 and that for the Hayashi disk the Mach number decreases at fixed eccentricity as a function of planar distance from the Sun. In the case M  1, CD = 2, regardless of K and R (Landau and Lifschitz, 1987). When M ∼ 1, the situation is more complicated since CD is a function of M. For R  103 , CD ∼ 2 when M > 1 and CD ∼ 0.5 when M < 1. For smaller R, CD ∝ R−1 at constant M. Therefore, it is prudent to estimate R so that we know what regime we are in when M ∼ 1. The Reynolds number for a spherical object of radius equal to the comet’s radius is given by R=

2ρg rc vr , χ

(15)

where χ is the gas viscosity. From Adachi et al. (1976) one has R=6

cs M , cm K

(16)

where  8RT πmH

The above equation shall be used in the next subsection.

cm =

2.1. Reynolds number, Mach number and Knudsen number

is the molecular speed and we used, in accordance with Adachi et al. (1976), χ = 13 ρg cm rc K. Hence  γπ M M R=6 (18) ≈ 4.44 , 8 K K

Adachi et al. (1976) discuss in some detail how the drag coefficient depends on three numbers, all of which are properties of the motion of the comet through the gas: the Reynolds number, R, which is a measure of the turbulence of the gas in the wake of the comet, the Mach number, M and the Knudsen number, K. The Mach number is the ratio of the relative velocity of the comet with respect to the gas to the sound speed in the gas, i.e., M = vr /cs . The Knudsen number is the ratio of the mean free path of the gas molecules, λ, to the size of the comet, i.e., K = λ/rc . Adachi et al. (1976) categorizer the force law according to the value of M, and further subdivide it according to the value of K. However, we can estimate what category the comet’s motion resides in by estimating the values of M, K and R from the properties of the gas and the comet.

(17)

where the ratio of specific heats γ = 7/5 was used for molecular hydrogen. The Knudsen number is given by K=

2mp 1.67 × 10−8 = , ρg rc πd 2 ρg rc

(19)

where mp is the mass of a proton and πd 2 is the cross-section of the hydrogen molecule for deflection due to an interaction with another molecule. Using the reference Hayashi disk and taking the size of a typical comet to be ∼2 km (Francis, 2005), we get, staying in the plane of the gas and for a nearly-circular

416

R. Brasser et al. / Icarus 191 (2007) 413–433

orbit, 

s R ≈ 1.57 × 10 e 1 AU 7

−3 

 rc . 2 km

(20)

The Reynolds number R > 102 provided e  0.05 and the Reynolds number is also larger for larger objects. In the regime M  1 and R  103 , CD = 2 (Adachi et al., 1976). When moving out of the plane, the density of the gas decreases—exponentially for the Hayshi disk—and thus the Knudsen number increases and the Reynolds number decreases. Once 1  R  102 the value of CD changes from about 0.5 to 2 as long as M  1. Furthermore, once the Knudsen number K  1, the drag coefficient is only a function of the Mach number—in fact, CD ∝ M−1 —so that the drag force FD ∝ vr (Phillips, 1975). For the Hayashi disk one can compute where K = 1 occurs in the plane, which yields a solution s = sc . Subsequently computing where K = 1 out of the plane and solving for z, one obtains    11 sc ln AU. zc = zs (21) 4 s For a comet size of about ∼2 km and for a minimum-mass density at 1 AU ρg0 = 1.4 × 10−9 g cm−3 , sc = 75 AU. Therefore, once z > zc then (Phillips, 1975) CD =

16 8 cm = KR 3 vr

(22)

because of slip acting on the comet and so 4π 2 r ρ g cm v r . (23) 3 c For orbits with very small eccentricity and inclination such √ that e < η and i < η, M = η = 0.0462s 1/4 for the reference Hayashi disk. This yields a Reynolds number 101  R  103 so that 0.5  CD  5. In summary, we decided to use the following formulation to compute CD when K < 1: FD =

• when M  1, CD = 2 independent of R; • when M < 1 and R > 103 then CD = 0.44 for small values of M and CD = 2 when M  1. We decided to use a quadratic profile in M so that CD = 0.44 + 1.56M2 ; • when R < 103 and for small Mach numbers CD = 24(1 + 0.15R0.687 )/R (Clift et al., 1978). Using a similar quadratic formulation in M yields CD = 2M2 + (1 − M2 )[24(1 + 0.15R0.687 )/R]; and when K  1 then ρg = 0. 2.2. Quasi-circular orbits Using the drag law in Eq. (10), when e, i and η are small, Adachi et al. (1976) show that  

a ˙ 5 + 4α 2 1 2 = −2Hge ζ (e, i, η) η + e + sin i , a 16 8

e ˙ = −2Hge ζ (e, i, η), e

˙

i = −Hge ζ (e, i, η), (24) sin i where  3 ρg (a) CD GM Hge = 8 ρc rc a     ρg0 1 km 1 g cm−3 = 0.0785CD rc ρc 1.4 × 10−9 g cm−3  α+1/2 1 AU × yr−1 , a  5e2 sin2 i vr ζ (e, i, η) = (25) , + + η2 = 8 2 vK and ρg (a) = ρg0 (a/1 AU)−α , where they use the semi-major axis a as a mean value of the position. The quantity Hge has a unit of frequency so that its reciprocal can be thought of as a fundamental timescale of the damping of the eccentricity. According to Adachi et al. (1976), Eqs. (24) give good results when the quantities e, sin i and η are of comparable magnitude. The reader should keep in mind that to derive the above equations, the orbit-averaging principle has been used and therefore should be interpreted with some caution. Averaging the differential equations describing the change in orbital elements over the orbit is fine as long as the change in elements over one orbit is small compared to the value of the elements themselves. Therefore, if an object has an orbit for which e/e (or i/i) is of order unity, the orbit-averaging principle cannot be used. The quasi-circular case allows one to examine some relevant timescales in the problem, which are discussed next. 2.2.1. Decay times Equations (24) allow one to define a decay time as Tx = x/ x ˙ where x is one of the elements a, e or i. The nature of Eqs. (24) does not allow easy separation of variables, so the timescales listed below are derived under special circumstances. In the equations below, the unit of time is years, with Hge given by Eq. (25) above. When i = 0 and η = 0 then  2 1 0.632 −3 = e , Ta = 5 Hge e3 Hge Te = 2e2 Ta ,

(26)

when e = 0 and η = 0 then √ csc3 i 5.657 csc3 i, Ta = 4 2 e = Hg Hge Ti = 2 sin2 i Ta .

(27)

Last, if e η and i η then Ta =

1 0.500 −2 = η , Hge 2Hge η2

Te = 2ηTa , Ti = 4ηTa .

(28)

Star clusters and comet cloud II

From the above equations the conclusion is that when e, i  η, the decay in semi-major axis is fastest due to orbital eccentricity, while when e, i ∼ O(η) the decay is mostly governed by η; the magnitude of the latter is typically O(10−3 ). In Fig. 1 the evolution of two test particles of radius ∼2 km in the minimum-mass Hayashi disk with the effect of gas drag is plotted. One particle has i = 10◦ and e = 0 while the other has e = 0.2 and i = 0. The top-left panel shows the evolution of a for the eccentric case—solid line—and the inclined case— dashed line. The eccentricity and the inclination evolution of both are plotted in the top-right and bottom-left panels. The bottom-right panel shows the evolution of q. Since the comet moves faster than the gas when it is at pericenter, the drag decreases the apocenter distance. At apocenter, the comet moves slower than the gas so that the drag increases the pericenter distance. The orbit will circularize at a distance between the two extremes, which depends on the power-law of the gas, since the latter determines how well angular momentum is conserved. This might be used as a mechanism to lift the pericenters of comets away from Jupiter and therefore prevent them from being ejected by this planet. From Fig. 1, one can see that q decreases in the circular, inclined case, while q increases in the eccentric, planar case. For orbits with small eccentricity and inclination, one can estimate which effect dominates. Using Eqs. (24) one can compute that  1

q ˙ = −2aHge ζ (e, i, η) −e + η + e2 + sin2 i − eη − e3 8  1 − e sin2 i . (29) 8

417

Assuming that both e  η and i  η, we have q ˙ > 0 for √ √

 

2 2 e(1 − e) 1 − e + e2 ≈ 2 2e + O e5/2 . (30) sin i < 1−e √ ˙ > 0. For most orbits sin i < 2 2e so that generally q The evolution of the orbit when it is near-parabolic is of interest given the nature of the problem to be explored in this paper. As such, the parabolic, planar case is examined next. 2.3. Near-parabolic orbits For near-parabolic orbits, all the drag occurs near pericenter, most of it within a distance q  r  2q, and one can use the parabolic case to determine the change in energy and angular momentum per passage through pericenter. While not strictly correct for bound orbits, it should be adequate for orbits with e  0.7. The derivation of the planar parabolic case has been done in Appendix A. The results for the change in angular momentum h and semi-major axis are 3/2   q p Zg B1 (α) h = 2 2GM (31) 1 AU and p

(1/a) = 8Zg B2 (α),

(32) p

where B1 (α) and B2 (α) are given in Appendix A and Zg , also given in Appendix A, is the parabolic equivalent of Hge and, although not strictly correct for the parabolic case, we used (1/a) = −(2/GM ) E. From Eq. (31) above, usually h ∼ 10−3 . Numerical simulations of near-parabolic orbits

Fig. 1. This figure is a summary of the action of gas drag on a comet’s orbital elements. The solid line shows an orbit with initial elements (a, e, i) = (1, 0.2, 0) while the dashed line shows an orbit with (a, e, i) = (1, 0, 11◦ ). The labels above each panel show what element is plotted versus time. Going clockwise from the top-left, these are a, e, q and i.

418

R. Brasser et al. / Icarus 191 (2007) 413–433

using SWIFT (Levison and Duncan, 1994), confirm the predictions of Eqs. (31) and (32), including the facts that (1/a) ∝ q −α and that h ∝ q 3/2−α . Therefore, one can treat each passage through pericenter as being an energy kick of constant magnitude, providing that the eccentricity stays large. The energy kicks for inclined orbits range about 0.8 of the above value for i ∼ 15◦ to about 0.1 of the above value when i ∼ 60◦ . Combined with the random energy kicks the comet receives from the planets with each passage through pericenter, one has the ingredients for a random-walk process in energy, which shall be discussed in more detail in Section 4.2. Instead, we digress for a moment to discuss the potential generated by the gas disk. 3. The disk potential During the course of this investigation, the authors became aware of the effect of the gravitational potential of the primordial solar nebula. Hence its effects are described here. In general solving Poisson’s equation to obtain the potential from a density distribution cannot be done analytically (Binney and Tremaine, 1987). Once one is outside the density distribution one can solve Laplace’s equation using spherical harmonics, but one can only use this as an approximation far away from the source to avoid having to use many terms in the series expansion and compute the corresponding moments. This implies that numerical integration of the density distribution is the only method with which to obtain the potential and the accelerations derived from it. For the Hayashi gas disk, the potential cannot be written in closed form and one has to resort to numerical integration of the density distribution by integrating the potential generated by a ring over all rings. The gravitational potential of a ring residing at (R, Z) in cylindrical coordinates and evaluated at (r, z) is given by (Binney and Tremaine, 1987) 2GMr Φr (R, Z, r, z) = −  π (r + R)2 + (z − Z)2    rR ×K 2 , (r + R)2 + (z − Z)2

(33)

where Mr = 2πRρ(R, Z) dR dz is the mass of the ring, ρ(R, Z) is given by Eq. (5) and K(k) is the elliptic integral of the first kind. The potential for a Hayashi gas disk situated from ri to rd evaluated at (r, z) is then rd ∞ Φh (r, z) = 4π

Φr (R, Z, r, z) dR dZ.

(34)

ri 0

The method employed for computing the above integral is described in Section 5. Far away from the gas disk, i.e., when r  rd or z  zs , the potential can be written in the form   J2 rd2 P2 (cos θ ) GMd 1− Φ(r, z) = − (35) , r r2 where Md is the mass of the disk, rd is the outer edge of the disk, P2 (cos θ ) is the second Legendre polynomial and cos θ =

z/r. The factor J2 rd2 is given by (e.g., Binney and Tremaine, 1987) 



2 ρ r  , z r  2 − 2z 2 r  dr  dφ  dz . J2 rd = (36) For the Hayashi disk, Eq. (36) becomes

 5/2  5/2  2r02 π 3/2 ri rd 2 2 J2 rd = ρ0g r0 z0 − Md 5 r0 r0  3  3  2 z ri rd − 0 − 3 r0 r0   2π 3/2 rd 5/2 ≈ ρ0g r04 z0 , 5Md r0

(37)

where r0 is the distance from the Sun corresponding to the value of ρ0g , ri is the inner edge of the disk, rd is the outer edge of the disk and z0 is the scale height at r0 . The latter approximate equality is only valid when rd  ri . 4. Timescale arguments In Brasser et al. (2006) an argument was made that comets will only be deposited in the OC if the timescale for lifting their pericenters away from the planets by the action of an external perturber in the form of gas and stars of an embedded cluster, tq , is shorter than the typical diffusion time in energy caused by the planets, td . However, with the presence of the gas, there is another timescale involved: the time it takes for the gas to increase the value of the pericenter away from the planet, as well as to decrease the semi-major axis. The first, the increase in q, has to be obtained from numerical experiments—although as shown later, it is closely related to the timescale to decrease the semi-major axis by the gas—while the second, the decay in a, can be computed analytically. Before venturing there, we digress for a moment to state the parameters we use throughout most of this work. 4.1. Parameter space Due to the nature of the problem, we were required to make several choices on what set of parameters to choose for the current study. These are: gas density, extent of the gas disk, comet density, comet size and the inclusion of Uranus and Neptune. Each will be discussed in detail. The drag is determined by the parameter ρ g CD (38) . rc ρc Given that the model of the disk automatically yields CD by the scheme given in Section 2.1 (which might depend on rc for small Reynolds numbers), that leaves three free parameters. Note that insofar as gas drag is concerned, reducing the comet’s radius by some factor is equivalent to increasing ρg by the same amount. The mean cometary density is a much better defined quantity than the mean nuclear radius, so fixing the comet’s density leaves two parameters. The comet’s density is set to ρc = 0.6 g cm−3 , which is based on data from Ag =

Star clusters and comet cloud II

Comet 81P/Wild 2 and 9P/Tempel 1 by Davidsson and Gutierrez (2006) and Davidsson et al. (2007). The comet’s radius might be obtained indirectly via the absolute magnitude, H , but this yields uncertain results. It has been decided to use a radius of 1.7 km, which is in rough agreement with using the data from 81P/Wild 2 and 9P/Tempel 1. The only free parameter left to explore is the gas density, ρg . The density profile of the gas is another parameter that might be tinkered with but instead we chose to use the standard Hayashi model, for which α = 11/4. Substituting the above of ρg , ρc and rc in Eq. (25) gives Hge = 7.70 × 10−2 CD (a/1 AU)−13/4 yr−1 . 4.2. Torquing and diffusion times In Section 2.3 we mentioned that near-parabolic comets will undergo a random walk in energy based on planetary energy kicks and the effect of gas drag near the pericenter of the orbit. We now discuss this in more detail. Following the work of van Woerkom (1948), the energy diffusion of an eccentric comet under the influence of the gas and the planets can be approximated by a Fokker–Planck equation. At the epoch t let ν(1/a, t) d(1/a) be the number of comets passing through pericenter per cometary period and having reciprocal semi-major axes in the range 1/a ∈ [1/a, 1/a + d(1/a)]. The Fokker–Planck equation for the problem is √   GM 3/2 ∂ν(x, t) 1 ∂ 2 ν(x, t) ∂ν(x, t) = x + D2 D1 , ∂t 2π ∂x 2 ∂x 2 (39) where x = 1/a, and ∞ Dn =

un h(u) du,

(40)

−∞

where h(u) is the distribution function in the energy change u (Everhart, 1968). Please see Appendix B for more detail. Formally, h(u) is defined as h(u) = N/ u, with N being the number of comets receiving an energy perturbation u in the range [u, u + u]. The Fokker–Planck equation (39) has been derived by van Woerkom (1948) in the case when D1 = 0, and has been studied and solved by Yabushita (1979, 1980) in the evolution of long-period comets. The case D1 = 0 corresponds to ignoring the far tails of h(u) and that there is no gas. Additionally, the Fokker–Planck approximation is only valid when the energy perturbations (1/a) 1/a and when two successive conjunctions between the comet and the planet are chaotic. The latter effect translates into there being a threshold value of in semi-major axis for which the process holds. This value of semi-major axis is computed in Appendix C. The presence of the gas makes the energy distribution function h(u) asymmetrical about zero, ignoring for now the far tails, which causes D1 = 0. In practice, however, for the minimum-mass Hayashi √ model and a comet size of about 2 km, the coefficient D1 ∼ D2 for i  5◦ and q in the region of the pre-LHB giant planets, i.e., 5 < q < 15 AU. This is obtained by evaluating equation (67) given in Appendix A and comparing it

419

with the rms energy changes caused by the planets in the same range in pericenter distance, q. Hence for near-parabolic orbits the gas plays an important role in the dynamics of the comets. The gas is often unable to bind hyperbolic comets back to the Solar System on their way out when being ejected by Jupiter or Saturn. However, since the term involving D1 is linear rather than quadratic—after n steps one has traveled a distance √ n units from the origin while in the quadratic term this is just n—the linear term generally dominates the motion. According to DQT87, the diffusion time in semi-major axis caused by planetary kicks is td = P (x)x 2 /D2 , where x = 1/a and P (x) is the period of the orbit. The time td is an estimate of the time it takes to change x by order of itself, provided that D2 x 2 . By the same token, one can then use a similar argument to say the corresponding timescale for a change in semi-major axis due to the gas drag is tdg = P (x)x/D1 . From the ratio of these two, one can compute the value of semimajor axis beyond which diffusion caused by planetary energy changes dominates over that of energy changes caused by the gas. Setting td = tdg gives D1 (41) , D2 and since, for our reference disk in the vicinity of the giant plan√ −1/2 ∼ ets for comets about 2 km, D1 ∼ D2 , we have ag ∼ D2 −1 D1 . For larger bodies, especially bodies the size of (90377) Sedna, this condition no longer holds and the motion is entirely dominated by that of planetary scattering. When a < ag the gas drag will win over the typical planetary perturbation because it enters linearly in the Fokker–Planck equation and an orbit with initial pericenter near ap and eccentricity e0 , will decay and circularize outside the planet’s orbit with a ∼ ap (1 + e0 ). However, once a > ag the planetary per−1/2 turbations dominate and the comet can diffuse to ad ∼ D2 , which corresponds to the reciprocal of the rms energy kick given by the distribution h(u), provided ad > ag —which is generally the case—from where it is either ejected or launched to the OC. A plot showing ag and ad as a function of the comet’s pericenter distance, q, is given in Fig. 2 for a minimum-mass Hayashi disk. The solid line with the ‘+’ symbols shows ag for i = 0, the dashed line with the ‘×’ symbols displays ag for i = 30◦ . The value of ad for i = 0 is displayed by the dashed line with ‘∗’ symbols while ad at i = 30◦ is plotted using the square symbols, , and the dotted line. The letters ‘JSUN’ indicate where the planets are on the q-axis. The value of ag was truncated at 10 AU for q close to Jupiter because for a  10 AU with q ∼ 5, the diffusion approximation is invalid, since the orbit is not sufficiently near-parabolic and the encounters are correlated. The value of ag ∝ ρ0g , so that doubling the density will double the value of ag . The reader should be aware that the minima in the values of ag correspond to local maxima of D2 rather than minima of D1 , and are thus caused by the planets and not by the gas. From Fig. 2 one can see that any comet with its pericenter close to Jupiter needs a ∼ 2aJ in order to be able to diffuse outward to a ∼ ad if i is small. For the other three giant planets, ag =

420

R. Brasser et al. / Icarus 191 (2007) 413–433

Fig. 2. A plot showing ag and ad as a function of the comet’s pericenter distance, q for a minimum-mass Hayashi disk. The solid line with the ‘+’ symbols shows ag for i = 0, the dashed line with the ‘×’ symbols displays ag for i = 30◦ . The value of ad for i = 0 is displayed by the dashed line with ‘∗’ symbols while ad at i = 30◦ is plotted using the square symbols, , and the dotted line. The letters JSUN correspond to the positions of the giant planets on the q-axis.

however, ad  2ap so that the comets either first need to be slung to an orbit with a ∼ ad using one close encounter well inside the Hill sphere, before they can diffuse outward without the gas dragging them back in and circularizing their orbits outside the orbit of the planet; or they have to be passed down to Jupiter. The danger with any comet under Jupiter’s control is that ad ∼ 200 AU (see, e.g., DQT87) which is much larger than the energy range for Oort cloud comets in the current Galactic environment, so that the chances of the comet landing in the OC are very small (about 3%; Dones et al., 2004). Fortunately, the fact that the comet may spend some time on a moderately eccentric orbit with q close to Jupiter implies the comet’s pericenter may be lifted away from Jupiter by the gas so that the other three planets—whose kicks are more gentle—can control the diffusion evolution and deposit the comet in the OC. Now that we investigated the role of the gas drag on the energy of the comet for near-parabolic orbits, we can examine the timescale to lift the pericenter away from the influence of a perturbing planet. 4.3. Lifting time of the pericenter The torquing time of the pericenter due to the gas, tqg , is defined here as the time it takes to increase q by a quarter of its original value. This can only be done for fairly eccentric orbits. The 25% increase is chosen because when q = 1.25aJ , where aJ is the semi-major axis of Jupiter, then D2 is a quarter of what it is at Jupiter, so the rms kicks are half as strong. From numerical experiments it turns out that, for planar orbits  11/4  a q tqg ∼ 5500 (42) yr ∼ 2tdg . 5 AU 50 AU

The value of tgq scales linearly with ρg0 and rc and ρc . For inclined orbits, one has two extreme cases: ω = 0 or 180◦ and ω = 90 or 270◦ . We have computed tgq for two values of inclination. These are tqgi ∼ 5.3tqg ,

i = 15◦ , ω = 0,

tqgi ∼ 7.3tqg ,

i = 30◦ , ω = 0,

tqgi ∼ 3.4tqg ,

i = 15◦ , ω = 90◦ ,

tqgi ∼ 6.0tqg ,

i = 30◦ , ω = 90◦ .

(43)

Using the above equations and the value of td computed numerically at several values of q, one can create figures for each of the planets where the values of tqg and td cross as a function of a. However, the diffusion limit is only valid when a < ad , i.e., when td > P (a). In practice, since for comets about 2 km 1/2 and the minimum-mass Hayashi nebula D1 ∼ D2 the semimajor axis where tgqi = td is approximately aqgi ∼

1 √ , 10 D2

(44)

for 15◦ < i < 30◦ . The value of D2 that should be used in the above equation is not the value for the planar case. The above equation implies that aqgi < ad , so that any comet with a > aqgi can diffuse to a = ad . Fig. 3 displays the values of td , tqg and the period P (a) as a function of semi-major axis when the value of q is in the vicinity of one of the major planets for planar orbits. Fig. 4 does the same for obits with i = 30◦ . Once again the minimum-mass Hayashi model and a comet size of 2 km was used. The labels reflect what each line represents. The comet will only be able

Star clusters and comet cloud II

421

Fig. 3. The plots show the values of P (a), td and tqg as a function of semi-major axis, for parabolic orbits with i = 0 and q locked to that of a giant planet, whose name is listed above each panel. The value of a where td and P (a) cross is the diffusion limit ad while the value where td and tqg cross corresponds to the value ag .

Fig. 4. Same as Fig. 3 except that here the orbits have i = 30◦ .

to diffuse to the OC when td < tqg . The crossing of td and tqg yields ag and the crossing of td and P (a) yields ad . For Jupiter the condition td < tqg is satisfied when i ∼ 0 so that those comets stay locked to Jupiter and may (quickly) dif-

fuse to a = ad and then either be ejected or launched to the OC, with ejection the more likely outcome. However, when i = 30◦ , the comets first need a very close encounter to be scattered to a ∼ 100 AU—meaning that the comets need an encounter

422

R. Brasser et al. / Icarus 191 (2007) 413–433

yielding (1/a) ∼ 1/aJ —before they can be fired to the OC. In the case of Saturn, it is only the low-inclination comets that can get out to the OC after having experienced a strong scattering event with (1/a) ∼ 1/aS , but objects with higher inclination will spiral in after receiving a large kick and settle outside the planet’s orbit because for inclinations 10◦ the kicks from Saturn are weaker than the effect of the gas (see Fig. 2). The reason Jupiter can send comets to the OC in the inclined case and for Saturn it is more difficult is because the cross-section for scattering a comet with (1/a) ∼ 1/ap is much larger for Jupiter than it is for Saturn. If Uranus and Neptune are included in the simulations, the probability of reaching the OC is even worse because even low-inclination comets cannot reach the OC and will settle on a circular orbit outside Neptune. Only if one of the latter two planets manages to launch the comet to an orbit with a  1000 AU in one shot—an unlikely event—can the comet diffuse to large values of a. However, the cross-sections for such strong scattering events by these planets are so small that the timescale for significant evolution of the comet’s orbital elements under the action of the gas is much shorter than the time it takes to experience a strong enough encounter. The presence of the gas hinders the scattering of the comets by the giant planets for not only does it lift the pericenter away from the scattering planet, but in addition by circularizing the orbit it reduces the encounter speed with the next planet and thereby reduces the total energy kick that next planet can give the comet. From the above simple argument it appears to be difficult to deposit km-sized comets in the OC if the density of the gas is the minimum-mass one and the gas disk extends far beyond Neptune, provided the density of the gas remains constant for timescales of ∼105 years. Usually the kicks from Jupiter and Saturn are too strong to deposit the comet in the OC while, for a gas disk extending far beyond Neptune, the comets cannot reach large values of a. Instead the comets are likely to be deposited outside of Neptune. The only way around this is to truncate the disk close to Saturn or Neptune or to reduce the density with time. Numerical simulations need to confirm this theory, and we turn to this next. 5. Numerical simulations and initial conditions The reader is referred to Brasser et al. (2006) for a complete description of the numerical methods employed for simulating the comets and the gas and stars of the embedded cluster. For this project the effect of gas drag and the gas potential was added to SWIFT’s RMVS3 routine (Levison and Duncan, 1994), as described below. A total of 2000 test particles was simulated in each run, which lasted up to 5 Myr. The model includes the gravitational effects of the Sun, Jupiter and Saturn (at 5.45 and 8.18 AU), and in half of the cases studied those of Uranus and Neptune (at 11.5 and 14.2 AU, as proposed by Tsiganis et al., 2005), together with the tidal field of the cluster and passing stars; the effect of the latter are computed by direct integration. The method of generating the initial parameters for the passing stars are discussed in Brasser et al. (2006). In the runs with only Jupiter and Saturn, the comets were distributed between 4 and 10 AU.

In the runs with Uranus and Neptune present, the comets were started between 10 and 16 AU, to allow for Jupiter and Saturn having cleared the regions surrounding them. The computations were performed with a time-step of 0.15 yr. This small step is necessary to adequately resolve the plane crossings of the comets through the gas at about 5 AU when they are on inclined orbits. The same time step was used for the Sun. The positions and velocities of the Sun, planets and test particles were written to disk every 10 kyr in the heliocentric frame, except for those of the Sun, which are written with respect to the cluster center. The orbits of the test particles were computed for up to 5 Myr, or until a particle became unbound and/or was further than some specific distance from the Sun (either a Plummer radius for low-density clusters or half a Plummer radius for high-density clusters), or it collided with a planet or came within 0.005 AU of the Sun. The effect of the gas drag is computed by first evaluating the density and the parameter η at the instantaneous position of the comet √ as well as the Kepler velocity, which is simply vK (r) = GM /r. Next the Mach number and the Knudsen number are evaluated. These we used to compute the Reynolds number and CD according to the recipe given in Section 2.1. This was then used together with the relative velocity between the gas and the comet to compute the acceleration produced by the gas drag. The effect of the gravitational potential of the gas disk is computed using Eqs. (33) and (34). First, a logarithmic grid is generated in cylindrical distance, r, and vertical distance, z, that extends far enough beyond the gas disk that Eq. (35) can be used to evaluate the potential and thus its derived accelerations correctly beyond the grid boundaries. After the grid was generated, the Hayashi potential was evaluated at the grid points by computing the integral in Eq. (34) using the CubPack++1 integration package. The radial and vertical accelerations were computed using numerical differentiation and stored in a similar grid, which is read in by SWIFT. The accelerations at any location are then obtained using bilinear interpolation on the grid, or by using Eq. (35) beyond the grid (typically when r  10rd ). In each of the runs the gas—and thus the gas potential— decays exponentially with an e-folding time of 2 Myr. Thus the force from the disk was rescaled by the appropriate factor. Only two central densities for the cluster are chosen: ρ0 = 103 and 105 M pc−3 . The first density acts fairly similar to the current Galactic environment in terms of torquing timescales, while the latter density is one that has much stronger torquing and readily produces objects like (90377) Sedna. The production runs are listed in Table 1 and the same runs have been done for both cluster densities. The existence of an inner edge of the solar gas disk at 6.2 AU can be argued because Jupiter is thought to have cleared out most of the gas inside it and another strip about 2 Hill radii outside its orbit. The outer edges of the disk can be argued to be the product of photoevaporation (Adams et al., 2004). Clarke (2007) found that for stars in the Orion nebula, only 20% of stars have disks larger than 50 AU. It 1 http://www.cs.kuleuven.ac.be/~nines/software/cubpack.

Star clusters and comet cloud II

Observations are needed to either confirm this statement or rule it out. These will be discussed in the next section.

Table 1 This table shows the various parameters chosen for our runs Run (103 , 105 ) 1, 6 2, 7 3, 8 4, 9 5, 10

423

Inner edge (AU)

Outer edge (AU)

Planets

1 6.2 6.2 6.2 6.2

100 10 100 10 15

JS JS JS JSUN JSUN

6. Results

The first column lists the run number corresponding to each cluster density. The second and third columns display the inner and outer edge of the solar nebula and the fourth column displays what planets are included in the simulation.

is widely believed the Sun’s disk was truncated at 50 or 30 AU (Gomes et al., 2005), which dynamically corresponds to an infinite disk if Saturn is at 8 AU (or Neptune is at 14 AU). Adams et al. (2004) find that in high-density clusters, the disks can be truncated at 20 AU in a few Myr. Truncating the disk at 10 or 15 AU (see Table 1) was done to explore dynamically interesting regions and is not far from Adams et al.’s (2004) results. The reader should be aware that (90377) Sedna itself has a radius of more than 1000 km. As such, the gas drag has almost no influence on its motion and its current orbit can be explained from the results of the previous paper. However, it remains to be seen whether or not objects a few km in size or smaller can end up on orbits like that of (90377) Sedna. If this is not the case, then one can predict some sort of size-sorting might have occurred in the early stages of the Solar System’s formation.

In this section the results of the simulations are presented. 6.1. Overview The parameters of the runs are given in Table 1. The first column states that the runs numbered 1–5 used the same cluster central density, while 6–10 used a different cluster density: runs 1–5 have a central cluster density of 103 M pc−3 while runs 6–10 have a density of 105 M pc−3 . The reader is also reminded that the gas disk surrounding the Sun decays exponentially with an e-folding time of 2 Myr. Fig. 5 shows a vs q at the end of the simulation for runs 1–5, i.e., after 5 Myr of evolution. The top-left panel refers to run 1 (disk truncated at 1 and 100 AU, using a minimum-mass model: mm1-100), the top-right panel is for run 2 (disk truncated at 6.2 and 10 AU, minimum-mass model: mm6.2-10), the middle-left panel refers to run 3 (disk truncated at 6.2 and 100 AU, minimum-mass model: mm6.2-100), the middle-right panel to run 4 (disk truncated at 6.2 and 10 AU with Uranus and Neptune present, minimum-mass model: mmun6.2-10) and run 5 is the bottom panel (disk truncated at 6.2 and 15 AU with Uranus and Neptune present, minimum-mass model: mmun6.2-15). The hor-

Fig. 5. This figure displays the endstates in a−q space for runs 1 to 5. All runs have a cluster central density of 103 M pc−3 . The headers above each caption are indicative of the parameters used. The horizontal lines correspond to the orbital radius of the planets used. See text for details.

424

R. Brasser et al. / Icarus 191 (2007) 413–433

Table 2 This table is a summary of the outcome of the runs that were performed Run AB (q < 4.5)

JS UN OC (4.5 < q < 10) (10 < q < 16) (q > 35)

Ejected Active

1 2 3 4 5 6 7 8 9 10

22.0 78.5 75.7 19.1 16.6 22.1 76.2 18.8 18.9 14.1

0.0 4.6 0.05 18.4 9.7 0.0 4.2 0.0 15.2 8.1

59.6 2.3 2.3 0.7 0.8 60.0 2.5 2.0 1.1 0.7

15.7 0.2 0.2 37.8 54.9 14.8 0.4 73.4 35.6 56.5

0.0 0.8 0.05 3.3 1.9 0.0 1.1 0.0 4.1 2.0

97.2 83.5 80.0 66.0 77.3 96.8 83.1 94.4 69.1 77.3

The fist column lists the run, the second to fifth columns list the fraction of material with their pericenters in the ranges provided while the sixth column states the amount of ejected material. The seventh column lists the fraction of surviving material (i.e., material not removed by ejection or collision with the planets).

izontal lines show the pericenter distances corresponding to the orbital distances of Jupiter and Saturn, as well as those of Uranus and Neptune where these are present. The panel for run 1 shows that the comets end up either on circular orbits inside of Jupiter, as co-orbitals of Jupiter and Saturn or on nearly-circular orbits outside of Saturn. No comets end up in the Oort cloud. Table 2 lists where the pericenter of the comets end up. The letters AB mean the comets are decoupled from Jupiter and are in the ‘asteroid belt’—even though this does not correspond with the asteroid belt we know today. The letters JS mean the comets are in the Jupiter–Saturn region; UN means the comets are in the Uranus–Neptune region and OC means the comets have ended up in the Oort cloud.

As one can see from Table 2, in the case where the disk is truncated at 1 and 100 AU on the inside and outside respectively, close to 60% of the comets end up on circular orbits inside of Jupiter, about 20% ends up in exterior mean motion resonances with Saturn. We found material in the resonances shown in the top panel of Fig. 6. The material that ends up inside of Jupiter, however, is not in resonance with this planet. Instead, the density of material peaks at a ∼ 1.5 AU and quickly drops to zero at a = 1.8 AU. In other words: the material is inside the current asteroid belt and in the region where the terrestrial planets are now. The typical path of material ending up in this region is that a gentle scattering by Jupiter places the comets on an orbit with q ∼ 2 AU. For this value of q the time to decrease the apocenter away from Jupiter’s orbit is much shorter than the time to receive a strong scattering by Jupiter, so that the orbit quickly circularizes and subsequently spirals in on a timescale given in Eq. (28). Moving on to run 2, the top-right panel of Fig. 5 shows a very different outcome. Since the inner edge of the disk is truncated at 6.2 AU, almost no material ends up inside of Jupiter. Having the outer edge truncated at 10 AU means that the material does not end up on (nearly) circular orbits outside of Saturn due to gas drag. Instead, the majority of the comets have q in the JSzone (78%) with almost all of these on eccentric orbits with q ∼ 9 AU in exterior resonances with Saturn. We have zoomed in on the region just outside of Saturn for runs 1 and 2—the result from run 3 is virtually the same as for run 1—and provided the results in Fig. 6. The exterior resonances with Saturn are labeled as such. It is clear that some resonances, such as the 2:3S and the 1:2S, are preferred over others (e.g., the 2:5S or the 3:5S).

Fig. 6. This figure shows q vs a for material trapped in exterior mean-motion resonances with Saturn. The top-panel is for run 1 (mm1-100) while the bottom panel is for run 2 (mm6-10). The dashed line corresponds to a circular orbit. Labels above a clump of dots show each resonance that is populated. Note that in the bottom panel, where the disk is truncated at 10 AU, almost all comets have q ∼ 9 AU.

Star clusters and comet cloud II

The 3:4S, 2:3S, 3:5S and 1:2S are the most heavily populated in the case the disk is truncated at 100 AU, while the 2:3S and the 1:2S are most heavily populated when the disk is truncated at 10 AU. It may seem surprising that comets in resonances with Saturn are not on circular orbits. In the case of the bottom panel, where the disk is truncated at 10 AU, the orbits are expected to be eccentric with q ∼ 9 AU as is observed, but in the top panel, where the disk is truncated much further out, a remnant eccentricity seems strange. However, the case of resonant trapping of planetesimals migrating due to gas drag has been considered by Weidenschilling and Davis (1985) where it was shown that during the trapping process an equilibrium can be achieved in which eccentricity excitation due to resonant interactions can be balanced by gas drag damping. An approximation for the remnant eccentricity at a resonance in the presence of aerodynamic drag comes from the method of Kary et al. (1993). Using a = aS ((j + p)/j )2/3 and a = aS [((j + p)/j )2/3 − 1] for any j/(j + p) exterior resonance and writing k = ((j + p)/j )1/3 , we have √ 2E( 3/2) 3 3ηk 2 e + e π 2  √ ηk 2 (k 2 − 1)  √ − 2E( 3/2) + K( 3/2) e π 3η2 k 2 (k 2 − 1) = 0. − (45) 2 Evaluating η at Saturn’s orbit gives

 0.771e3 + 9.138 × 10−3 ke2 − 8.879 × 10−3 k 3 k 2 − 1 e

 − 5.567 × 10−5 k 4 k 2 − 1 = 0. (46) For the 1:2S resonance, solving Eq. (46) gives eeq = 0.113, which is very close to the observed value of about 0.105, and the approximation improves for resonances closer to Saturn. In all cases, the eccentricities found in our simulations agree very well with the theoretical predictions. It should be noted that there exists a range in sizes rc,m < rc rp of the comets for resonance trapping to occur for orbits with e η. We follow Weidenschilling and Davis (1985) who analyze this trapping for (j + 1)/j resonances. The minimum size a comet can have to be trapped is then     ρg ap 1 g cm−3 rc,m = 2.30 × 105 ρc 1 AU 1.4 × 10−9 g cm−3  −1 √ η mp × (47) km, 7/6 1/3 m⊕ C(α)j (j + 1) where j is the order of the resonance and C(α) is a function of the semi-major axis ratio, α, between the resonance and the planet. For this problem one can take α = [j/(j + 1)]2/3 since the precession of  is slow compared to the orbital period. The function C(α) is a function of the well-known Laplace coef(j ) ficients bs (α) and is given by (Weidenschilling and Davis, 1985) (j )

(j )

C(α) = (1 + 2j )b1/2 (α) + α

db1/2 (α) dα

−

δ1j , α2

(48)

425

where δ1j is a Kronecker delta. The Laplace coefficients are given by   −s j Γ (s − j )

1 + α2 (j ) bs (α) = 2 (49) 1 − α 2 Ps−1 , Γ (s) 1 − α2 μ

where Γ (x) is the Euler gamma function and Pν (x) is an associated Legendre function. For our standard Hayashi disk, given by Eqs. (5) and (9), Eq. (47) becomes     ρg0 ap −3/2 1 g cm−3 rc,m = 1060 ρc 1 AU 1.4 × 10−9 g cm−3  −1 1/2 mp j × km. (50) C(α)(j + 1)2 m⊕ Evaluating this for the 2:3S resonance with j = 2, mp = 95m⊕ , ap = 8.18 AU, ρc = 0.6 g cm−3 and ρg0 = 1.4 × 10−9 g cm−3 yields rc,m = 24 m. The reader should keep in mind that this size threshold only applies for comets approaching the resonance on nearly-circular orbits (Weidenschilling and Davis, 1985), and that this threshold comet size cannot be applied for cases with eccentricities e  η. What is important to note is the product μC(α) in the denominator of Eq. (47). The function C(α) → ∞ when α → 1 as one can see from Eq. (49), implying that a less massive planet will trap material at mean-motion resonances that lie closer to its orbit than a more massive planet. However, this argument breaks down once a − ap is inside the region of chaotic motion surrounding the planet. From Duncan et al. (1989) for a circular planetary orbit a − ap ∼ 0.04(mp /m⊕ )2/7 ap AU, i.e.,     mp 2/7 j + 1 2/3 (51) − 1 = 0.04 . j m⊕ One can solve Eq. (51) for j given mp and then compute both the closest stable resonance to the planet and the corresponding minimum size for trapping using Eq. (50). For Saturn the solutions are j = 4 and using the values of densities and planet mass above, yields rc,m = 8 m. For the 1:2S resonance, however, rc,m = 230 m, so that material spiraling in towards Saturn on orbits with e η and rc > 230 m will mostly get trapped in the 1:2S resonance and any material with rc < 8 m will reach Saturn. Hence the trapping range in comet size for intermediate resonances is small. From Fig. 6 one can see no material is trapped in resonances closer to Saturn than the 4:5S, confirming the above result. To test the above results, and to generate a possible outcome for the material inside of Jupiter from run 1 (mm1-100), we repeated the same run and placed Mars-sized embryos 12 Hill radii apart starting at 1 AU and ending at 3.4 AU. The result is that some material gets trapped in mean-motion resonances with the embryos, but a fair fraction of the comets reach 1 AU. However, when the embryos were restricted to less than 1.8 AU, i.e., far from Jupiter, then material that approached the outermost embryo was not excited and subsequently scattered inward, but settled into high-order exterior mean-motion resonances with said embryo, such as the 9:8, 12:11, 15:14 and 21:20. When we used Eq. (45) to compute the equilibrium eccentricities, we found excellent agreement with observations

426

R. Brasser et al. / Icarus 191 (2007) 413–433

Fig. 7. Various trajectories are plotted in the a−q plane of comets ending up in exterior mean-motion resonances with Saturn. The top two panels are for a disk truncated at 10 AU while the bottom panels are for a disk truncated at 100 AU. Once again the horizontal lines correspond to the orbital radii of the planets while the upward-sloping line corresponds to a circular orbit. See text for details.

once again and also with Kary et al. (1993). We refer the reader to Kary et al. (1993) for a more detailed study of the resonance trapping, which goes beyond the scope of this paper. The route a comet takes to land in the resonance is caused by a combination of scattering and drag. We have provided examples of different routes in Fig. 7. The top two panels show possible routes in the (q, a) plane when the disk is truncated at 10 AU, while the bottom two panels are for routes when the disk is truncated at 100 AU. Since our output interval is 10 kyr, a lot of the fine details are lost, but at least one can get some understanding of the underlying mechanism. The topleft panel shows a comet that gets scattered to a ∼ 10.7 with e ∼ 0.3. The comet is caught in the 2:3S resonance and the orbit is circularized at nearly constant semi-major axis until the eccentricity reaches the equilibrium value. This is an example of what we dubbed a ‘vertical case.’ The top-right panel shows a completely different story. The comet is locked in an interplay between Saturn and the gas. During this episode of being kicked and dragged, q increases ever so slightly, until it reaches a point where the kicks from Saturn are not strong enough to overcome the action from the gas. The comet then spirals inward at nearly constant q until it hits the 1:2S resonance, and gets trapped. We dub this type of motion a ‘horizontal case.’ The bottom-left panel is a case where the gas disk is truncated at 100 AU. The comet rattles around until it reaches an orbit with a ∼ 14 and q ∼ 9. At this point it is sufficiently eccentric and decoupled from Saturn that gas drag wins over kicks produced by Saturn. The orbit circularizes—both a and q increase—until e ∼ 0 and a ∼ 11. At this point the comet is close to the 2:3S resonance

(a = 10.72 AU) and gets trapped in it. Since the equilibrium eccentricity in the 2:3S resonance is about 0.07, the eccentricity increases at fixed semi-major axis until it reaches the equilibrium value. Last, the bottom-right panel shows another vertical case, but unlike the case in the top-left panel, the comet is circularized until q ∼ 11.5 AU while in the previous case, q ∼ 9. We can now focus our attention back to Fig. 5. The middleleft panel (mm6.2-100) of Fig. 5 once again shows that no material has ended up in the Oort cloud. This case is similar to the top-left panel, with the exception that now instead of the majority of the comets being inside of Jupiter, this region is no longer available so they end up on orbits outside of Saturn, with again most of them in a resonance with the latter planet. The end results are completely different once Uranus and Neptune are present. The existence of these planets prevents the comets from settling onto orbits in exterior resonances with Saturn and since we truncated the disk at 10 AU, the pericenters of the comets can be lifted away from Jupiter and kept at Uranus and Neptune. The result is that about 40% of the comets are in the UN–SD with 20% in the JS–SD and about 3% is in the OC, much more than in the previous cases. The material in the UN–SD will be followed in a later simulation to determine how much of that material eventually ends up in the OC. It should be noted that only about 66% of the comets are still active. The remainder has been ejected or collided with one of the giant planets. This number is also much higher than in the previous three cases, as is shown by the number of remaining comets in Table 2.

Star clusters and comet cloud II

427

Fig. 8. The endstates for runs 6 through 10 are displayed in a manner similar to Fig. 5.

A similar story applies for the bottom panel of Fig. 5, where the gas disk is now truncated at 15 AU, i.e., just outside of Neptune’s orbit. The result of this is that material in the UN–SD has had the pericenters lifted to just outside of Neptune’s orbit, where their evolution is slow: 55% of the comets are in the UN–SD. A large fraction of these are Uranus Trojans, while a smaller number are Neptune Trojans. The scattered objects have 14.2 < q < 15.5 AU and 15 < a < 25 AU. Fig. 8 depicts the endstates of runs 6–10, which have ρ0 = 105 M pc−3 . Apart from the OC being closer in because the cluster is denser, the end results are statistically the same as in the case of Fig. 5 (see Table 2). The results presented above are for comets the size of 1.7 km. This appears to be a ‘typical’ radius, as has been explained before. The above results also indicate that if the disk was truncated at a large distance from the planets (100 AU in our case) then we expect virtually none of this size comets to end up in the OC and instead expect them to end up on orbits outside of Saturn in cases in which Uranus and Neptune are not yet present. This material could then be used to grow Uranus and Neptune or, if those planets are already there, end up on nearly-circular orbits outside of Neptune and affect its migration. However, the deceleration due to drag is inversely proportional to the comet’s radius. Therefore, one would think that the presence of the gas acts as a size-sorting mechanism. If this is truly the case then simulations using different-sized comets should provide the details. This is discussed in the next subsection.

6.2. Size sorting Until now, we have only examined the results of simulations performed with different disk sizes (and therefore disk masses), but have not examined the effect of the sizes of the comets. In the case where the effect of the gas is pure aerodynamic drag, resizing the comets is equivalent to rescaling the disk density, but since the disk also acts as a perturbing potential, this simple relation no longer holds. We therefore chose to perform a set of simulations for comets with different sizes, keeping all the other parameters fixed. Thus we chose to use a standard disk ranging from 6.2 to 100 AU, a cluster central density of ρ0 = 103 M pc−3 and we included the influence of Uranus and Neptune. These simulations were run for 3.5 Myr, in order to better compare the end results with those of Brasser et al. (2006) (which had no gas present). The summary of the results is given in Table 3 and Fig. 9, which are snapshots in a–q space of the endstates of all the runs except that with rc = 300 km. As can be seen, once rc > 10 km or so, the results appear not very different. This is also apparent in Table 3 since the efficiency does not increase significantly for comet radii larger than about rc = 30 km, suggesting a change in behavior for cometary radii in this range. As has been explained in Section 4 and proven in the previous subsection, the comets need to be scattered to a  ag in order to end up in the Oort cloud. It was mentioned in Section 4 that this is an unlikely event for 2 km comets and as can be seen from the top-left and middle-left panels of Figs. 5 and 8,

428

R. Brasser et al. / Icarus 191 (2007) 413–433

Fig. 9. Snapshots in the a–q plane of the endstates at 3.5 Myr for runs with similar disks and cluster densities but different comet sizes. The labels above each panel denote the size of the comet used. As can be seen, once rc > 10 km or so, the end results appear very similar. Table 3 This table is similar to Table 2, with the exception that the first column now lists the radius of the comets used during the simulation Comet size (km)

AB

JS

UN

OC

Ejected

Kozai (q > 16, i > 40◦ , a < 103 )

Active

1.7 5 10 30 100 300 1000

0.1 0.5 0.6 0.6 0.3 0.4 0.4

5.4 4.3 4.4 3.8 7.7 8.1 7.8

88.7 76.2 54.8 33.6 25.4 23.9 22.6

0.0 1.2 4.2 8.1 10.7 10.8 11.8

1.3 8.1 25.1 45.1 46.3 47.4 53.9

0.0 0.0 0.05 0.2 0.05 0.5 0.2

92.6 84.6 68.6 50.1 49.3 49.1 47.1

The seventh column lists how many particles were affected by the Kozai resonance induced by the disk and had the parameters in the range shown at the top of the column. The cluster density and gas disk used are the same for all runs in this table. See text for details.

no comets ended up in the Oort cloud, confirming our theory. However, D1 ∝ rc−1 or, rather, ag ∝ rc−1 , and larger comets can reach the OC easier. This implies there is a threshold value of rc for which the value of ag is such that the planets can comfortably scatter material to the region where the comets can enter the OC. Looking at the results of Table 3, the value rc = 10 km seems to be a turning point, where the effects from the gas and the planets are about equal. In the case where rc = 30 km, about half of the comets are ejected, showing that the planetary perturbations largely overwhelm than the gas drag. In fact, tgq ∝ rc so that increasing the comet’s size by a certain factor increases

the time to drag it back by the same amount, or, in other words, it decreases ag . One can project this onto Fig. 4 to determine what size comets are needed for ag in the cases of Jupiter and Saturn to be of order 2ap . In the case of Jupiter, it turns out to be rc ∼ 10 km while for Saturn it is rc ∼ 30 km, in excellent agreement with simulations. In addition, the efficiency of depositing material in the OC seems to asymptote to 12%. From Brasser et al. (2006), the typical efficiency for a cluster of the same density is also roughly 12%, so that the current results for large objects are similar. This proves two things: first, the influence of Uranus and Neptune has little to no effect on the outcome, proving that the main contributors to the OC on this timescale are Jupiter and Saturn, as expected; and second, that for large enough bodies (larger than about 30 km), the gas drag plays virtually no role in the motion, as was expected too. There is, however, an important difference between the end results of these runs and those of Brasser et al. (2006). The disk potential induces the Kozai mechanism which can potentially increase the inclination of the object caught in it and raise its pericenter. We did not see this happening in our main runs summarized in Table 2 because the gas drag is too strong to allow this subtle effect to show. However, for the larger-sized comets, this effect was clearly visible. The natural question to ask is the following: how many comets end up on inclined, quasi-circular orbits with their pericenters well outside the influence of Neptune? We decided to use the criteria that q > 16 and i > 40◦ and a < 1000 AU. Some objects that end up in the OC also satisfy

Star clusters and comet cloud II

429

Fig. 10. This figure shows the effect of the Kozai mechanism induced by the gravitational potential of the gas disk on a comet. The two left panels display a and q as a function of time, while the two right panels show q vs ω. The libration of ω around 90◦ and 270◦ is particular, and common, of the Kozai mechanism. In the q vs ω panels, the bullets show the libration phase corresponding to the oscillations in q in the left-side panels while the ‘+’ symbols show are during the scattering phase. The top orbit has q ∼ aN so that it is not saved, but the bottom orbit has q > aN so that it is saved from close encounters.

those criteria, but naturally these are excluded from the sample. The fraction of objects that satisfy these criteria are listed under ‘Kozai’ in Table 3, while Fig. 10 is an example of two objects exhibiting this kind of behavior, with ω librating around ±90◦ . In any case, the presence of the primordial solar nebula acts as a size-sorting mechanism: for a large disk, small comets will be trapped in mean motion resonances with Saturn (or other planets, if they are present) and will not reach the Oort cloud, while large objects are unaffected by the drag and will reach the OC as in Brasser et al. (2006). This has implications for objects like (90377) Sedna and 2000 CR105 . In Brasser et al. (2006) we demonstrated that the orbits of such objects can be explained to have arisen from the time the Sun was in a dense (enough) star cluster. These objects are large enough to, in this model, have been unaffected by the gas drag. Therefore, if our model is representative of the actual events that occurred at that time, we do not expect many objects on (90377) Sedna-like orbits because the presence of the gas prevents them from attaining a largeenough semi-major axis. Instead, the smaller bodies end up as fuel for the formation of Uranus and Neptune or are the catalyst to start Neptune’s eventual migration. These bodies will subsequently be scattered by the giant planets and a small fraction (3%, Dones et al., 2004) will land in the current OC. One last thing the reader should keep in mind is that we neglected physical collisions between the comets. For large bodies this turns out to be not a problem but for small bodies the situation is not so clear: Stern and Weissman (2001) argue collisions

play a major role in OC formation while Charnoz and Morbidelli (2003) argue the opposite. Given that this is a simplified model, however, we believe that ignoring physical collisions is satisfactory. 7. Summary and conclusions Computer simulations have been performed of OC formation from comets in the Jupiter–Saturn region of the pre-LHB Solar System while the Sun was in an embedded cluster and with the additional inclusion of the dynamical effects of the primordial solar nebula. The gas disk was assumed to be initially a minimum-mass Hayashi disk with the density decaying exponentially with an e-folding time of 2 Myr. Our simulations showed that: • For disks where the outer edge is far (5 AU) from the outermost planet (either Saturn or Neptune), no OC is formed with km-sized bodies and the comets end up in exterior mean-motion resonances with the outermost planet. • When the inner edge of the disk is much closer to the Sun than Jupiter’s orbit, the majority of the comets (∼60%) end up on circular orbits inside of Jupiter and ∼20% end up beyond the outer planet. • When the outer edge of the disk is truncated a few AU beyond the orbit of the outermost planet, an OC will form but

430

R. Brasser et al. / Icarus 191 (2007) 413–433

the efficiency is much lower than in the case with no gas present. • Increasing the size of the comets decreases the effect of the gas drag. The effect of the gas drag becomes comparable to the scattering of Jupiter and Saturn for comets the size of 10 and 30 km, respectively. Therefore, the presence of the gas acts as a size-sorting mechanism with large bodies reaching the OC while small ones remain trapped in the planetary region. • A byproduct of some of our simulations are endresults with a substantial fraction of the comets in the Uranus–Neptune scattered disk. A study of this second phase of trapped material being scattered out is planned for the future. Acknowledgments We are indebted to Doug McNeil for making available to us his code for computing the Hayashi potential and the derived accelerations. R.B. and M.J.D. thank H.F.L. and Alessandro Morbidelli who made them aware of the effect of the gravitational potential of the primordial solar nebula while M.J.D. was on a visit to Obs. Nice, France. R.B. and M.J.D. thank NSERC of Canada for funding this research. R.B. also acknowledges funding from CITA through their National Fellowship program. H.F.L. acknowledges funding from NASA’s ORIGINS and PPS programs. Appendix A During the formation of the OC, most comets reach nearparabolic orbits before the tides and stars can lift their pericenters out of the planetary region. During this stage, the effect of the gas drag occurs close to pericenter and can be modeled fairly well as an energy kick applied at each pericentric passage. What is not obvious is that for the Hayashi model (Hayashi et al., 1985), the angular momentum averaged over the orbit is approximately conserved in the planar case. This phenomenon is computed here for a parabolic orbit, since this is a good approximation for comets in the primordial scattered disk on their way to the OC. Since the inclinations are small, typically less than 30◦ , the planar case still has some validity. For any Kepler orbit, the relation between the true anomaly, f and the distance to the Sun, r is r f˙ = h, p r= , 1 + e cos f where h is the specific angular momentum and p is the semi-latus rectum. The relation between h and p is that p = h2 /GM . For a parabola, e = 1 and p = 2q, so one has  r 2 f˙ = 2GM q, (52) 1 r = q sec2 f. (53) 2 Substituting Eq. (57) into Eq. (56) gives  2GM 4 1 dt, sec f df = (54) 2 q3 2

which can be integrated to give  2GM 1 1 2 (t − t0 ) = 2 tan f + tan3 f, 2 3 2 q3

(55)

and relates f to the time t . In a Cartesian frame where the Sun is at the origin and the pericenter is at (q, 0) the parabola is represented as x = r cos f = 2q − r,  r − 1. y = r sin f = 2q q

(56)

The components of the velocity of the comet are √  x 2GM q r ˙ x˙ = r˙ − y f = − − 1, r r q √ 2GM q y . y˙ = r˙ + x f˙ = (57) r r The circular velocity of the gas needs to be computed at the same coordinates as the velocity of the comet, i.e., at the same place in the (r, f ) plane. This yields    GM r y GM x˙g = − = −2q − 1, r r q r3     GM 2q x GM = − 1 . y˙g = (58) r r r r3 v − vg ), we need to comSince the drag force is F ∝ | v − vg |( pute x˙ − x˙g and y˙ − y˙g . These are      GM r 2q 2q x˙ − x˙g = −1 − , r q r r     GM 2q 2q 1+ − . y˙ − y˙g = (59) r r r The magnitude of the relative velocity is       GM 8q vr =  3− . r r

(60)

It is the above quantity that makes the integrals that appear below non-trivial. The dissipative drag force per unit mass is of the form   3 CD ρ0g r −α  A=− (61) vr vr , 8 rc ρc r0 where ρg0 is the density of the gas at r0 . The subsequent torque per unit mass on the comet is given by h˙ = r × A and is in the z-direction. Using the substitution r = q ξ 2 , the resulting torque is given by √ √ 2) 3ξ − 2 2 (ξ − p h˙ = GM Zg (62) , ξ 2α+3/2

Star clusters and comet cloud II

√ √    1 1√ 227 6 215 3 2 − RF 0, 1, − + 288 192 2 3 ≈ −0.0358846, (69) 

where p Zg

431

  3 ρg0 CD r0 α = 8 ρc rc q     ρg0 1 km 1 g cm−3 = 0.0785CD rc ρc 1.4 × 10−9 g cm−3 α  α  r0 1 AU AU−1 . × (63) q 1 AU

The total change in specific angular momentum per pericenter passage is therefore  ∞ 3 ξ q3 h = 4 h˙  dξ ξ 2 − 1 2GM 1





= 2 2GM

q 1 AU

3/2

p

Zg B1 (α),

(64)

where

∞ (ξ − √2) 3ξ − 2√2  B1 (α) = d ξ. ξ 2α−3/2 ξ 2 − 1

(65)

1

√ From Eq. (65) it is easy to see that h˙ = 0 when ξ = 2, i.e., when r = 2q. The integral diverges when α  3/2, meaning that if α  3/2, h is infinite. The change in specific energy of the orbit can be computed  as the orbit-averaged value of d(1/a)/dt = −2( v · A)/GM . Hence  √ √ ( 2ξ − 1) 3ξ − 2 2 d(1/a) 2GM p (66) Zg =2 dt q3 ξ 2α+9/2 p

so that, for Zg given by Eq. (63) (1/a) = 8Zg B2 (α) AU−1 , p

(67)

where

∞ (√2ξ − 1) 3ξ − 2√2  B2 (α) = d ξ. ξ 2α+3/2 ξ 2 − 1

 B2

11 4



√   212861 3 1 1√ 2 = RD 0, 1, − 21233664 2 3 √ √   310243 6 310243 3 + − 1769472 1176948   √ 2√ 1 1 2, 1 − 2 × RJ 0, 1, − 2 3 3 √ √   310243 3 1625851 6 + − 393216 2949120   1 1√ × RF 0, 1, − 2 2 3 ≈ 0.183456,

(70)

where RF (x, y, z), RD (x, y, z) and RJ (x, y, z, w) are Carlson’s elliptic integrals of the first, second and third kind respectively. One can convert the above equations to a relation involving the classic elliptic integrals by using

 RF 0, 1, 1 − k 2 = K(k),

 3 3 E(k) − 2 K(k), RD 0, 1, 1 − k 2 = 2 (71) 2 k (1 − k ) k

 3

 RJ 0, 1, 1 − k 2 , 1 − n = Π(n, k) − K(k) . (72) n The results agree numerically when using the substitutions above. Since B2 (11/4) > 0, (1/a) > 0 and a shrinks, as expected. Evaluating Eq. (67) at 5 AU for our reference Hayashi disk and a comet of size 1.7 km, we have (1/a) = 2.704 × 10−3 AU−1 . Numerical experiments with SWIFT give (1/a) = 2.646 × 10−3 AU−1 . For the angular momen√ tum, evaluation of Eq. (64) yields h = 2.068 × 10−3 GM while√numerical integration with SWIFT yields h = 2.711 × 10−3 GM . If the disk is truncated, the integrals in Eqs. (65) and (68) should be over the range the disk where the density is non-zero rather than to infinity.

(68)

1

Both of the integrals in Eqs. (65) and (68) are cubic elliptic integrals of the third kind. When trying to reduce these to Legendre’s standard integrals, one needs to evaluate limits that end in singularities, so that Carlson’s integrals (e.g., Carlson, 1999) are preferred. In addition, these latter integrals do not suffer from the restrictions that the Legendre integrals have. For the standard Hayashi model, 2α + 3/2 = 7 and 2α − 3/2 = 4. Using Carlson’s (1999) relations we have √     137 3 11 1 1√ B1 2 =− RD 0, 1, − 4 10368 2 3 √ √   215 3 215 6 + − 576 864   2√ 1 1√ × RJ 0, 1, − 2, 1 − 2 2 3 3

Appendix B. Energy diffusion In this section, we follow van Woerkom (1948) to derive the diffusion equation. At the epoch t let ν(1/a, t) d(1/a) be the number of comets passing through pericenter per year and having reciprocal semimajor axes in the range 1/a ∈ [1/a, 1/a + d(1/a)]. If N (1/a, t) is the total number of comets with the same value of 1/a then (van Woerkom, 1948) 2π N (1/a, t) = P (a) ν(1/a, t) = √ a 3/2 ν(1/a, t). GM

(73)

Define a function h(u) to be the normalized probability density function of the change in energy, u, per pericenter passage (Everhart, 1968). In other words h(u) = N/ u where N is the number of comets receiving a change in 1/a equal to u in the range [u, u + u]. For large-enough eccentricities the

432

R. Brasser et al. / Icarus 191 (2007) 413–433

function h(u) is independent of the semi-major axis, since the orbit is sufficiently near-parabolic that any subsequent pericenter passages keep the pericenter distance q approximately fixed. The fact that q is nearly fixed can be proven by conservation of the Tisserand parameter in the case (1/a) 1/q and, for large semi-major axis, using Q ≈ 2a. For the above equation to be valid, two successive conjunctions between the comet and the planet need to be chaotic (see Appendix C). During one year, the number of comets with a definite value of 1/a will increase by (van Woerkom, 1948) ∞ ν(1/a − u, t)h(u) du,

(74)

−∞

while those comets which originally had this value of 1/a will be perturbed to other orbits. However, since the above theory is only an approximation, the values ±∞ will suffice. The total number of comets with reciprocal semi-major axis 1/a will therefore increase by ∂N = ∂t

∞ ν(1/a − u, t)h(u) du − ν(1/a, t).

(75)

−∞

The integral-differential equation (75) can be reduced to a Fokker–Planck problem by expanding the quantity ν(1/a − u, t) around u = 0, i.e., assuming that the energy kicks u are √ small compared to 1/a. Using N (1/a, t) = (2π/ GM ) × a 3/2 ν(1/a, t) one has √

∂ν(1/a, t) 2π a 3/2 ∂t GM ∞  ∞ n n u ∂ ν(1/a, t) h(u) du − ν(1/a, t). = n! ∂un

and the constant c1 is determined by the initial conditions. The above two equations cannot be solved analytically and therefore one needs to resort to numerical methods with acceptable boundary conditions, which is, however, beyond the scope of this paper. Appendix C. Limiting semi-major axis for energy diffusion The previous appendix modeling the diffusion of comets under the influence of one or more of the giant planets is only applicable to high-eccentricity orbits with the pericenter close to the orbit of a planet. In essence, the diffusion approximation is only correct when the encounters between the planet and the comet are uncorrelated; in other words, the planet’s phase at the time of perihelion passage of the comet varies chaotically because of perturbations by the planet on the comet’s orbit with every passage. The limiting orbital energy or semi-major axis at which this phenomenon occurs is derived here. Fernandez (2005) explains that for very eccentric comets, the binding energy, x, is small compared to xp so that the perturbations per pericenter passage, x can be a significant fraction √ of x. The period of the comet is given by P = 2πx −3/2 / GM∗ and so P 3 x (81) = . P 2 x The changes in the period, or binding energy, become chaotic when P ∼ 12 Pp —see, e.g., Duncan et al. (1989)—and is illustrated in Fig. 11 in the case the planet’s orbit (in this case Jupiter), is circular. When this concept is introduced in Eq. (81),

(76)

−∞ n=0

Assuming that the sum and integral signs can be exchanged, ∂ ∂ truncating the series when n = 2, and using ∂u = − ∂(1/a) and x = 1/a √   GM 3/2 ∂ν(x, t) 1 ∂ 2 ν(x, t) ∂ν(x, t) = x + D2 D1 , ∂t 2π ∂x 2 ∂x 2 (77) where ∞ Dn = (78) un h(u) du. −∞

Strictly speaking, the limits of the integration are not ±∞ but are restricted by the radius of the planet. In addition, it means taking the integration to the far tails of h(u), where the diffusion approximation breaks down. The Fokker–Planck equation (77) can be solved by separation of variables so that ν(x, t) = μ(t)ξ(x) where μ(t) and ξ(x) satisfy dμ(t) = c1 μ(t), dt c1 ξ(x) d2 ξ(x) D1 dξ(x) + = 0, − 2 D2 dx dx D2 x 3/2

(79) (80)

Fig. 11. This figure is a sketch depicting Jupiter on a circular orbit of period PJ and a comet on a period N × PJ . Initially, the conjunction happens at J1 –C1 –S. If the comet’s orbit were to remain unperturbed, the same conjunction would occur after time N ×PJ as J2 –C2 –S. However, if the comet’s period is perturbed by 12 PJ , then the situation at the next pericenter passage will be J2 –S–C2 , i.e., Jupiter is now on the other side of the Sun. Figure courtesy of Fernandez (2005).

Star clusters and comet cloud II

Table 4 This table is a summary of the semi-major axis beyond which a comet with pericenter locked to a planet is said to be diffusing Planet

ap (AU)

10μ/ap

ac (AU)

Jupiter Saturn Uranus Neptune

5.45 8.18 11.5 14.2

1.75 × 10−3 3.50 × 10−4 3.80 × 10−5 3.63 × 10−5

23 55 164 189

The first column lists the planet’s name, the second column is the planet’s semimajor axis in the pre-LHB model, the third column lists a typical value of the rms energy change per pericenter passage. Last, the fourth column shows the critical semi-major axis beyond which the comet is said to be diffusing.

the critical semi-major axis becomes   2 P 2/5 ac  . 3 x

(82)

Given that for most perihelion passages, x ∼ 10μ(1/ap ) (DQT87) we have ac ∼ (30μ)−2/5 ap ,

(83)

where μ = mp /M∗ . Table 4 lists the various values of μ and ac for each planet on the pre-LHB configuration orbits. References Adachi, I., Hayashi, C., Nakazawa, K., 1976. The gas drag effect on the elliptic motion of a solid body in the primordial solar nebula. Prog. Theor. Phys. 56, 1756–1771. Adams, F.C., Hollenbach, D., Laughlin, G., Gorti, U., 2004. Photoevaporation of circumstellar disks due to external far-ultraviolet radiation in stellar aggregates. Astrophys. J. 611, 360–379. Binney, J., Tremaine, S., 1987. Galactic Dynamics. Princeton Univ. Press, Princeton, NJ. Brasser, R., Duncan, M., Levison, H., 2006. Embedded star clusters and the formation of the Oort cloud. Icarus 184, 59–82. Carlson, B.C., 1999. Toward symbolic integration of elliptic integrals. J. Symb. Comp. 28, 739–753. Charnoz, S., Morbidelli, A., 2003. Coupling dynamical and collisional evolution of small bodies: An application to the early ejection of planetesimals from the Jupiter–Saturn region. Icarus 166, 141–156. Clarke, C.J., 2007. The photoevaporation of discs around young stars in massive clusters. Mon. Not. R. Astron. Soc. 376, 1350–1356. Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops and Particles. Academic Press, New York. Davidsson, B.J.R., Gutierrez, P.J., 2006. Non-gravitational force modeling of Comet 81P/Wild 2. Icarus 180, 224–242. Davidsson, B.J.R., Gutierrez, P.J., Rickman, H., 2007. Nucleus properties of Comet 9P/Tempel 1 estimated from non-gravitational force modeling. Icarus 187, 306–320. Dones, L., Weissman, P.R., Levison, H.F., Duncan, M.J., 2004. Oort cloud formation and dynamics. In: Johnstone, D., Adams, F.C., Lin, D.N.C., Neufeld, D.A., Ostriker, E.C. (Eds.), Star Formation in the Interstellar Medium, vol. 324. ASPC, San Francisco, p. 371.

433

Duncan, M., Quinn, T., Tremaine, S., 1987. The formation and extent of the Solar System comet cloud. Astron. J. 94, 1330–1338. Duncan, M., Quinn, T., Tremaine, S., 1989. The long-term evolution of orbits in the Solar System—A mapping approach. Icarus 82, 402–418. Everhart, E., 1968. Change in total energy of comets passing through the Solar System. Astron. J. 73, 1039–1052. Fernandez, J.A., 2005. Comets, Nature, Dynamics, Origin, and Their Cosmogonical Relevance. Springer-Verlag, New York. Francis, P.J., 2005. The demographic of long-period comets. Astrophys. J. 635, 1348–1361. Gladman, B., Duncan, M., 1990. On the fates of minor bodies in the outer Solar System. Astron. J. 100, 1680–1693. Gomes, R., Levison, H.F., Tsiganis, K., Morbidelli, A., 2005. Origin of the cataclysmic late heavy bombardment period of the terrestrial planets. Nature 435, 466–469. Hayashi, C., Nakozawa, K., Nakagawa, Y., 1985. Formation of the Solar System. In: Black, D.C., Matthews, M.S. (Eds.), Protostars and Planets II. Univ. of Arizona Press, Tucson, AZ. Higuchi, A., Kokubo, E., Mukai, T., 2002. Cometary dynamics: Migration due to gas drag and scattering by protoplanets. In: Warmbein, B. (Ed.), Proc. ACM 2002. ESA Publications Division, Noordwijk, Netherlands. Ida, S., Makino, J., 1993. Scattering of planetesimals by a protoplanet— Slowing down of runaway growth. Icarus 106, 210–227. Kary, D.M., Lissauer, J.J., Greenzweig, Y., 1993. Nebular gas drag and planetary accretion. Icarus 106, 288–307. Kiang, T., 1962. The effect of a resisting medium on elliptical orbits. Mon. Not. R. Astron. Soc. 123, 359–382. Lada, C.J., Lada, E.A., 2003. Embedded clusters in molecular clouds. Annu. Rev. Astron. Astrophys. 41, 57–115. Landau, L.D., Lifschitz, E.M., 1987. Course of Theoretical Physics, vol. 6: Fluid Mechanics, second ed. Pergamon, New York. Levison, H.F., Duncan, M.J., 1994. The long-term dynamical behavior of shortperiod comets. Icarus 108, 18–36. Oort, J.H., 1950. The structure of the cloud of comets surrounding the Solar System and a hypothesis concerning its origin. Bull. Astron. Inst. Neth. 11, 91–110. Phillips, W.F., 1975. Drag on a small sphere moving through a gas. Phys. Fluids 18, 1089–1093. Stern, S.A., Weissman, P.R., 2001. Rapid collisional evolution of comets during the formation of the Oort cloud. Nature 409, 589–591. Tanaka, H., Ida, S., 1996. Distribution of planetesimals around a protoplanet in the nebula gas. Icarus 120, 371–386. Tsiganis, K., Gomes, R., Morbidelli, A., Levison, H.F., 2005. Origin of the orbital architecture of the giant planets in the Solar System. Nature 402, 635–638. Weidenschilling, S.J., 1977. Aerodynamics of solid bodies in the solar nebula. Mon. Not. R. Astron. Soc. 180, 57–70. Weidenschilling, S.J., Davis, D.R., 1985. Orbital resonances in the solar nebula: Implications for planetary accretion. Icarus 62, 16–29. van Woerkom, A.J.J., 1948. On the origin of comets. Bull. Astron. Inst. Neth. 10, 445–472. Yabushita, S., 1979. Statistical tests of distribution of perihelion points of longperiod comets. Mon. Not. R. Astron. Soc. 189, 45–56. Yabushita, S., 1980. On exact solutions of diffusion equation in cometary dynamics. Astron. Astrophys. 85, 77–79.