A&A 516, A72 (2010)

Oort cloud formation at various Galactic distances

R. Brasser¹ - A. Higuchi² - N. Kaib³

1 - Dept. Cassiopée, University of Nice - Sophia Antipolis, CNRS, Observatoire de la Côte d'Azur, Nice, France
2 - Dept. of Earth and Planetary Science, Tokyo Institute of Technology, Ookayama, Meguru-ku, Tokyo, Japan
3 - Astronomy Department, University of Washington, Seattle, WA, USA

Received 17 February 2010 / Accepted 17 March 2010

Abstract
In this study we present the results from numerical simulations of the formation of the Oort comet cloud where we positioned the Sun in various parts of the disc of the Galaxy, starting at 2 kpc up to 20 kpc from the Galactic centre. All simulations were run for 4 Gyr. We report that the final trapping efficiency of comets in the Oort cloud is approximately 4% and is almost independent of the solar distance from the Galactic centre. This efficiency is not enough to explain the flux of long-period comets and be consistent with the mass of the protoplanetary disc. In addition, the population ratio between the Oort clouds and their corresponding scattered discs is at least two orders of magnitude lower than observed today. Solar migration from the inner regions of the Galaxy - where the initial trapping efficiency is higher - to farther regions - where the retention is higher - does not add enough of an effect to increase the efficiency to the necessary value of approximately 20%, so that other mechanisms for forming the Oort cloud need to be investigated.

Key words: comets: general ; Oort cloud

1 Introduction

Ever since Oort (1950) published his results on the orbital energy distribution of long-period and new comets nearly sixty years ago, astronomers have studied the formation and evolution of this quasi-spherical cloud of comets surrounding the Sun. Most studies to date have concentrated on its formation and evolution in our current Solar System and with the Sun at its current distance from the Galactic centre on a circular, coplanar orbit. One such study was performed by Dones et al. (2004) - based on the model of Levison et al. (2001) - to simulate the formation of the Oort cloud in the current Galactic environment, starting with test particles on nearly circular, nearly planar orbits with the giant planets on their current orbits. After 4Gyr their results have shown that only 2.5% of comets ended up in the outer cloud (semi-major axis a >20 000 AU; Hills 1981), down from about 5% around 1Gyr, while the overall efficiency of depositing material in the cloud after 4Gyr is about 5% (down from 7.5% at 1Gyr). The inner and outer clouds contain roughly equal mass. Kaib & Quinn (2008) have performed similar simulations to those of Dones et al. (2004) and report only 1% of their comets winding up in the outer Oort cloud after 4.5 Gyr. These low efficiencies are difficult to reconcile with the estimated mass in solids in the primordial solar nebula and planet migration.

The problems associated with the above studies were highlighted by the results of Kaib & Quinn (2009), who conclude that most long-period comets (LPCs) could originate in the inner Oort cloud rather than the outer Oort cloud because of a newly-discovered dynamical path connecting the two. In the most extreme case, their results indicate that the outer Oort cloud could be devoid of comets and that all the LPCs originate in the inner cloud. This result requires a new approach to Oort cloud formation and evolution because the aim has shifted from depositing most objects in the outer Oort cloud to building a more substantial inner Oort cloud.

Depositing these comets is a problem in itself. Kaib & Quinn (2009) claim that the Oort cloud contains approximately  kilometre-sized comets, although this number is fairly uncertain because it is based on the observed LPC flux (Francis 2005); the number of comets in the cloud probably varies between . Since it is likely that the Oort cloud formed in two stages (Brasser 2008) and that there were approximately 10¹² kilometre-sized bodies available to deposit in the cloud during the second stage (Morbidelli et al. 2009), the efficiency of depositing material in the cloud in the current environment of the Sun is almost an order of magnitude too low if the cloud was generated solely during the second stage.

Star clusters are able to increase the trapping efficiency of comets in the inner Oort cloud during the first stage (Brasser et al. 2006; Kaib & Quin 2008) but then aerodynamic gas drag prevents kilometre-sized bodies from entering the cloud (Brasser et al. 2007) so that in the most extreme case the first stage did not make any contribution to the cloud. Even if gas drag is not a problem, the star clusters that are able to produce 90 377 Sedna have an Oort cloud that is devoid of comets farther than 5000 AU from the Sun, and the mechanism of Kaib & Quinn (2009) does not work at shorter distances. The transfer of comets with a<5000 AU to orbits with a>5000 AU is not efficient enough to supply the number of comets of Kaib & Quinn (2009) unless the mass of the disc between 5 and 15 AU during the first stage exceeded 150 (Brasser 2008). In summary, the dilemma of reconciling the low trapping efficiency with the flux of LPCs and the existence of Sedna still needs to be solved.

Here we make an attempt to resolve some of these issues by changing one parameter that has been kept constant in the previously mentioned studies: the strength of the Galactic tides. This is accomplished by varying the location of the Sun in the Galaxy, and thus corresponds to a change in the local density. This variation in distance is justified by evidence that the Sun migrated because of encounters with spiral arms, and it could have been scattered in and out many times yet remained on a nearly circular orbit (Sellwood & Binney 2002). The same study shows that the current orbit of the Sun is consistent with an origin between 4 kpc and 12 kpc from the centre of the Galaxy. In all likelihood the Sun's birthplace cannot be traced because of its chaotic dynamics within the Galactic disc, but perhaps some information can be gained from the Oort cloud. In addition, we can investigate whether the problem of trapping efficiency and LPC flux can be solved if the Sun had resided in different parts of the Galaxy.

One question that demands an answer up front is how this study is different from those in Brasser et al. (2006) and Kaib & Quinn (2008), who studied the formation of the Oort cloud while the Sun was in a star cluster, because the latter studies were able to deposit a much higher percentage of comets in the Oort cloud than the standard model of Dones et al. (2004). The problem with both of these studies was highlighted by Brasser et al. (2007), who showed that aerodynamic gas drag from the solar nebula prevented small comets from reaching the Oort cloud, so that it is possible that all the comets must have arrived later. This second phase of Oort cloud formation happened approximately 600 Myr after the cluster stage during what is called the late heavy bombardment - a global instability in the outer solar system that triggered an intense bombardment of the terrestrial planets (Tsiganis et al. 2005; Gomes et al. 2005) - and if the Sun was in its current Galactic location, we are back to square one. However, if the Sun was in a different part of the Galaxy, preferably closer to the Galactic centre than it is today, then the possibility exists of trapping more comets in the cloud. This scenario is what we aim to investigate here.

The outline for this paper is as follows. In Sect. 2 theoretical considerations are discussed and some predictions are made regarding the inner edge and extent of Oort clouds surrounding solar type stars in other parts of the Galaxy; trapping efficiencies and erosion timescales are also discussed. In Sect. 3 we describe our numerical methods. Section 4 shows the results of the numerical simulations. In Sect. 5 we discuss some of the results, and in the final section conclusions are drawn.

2 Theoretical considerations

The most widely used model for computing the force on a comet caused by the Galactic tidal field is the one proposed by Heisler & Tremaine (1986). In a rotating co-ordinate system where the -axis points away from the Galactic centre, the -axis is towards the south Galactic pole, and the -axis is perpendicular to these according to the right-hand rule (Levison et al. 2001), the force is given by

where A and B are the Oort constants, G is the gravitational constant, and the density of material in the Galaxy in the neighbourhood of the star. Equation (1) can be rewritten as (Levison et al. 2001)

where is the angular speed of Galactic rotation measured at a distance R from the Galactic centre, and (Matese & Whitmire 1996) is a measure of the gradient of the Galactic rotation curve. Levison et al. (2001) argue that is a small parameter and that it can be ignored. They also chose A=-B = 13 km s^-1 kpc^-1, hence  km s^-1 kpc^-1, the same as Heisler & Tremaine (1986). For the purpose of this paper, we usually have so we keep the parameter .

In the inertial frame centred on the star, (x,y,z), the acceleration a comet experiences from Galactic tides becomes

At the solar orbit, R = 8 kpc,  km s^-1 kpc^-1, and  km s^-1 kpc^-1 (Levison et al. 2001), with   pc^-3 (Holmberg & Flynn 2000); therefore the z-component of the galactic tide is approximately an order of magnitude stronger than the radial components in the solar neighbourhood, which is the reason the radial components have been ignored in many of the previous studies (e.g. Heisler & Tremaine 1986; Matese & Whitmire 1996). Before getting back to the influence of the tides on the comets, we first discuss the potential of the Galaxy and how the Oort constants and local matter density are derived from it.

2.1 Galactic potential and tidal field

To study the influence of the Galactic tides on the Oort cloud at various locations in the Galaxy, a model is needed for the potential of the Galaxy from which one can derive the Oort constants A and B and the matter density , while other parameters are all derived from these three. We have chosen the model for the Galactic potential described by Flynn et al. (1996), with some of the parameters modified to coincide with the most current value (2.66 kpc) of the scale length of the density in the disc (Chris Flynn, 2008; personal communication). The equations employed in Hernquist (1993) are equally viable, but they involve the exponential integral function, , for the potential of the halo, while for the disc and bulge the potential cannot be written in closed form. We therefore chose not to adopt his equations at this time.

The potential consists of a component due to the bulge/stellar halo and inner core , the disc , and the dark halo . The core potential consists of two Plummer spheres, the disc is made up of of three Miyamoto-Nagai components (Miyamoto & Nagai 1975), and the dark halo has a spherical, logarithmic form. From Flynn et al. (1996),

where and are the core radii and masses of the bulge/stellar halo term and inner core term, is the core radius of the halo, and is the velocity at large radius relative to ; R²=x²+y² and are the mass components of the disc; and a_n and b are related to the scale length and scale height of the disc, which are respectively the distance in the midplane and vertical direction over which the mass density falls off by a factor of e. By writing , one has  . An overview of all the constants is given in Table 1.

Table 1:   Parameters used in the potential of the Galaxy.

From the relation between the potential and matter density, we find for the above potentials (Binney & Tremaine 1987)

The density in the midplane is obtained by evaluating the above equations with z=0. At R=8 kpc,   pc^-3, close to the current value of 0.1  pc^-3 (Holmberg & Flynn 2000). The surface density in the disc at the same distance is   pc^-2, corresponding to a scale height of b=283 pc (we used b=300pc in Table 1). The best fit for the midplane density in the range  kpc is an exponential with a scale length of 3.29 kpc:   pc^-3. The rotational velocity is given by and is 229 km s^-1 at R=8 kpc. Figure 1 displays the rotation velocity in the midplane as a function of Galactic distance. The solid line is the total speed, the dashed line the contribution due to the bulge and inner core, the dashed-dot line is the contribution by the disc, and the dotted line accounts for the contribution by the dark halo. Similarly, Fig. 2 plots the midplane density as a function of Galactic distance. The solid line represents the total density, the dashed line the contribution by the bulge and inner core, the dashed-dot line from the disc and the dotted line from the dark halo.

The Oort constants A and B are related to the potential through (Binney & Tremaine 1987)

The Oort constants, as well as and |A+B|, are plotted in Fig. 3. The top panel shows A (solid line) and |B|(dashed line). The middle panel is and the lowest panel shows |A+B|.

From all these relations, one can compute the strength of the vertical tide compared to the planar tides, T_Z, i.e. the ratio of vs. and , which are obtained from Eq. (3). The value of T_Z is calculated by averaging the x and y tides over one Galactic rotation and making an approximation that x, y, and z have comparable magnitudes. We have

At 8 kpc, T_Z=11.28. If were zero, then T_Z=10.04. Heisler & Tremaine (1986) find T_z=15, but they use  pc^-3, which is almost twice as large as ours. Figure 4 plots T_Z as a function of Galactic distance for both symmetric and asymmetric tides ( and nonzero, respectively).

Now that we have a model for computing the tidal forces, we now turn to the passing Galactic field stars in the next section.

2.2 Passing stars

The issue of computing the action of passing stars on the Oort cloud contains several complexities, such as the correlation between age and velocity dispersion and the abundance of each spectral type as a function of Galactic distance. The latter is difficult to determine, while some attempts have been made to determine the first.

Lewis & Freeman (1989) have observed the velocity dispersion of the old disc as a function of the distance in the plane from the Galactic centre. They were able to fit an exponential curve through their data for the radial and tangential velocity dispersions ( and , respectively) with two different kinematic scale lengths. The best fit through their data is

While this gives a 1D velocity dispersion at R=8 kpc of 24.6 km s^-1, similar to what is used by Heisler et al. (1987) and found by Garcia-Sanchez et al. (2001) for late spectral types, the 1D dispersion at 20 kpc is only  km s^-1, which seems very low and is much lower than found numerically by Rodionov & Orlov (2008). Since the Galactic rotation curve is most likely flat between 9 and 15 kpc (Demers & Battinelli 2007), the model specified in Eqs. (4) seems a good representation of the reality. Since these models yield an exponential decay of the matter density in the disc, there is no a priori reason to distrust the exponential fits given in Eqs. (8) above.

However, as argued by Evans & Collett (1993), the ratio of for the fit Lewis & Freeman (1989) used is incorrect. Instead, Evans & Collett (1993) suggest a fit , with b = 4.5 kpc and  km s^-1. Unfortunately, Evans & Collett (1993) do not give a fit for . Attempting a similar fit to the data by Lewis & Freeman (1989) for yields poor results. The best fit we could find was a simple power law with index -1.11 for a broad range of Galactic distance. Bienaymé & Séchaud (1997) obtain a very long kinematic scale length in the solar neighbourhood, usually longer than 20 kpc, which they attribute to them having used more blue stars, whose kinematics are similar to that of gas, implying that the scale length in these regions is not well known. As a result, a power-law fit is probably okay. In the top two panels of Fig. 5 we plot and of Lewis & Freeman (1989) with the best fits mentioned earlier.

If we assume that the above fits are representative of the reality in the range  kpc, we can compute the 1D isotropic velocity dispersion as . Taking , which is a reasonable assumption, and calculating the resulting values of in the range we are interested in yields a best fit of  km s^-1. The value of obtained this way from the data of Lewis & Freeman (1989) and its best fit are depicted in the bottom panel of Fig. 5, and yields  km s^-1 at 20 kpc, in reasonable agreement with the results of Rodionov & Orlov (2008).

The way the passing stars are implemented in our simulations follows. We assume that the above change in velocity dispersion as a function of Galactic distance is the same for all spectral types. In other words, we assume that the relative increase or decrease in velocity dispersion with Galactic distance is the same across all spectral types. For example, the local 1D velocity dispersion is about 25 km s^-1 for M-dwarfs (Garcia-Sanchez et al. 2001). At 14 kpc this would be  km s^-1. For B-stars the velocity dispersion is about 8.5 km s^-1 at 8 kpc, which becomes about 5 km s^-1 at 14 kpc by this method.

Now that we have an idea of how to model the passing Galactic field stars, in the next section we look at some relevant time scales that naturally enter into the problem.

2.3 Time scales and relevant quantities

The tidal torquing time of the comet's pericentre due to the Galactic tide is given in Duncan et al. (1987) as . Since the planar tides are not always small compared to the vertical tide, we cannot compute using only the vertical tide. From Fouchard et al. (2005), who adds the contribution of the radial tides, we compute t_q, in Myr, as

where we used and and

Here is the correction for including the radial tides. The assumptions on and are fair since the comets tend to cluster around and an inclination of the planetary system's plane with respect to the Galactic plane of 60 is typical as well. Only when does one need to take its value into account. However, the chances of finding a planetary system with such an inclination is rare if one makes the fair assumption that the systems are isotropically distributed. It turns out that for  kpc, which implies that the planar tides do not contribute much to lifting the comets away from scattering planets, even though their strength is close to the vertical tide for R > 16 kpc. In what follows, we use t_q to estimate where the inner edge of the Oort cloud is in the limiting cases where the planets are very heavy or very light.

For heavy planets the typical energy change a comet receives per pericentre passage is much greater than the binding energy of a comet in the cloud, so that only a low percentage of comets scattered by such planets will land in the cloud. The semi-major axis at which this becomes possible is computed from setting the tidal torquing time equal to the orbital period (Duncan et al. 1987). The corresponding minimum energy to which comets will be launched based on this time scale, using  AU and q=(9/8)a_p (Brasser & Duncan 2008), is

For the potential that we used, a best fit through yields the following relation for the minimum semi-major axis at which the tides can lift the comets into the cloud:

where R is measured in kpc and the fit has errors less than 3%. Using a_p = 5 AU at 2 kpc, the result is AU, at 8 kpc we have , AU and AU at 20 kpc.

For light planets, the comets diffuse towards the Oort cloud rather than being launched and are lifted at a semi-major axis corresponding to where the diffusion time equals the tidal torquing time i.e. t_q=t_d (Duncan et al. 1987). Here the diffusion time t_d is given by Duncan et al. (1987) as

where is the mass of the scattering planet, and a number measuring the root-mean-square energy change a comet experiences per pericentre passage. The latter is defined through . Duncan et al. (1987), Tremaine (1993), and others traditionally used , which is correct for most purposes and for ecliptic inclinations less than 30. Once again, when setting AU and taking q = (9/8) a_p we have

Repeating the same exercise as before, the best fit to the factor yields

Using Neptune with , a_p = 30 AU, and , we have AU at 2 kpc, AU at 8 kpc, and AU at 20 kpc (which is probably not correct because the launching formula with the same mass and semi-major axis gives AU). The values quoted just now are correct to within factors of two because of the value of . As can be seen, the semi-major axis at which comets get lifted into the cloud in the diffusion case is a much steeper function of R than in the launching case. In fact, the comets will diffuse towards the cloud if (Brasser & Duncan 2008)

At R=20 kpc, for a_p=30 AU, while for a_p=5 AU. At 2 kpc, setting a_p=30 AU yields and at a_p=5 AU. Equation (16) should be compared with the condition that the halflife of the comets scattered by a planet of mass be less than the age of the stellar system, t_*. This is a slightly stronger criterion than the condition t_d < t_* used in Tremaine (1993) and Brasser & Duncan (2008). Since (Brasser & Duncan 2008), t_1/2 < t_* yields

Taking t_* = 5 Gyr, a_p = 30 AU, and , Eq. (17) yields , so that if one wants to build the Oort cloud efficiently at R=20 kpc, a planet mass comparable to Uranus ( ) and Neptune ( ) is needed at a distance of 30 AU or more, and it takes approximately 5 Gyr for half the comets to be eliminated from the planet's region (assuming there is just the one planet; with multiple planets the process will be much faster).

Now that we know the time scales involved for massive and light planets to scatter the comets and in which part of the cloud these comets are most likely to end up, we inspect a few other time scales that naturally enter the problem. The first one is the half life for a comet in the cloud subjected to stellar perturbations, which is a measure of how long the comets can survive in the cloud, of how much the cloud wears down over time, and by proxy the total trapping efficiency of comets in the cloud. The second is the Kozai cycle time, which is a measure of how long the comet can stay in the cloud, before its pericentre reaches the planets again.

Bailey (1986), following Hut & Tremaine (1985), computes the half life for a comet in the Sun's Oort cloud subjected to perturbations of passing stars, which is given by (Brasser et al. 2008)

From Garcia-Sanchez et al. (2001), we have  pc^-3. Assuming the same scale length applies to the density of stars than the rest of the Galaxy, then pc^-3. Also from Garcia-Sanchez et al. (2001),

Using previous fitting relations and Bailey's (1986) formula, we get

which, as one can see, increases exponentially with R for fixed a. The above result is a bit misleading since the behaviour of a with R is not taken into account.

There is another times cale in this problem: the Kozai period of the comets, induced by the perturbations of the Galactic tide. The Kozai time measures how long it takes for the argument of pericentre of the comet to revolve one cycle (or, if the latter is librating, to perform one oscillation), and because the oscillations of the eccentricity and inclination of the comet are coupled to the evolution of the argument of pericentre (Kozai 1962), it is also the time it takes for the pericentre of the comet to perform one cycle and re-enter the planetary region. This time scale is a measure of how long a comet can be trapped in the cloud before it reaches the planetary region once again, provided it has not been stripped or has its trajectory reset by a passing star (Weissman 1980), and the stripping occurs on a time scale t_1/2^B. Once the pericentre of the comet re-enters the planetary region, it can be ejected by the planets.

Formally the Kozai period only exists when the radial tides do not, because then the problem is integrable. However, if the strength of radial tides is smaller than the vertical tide, then the Kozai period can be computed approximately if the semi-major axis of the comet's orbit is not too large. Following Higuchi et al. (2007), the Kozai period becomes

where n is the mean motion of the comet, are constants that depend on the orbit's energy and angular momentum, and is Legendre's elliptic integral of the first kind. For nearly-radial orbits, one can expand in powers of q/a. Doing this for or yields

while for or 270

Averaging between both these extremes, the period over the range becomes

Now that we have an understanding of the time scales involved and the parts of the cloud the comets end up as a function of the position and mass of the planet, we are ready to discuss the numerical simulations.

3 Method

We used a new code to perform the numerical simulations called SCATR (Symplectically-Corrected Adaptive Timestepping Routine) (Kaib et al. 2010), which is based on SWIFT's RMVS3 (Levison & Duncan 1994). SCATR is a modified RMVS3 to include passing stars and the influence of the Galactic tide (Levison et al. 2001; Dones et al. 2004). The major advantage of SCATR over RMVS3 is that it is able to use much longer time steps at large distances from the Sun.

For our simulations we used initial conditions similar to Dones et al. (2004) and Kaib & Quinn (2008). The giant planets are on their current orbits, after their instability at the time of the late heavy bombardment (Tsiganis et al. 2005; Gomes et al. 2005), and we used 2500 test particles to represent the comets, originally situated between 4 and 32 AU. We placed the Sun on orbits at 2, 4, 6, 8, 14 and 20 kpc from the Galactic centre, and the simulations were run for 4 Gyr. The criteria for removing comets were as follows: either they were hyperbolic or, if bound, farther than a certain distance from the Sun, which is computed below. A comet is deemed to be in the Oort cloud when its pericentre distance q>40 AU, its semi-major axis a>1000 AU, and its eccentricity e<1. All simulations for this study were run on the CRIMSON Beowulf cluster at OCA.

Table 2:   Stellar data from Garcia-Sanchez et al. (2001).

The method of generating passing stars is based on a mixture of methods by Heisler et al. (1987) and by Rickman et al. (2008). Table 2, which is taken from Garcia-Sanchez et al. (2001), shows the values of the mass m_i, number density n_i, and velocity dispersion for a large class of main-sequence stars in the solar neighbourhood. The velocity dispersion is taken to be the 1-dimensional isotropic value and is obtained by dividing the values in Garcia-Sanchez et al. (2001) by . The outline follows.

• Calculate the time of the mth stellar encounter t_m, according to the method of Heisler et al. (1987). The probability that the the next encounter happens at time t is given by , where f is the mean encounter frequency of stars with impact parameters relative to the Sun shorter than some chosen maximum impact parameter . From Heisler et al. (1987) we have

  (24)

  
  
  where the sum is taken over stellar spectral types. Since both and n_i change with Galactic distance, R, the value of needs to be adjusted to ensure that the passing stars are entered much farther away than the typical maximum apocentre of a comet's orbit (2 pc at the Sun's current location).
• Next we need to know the spectral type of the star. Again we use the method outlined in Heisler et al. (1987). Define and for any . We then choose a uniform random number between 0 and 1, . The star then has type i if . Rickman et al. (2008) base the star's type on the encounter frequency of Garcia-Sanchez et al. (2001). While this method is as good as that of Heisler et al. (1987) we do not know the encounter frequencies at other Galactic distances, while we are able to compute estimates of n_i and (and therefore ) as outlined earlier.
• Next we choose the star's entry point on a sphere of radius . The position vector is generated as

  =

  (25)

  
  
  where and are chosen randomly with uniform distribution in the intervals and .
• The next step is to choose the magnitude of the velocity. For this we follow the outline of Rickman et al. (2008). The peculiar velocity ( ) of a star with respect to its local standard of rest is calculated from a spherical Maxwellian distribution as follows. We compute three random numbers that have a Gaussian probability distribution with a mean of 0 and standard deviation of 1. We use the formula

  (26)

  
  
  with l=5 where are random numbers between 0 and 1. The peculiar velocity of the star is then . The magnitude of the star's heliocentric velocity is subsequently calculated by combining the vector with the Sun's apex velocity . It should be noted that values for only exist for the solar neighbourhood (Table 2). We return to this issue shortly. We assume the two vectors are randomly oriented so that we have

  (27)

  
  
  where is uniformly chosen between -1 and 1. Since the flux is dependent on velocity, the latter velocity needs to be modified. Compute a high velocity for each spectral type, called V_0i. Generate a random number . Then, if , keep the value of V_i, otherwise recompute it according to the above until the condition for is fulfilled. We chose .
• Last, we need to know the components of the velocity vector of the star. Define v₃ as the component anti-parallel to the radius vector and v₁ and v₂ as the components orthogonal to it. Then we have where is a random number between 0 and 1, and and , where . The components of the velocity are then
  
  
  where the same and that were used for the radius vector need to be inserted.

The above description ensures that the stellar impact parameter b is uniform in b². Using this method at the Sun's radius with the values of given in Table 1 yields an average velocity for all stars of  km s^-1 with a standard deviation  km s^-1, the same as in Rickman et al. (2008). However, as mentioned earlier, the value of is not known at other locations (apart from the Sun not physically being there). There are therefore three natural choices for : we could set it to zero, the velocity dispersion for G0 stars , or . The Sun's total velocity is comparable to the velocity dispersion of G0 stars , while a more typical G0 star has . The choice , though it was used by Heisler et al. (1987), seems artificial, so it was rejected. The choice seems more logical, but it would mean copying the Sun's motion, which is not typical. Therefore, we decided to set . For a G0 star at 8 kpc this increased the mean encounter speed to  km s^-1 with a dispersion  km s^-1, and we used like Garcia-Sanchez et al. (2001). The reader may note that the increased velocity of the stars from the method employed by Rickman et al. (2008) will result in weaker perturbations on the Oort clouds. However, the increased velocity of the stars also results in an increased flux (of about 30% at 8 kpc), and from numerical experiments at 8 kpc it was found that two effects more or less cancel each other out.

To estimate the value of , we proceed along the lines of Antonov & Latyshev (1972) by studying the zero-velocity curves around the Sun. According to Antonov & Latyshev (1972), the Jacobi constant for the system consisting of a comet orbiting the Sun in the co-ordinate system erected in Sect. 2 is

where C is the ``third Oort constant'', given by

and is a measure of the vertical frequency of the Sun's motion. After computing it, a best fit for  kpc is an exponential with scale length 6.2 kpc:  km s^-1 kpc^-1 with errors less than 5%. At 8 kpc it evaluates to C = 68.3 km s^-1 kpc^-1, in good agreement with current observations. According to Antonov & Latyshev (1972), the critical value of beyond which no bound orbits exist is . The lengths of the semi-axes of the tri-axial ellipsoid at the critical value of then is

where . We then set the size of the sphere around the Sun to . A plot of x_c, y_c, and z_c as a function of Galactic distance is shown in Fig. 6.

4 Results of numerical simulations

In this section the results of our numerical simulations are described. We are interested in the following issues: the internal structure of the cloud, the formation efficiency, the formation duration, and the extent and size of the cloud.

We first focus on the size and extent of the cloud. In Fig. 7 we plot the semi-major axis of the inner edge of the clouds as a function of Galactic distance. The line denotes the analytical result of Eq. (15) where we used a_p=25 AU,   , and . The bullets denote the value of semi-major axis inside of which 5% of the comets reside. The upper and lower marks of the error bars denote the same for 1% and 10%. As one can see, the analytical result follows the bullets rather well. This leads us to conclude that, generally, the analytic result can be used. The next figure, Fig. 8, plots the cumulative distribution in semi-major axis of Oort clouds at 2 (leftmost line), 8, 14, and 20 kpc (rightmost line). These distributions show the size and extent of the clouds. The distributions shift towards higher values of semi-major axis with increasing Galactic distance, as expected. Additionally, the extent of the cloud decreases with increasing distance. At 2 kpc the ratio between the inner and outer edges is approximately 30, while at 20 kpc this ratio is close to 8. Thus, closer to the Galactic centre, the comets are spread over a winder range of binding energies.

We now turn to the matter of the formation efficiency of the cloud. We depict this efficiency at various Galactic distances in Fig. 9. The top panel depicts the efficiency of depositing the comets in the Oort cloud after 4 Gyr, while the bottom panel shows the efficiency of depositing comets assuming that no subsequent decay takes place; i.e., all the comets that make it to the cloud are deposited there at the same time, and there is no subsequent erosion caused by the passing stars and the Galactic tide. The bottom panel should only be used as an indicator of a theoretical upper limit to the efficiency. It is surprising that the top panel shows small deviations in the overall efficiency as a function of Galactic distance. The difference is more pronounced in the lower panel but this is to be expected: the density is higher closer to the Galactic centre so that comets are trapped closer to the Sun and thus, in principle, have a higher efficiency. That the efficiency after 4 Gyr is more or less independent of Galactic distance implies that there is more wearing down of the cloud in the inner parts of the Galaxy, even though the comets are closer to the Sun. The low trapping efficiencies at short Galactic distance do not bode well for explaining the LPC flux. We return to this in the next section.

A good measure of the internal structure and distribution of the cloud is to calculate its density profile, n(a), which measures the volume density of comets with semi-major axis between a and , i.e., . For most clouds n(a)is a power law, i.e., . Figure 10 depicts a few such density profiles for clouds at various Galactic distances: the top-left panel is for the cloud at 2 kpc, the top-right panel is at 8 kpc, the bottom-left panel refers to 14 kpc, and the last panel refers to the cloud at 20 kpc. The line shows the best fit, all of which are power laws. The index ranges from -2.7 at 14 kpc to -3.3 for the outer parts of the cloud at 20 kpc. The inner part at 20 kpc can be fit by . Dones et al. (2004) has already reported that the inner parts of the Oort cloud have a flat density distribution, which is clearly visible in the last panel.

It is interesting to compare this result with earlier simulations. Dones et al. (2004) finds for the whole Oort cloud, though they also employ a two-stage power law that gives a steeper slope in the outer parts (-3.7). For the cloud at 8 kpc, we have , similar to what Dones et al. (2004) find. This is to be expected since the simulations presented here are identical to those of Dones et al. (2004), so it is encouraging that we get the same result, and it proves that the slightly different values of the Oort constants A and B, as well as the increased stellar encounter speeds, make little difference for the final result.

Another measure to identify the inner structure of the cloud is to plot the mean eccentricity and ecliptic inclination , of objects in certain parts of the cloud. These quantities have been depicted in Fig. 11, which shows the values of and at the 10th (bullets), 50th (open circles) and 90th (triangles) percentiles in a as a function of Galactic distance. The straight lines at and are values for a system that is isotropic and in virial equilibrium. The data was gathered as follows: comets that are in the cloud were sorted according to semi-major axis, which yields the cumulative distribution, F(a). For the data points corresponding to the 10th percentile, we averaged the eccentricity or inclination of comets with semi-major axes between F(0.05) and F(0.15). We performed similar methods for the other data points. The inner regions and around the median distance the clouds are always more eccentric than the virial value, while the comets are more circular when farther than the median distance. Regarding the inclinations, the inner parts are always prograde, similar to what Dones et al. (2004) and Kaib & Quinn (2008) find. Around the median, the clouds are as good as isotropic, but at large distances from the Sun the clouds are predominantly prograde yet again.

Some of the effects depicted in Fig. 11 are also visible when plotting the comets in a-q space at the end of the simulations. Figure 12 shows a sample of Oort clouds at the end of the simulations. The top-left panel depicts the cloud located at 2 kpc from the Galactic centre, the top-right panel is at 8 kpc, the bottom-left panel refers to 14 kpc, and the bottom-right panel refers to 20 kpc. The structures are different at various distances. Two effects are clearly visible. The first is that in each panel where the distance to the Galactic centre is farther than the previous one, the distance of the cloud to the Sun increases, as we showed earlier. The second effect pertains to the distribution in q: these values are systematically higher in proportion to the semi-major axis at the outer edges of the cloud, while the values of q span a much wider range around the median. Of course, near the inner edge, the range in q is limited because of lower values of the semi-major axis and because the comets have not had the time to complete one Kozai cycle.

One last question one could ask is the formation time of the cloud, especially when comparing the formation time with the orbital period of the Sun around the Galactic centre. One could envision that, if the formation time is comparable to the orbital period of the Sun, the cloud could have encountered several spiral arms, depending on how close the Sun orbits the co-rotation radius, and these encounters might have caused the Sun to migrate (Sellwood & Binney 2002). In Fig. 13 we plot the cumulative time for comets the reach the cloud, for clouds at 2, 8, 14, and 20 kpc (left to right, respectively). As one can see, the formation time increases with increasing Galactic distance, which is caused by the planets having to scatter the comets onto larger orbits before the tides can lift them into the cloud. The fraction of comets deposited into the cloud scales approximately logarithmic with time. The median formation time of the cloud ranges from 50 Myr at 2 kpc to 150 Myr at 20 kpc, which is shorter than or comparable to the orbital period of the Sun at most Galactic distances. At 8 kpc, over 70% of the comets that survive after 4 Gyr in the cloud are already in the cloud after one Galactic orbit so that during the formation process the Sun can probably be considered to be on a fixed orbit. It is only during the subsequent evolution of the cloud after the formation process that the migration of the Sun becomes important. We discuss this in the next section.

5 Discussion and implications

The results presented in the previous section are to be expected based on earlier work of Dones et al. (2004) and Kaib & Quinn (2008), as well as the theory presented in section 2. Nontheless, it contains a few surprises that we believe should be discussed. However, first there are a few other points that we want to clarify.

The first issue is that the formation efficiency is more or less the same for most Galactic distances: 3% 1%. The second issue deals with solar migration within the Galaxy, and needs further discussion. Sellwood & Binney (2002) have demonstrated that solar encounters with spiral arms can force the Sun to migrate significantly while keeping its orbit nearly circular. The current position of the Sun is consistent with a formation between 4 and 12 kpc. In the simulations presented here we kept the Sun at the same distance from the Galactic centre for 4 Gyr. Since it is likely that the Sun has migrated during that time, the results of these simulations should be viewed as a proof of concept rather than as indicative of the current Oort cloud and as representative of its population. A way that migration will affect the structure of the Oort cloud can be understood from the following example in which we assume that the Oort cloud formed in two stages as outlined in Brasser (2008). Suppose that during the time of the late heavy bombardment - a phase of intense bombardment of the terrestrial planets, most likely caused by a dynamical instability of the outer planets approximately 600 Myr after the formation of the giant planets (Tsiganis et al. 2005; Gomes et al. 2005) - the Sun was closer to the Galactic centre than it is today, at 4 kpc or so. During the time of the late heavy bombardment, a lot of material very quickly becomes available for the giant planets to scatter (Gomes et al. 2005). From the results presented in the previous section, it appears that most comets are in the Oort cloud after nearly a hundred million years and that there is subsequently an exponentially-decreasing flux of comets reaching the cloud until the end of the simulations. This can be understood from Fig. 13: the cumulative distribution is more or less linear as a function of , so that a histogram of number of comets deposited vs is more or less flat, thus indicating a logarithmic increase in the number of comets being deposited in the cloud. In the inner parts of the Galaxy the flux of stellar encounters and short Kozai cycle time cause the Oort cloud to erode very quickly. However, the majority of comets can be saved if the Sun were to migrate outwards shortly after the late heavy bombardment so that the retention of comets will be higher than closer to the Galactic centre. This would help the scenario of Kaib & Quinn (2009), because when the Sun is at 4 kpc most of the comets get deposited in the inner cloud, where they argue most of the LPCs could originate from. In other words, the potential exists for creating a more massive inner Oort cloud than in the scenario of Dones et al. (2004) and Kaib & Quinn (2008) if the circumstances are just lucky enough. Alternatively, if the Sun formed much farther out, all of the comets would have been deposited in the outer cloud and the problem of the flux of LPCs need not bother with a possible inner Oort cloud population. Of course, in that case the low formation efficiency makes it difficult to reconcile the LPC flux with the remaining mass of the primordial solar nebula, which became fully destabilised during the late heavy bombardment.

Unfortunately, it is very unlikely that solar migration will save the day if the following is correct. Kaib & Quinn (2009) argue that there are approximately  km-sized comets in the whole Oort cloud. Morbidelli et al. (2009) state that there are approximately 10¹² km-sized bodies in the trans-Neptunian disc before the late heavy bombardment. To get the required Oort cloud population, we need a trapping efficiency of approximately 20% after 4 Gyr, which is much higher than attained in any of the simulations presented in this work. One could argue that the first Oort cloud, which was formed when the Sun was still in its birth cluster, could supply the necessary flux of LPCs in the manner explained by Kaib & Quinn (2009). The formation efficiency of that first Oort cloud is approximately 10% (Brasser et al. 2006; Kaib & Quinn 2008), which might be just enough to explain the LPC flux. However, if 90377 Sedna formed by the same mechanism, then most of the comets in this first cloud are inside 5000 AU, where the Kaib & Quinn (2009) mechanism cannot operate. Additionally, the gas drag from the solar nebula prevents km-sized comets from being scattered in the early age of the Solar System (Brasser et al. 2007), so that the fist Oort cloud would only consist of large bodies or would be altogether empty. Indeed, no Sedna analogues have been found thus far (Schwamb et al. 2009), which begs the question of whether or not it is a fluke, i.e., a 2- or 3 event where an unlikely close stellar passage with the Sun occurred while it was still in its birth cluster. Last, the observed population ratio between the Oort cloud and the scattered disc (defined as comets having  AU and a > 40 AU) is approximately 1000:1 (Duncan & Levison 1997; Levison & Duncan 1997; Levison et al 2008), while simulations tend to yield a 5:1 ratio (Kaib & Quinn 2008) to 20:1 (Dones et al. 2004). We find a population ratio between 15-25:1 for our simulations presented here, with no clear correlation between this ratio and Galactic distance. Unfortunately it is much less than the observed value and is a severe shortcoming of these simulations (Dones et al. 2004; Kaib & Quinn 2008). Apparently changing the strength of the tides does not alter this ratio.

We should emphasise, though, that the above numbers are fraught with large uncertainties. The number of comets in the cloud is inferred from the LPC flux, which currently stands at approximately 2-3 comets per year with q<4 AU (Francis 2005). An error of 1 comet per year on either side increases the number of comets in the cloud appropriately, hence the necessary formation efficiency, with a lower LPC flux better supporting our models. Another uncertain parameter is the number of objects in the scattered disc. Even though dynamical models based on the flux of Jupiter-family comets yield approximately  km-sized bodies, the precise number depends on the size distribution. For the Kuiper belt, the best fit is a two-stage power law (Bernstein et al. 2004, Fraser & Kavelaars 2009; Fuentes & Holman 2008; Fuentes et al. 2009) with a differential size distribution where for objects with diameters  km, and for smaller objects. Taking and using the fact that there are approximately 10⁵ objects in the scattered disc with D > 50 km (Levison et al. 2008) gives us approximately  km-sized bodies in the scattered disc; this can be pushed to objects when using . The latter number is still two orders of magnitude lower than the inferred population of the Oort cloud. However, the parameter with the largest uncertainty is the conversion from absolute magnitude to diameter. For scattered disc objects, this involves the albedo but for comets the conversion is much poorer because the coma hides the nucleus. Therefore, we may be mislead thinking that the flux of LPCs of 2-3 comets per year is for km-sized objects while these bodies might, in fact, be a lot smaller. Only if we underestimate the sizes of Oort cloud comets by at least an order of magnitude do the two population ratios begin to agree, but this seems unlikely. Therefore, in light of the above problems, we tend to rule out the formation of the Oort cloud by this mechanism. The only alternative mechanism that can explain the high population ratios we can think of is the trapping of comets while the Sun was in its birth cluster, provided the other stars could scatter the comets out and not be hindered by gas drag. The trapping scenario is being explored (Levison et al. 2010).

6 Summary and conclusions

We performed numerical simulations of the formation and evolution of the Oort cloud with the Sun orbiting the Galactic centre at a distance ranging from 2 kpc to 20 kpc. We used an up-to-date model for the Galactic potential and computed the solar orbit, matter density and other parameters from it that are needed to calculate the Galactic perturbations on the comets in the cloud. We used the observational results of Lewis & Freeman (1989) to compute the 1D velocity dispersion of stars at various Galactic distances, from which we constructed a method of computing the action of these passing stars on the Oort cloud in our integrator called SCATR (Kaib et al. 2010) based on methods by Heisler et al. (1987) and Rickman et al. (2008).

We have demonstrated that the Oort cloud can readily form in the other parts of the Galaxy. For the clouds that were generated, we established that the inner edge of the cloud - defined as the 5% mark in the cumulative semi-major axis distribution - fits approximately with the theory of Eq. (15). The clouds tend to have a flattened, very eccentric inner part, which is caused by the tides and stars torquing the comets very slowly, and a mostly isotropic middle and outer part, with mean eccentricity fairly close to the virial value of 2^-1/2 and mean ecliptic inclination close to 90. However, the eccentricities close to the outer edge of the cloud tend to be below the virial value, and the ecliptic inclinations favour prograde motion.

The efficiency of trapping comets in the cloud after 4 Gyr ranges from 2% at 2 kpc to 4.5% at 8 kpc and is more or less independent of the distance to the Galactic centre. The difference in efficiency between a cloud at 2 kpc and one at 8 kpc is caused by the stronger erosion of the cloud closer to the Galactic centre, owing to the higher density of stars in these regions and the shorter Kozai cycle. Therefore, the efficiency cannot be used to trace a possible migration of the Sun in the Galaxy. We conclude that this efficiency is not high enough to explain the current flux of LPCs, even if the Sun was originally closer to the Galactic centre and subsequently migrated outwards. In addition, the population ratio between the Oort cloud and the scattered disc is lower than the observed value by approximately two orders of magnitude. Therefore, assuming that gas drag prevents the formation of an Oort cloud during the cluster stage, the only viable alternative to explaining the number of comets in the Oort cloud is through the trapping of interstellar comets while the Sun was in its putative birth cluster.

Acknowledgements

We wish to thank Chris Flynn for very helpful discussions of the Galactic potential and the stellar velocity dispersion during the early stages of this study. We also thank Alessandro Morbidelli and an anonymous reviewer for helpful comments. R.B. thanks Germany's Helmholtz Alliance for financial support through their `Planetary evolution and Life' programme. A.H. thanks the Japan Society for the Promotion of Science for supporting this work. N.A.K. is financed by a NASA Earth and Space Science Fellowship.

References

All Tables

Table 1:   Parameters used in the potential of the Galaxy.

Table 2:   Stellar data from Garcia-Sanchez et al. (2001).

All Figures

---

Copyright ESO 2010