Direct N-body simulation of star clusters or other stellar systems has proven itself an extremely powerful tool to study the structure and evolution of stellar systems, since the pioneering work by von Hörner ([1960,1963]). The only known way to do experiments on stellar systems is to construct their models in computers, because we cannot do laboratory experiments on stellar systems,
Of course, N-body simulation is not the only way to construct computer models of star clusters. One could use Monte-Carlo (Henon , for recent development see Giersz ) or direct integration of Fokker-Planck (FP) equations (Cohn , Takahashi , Drukier et al. ). In particular, recent advance in the treatment of the two-dimensional [f(E,J)] FP equation made it possible to use FP code to the study of the evaporation of clusters in the tidal field of the parent galaxy.
The main advantage of these methods over direct N-body simulation is the calculation cost. While it is still out of reach to perform the direct simulation of, say, a typical globular cluster with , the FP equation is valid in the limit of . Therefore, when the assumption of large-N limit is good, methods based on the FP approximation can deliver reliable results for the calculation cost orders of magnitudes smaller than that of direct N-body calculation. In principle, one can do the N-body simulation with small N and then try to extrapolate that result to larger N. However, this problem, which is sometimes called ``the scaling problem'' has proven itself a complex and difficult problem(Heggie et al. ). The main reason of the difficulty is that there are many physical processes of which timescales depend on the number of particles in different ways, and that there is no way to adjust all relevant timescales consistently.
On the other hand, methods based on FP approximation do have their own limitations. For example, it's pretty difficult to extend the formulation beyond the spherical symmetry. Moreover, they also suffer the problem of the time scaling, only difference is that they are affected from the opposite direction. For example, when we want to include the dynamical effect of the slowly varying potential, FP approximation goes into trouble. Consider the tidal stripping of globular clusters with non-circular orbit. The orbital timescale is of the order of years. The orbital timescale of the stars around the tidal boundary of the cluster is also years. On the other hand, the half-mass relaxation time is of the order of years. Therefore, the timestep to integrate the FP equation is at the largest years, and in practice much smaller than that. Clearly, the assumption that the dynamical timescale is smaller than the Fokker-Plank timestep is broken, and FP calculations tends to grossly overestimate the escaper rate.
As beautifully demonstrated by Takahashi and Portegies Zwart (), one can incorporate the escaping timescale into FP calculation. However, in order to do so the model parameter must be adjusted so that the agreement with reference N-body calculations is achieved. Moreover, this agreement has been so far achieved only with N-body simulations with simple static ``tidal cutoff'', where no tidal field is imposed and stars outside the tidal radius are just removed. Whether or not FP calculations can reproduce the result of N-body calculations with tidal fields is an open question.
The ultimate solution for the scaling problem is to perform N-body simulations with sufficiently large number of particles. In the first half of this paper, we describe our effort in that direction, the GRAPE project to develop and use special-purpose computers for N-body simulations. In the last half of this paper, we'd like to discuss briefly our recent work on N-body simulations of galactic centers with massive black holes, where again the scaling problem shows up.