Scattering of Planetesimals by a Planet: Formation of Comet Cloud Candidates

We have investigated the first dynamical stage of comet cloud formation, the scattering of planetesimals by a planet. The orbits of planetesimals were calculated using circular restricted three-body formalism. We obtained the probabilities of the following results of scattering as functions of the orbital parameters of the planets and planetesimals: (1) collision with the planet, (2) escape from the planetary system, and (3) candidacy as a member of the comet cloud (planetesimals with large semimajor axes). We also derived simple empirical formulae for these probabilities that are accurate enough for order-of-magnitude estimation. We found that a planetesimal with an initial eccentricity of e ≳ 0.4 can escape from the planetary system or be a candidate for an element of the comet cloud due to scattering by a planet. As the energy range of the comet cloud is narrow, the probability of any planet producing escapers is always much higher than that of producing candidates. Using the probabilities and assuming a distribution of planetesimals, we obtained the efficiencies of collision, escape, and candidacy for a given planet. We applied the results to the solar system and found that, among the four giant planets, Jupiter is the planet most responsible for producing candidate elements of the Oort Cloud, as long as the inclination of planetesimals is constant or proportional to the reduced Hill radius of each planet.


INTRODUCTION
The Oort Cloud is a spherical comet reservoir surrounding the solar system. It is generally accepted that it consists of more than10 12 comets, and its size is on the order of 10 4 -10 5 AU (e.g., Weissman 1990; Dones et al. 2004). The existence of this comet reservoir was first proposed by Oort (1950). He suggested that planets scattered small bodies in the planetary region outward and that passing stars raised their perihelia out of the planetary region. There is general agreement that the Oort Cloud comets are the residual planetesimals of planet formation.
The standard scenario of the Oort Cloud formation consists of two dynamical stages: (1) giant planets raising the aphelia of planetesimals to the outer region of the solar system and (2) the Galactic tide, passing stars, and giant molecular clouds pulling their perihelia out of the planetary region. The first dynamical stage has been studied analytically by several authors (e.g., Safronov 1972;Weidenschilling 1975;Fernández 1978). Safronov (1972) estimated the ejection rate of planetesimals by the four giant planets and the relative importance of the planets for ejection and comet cloud formation. Using a Monte Carlo, Ö pik-type code, Fernández (1978) calculated the probability of a planetesimal colliding with a planet, being ejected from the solar system, and having a near-parabolic orbit with close encounters with planets. Both Safronov (1972) and Fernández (1978) concluded that Jupiter and Saturn hardly contribute to the formation of the Oort Cloud, since their ejection rates are too high because of their large masses, and that Uranus and Neptune may have had important roles in the Oort Cloud formation.
The first direct numerical calculation of the overall formation of the Oort Cloud was performed by Duncan et al. (1987). Their calculation included the gravity of the four giant planets and passing stars, and the Galactic tide. They concluded that the density profile between 3000 and 50,000 AU is roughly a power law proportional to r À3.5 , where r is the heliocentric distance, and that the inner Oort Cloud (semimajor axis <20,000 AU) contains roughly 5 times as many comets as the outer Oort Cloud. Dones et al. (2006) performed a similar calculation but with small initial eccentricities of planetesimals, which may be more realistic initial conditions. Their calculations showed that planetesimals typical of those that form the Oort Cloud are originally from the Uranus-Neptune region and are given the last scattering by Saturn. Both papers deal only with the solar system and do not focus on the general properties of planet-planetesimal scattering. Therefore, their results cannot be applied to other planetary systems. Tremaine (1993) first considered the formation of comet clouds in other planetary systems. Using the results of previous studies (e.g., Duncan et al. 1987;Heisler & Tremaine 1986), he derived a condition for comet cloud formation as a function of the mass (mp) and semimajor axis of a planet (a p ). However, this condition is based on the previous calculations, which are applicable only to the solar system.
In order to construct a general theory of comet cloud formation applicable to general planetary systems, it is necessary to clarify the elementary processes with parameters other than those for the solar system. In this paper we investigate the first dynamical stage of comet cloud formation using numerical calculations. In the late stage of planet formation, the scattering of planetesimals by a planet results in one of four events: (1) collision with the planet, (2) escape from the planetary system, (3) survival in the planetary system, and (4) a fall onto the central star. We call the planetesimals with fates 1, 2, and 4 ''colliders,'' ''escapers,'' and ''fallers,'' respectively. The fate 3 planetesimals with large semimajor axes, for example, larger than 3000 AU, are candidates to be elements of the comet cloud (hereafter ''candidates''). We calculate the probabilities of producing colliders, escapers, and candidates by orbital integration using circular restricted three-body formalism. We integrate the orbits of the planetesimals for one Kepler period to evaluate the probability for one encounter. The ratios of these fates depend on the parameters of the planetary system, i.e., the orbital elements of the planetesimals and the orbital elements and mass of the relevant planet. We investigate the dependence of the ratios on these parameters. We also present simple and useful empirical fitting formulae for these probabilities.
The outline of this paper is as follows: We describe the numerical method, model, and initial conditions in x 2. In x 3 we present the results for the probabilities of colliders, escapers, and candidates. We derive empirical fitting formulae of the probabilities in x 4. In x 5, using the probabilities, we estimate the efficiencies of planets for producing colliders, escapers, and candidates and apply the results to the solar system. Section 6 is devoted to a summary and discussion.

Model and Initial Conditions
We consider a star-planet-planetesimal three-body problem in which a planetesimal is treated as a massless particle. We consider planets of a p ¼ 1, 5, 10, and 30 AU and m p ¼ 0:1m J , 0.3m J , 1m J , 3m J , and 10m J , where m J 0:001m Ã and m Ã is the mass of the star, which is set as m Ã ¼ 1 M . The orbit of the planet is assumed to be circular. We set the densities of the planet and star as 1 g cm À3 .
We set up an initial planetesimal disk that consists of planetesimals and a planet. All the planetesimals have the same eccentricity e and inclination i. The range of e is from 0.1 to 0.9 with an interval of 0.1, and the range of i is from 0 to 0.1 rad with an interval of 0.01. The semimajor axes of the planetesimals a are uniformly distributed over the planetesimal disk; in other words, the surface number density of planetesimals is proportional to a À1 . We call the region of a where the coplanar orbits of the planet and a planetesimal cross the ''orbit-crossing region.'' The inner and outer edges of the orbit-crossing region for unperturbed orbits are a min ¼ a p /(1 þ e) and a max ¼ a p /(1 À e), respectively, where a p is the semimajor axis of the planet. We distribute planetesimals in a range wider than the orbit-crossing region, since gravitational focusing of the planet is effective. The argument of the perihelion !, the longitude of the ascending node , and the mean anomaly l of the planetesimals are distributed randomly. In the disk, a ring in a with a width of 1 AU contains 10 6 or 10 7 planetesimals. The initial disk parameters are summarized in Table 1.

Integration Method
The orbits of the planetesimals are integrated numerically using the fourth-order Hermite scheme (Makino & Aarseth 1992) with a hierarchical time step (Makino 1991). The equation of motion for a planetesimal is where G is the gravitational constant, r p is the heliocentric position of the planet, and r is the heliocentric position of the planetesimal. The last term of the right-hand side represents the indirect term. We calculate the orbits of all the planetesimals for 1 Kepler period (T K ). During the orbit integration, if the separation between the planet and a planetesimal becomes smaller than the radius of the planet R p or the heliocentric distance of the planetesimal becomes smaller than the radius of the central star R Ã , the planetesimal is counted as a collider in the former and a faller in the latter. Orbital elements of planetesimals except colliders and fallers are checked after T K . If the perihelion distance of the planetesimal is smaller than R Ã , it is also counted as a faller. If the eccentricity of the planetesimal is larger than 1, it is counted as an escaper. A planetesimal with a > a can is counted as a candidate, where a can is the minimum semimajor axis of a candidate for inclusion in the comet cloud. If the separation between the planet and a planetesimal after T K is smaller than the Hill radius of the planet, r H ¼ a p (m p /3m Ã ) 1 = 3 , we discard the planetesimal. This is just because its orbit changes a great deal due to strong interaction with the planet.

Definitions of Probability and Efficiency
We denote the probabilities of producing colliders, escapers, and candidates per T K as P col , P esc , and P can . Using these probabilities P, we define the efficiencies K of a planet. Efficiencies K col , K esc , and K can represent the expected numbers of colliders, escapers, and candidates per unit time and are defined as where n s is the surface number density of planetesimals, given as n s ¼ n 1 a , where n 1 is the reference surface number density  at a ¼ 1 AU, is the power-law index of the radial distribution, and a in and a out are the inner and outer edges of the planetesimal disk. In the present paper we do not show or discuss the probability and efficiency of becoming a faller.

PROBABILITIES
We numerically evaluate P and investigate its dependences on the parameters of the planetesimals and planets: a, e, i, a p , and m p . We show the results in order of P col , P esc , and P can . We refer to ''case 0'' as the standard case. For each P we first present the results of the standard case and the dependences on a and e. Next we show the dependences of P on the other parameters. The details of the parameter dependences are discussed in x 4. Figure 1 shows P col against a for e ¼ 0:5 in the standard case. In this case a min ¼ 3:3 AU and a max ¼ 10 AU. Collision takes place over the orbit-crossing region. The probability P col has two peaks, one at each end of the region. The relative velocity between a planetesimal and the planet, v r , explains this feature. In the two-body approximation, the collisional cross section of the planet is given by

Standard Case
where V esc is the surface escape velocity from the planet, given by For the parameter range of this study, v r < V esc and gravitational focusing is effective (e.g., Kokubo & Ida 1996). Thus, the first term on the right-hand side of equation (3) can be neglected, and col is proportional to v À2 r . The number of collisions per T K is proportional to col . The relative velocity for the unperturbed orbit of a planetesimal is given as where v p is the Kepler velocity of the planet (e.g., Bertotti et al. 2003). Figure 2 shows v r against a over the orbit-crossing region for i ¼ 0:05. At both ends of the orbit-crossing region, v r is smaller, and thus, P col is larger. The real relative velocity may be slightly different from equation (4) because of planetary perturbation. However, the essential behavior of the real relative velocity is almost the same as that of equation (4). Figure 1 also shows P col for e ¼ 0:7. In this case a min ¼ 2:9 AU and a max ¼ 16:7 AU. The overall shape of P col is similar to that for e ¼ 0:5. The region for collision expands and P col decreases in its entirety as e increases.

Parameter Dependences
We compare P col for different values of i, a p , and m p in Figure 3. The overall shape of P col does not change with different i, a p , and m p . Figure 3a compares P col for i ¼ 0:05 and 0.07. The value of P col decreases with increasing i. For i ¼ 0:07, P col is about two-thirds of P col for i ¼ 0:05. Figure 3b shows P col for a p ¼ 5 and10 AU. The orbit-crossing region shifts outward and expands with a p , while P col decreases with increasing a p . For a p ¼ 10 AU, P col is about half of P col for a p ¼ 5 AU at the same a/a p . Figure 3c shows P col for m p ¼ 1m J and 0.3m J . The value of P col increases with m p . For m p ¼ 0:3m J , P col is about onefifth of P col for m p ¼ 1m J . Figure 4 shows P esc against a for e ¼ 0:7 in the standard case, in which a min ¼ 2:9 AU and a max ¼ 16:7 AU. Escapers appear over most of the orbit-crossing region. In this region P esc increases with a and suddenly drops before the outer edge of the region. This increase is also explained by v r . The gravitational radius (impact parameter for 90 deflection) of a planet is given by r g ¼ Gm p /v 2 r , which means that planetesimals with smaller v r are easily scattered at large angles. The increase in P esc with a may reflect that v r decreases with increasing a for a k a p (Fig. 2). Figure 4 also shows P esc for e ¼ 0:8, where a min ¼ 2:8 AU and a max ¼ 25 AU. The overall shape of P esc is similar to that for e ¼ 0:7. The value of P esc decreases with increasing e at constant a, although the maximum value of P esc increases with e. The maximum value of P esc for e ¼ 0:8 is about twice that for e ¼ 0:7.

Standard Case
We find that there are no escapers for e 0:3. This is explained by a simple fly-by theory (e.g., Madonna 1997) based on the two-body approximation. By using v r and v p , the velocity of a planetesimal is written as There is also  the minimum relative velocity v min r for escape derived from the fly-by theory. A planetesimal gains the largest additional velocity if it is scattered in the direction of the orbital motion of the planet. In this case v min r is given by v min Fernández 1978). We plot v min r also in Figure 2. In the region where v r is smaller than v min r , no planetesimals can escape. To satisfy the condition v r > v min r , we need e k 0:4. Thus, planetesimals initially with e P 0:4 cannot escape. Even for e k 0:4, the condition is not satisfied near the end of the orbit-crossing region, and this explains the sudden drop in P esc before the end of the region. However, in the case of large m p , the simple fly-by theory is no longer valid. In case 16 there exist escapers even for e ' 0.

Parameter Dependences
In Figure 5 we compare P esc for different values of i, a p , and m p . The overall shape of P esc , the increase and sudden drop against a, does not change with different i, a p , and m p . Figure 5a compares P esc for i ¼ 0:05 and 0.07. Similarly to P col , P esc decreases with increasing i. The probability P col for i ¼ 0:07 is about two-thirds of P col for i ¼ 0:05. In Figure 5b we show P esc for a p ¼ 5 and 10 AU. The orbit-crossing region shifts outward and expands with a p , while P esc at the peaks is almost the same, which may imply that P esc is scaled by a/a p . Figure 5c compares P esc for m p ¼ 1m J and 0.3m J . Compared to P col , P esc increases steeply with m p . The probability P esc for m p ¼ 0:3m J is about 1/10 of P esc for m p ¼ 1m J .

Candidacy To Be a Member of a Comet Cloud
3.3.1. Standard Case Figure 6 shows P can against a for e ¼ 0:7 and in the standard case, in which a min ¼ 2:9 AU and a max ¼ 16:7 AU. We use a can ¼ 3000 AU, which corresponds to the inner edge of the inner Oort Cloud (Dones et al. 2006). Candidates appear almost over the orbit-crossing region. In this region P can increases with a and suddenly drops before the end of the region. This behavior is similar to that of P esc because large scattering is required for producing candidates, as well as escapers. However, the dependence on a is stronger than that of P esc . Figure 6 also shows P can for e ¼ 0:8. In this case a min ¼ 2:8 AU and a max ¼ 25 AU. Similarly to P esc , the value of P can decreases with increasing e at constant a, and the maximum value of P can increases with e. The maximum value for e ¼ 0:8 is about 3 times that for e ¼ 0:7. There are no candidates for e 0:3. It should be noted that P can is always much smaller than P esc .

Parameter Dependences
In Figure 7 we show P can for a variety of parameters including i, a p , m p , and a can . The overall shape of P can does not change with i, a p , m p , and a can . Figure 7a compares P can for i ¼ 0:05 and 0.07. The probability P can for i ¼ 0:07 is about two-thirds of P can for i ¼ 0:05. The decrease in P with increasing i is a feature common to P col , P esc , and P can . Figure 7b compares P can for a p ¼ 5 and 10 AU. The orbit-crossing region shifts outward and expands with a p , and the maximum value of P can increases with a p . The maximum value of P can for a p ¼ 10 AU is about twice that for a p ¼ 5 AU. Figure 7c compares P can for m p ¼ 1m J and 0.3m J . The probability P can for m p ¼ 0:3m J is about 1/10 of P can for m p ¼ 1m J . The value of P can increases with m p in a manner similar to P esc .  We plot P can for a can ¼ 3000 and 5000 AU in Figure 7d . The value of P can decreases with increasing a can . For a can ¼ 5000 AU, the probability P can is about half of P can for a can ¼ 3000 AU. Figures 8, 9, and 10 show P against a for i ¼ 0, i.e., cases 1, 13, and 19. Each panel compares P for different values of e, a p , and m p . In all panels the qualitative features of P are almost the same as for i 6 ¼ 0, but the values are higher. For example, in case 1, P col , P esc , and P can for i ¼ 0 are typically $100, $10, and $5 times higher than in the standard case, respectively. The parameter dependences of P, except on a can , are weaker than those in all cases for i 6 ¼ 0.

EMPIRICAL FITS FOR PROBABILITIES
We derive simple empirical formulae for P using the results of numerical integration, which is useful in estimating P for general planetary systems. For the empirical formulae for P fit , we use simple power-law fitting in a, e, i, a p , m p , and a can , where f is a numerical factor and , , , , , and represent power-law indices. The last term a can is only for candidacy. This assumption works well to reproduce the numerical results. By comparing P fit with P, we estimate the values of f, , , , , , and for collision, escape, and candidacy, and those for i ¼ 0 empirically in simple figures. We adopt integers or simple fractions for their power-law indices. The empirical fits are sufficiently accurate for order-of-magnitude estimation of P. We plot P fit in Figures 1 and 3-10 over the orbit-crossing region.

Collision with a Planet
By comparing equation (5) with P col we obtain the empirical formula where we approximate P col as being constant in a, neglecting peaks at the ends of the orbit-crossing region. Although the fit of P Bt col does not reproduce the peaks at both ends of the orbitcrossing region, it approximates the average value of P col for 0:1 e 0:9. Equation (6) is valid in all cases for i 6 ¼ 0 in Table 1.

Escape from a Planetary System
We consider the cases for e ! 0:5 and m p 3m J , excluding cases 23 and 24, in which the simple fly-by theory is not valid. We replace e with (1Àe) in equation (5) because P esc strongly depends on the perihelion distance q ¼ a (1 À e) rather than on e (e.g., Duncan et al. 1987).
By comparing equation (5) with P esc , we obtain P Bt esc $ 4 ; 10 À6 a a p

3
(1 À e) sin i À1 m p m J 2 : Comparing P esc and equation (7) we find that the overall feature of P esc for i 6 ¼ 0 is approximated by equation (7), except for  e ¼ 0:9 in case 2 and for e ¼ 0:9 in case 22. In these exceptional cases P esc is roughly given by equation (10).
Note that P Bt esc reproduces only the increase in P esc with a and overestimates P esc for a after the rapid decrease around a max . The overestimation is due to the difference between the orbitcrossing region and the escape region. If these two regions are completely the same, there is almost no overestimation. However, the outer edge of the escape region, where v r > v min r , is smaller than the outer edge of the orbit-crossing region, a max . Since this difference is larger for smaller e, the overestimation of P esc by P Bt esc is larger for smaller e. We find that P Bt esc agrees with the increase of P esc with a for e ! 0:5.
For a planetesimal with a 3 a p , in other words, with e % 1, the energy change of planetesimals depends on a and e only through q (Duncan et al. 1987;Dones et al. 1996). Hence, under such conditions, equation (7) is no longer valid.

Candidacy To Be a Member of a Comet Cloud
The parameter ranges of e and m p we consider here are e ! 0:4 and m p 3m J , and we replace e with (1Àe) in equation (5) for the same reason as for escape. By comparing equation (5) with P can , we obtain P Bt can : P Bt can $ 1:2 ; 10 À5 a a p

5
(1 À e) 2 sin i À1 m p m J 2 a can a p À1 : Equation (8) is in good agreement with the increase in P can for a can k 10a and for e ! 0:6. Similarly to P Bt esc , P Bt can ex-presses only the increase in P can with a but not the decrease around a max . We find that the behavior of P can for e ¼ 0:9 in cases 2, 3, and 22 is expressed by equation (11) rather than equation (8).
In previous studies, many authors mentioned Jupiter's inefficiency in forming candidates because of its high mass; however, we find that the ratio of P Bt can to P Bt esc does not depend on m p from equations (7) and (8). For a planetesimal with a 3 a p or e % 1, equation (8) is not valid for the same reason as for escape (Duncan et al. 1987;Dones et al. 1996).

Two-dimensional Cases
We consider P in cases for i ¼ 0. By comparing P col , P esc , P can , and P fit , we obtain P Bt 2D : P Bt 2D;esc $ 2 ; 10 À3 a a p

7=2
(1 À e) 3=2 m p m J a can a p À1 : As seen in Figures 8-10, equations (9)-(11) reproduce P well. We find that the ratio of P Bt 2D;can to P Bt 2D;esc also does not depend on m p from equations (10) and (11).

Geometric Interpretation
In this section we try to understand the relation between P fit and P Bt 2D using a geometric interpretation. We assume that the probability P fit is approximated by the ratio of cross sections for collision, escape, and candidacy to the total cross section: In i 6 ¼ 0 cases we define the total cross section as the cylindrical area at a p through which planetesimals pass, which is estimated as Using equation (12)  : In i ¼ 0 cases P fit is approximately given by the ratio of twice the impact parameter b to the length of the orbit of the planet:  Using equation (17) and P Bt 2D obtained in x 4.4, we obtain b: b col $ 2:5 ; 10 À3 e À1 a p AU 1=2 m p m J

7=2
(1 À e) 3=2 a p m p m J a can a p À1 : The power-law index of b col for m p , col ¼ 2/3, reflects that gravitational focusing is effective (e.g., Kokubo & Ida 1996). The dependences of b esc and b can on m p are the same as that of r g on m p .
From equations (14)-(16) and (18)-(20) we find the relations $ b 2 for collision and escape and $ 2b esc b can for candidacy. These relations suggest that b col , b esc , and b can can be interpreted as the effective maximum impact parameters. The cross sections for collision and escape, col and esc , are roughly given by circular areas with radii of b col and b esc , respectively. The cross section of candidacy is consistently approximated as a narrow ring area with a width b can and a radius b esc just outside of esc . Note that the actual shape of is not a circle or a ring but rather a complicated figure. However, the geometric interpretation with the effective impact parameters is almost consistent with the numerical results in both the i 6 ¼ 0 and i ¼ 0 cases.
In equation (12) we assume total includes . However, the width of the total cross section a p sin i decreases with i and becomes smaller than the width of (2b) for small i. Then the behavior of P for i 6 ¼ 0 becomes similar to that for P Bt 2D rather than for P fit . We can roughly estimate the parameter ranges in which P Bt col and P Bt esc are valid from the condition b a p sin i, which leads to where we use b ¼ / ð Þ 1 = 2 , which is about 50% smaller than the value of b of equation (18) or (19). The criterion for P Bt can to be valid is b esc þ b can < a p sin i. This criterion reduces to the same criterion as for P Bt esc (eq. [22]) because b can Tb esc . However, P Bt can is more sensitive to the criterion than P esc because can is the area outside of esc . These criteria are all consistent with the results of the numerical calculations, as we mention in xx 4.1-4.3.

EFFICIENCIES
We calculate the efficiencies of collision, escape, and candidacy from equation (2) and the numerical results for P. To obtain realistic efficiencies we need to know the orbital distribution of planetesimals around a planet during or after planet formation. However, this distribution is uncertain so far. In the present paper, as the first step, we adopt simple disk models: (1) flat disk ( ¼ 0), (2) standard disk ( ¼ À3/2), and (3) standard disk with a solar system planet. In model 3 we consider the four giant planets of the solar system. We use n 1 ¼ 1 in the disk models for simplicity. We also obtain empirical fits of the efficiency K fit from equation (2) and P fit . The empirical fits K fit are shown with K in Figures 11-13. We plot K Bt col for 0:1 e 0:9, K Bt esc for e ! 0:5, and K Bt can for e ! 0:4. Figure 11 shows K against e in the standard case ¼ 0, a in ¼ 0, and a out ¼ 1. The collision probability K col decreases gradually with increasing e. On the other hand, for e ! 0:4, K esc and K can increase monotonically with e. As we have already seen in xx 4.2 and 4.3, K esc and K can are positive only for e ! 0:4. The relative magnitude between K col and K esc or K can varies with e. Only K col has a positive value for e 0:3. For e ! 0:5, K esc is larger than both K col and K can , and K can exceeds K col for e ! 0:8. The dependences of K on i, m p , and a can are independent of and the same as for P.

Flat Disk
The empirical formulae K fit integrated over the orbit-crossing region for ¼ 0 are where we neglect the terms of (1 þ e) with power-law indices smaller than À 1 2 . These empirical fits are also plotted in Figure 11. The fit of K Bt col shows good agreement with K col . The differences are within $10% for all e. For e ¼ 0:5, K Bt esc agrees with K esc within a factor of $3. For e ! 0:6, K Bt esc agrees with K esc within $70% error. The fit of K Bt can agrees with K can within a factor of $3 for e ¼ 0:4 and within $50% error for e ! 0:5. These errors of K Bt esc and K Bt can decrease with increasing e because the regions where escapers and candidates appear overlap well with the orbit-crossing region for large e, as shown in Figures 4 and 6.
We also compare K and K fit in other cases and find that K is well approximated by K fit . Their errors are typically about the same as for the standard case. In the worst cases K fit agrees with K within $50% for collision for 0:1 e 0:9 and within a factor of $3 for escape for e ! 0:5, except for case 16. In the worst cases for candidacy K Bt can agrees with K within a factor of $5 for e ¼ 0:4 and $2 for e ! 0:5, except for case 16.
The relative magnitude of K also varies with a p and m p . Using equations (23)-(25) we can derive the dependences of the ratios as K esc /K col / a p m 2 = 3 p for e ! 0:4 and K can /K col / a 2 p m 2 = 3 p a À1 can for   e ! 0:5. These relations imply that planets with large semimajor axes and/or large mass produce escapers and candidates effectively compared to the production of colliders. Figure 12 shows K against e in the standard case and for ¼ À3/2, which is the standard value for protoplanetary disks (Hayashi 1981). The integration range is from a in ¼ 0 to a out ¼ 1.

Standard Disk
The dependences of K col , K esc , and K can on e are stronger than for the flat disk model. The relations among K col , K esc , and K can are almost the same as for the flat disk model. The efficiencies K esc and K can exceed K col at e ¼ 0:5 and e ¼ 0:9, respectively.
The empirical formulae K fit for ¼ À3/2 are K Bt col $ 9 ; 10 À6 e À1 sin À1 i : These empirical fits are plotted in Figure 12.
The fit of K Bt col agrees well with K col within $10% error for e 0:6. For e ! 0:7 the errors are larger and within $40%. The fit of K Bt esc agrees with K esc within a factor of $2 for e ¼ 0:5 and within $30% error for e ! 0:5. The fit of K Bt can agrees with K can within a factor of $2 for e ¼ 0:4. For e ! 0:5, K Bt can and K can show good agreement within $15% error. Compared with the flat disk model, the errors between K and K fit of escape and candidacy are small. This is because the contributions of P at large a, where the differences between P and P fit are large, are weakened for small .
The efficiencies in other cases are also approximated by K fit , with errors similar to the standard case. In the worst cases K fit agrees with K within $50% error for collision for 0:1 ! e ! 0:9 and within a factor of $3 for escape for e ! 0:5. In the worst cases for candidacy, K Bt can expresses K can within a factor of $4 for e ¼ 0:4 and $2 for e ! 0:5.

Application to the Solar System
We apply the results of P to the solar system and compare the K-value of the four giant planets, Jupiter, Saturn, Uranus, and Neptune, using the empirical fits K Bt col , K Bt esc , and K Bt can . We use the present values of a p and m p and assume the eccentricities of the planets as e p ¼ 0. We adopt ¼ À3/2, which corresponds to the standard protoplanetary disk for the solar system. The integration range is from a in ¼ 0 to a out ¼ 1. Figure 13 shows K Bt col , K Bt esc , and K Bt can for the giant planets against e for i ¼ 0:05 and a can ¼ 3000 AU. Jupiter always has the highest K fit because of having the largest m p . Among the four planets, the inner planets have higher K Bt col , and the massive planets have higher K Bt can . For escape, Uranus and Neptune have almost the same values of K Bt esc . The relative magnitudes between K Bt col and K Bt esc and those between K Bt esc and K Bt can do not vary with the planets. All K Bt esc values exceed K Bt col and K Bt can for e ! 0:5. Only the relative magnitudes between K Bt col and K Bt can vary slightly with the planets. For Jupiter, Saturn, and Uranus K Bt can exceeds K Bt col for e ¼ 0:9. For Neptune K Bt can exceeds K Bt col for e ! 0:8.
Next we consider the case in which i is proportional to the reduced Hill radius of the planet: h ¼ r H /a p . For each planet we set the inclination of the planetesimals to i ¼ h. This application reflects that planetesimals around a planet are excited to the degree that i / h (Ida & Makino 1993). Figure 14 shows K Bt col , K Bt esc , and K Bt can for the giant planets against e for i ¼ h and a can ¼ 3000 AU. The relative importance of the planets in each K fit and the relation among K fit values for a planet are almost the same as for i ¼ 0:05. Compared to the case for i ¼ 0:05, K fit slightly decreases for Jupiter because for Jupiter i > 0:05. For other planets, as i < 0:05, K fit increases. In this model Jupiter still has the highest K fit , despite the disadvantage of having the highest i among the four planets.
In the real solar system planets have finite e p ; thus, the results of circular restricted three-body formalism may not be directly applicable. However, the e p values of the four giant planets are small (e p ' 0:05) and do not make any large differences in K. We perform calculations in the standard case with e p ¼ 0:05 and find that the difference between K-values for e p ¼ 0 and 0.05 is less than 10% for collision and less than 3% for escape. For candidacy the difference is typically within $5%, with a maximum difference of $25%.

SUMMARY AND DISCUSSION
We performed numerical calculations of the first dynamical stage of comet cloud formation, scattering of planetesimals by a planet. The orbital evolution of planetesimals was investigated using circular restricted three-body formalism. We considered planets of m p ¼ 0:1 10 ð Þm J and a p ¼ 1 30 AU and planetesimals of e ¼ 0:1 0:9 and i ¼ 0 0:1.
We obtained the probabilities P of collision with a planet, escape from a planetary system, and candidacy for inclusion in a comet cloud for a single encounter as functions of orbital parameters of planets and planetesimals. We found that a planetesimal with an initial eccentricity of e k 0:4 can escape from the planetary system or be a candidate for the comet cloud due to scattering by a planet. The probability of any planet producing escapers is always much higher than that of producing candidates, since the energy range of the comet cloud is narrow. Furthermore, the production ratio of candidates to escapers is independent of m p . We also derived simple empirical formulae for these probabilities that are accurate enough to use for orderof-magnitude estimation. Using the probabilities and assuming the distribution of planetesimals, we obtained the efficiencies K of planets for collision, escape, and candidacy.
We applied the results to the giant planets in the solar system and the standard disk model for solar system formation. We found that among the four giant planets, Jupiter is most responsible for producing candidates for elements of the Oort Cloud insofar as the inclination of planetesimals is constant or proportional to the reduced Hill radius of each planet.
Simulations in Dones et al. (2006) showed that the typical comet in the Oort Cloud is a planetesimal originally from the Uranus-Neptune region, placed in the Oort Cloud by Saturn. Dones et al. (2006) calculated the efficiency for each planet, which is the ratio of the number of planetesimals remaining in the Oort Cloud to that of planetesimals that had close encounters with the planet after 4 Gyr. The efficiencies for Jupiter and Saturn are $2% and are about 1/10 of those for Uranus or Neptune. They concluded that Jupiter and Saturn eject from the solar system many planetesimals that had close encounters with them and that their efficiencies for populating the Oort Cloud are low compared to those for Uranus and Neptune. On the other hand, what we evaluate as P can is the probability of the formation of candidates from among all planetesimals with a certain e and i and on crossing orbits with a planet, taking into account not only close encounters with the planet but also distant encounters. Our results show that under the same conditions for planetesimals, Jupiter has the highest P can among the four planets, or the highest potential ability to form candidates.
In this paper we investigated the elementary process of scattering of planetesimals by a planet and applied the results to simple planetesimal disk models. In order to construct a more realistic formation scenario of a comet cloud, we have to clarify the number and orbital distributions of planetesimals around planets during and just after planet formation. The distributions are affected by the structure of planetary systems, the interactions among planetesimals, and the existence of the gas disk. By applying our results to the realistic distributions of planetesimals, we will be able to discuss a more realistic scenario of comet cloud formation.
This work was partially supported by the Ministry of Education, Culture, Sports, Science, and Technology, Japan, the 21st Century COE Program ''Origin and Evolution of Planetary Systems'' and a Grant-in-Aid for Scientific Research on Priority Areas, ''Development of Extrasolar Planetary Science.'' A. H. is supported by a Japan Society for the Promotion of Science Research Fellowship for Young Scientists. We wish to thank E. I. Chiang for his useful comments. We would also like to thank the anonymous referee for valuable comments and suggestions.