Article Views Abstract Full Text (HTML) Full Text (PDF) Article Tools Save to My Folders Download Citation Alert Me When Article is Cited Post to CiteULike E-mail This Page Submit an E-Letter Commercial Reprints and Permissions View PubMed Citation Related Content Similar Articles In: Science Magazine Web of Science PubMed Search Google Scholar for: Articles by Hut, P. Articles by Makino, J. Find Citing Articles in: Web of Science (21) CrossRef Google Scholar Scopus My Science My Folders My Alerts My Saved Searches Sign In More Information More in Collections Astronomy Related Jobs from ScienceCareers Astronomy

Science 22 January 1999: Vol. 283. no. 5401, pp. 501 - 505 DOI: 10.1126/science.283.5401.501 Prev | Table of Contents | Next Review Astrophysics on the GRAPE Family of Special-Purpose Computers Piet Hut and Jun Makino The GRAPE-4, the world's fastest computer in 1995-1997, has produced some major scientific results through a wide diversity of large-scale simulations in astrophysics. Applications have included planetary formation, the evolution of star clusters and galactic nuclei, and the formation of galaxies and clusters of galaxies. P. Hut is at the Institute for Advanced Study, Princeton, NJ 08540, USA. J. Makino is in the Department of General Systems Study, College of Arts and Sciences, University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan. Computational physics has emerged as a third branch of physics, grafted onto the traditional pair of theoretical and experimental physics. At first computer use seemed to be a straightforward offshoot of theoretical physics, providing solutions to sets of differential equations too complicated to solve by hand. But soon the quantitative improvement in speed yielded a qualitative shift in the nature of these computations. Rather than asking particular questions about a model system, we now use computers more often to model the whole system directly. Answers to relevant questions are then extracted only after a full simulation has been completed. The data analysis following such a virtual lab experiment is carried out by the computational physicist in much the same way as it would be done by an experimenter or observer analyzing data from a real experiment or observation. Recent increases in computer speed are more modest than what could be expected purely from the ongoing miniaturization of computer chips. Because the number of transistors on a single chips doubles every 1.5 years, a chip now contains a hundred times more transistors than it did 10 years ago. With a clock speed increase of more than a factor of 10, one might have expected a speed increase of more than a factor of 1000 over the last decade. However, the actual speed increase of a typical computer chip has been at most a factor of 100, lagging far behind theoretical expectations. The reason for this relatively poor performance lies in the costs resulting from the growing complexity of a general-purpose chip. Designing a chip for only one specific purpose yields a rapidly growing payoff. Therefore, the time seems ripe to explore which types of calculations can be realized directly in hardware, in the form of special-purpose computers, rather than run in software on general-purpose computers. One such project has resulted in the GRAPE (short for Gravity Pipe) family of special-purpose hardware, designed and built by a group of astrophysicists at the University of Tokyo (1). Like a graphics accelerator speeding up graphics calculations on a workstation without changing the software running on that workstation, the GRAPE acts as a Newtonian force accelerator in the form of an attached piece of hardware. In a large-scale gravitational N-body calculation, where N is the number of particles, almost all instructions of the corresponding computer program are thus performed on a standard workstation, while only the gravitational force calculations, in the innermost loop, are replaced by a function call to the special-purpose hardware. Specifically, the force integration and particle pushing are all done on the host computer, and only the interparticle force calculations are done on the GRAPE (Fig. 1). This may seem problematic, given that the intrinsic speed of the GRAPE is faster than that of the host computer (an ordinary workstation) by a factor of 10,000. However, the interparticle calculations require a computer processing power that scales with N2, whereas all other actions on the host scale only in proportion to N. Therefore, each doubling of the number of particles doubles the work load on the GRAPE relative to that on the workstation. In this way, no matter how slow the workstation, it will be able to keep up with the GRAPE for large enough values of N. Fig. 1. The GRAPE-4 hardware. [View Larger Version of this Image (89K GIF file)]

For some applications, more efficient algorithms have been devised that require the computation of a number of interparticle force calculations that scales with N log N, rather than N². It turns out that even these methods can still be efficiently run on the GRAPE (2); although the asymptotic scaling advantage is not very large in that case, the overall coefficient in the scaling relation turns out to favor the use of the GRAPE. Some versions of the GRAPE (Table 1) allow arbitrary force implementations for applications such as molecular dynamics. For example, the MDGRAPE has been used to study the structure of protein molecules (3). However, most GRAPEs have been used to study astrophysical problems such as those reviewed below.

Table 1. Summary of GRAPE hardware (1). Speed is expressed in flops (floating-point operations per second). Machine Year Peak Speed Notes Limited-precision data path GRAPE-1 1989 240 Mflops Concept system GRAPE-3 1991 15 Gflops 48 custom chips, 10-MHz clock GRAPE-5 (1999) ~1 Tflop Under development Full-precision data path GRAPE-2 1990 40 Mflops IEEE precision, commercial chips HARP-1 1993 180 Mflops Force and its time derivative GRAPE-4 1995 1.1 Tflops  1692 custom chips, 32-MHz clock GRAPE-6 (2000) ~200 Tflops  Under development Arbitrary force law    GRAPE-2A 1992 180 Mflops Force look-up table MDGRAPE 1995 4 Gflops  Custom chip with force look-up table MDGRAPE-2 (2000) ~100 Tflops  Under development

Star Cluster Evolution

A globular star cluster (Fig. 2) (4) typically contains about a million stars, packed together much closer than the stars in the solar neighborhood. Such a cluster describes a wide orbit around the parent galaxy, well separated from the stars in that galaxy. If we lived in the core of a dense globular cluster, the brightest stars would appear as bright as the full moon, which would make them too bright to look at directly, given their point-like nature. Optical astronomy of anything but the nearby stars in the same globular cluster would be rather difficult in such a situation.

---

Fig. 2. (left). A globular cluster, M15, as seen by the Hubble Space Telescope (4).
Fig. 3. (right). The evolution of the central density _c. Thirty time units correspond roughly to one initial half-mass relaxation time. Curves for different values of N are vertically shifted by 3 units. [View Larger Versions of these Images (1 + 1K GIF file)]

---

In the simplest approximation, we can study a globular cluster as a collection of point masses, reducing the problem to the gravitational N-body problem, which was solved by Newton for N = 2 but was only studied in detail for N > 2 when computers became available. Any localized distribution of particles will tend to become spherical, as a result of forgetting the initial conditions, on a two-body relaxation time scale:

(1)

where the crossing time t_cr is a measure of the time it takes for a typical star to move across the cluster.

Heat is transported through the cluster, as a consequence of many two-body encounters, on the time scale t_rel. On longer time scales, any self-gravitating star system is unstable. Because the system tends to relax toward a Maxwellian velocity distribution, there are always some stars that acquire a velocity that exceeds the escape velocity, after which they are lost from the system. Other stars tend to congregate in the central regions, which grow denser at an ever-increasing rate, because higher density implies more frequent encounters and hence a faster two-body relaxation.

This runaway redistribution of energy and mass leads to a phenomenon called gravothermal collapse, often called core collapse, which takes place on a time scale t_cc ~ 10t_rel. Core collapse was hinted at in numerical simulations (5) in the 1960s and was verified through direct N-body simulations (6) and modeled by semianalytic methods (7) in the 1970s. Core collapse is a fundamental feature of long-term stellar-dynamical evolution, showing the instability that results from the negative specific heat inherent in self-gravitating systems. During core collapse, at first the system can be modeled as passing through a series of self-gravitating equilibrium models exhibiting a maximum entropy for a finite central concentration. Once this maximum is passed, subsequent evolution will increase the entropy, and the structure of the star cluster is forced to deviate from that of an equilibrium model.

Even in an idealized system of self-gravitating point particles, core collapse will be halted before an infinite central density is reached. When the central density is high enough, occasional close encounters between three unrelated particles will form bound pairs (binary stars in the case of star clusters), with the third particle carrying off the excess kinetic energy required to leave the other two particles bound. Subsequent encounters between such pairs and other single particles tend to increase the binding energy of these pairs, which leads to a heating of the surrounding system of single particles.

When enough pairs have been formed in this way, the resulting energy production will reverse the process of core collapse. After reaching a minimum radius and a maximum density, the core region will expand again. Core collapse, when it threatens to occur by the collective effects of two-body relaxation, can thus be narrowly averted by a handful of crucial three-body or four-body reactions in the dense core of a nearly collapsed cluster. What will happen next depends on the total number N of particles in the system. If this number is sufficiently small, N  10,000, the whole system will slowly and steadily expand. In this case, the steady energy production in three-body encounters in the center of the system is in steady-state equilibrium with the continuous loss of energy through the outskirts of the system.

If the total number of particles exceeds 10⁴, however, a different behavior emerges. The more particles there are in the system, the higher the central density must become to halt core collapse. As a result, the postcollapse phase features a short relaxation time in the center of the cluster, shorter than the relaxation time in the outer regions, where most of the particles can be found. From the point of view of the inner core dynamics, the bulk of the mass further out seems almost frozen. It is this discrepancy in time scales that can cause the inner core to become "impatient" and to revert to a local collapse, triggered by the slightest fluctuation in the direction of the energy flow produced by stochastic three- and four-body interactions.

What happens then is that about 1% of the inner particles will go into a coherent collapse, locally reminiscent of the original core collapse. As before, bound pairs of particles spring into action, generate energy, and manage to reverse the collapse in the nick of time, preventing an infinite central density from building up. This process repeats itself, leading to irregular oscillations of the core of the cluster.

The existence of these oscillations was unknown until 1983, when they were first found in approximate simulations (8). Dubbed "gravothermal oscillations," they were subsequently analyzed in detail with semi-analytic methods (9). Their occurrence was confirmed in a variety of approximate numerical simulations (10) and was shown to correspond to low-dimensional chaos for large values of N (11). Direct verification of the existence of these oscillations was attempted, using the fastest supercomputers available, but these attempts were unsuccessful (12, 13).

These oscillations were difficult to confirm through direct N-body simulations because it has not been possible to model star cluster evolution with more than 10,000 particles until the advent of the GRAPE-4. This may seem surprising, given the fact that cosmological simulations now routinely handle up to 1000 million particles. The main difference between the two types of calculations lies in the higher accuracy required for star cluster simulations, together with the much larger number of time steps required, relative to cosmological simulations. Following the gravothermal collapse requires a very accurate integration of the equations of motion. The required accuracy is difficult to achieve using approximate schemes like tree codes (14); therefore, traditional direct-summation schemes must be used. Even on supercomputers, the maximum particle number is thus limited to about 10,000. In the case of star clusters, each particle stands for an individual star and thus has a direct physical meaning. In the case of a cosmological simulation, each galaxy is represented by a relatively small number of particles that sample the distribution of stars in phase space. Each particle thus represents the average behavior of many millions of stars. The time steps used can therefore be much larger than would be the case if we were to follow the close encounters of individual stars.

The existence of gravothermal oscillations was proven when they were seen in a direct N-body simulation on the GRAPE-4 that was able to incorporate N values beyond N = 10,000 (15) (Figs. 3 and 4). After core collapse, the fluctuations in central density grow with increasing N values (Fig. 3). For the largest N values displayed, the typical behavior of core oscillations emerges, with its deep and long-lasting troughs punctuated with brief interludes of high core density. For smaller N values, some oscillatory behavior seems to be present but less pronounced. The results of the central density evolution, while suggestive, do not answer the question of the existence of gravothermal oscillations.

---

Fig. 4. Changes in central density  and central velocity dispersion v_c² for a simulation with 32,000 particles. Each data point presents a time average, obtained by averaging  and v_c² over 80 snapshots. Arrows indicate the direction of evolution. [View Larger Version of this Image (21K GIF file)]

---

Figure 4 provides the proof of the gravothermal nature of these oscillations (16). The thermodynamic cycle exhibited by the central density and "temperature" (as measured by the velocity dispersion) is traversed in the opposite direction from that of a Carnot engine: The decompression stage takes place at a lower temperature than the compression stage. This is a reflection of the negative heat capacity of self-gravitating systems: Compression leads to a temperature increase resulting in more heat loss and hence more compression, with the opposite effects holding during decompression. The period of decompression finishes when the core expands beyond the central isothermal area (17).

Black Hole Spiral-In

When two galaxies collide, they are likely to stick together if their relative speed is not too high. Within a few crossing times, the transient ripples and distortions will be smoothed out, and the resulting single galaxy will settle down into a new equilibrium configuration. While all this is going on, the dense nuclei of both galaxies will spiral in as a result of dynamical friction in the central regions of the collision. Finally they, too, merge to form a single dense nucleus.

Recently, many black holes at the centers of galaxies have been detected, with masses spanning a range from 1 million to 1000 million solar masses, up to 0.1% of the mass of the parent galaxy (18). Many, if not most, galaxies harbor such massive black holes in their center. It is far from clear what happens when two such galaxies collide. At first, the two black holes will keep circling each other within the single, newly formed dense nucleus. Although dynamical friction tends to let them spiral in rapidly at first, this process becomes considerably less efficient by the time the amount of mass in stars between the two holes becomes smaller than the mass of the holes themselves. The stars that initially tend to be most efficient in providing a braking mechanism are scattered into different orbits. As a result, the system may reach a stagnation point in which little further dynamical friction occurs.

The prediction of this stagnation process was made almost 20 years ago (19), and since then many attempts have been made to check this prediction quantitatively, using large-scale N-body calculations. Until the advent of the first GRAPEs, this problem was completely intractable, even on the largest supercomputers available. One reason that the GRAPE computers are suitable for this type of problem is the intrinsically high dimensionality of the problem. With two black holes in an eccentric orbit around each other, there is no symmetry in either configuration space or velocity space. As a result, the stellar dynamics problem is truly six-dimensional when seen as a fluid flow in phase space.

In contrast, modeling a globular cluster is often done by assuming spherical symmetry, which leaves only one spatial dimension (radial) and two velocity dimensions (radial and tangential) to worry about. In practice, further simplifications have often been made in which the distribution function of the stars is assumed to be dependent only on energy, or sometimes on energy and angular momentum. Fokker-Planck methods have therefore been useful in modeling globular clusters, especially during the core collapse phase. After core collapse, during the reexpansion phase, the effects of binaries must be taken into account, an extremely granular process that defies the main Fokker-Planck assumptions of smoothness of the distribution function. Even so, it has been useful to compare the full N-body calculations in the postcollapse domain with approximate Fokker-Planck treatments. However, a Fokker-Planck treatment of a six-dimensional system is completely impractical from the outset.

The first attempts to use the GRAPE to tackle this problem were made in 1990 (20) with the GRAPE-2, followed by more recent attempts (21) on the GRAPE-4. Three important conclusions have emerged from these studies:

1) When two identical galaxies, each harboring a central black hole, merge, they will produce a merger remnant with a ratio of core radius r_c to half-mass radius r_h that is comparable to that of the original galaxies. In contrast, galaxies without black holes tend to produce merger remnants in which r_c/r_h is smaller than in the original galaxies. In the former case, r_c/r_h ~ M_BH/M_tot, where M_BH is the mass of the central black hole and M_tot is the mass of the whole galaxy.

2) This "core," formed around the black hole binary after the merging of the two galaxies, does not have a completely flat density distribution in the center. In fact, it looks more like the "weak cusps" observed in many galaxies by the Hubble Space Telescope (22). The formation mechanism of this cusp is not well understood.

3) Whether or not a black hole binary, lurking in the core of a merger remnant, has had time to spiral in within the current age of the universe, and under what circumstances, is still largely an open question. We expect the continuum limit to be reached for N  10⁷. These calculations will only be feasible with the GRAPE-6 (Table 1).

Formation and Evolution Processes, from Planets to Galaxy Clusters

We next discuss how the GRAPE computers have been used to study the origin of structure in the universe, from very small scales to the largest scales that can be observed. On the small end, the coagulation of grains and boulders to form planets has been modeled in an effort to understand the formation process of our own planetary system, as well as that around nearby stars. Increasing the length scale of interest by a factor of 1000 million, we discuss the formation of galaxies. Multiplying the size by another factor of 1000, we reach the scale at which rich clusters of galaxies evolve.

Planet formation. After the sun was formed, some matter of the protosolar nebula was left in a disk around the sun. Grains that condensed out of the original gas coagulated through collisions to form larger and larger particles, the size of pebbles, boulders, and larger protoplanetary bodies. To model this process in detail has turned out to be difficult, because significant evolution takes place on a time scale larger than a crossing time by a factor of 1 million or more.

The main stumbling block has been the need to simultaneously model the presence of a wide variety of particle sizes, or equivalently, masses. A little more than 10 years ago, it was realized that dynamical friction plays an essential role in planetary formation (23, 24). This process forces more massive particles to have smaller random velocity, which effectively increases their collision cross section. Thus, massive particles can grow much more rapidly than less massive particles.

Kokubo and Ida (25) used the GRAPE-4 to model this type of growth of planetesimals, under the assumptions that the accretion was perfect (that is, the collisions were totally inelastic) and that there was no gas left in the system to cause nongravitational drag on the particles. They found the mass distribution to relax quickly to a continuous power-law mass distribution with dN/dm  m^-2.5, where N is the cumulative number of bodies, independent of the initial mass distribution. [This result was later derived analytically (26).] Their most interesting result was that the heaviest body would subsequently detach from the continuous power-law distribution, featuring a much more rapid growth in mass, called runaway growth, that could lead to the formation of a planet.

Kokubo and Ida (27) again used the GRAPE-4 to study the later stages of planet formation on a more global scale. The earlier local runaway studies, leading to the formation of a single protoplanet, give rise to multiple protoplanet formation when a large fraction of the protoplanetary disk is modeled. They found that such protoplanets are formed and keep growing independently, provided their orbital separations are wide enough. After a while, the growth rate of these protoplanets slows down, because their gravitational perturbations increase the random motion of the swarm of planetesimals in which they are embedded. A continuous mass distribution of relatively light planetesimals can thus coexist with a small number of large protoplanets for millions of years.

Galaxy formation. To study the formation of a single galaxy, it is important to model its environment, out to large distances, given the long-range character of the gravitational force, which through tidal effects influences the angular momentum distribution within the contracting gas clouds destined to form galaxies. In addition, it is essential to model the gasdynamical effects that influence the early phases of galaxy formation. Although the GRAPE has been designed primarily for stellar-dynamical computations, it has proven to be flexible in accommodating deviations from an inverse square law. A key property of the GRAPE hardware is that it uses the interparticle distances, which are computed in order to calculate the pairwise gravitational forces, to construct for each particle a list of neighboring particles that reside within a prescribed distance.

Using this neighbor list, hydrodynamical simulations can be run on the front-end workstation. The prime example here is smoothed particle hydrodynamics (28). Examples of these types of simulations include the formation of galaxies (29), the physical origin of Lyman- and metal line absorption systems (30), the structure of galaxy clusters (31), and the fragmentation of molecular clouds (32).

Simulations of galaxy formation have demonstrated that structure, kinematics, and chemical evolution of model galaxies that form in hierarchical clustering scenarios agree with corresponding properties of observed galaxy populations (33). The major shortcoming is that simulated galaxies are too concentrated. This is usually referred to as the angular momentum problem (29) and suggests that efficient feedback due to late stages of stellar evolution (for example, winds and supernovae) is needed for a successful galaxy formation model.

Simulations of damped Lyman- absorption systems demonstrated that nonequilibrium dynamics can easily explain the apparent discrepancy between the observed high velocity of low ionization lines and the relatively small circular velocity predicted by hierarchical models of structure formation (30). The evidence that damped Lyman- absorbers at high redshift are related to large rapidly rotating disks, which would disagree with the hierarchical clustering hypothesis, is thus not compelling (34).

Galaxy cluster evolution. Galaxy formation is an extremely slow process, starting sometime within the first 1000 million years after the Big Bang and continuing today. Most galaxies are formed in isolation or in small groups, but some galaxies are formed in much richer groups called clusters of galaxies, or even superclusters of galaxies. The typical properties of galaxies formed in such clusters differ from those of galaxies formed elsewhere. For example, most galaxies in clusters are elliptical, whereas most field galaxies are spiral (35).

To what extent do these differences reflect the different formation history of the galaxies, as they may have been affected, for example, by the much higher matter density in the sites where rich clusters of galaxies were born? And to what extent do the differences reflect later modifications to the galaxies resulting from the different dynamical environment of a rich cluster? In attempts to resolve this nature-versus-nurture debate, the GRAPE has been used to model the internal evolution of a rich galaxy cluster.

Apart from the calculations by Bartelmann and Steinmetz (31), already mentioned in the previous section, earlier work by Funato et al. (36) simulated the evolution of clusters of galaxies containing 32 to 128 galaxies. What they found is that "passive" evolution of galaxies, caused by mutual encounters as well as by the influence of the tidal field of the parent cluster, alters the mass and size of individual galaxies. In particular, they found that passive evolution leads to a distribution of masses with M()  ⁴, where  is the internal velocity dispersion of the stars within a galaxy.

To understand the detailed mechanism of this passive evolution, Funato and Makino (37) used the GRAPE-4 to study a large number of encounters between two isolated galaxies, in order to determine how the resulting changes of mass and binding energy depend on the models used for the galaxies and on the parameters describing the type of encounter. They then estimated the cumulative effect of encounters in the setting of a rich cluster of galaxies. They again found that the mass distribution of galaxies, for the mass M of a galaxy as a function of its velocity dispersion , tends to approach M()  ⁴. This resembles the observational Faber-Jackson relation, the empirical result that the luminosity of a galaxy L()  ⁴ for elliptical galaxies. Note that the remnants of collisions between galaxies typically resemble elliptical galaxies, even if the progenitors were spiral galaxies or other types of galaxies. Because it is also reasonable to assume that M  L, this agreement with observations suggests that the encounters of galaxies play an important role in the evolution of galaxies in a cluster of galaxies.

REFERENCES AND NOTES

1. D. Sugimoto, et al., Nature 345, 33 (1990) , J. Makino and M. Taiji, Scientific Simulations with Special-Purpose Computers--The GRAPE Systems (Wiley, Chichester, UK, 1998).
2. J. Makino, Publ. Astron. Soc. Jpn. 43, 621 (1991); E. Athanassoula, A. Bosma, J.-C. Lambert, J. Makino, Mon. Not. R. Astron. Soc. 293, 369 (1998).
3. Y. Komeiji, J. Comp. Chem. 18, 1546 (1997).
4. P. Guhathakurta, B. Yanny, D. P. Schneider, J. N. Bahcall, Astron. J. 111, 267 (1996).
5. M. Hénon, Ann. Astrophys. 24, 369 (1961).
6. S. J. Aarseth, M. Hénon, M. R. Wielen, Astron. Astrophys. 37, 183 (1974).
7. D. Lynden-Bell and P. P. Eggleton, Mon. Not. R. Astron. Soc. 191, 483 (1980).
8. D. Sugimoto and E. Bettwieser, ibid. 204, 19 (1983); E. Bettwieser and D. Sugimoto, ibid. 208, 493 (1984).
9. J. Goodman, Astrophys. J. 313, 576 (1987).
10. H. Cohn, P. Hut, M. Wise, ibid. 342, 814 (1989). J. L. Breeden, H. N. Cohn, P. Hut, ibid. 421, 195 (1994); R. Spurzem and P. D. Louis, in Structure and Dynamics of Globular Clusters, S. G. Djorgovski and G. Meylan, Eds. (ASP Conf. Ser. Vol. 50, Astronomical Society of the Pacific, San Francisco, 1993), pp. 135-136.
11. J. L. Breeden and H. N. Cohn, Astrophys. J. 448, 672 (1995).
12. J. Makino and D. Sugimoto, Publ. Astron. Soc. Jpn. 39, 589 (1987).
13. R. Spurzem and S. J. Aarseth, Mon. Not. R. Astron. Soc. 282, 19 (1996).
14. J. Barnes and P. Hut, Nature 324, 446 (1986) .
15. J. Makino, in Dynamical Evolution of Star Clusters, P. Hut and J. Makino, Eds. (IAU Symp. 174, Kluwer, Dordrecht, Netherlands, 1996), pp. 151-160.
16. ___, Astrophys. J. 471, 796 (1997).
17. More realistic treatments of star clusters, where the evolution of individual stars is modeled, are under way [ T. Fukushige and D. C. Heggie, Mon. Not. R. Astron. Soc. 276, 206 (1995); S. Portegies Zwart, P. Hut, J. Makino, S. L. W. McMillan, Astron. Astrophys., in press; S. J. Aarseth and D. C. Heggie, Mon. Not. R. Astron. Soc. 297, 794 (1998).
18. J. Kormendy and D. Richstone, Annu. Rev. Astron. Astrophys. 33, 581 (1995).
19. M. C. Begelman, R. D. Blandford, M. J. Rees, Nature 287, 307 (1980) .
20. T. Ebisuzaki, J. Makino, S. K. Okumura, ibid. 352, 212 (1991) .
21. J. Makino and T. Ebisuzaki, Astrophys. J. 465, 527 (1996); J. Makino, ibid. 478, 58 (1997).
22. T. R. Lauer, et al., Astron. J. 110, 2622 (1995).
23. R. Greenberg, J. F. Wacker, C. R. Chapman, W. K. Hartmann, Icarus 35, 1 (1978).
24. G. W. Wetherill and G. R. Stewart, ibid. 77, 330 (1989).
25. E. Kokubo and S. Ida, ibid. 123, 180 (1996).
26. J. Makino, T. Fukushige, Y. Funato, E. Kokubo, New Astron., in press.
27. E. Kokubo and S. Ida, Icarus 131, 171 (1998).
28. M. Umemura, et al., Publ. Astron. Soc. Jpn. 45, 311 (1993); M. Steinmetz, Mon. Not. R. Astron. Soc. 278, 1005 (1996).
29. M. Navarro and M. Steinmetz, Astrophys. J. 478, 13 (1997).
30. M. Haehnelt, M. Steinmetz, M. Rauch, ibid. 465, L95 (1995).
31. M. Bartelmann and M. Steinmetz, Mon. Not. R. Astron. Soc. 283, 431 (1996); T. Suginohara and J. P. Ostriker, Astrophys. J. 507, 16 (1998).
32. R. S. Klessen, A. Burkert, M. Bate, Astrophys. J. 501, L205 (1998).
33. M. Steinmetz and E. Müller, Mon. Not. R. Astron. Soc. 276, 549 (1995).
34. M. Haehnelt, M. Steinmetz, M. Rauch, Astrophys. J. 495, 647 (1998).
35. J. Binney and M. Merrifield, Galactic Astronomy (Princeton Univ. Press, Princeton, NJ, 1998).
36. Y. Funato, J. Makino, T. Ebisuzaki, Publ. Astron. Soc. Jpn. 45, 289 (1993).
37. Y. Funato and J. Makino, in preparation.
38. We thank Y. Funato, E. Kokubo, and M. Steinmetz for discussing their GRAPE results with us. We also thank an anonymous referee for helpful comments.