Tidal mass loss in star clusters and treatment of escapers in FokkerPlanck models
Abstract
This paper presents a new scheme to treat escaping stars in the orbitaveraged FokkerPlanck models of globular star clusters in a galactic tidal field. The existence of a large number of potential escapers, which have energies above the escape energy but are still within the tidal radius, is taken into account in the models. The models allow potential escapers to experience gravitational scatterings before they leave clusters and thus some of them may lose enough energy to be bound again. It is shown that the mass evolution of the FokkerPlanck models are in good agreement with that of body models including the full tidalforce field. The massloss time does not simply scale with the relaxation time due to the existence of potential escapers; it increases with the number of stars more slowly than the relaxation time, though it tends to be proportional to the relaxation time in the limit of a weak tidal field. The FokkerPlanck models include two parameters, the coefficient in the Coulomb logarithm and the coefficient controlling the efficiency of the mass loss. The values of these parameters are determined by comparing the FokkerPlanck models with the body models. It is found that the parameter set works well for both singlemass and multimass clusters, but that the parameter set is another possible choice for multimass clusters.
keywords:
stellar dynamics – globular clusters: general – galaxies: star clusters: general – methods: numerical1 Introduction
The numerical integration scheme of the orbitaveraged FokkerPlanck (FP) equation developed by Cohn (1979) has been one of the most useful tools for simulating the dynamical evolution of globular star clusters. In addition to twobody relaxation, many physical processes have been incorporated into FP models to achieve realistic modelling of the globular cluster evolution; these processes include tidal cutoff, binary heating, disc and bulge shocks, mass loss via stellar evolution, etc. (see Shin, Kim & Takahashi 2008 for a recent example of detailed FP modelling).
In this paper we consider the dynamical evolution of globular clusters in a steady galactic tidal field. Our main purpose is to investigate what boundary condition can give a better description of escape of stars from clusters in the tidal field. This study has been motivated by the studies by Fukushige & Heggie (2000) and Baumgardt (2001).
Fukushige & Heggie (2000) found that a large fraction of stars with energies above the escape energy (i.e. potential escapers) take much longer escape time than the dynamical time. Until their study it had been generally thought that the escape timescale is of the order of the dynamical time and that the massloss times of the clusters essentially scale with the relaxation time, which is much longer than the dynamical time. The findings of Fukushige & Heggie (2000) indicate that this simple scaling may be spoiled by potential escapers with long escape times.
In fact Baumgardt (2001) performed body simulations and showed that the massloss times (lifetimes) of clusters do not scale with the relaxation time but scale with . He concluded that the reason is that some of potential escapers are scattered back to lower energies before they leave the cluster. More recently Tanikawa & Fukushige (2005) showed that the dependence on the relaxation time changes with the strength of the tidal filed. These two studies have revealed that the behavior of potential escapers greatly influences the rate of mass loss from clusters in the tidal field.
The effects of long escape times and rescattering of potential escapers have never been considered in previous FP models in the literature, but it was assumed that escapers leave a cluster on the dynamical timescale, as is described in detail in Section 2. Since the effect of the galactic tidal field is essentially important to the cluster evolution, it is necessary to find a way to include the effect into FP models as precisely as possible.
We should mention that Takahashi & Portegies Zwart (1998, 2000) compared FP and body models of star clusters in the tidal field and found good agreement between these two theoretical models over a wide range of initial conditions. They showed that the use of anisotropic FP models with the apocentre escape criterion (Takahashi, Lee & Inagaki, 1997) and the dynamicaltime removal of escapers (Lee & Ostriker, 1987) is necessary to obtain such good agreement. However, note that in their body models the tidal force field is not included but the tidal cutoff is applied. Takahashi & Portegies Zwart (2000) confirmed that the difference between tidal cutoff and selfconsistent tidal field body models is small for a particular set of initial conditions, but did not do systematic investigations on this problem.
In this study we have devised a new scheme to treat escapers in FP models. The scheme defines a region of potential escapers in phase space and allows them to be scattered again. Comparing the results of FP models calculated with the new scheme with the results of body models, we examine the accuracy of the FP models.
2 Fokker–Planck models of star clusters in a steady tidal field
2.1 Basic assumptions
The orbitaveraged FP equation is derived under the assumption of spherical symmetry of star clusters (Cohn, 1979). Therefore the tidal field, which is not spherically symmetric, cannot be directly incorporated into orbitaveraged FP models. In FP models the effect of the tidal field is taken into account by imposing a tidal cutoff radius on the cluster, which is treated as an isolated system in other respects. Under these assumptions the distribution function of stars at time depends only on the energy of a star per unit mass, , and the angular momentum per unit mass, .
2.2 Classical treatments of escapers
First we summarise classical treatments of escapers used in FP models of previous studies.
2.2.1 Escape criteria in phase space
In previous studies, two kinds of criteria were adopted to define an escape region in space:

Energy criterion
(1) 
Apocentre criterion
(2)
where is the cluster mass and is the apocentre radius of a star having energy and angular momentum . It is assumed that a star is destined to escape once it enters into the escape region.
The apocentre criterion (Takahashi et al., 1997) is considered to be more realistic, at least as long as the tidal field is modelled as a radial cutoff, and in fact gives better agreement between FP and body models (Takahashi & Portegies Zwart, 1998, 2000). For isotropic FP models, where the distribution function does not depend on , only the energy criterion can be applied (e.g. Lee & Ostriker 1987).
2.2.2 Removal of escapers
In previous studies stars in the escape region are assumed to leave the cluster inevitably, as mentioned above. It is also assumed that the time required for this travel is of the order of the dynamical time at the tidal radius. Considering this travel time, Lee & Ostriker (1987) applied the following equation to the distribution function in the escape region:
(3) 
where is a dimensionless constant determining the efficiency of escape (see also Lee, Fahlman & Richer 1991). The timescale is an orbital timescale at the tidal radius defined by
(4) 
where is the mean mass density within the tidal radius.
Since the dynamical time is generally much smaller than the relaxation time in globular clusters, we may assume that escapers leave the cluster immediately after they enter into the escape region, when we are interested only in the evolution on the relaxation timescale. This assumption leads to the boundary condition
(5) 
on the tidal boundary (e.g. Chernoff & Weinberg 1990).
2.3 A new treatment of escapers
The boundary condition of equation (3) takes account of the fact that stars satisfying the escape criterion, i.e. potential escapers, need time to actually leave the cluster. However the effect of rescattering of potential escapers is not considered there. Here we propose a new scheme in which the rescattering effect is taken into account.
First we summarise basic assumptions and equations. Suppose that the cluster is on a circular orbit, with radius and angular velocity , round the centre of a spherical galaxy. We consider the motion of a star in the rotating coordinate system moving with the cluster; the origin is at the cluster centre, the axis points to the galactic centre, and the axis is in the cluster orbital plane. If the cluster and the galaxy are treated as point masses and () and the size of the cluster is much smaller than , there exists a conserved quantity known as the Jacobi integral given by
(6) 
(cf. Spitzer 1987, Chapt. 5). Here is the velocity of the star measured in the rotating frame, is the distance from the star to the cluster centre, and the angular velocity is given by
(7) 
The third term on the rightside in equation (6) is a combination of the centrifugal and tidal potentials.
The effective potential is defined as
(8) 
A contour plot of is shown, e.g., in Fig. 5.1 of Spitzer (1987). The effective potential has the saddle points at , where
(9) 
and
(10) 
The equipotential surface passing through these saddle points intersects with the axis at , where
(11) 
The necessary condition for escape of a star from the cluster is given by
(12) 
Note that equations (10), (11), and (12) are valid for any spherical galactic potential. Fukushige & Heggie (2000) found that the timescale for escape of stars with varies as
(13) 
With this relation in mind we have devised a new scheme to follow the evolution of potential escapers. In this scheme the evolution of the distribution function for potential escapers is described by
(14) 
where the first term on the rightside is the FP collision term and the second term represents mass loss due to escape. Here the escape timescale is given by
(15) 
where is a dimensionless numerical constant. It should be noted that energy , not the Jacobi integral , is used in equations (14) and (15). Energy does not include the centrifugal and tidal potentials. Despite this difference, we use the same critical value of energy
(16) 
where the tidal radius is identified with . One might think that using equation (15) with equation (16) is too crude an approximation, but it brings good agreement between FP and body models as is shown in Section 3.
The most important difference between equations (3) and (14) is that the latter includes the collision term. Thus equation (14) allows potential escapers to be scattered back to lower energies. The effect of mass loss is included in both equations in a similar way, though the functional forms of the escape timescale are different.
2.4 The FokkerPlanck code
The FP code used in the present study is essentially the same as that used by Takahashi & Portegies Zwart (2000), but adopts the new scheme for treating escapers described above. The code calculates the evolution of the distribution function . Unlike Takahashi & Portegies Zwart (2000), stellar evolution is not considered in the models presented in this paper. Instead the effect of heating by threebody binaries is considered in the manner described in Takahashi (1997).
For all the models presented in the present paper, 201 energy mesh points, 51 angularmomentum mesh points, and 101 radial mesh points are used. The meshes are constructed as described in Takahashi (1995). When calculating the evolution of multimass clusters, 10 discrete masscomponents are used to represent a continuous mass function.
Our FP models have two free parameters: one is in equation (15) and the other is in the Coulomb logarithm appearing in the FP collision term. How the value of is determined is described in Section 3. We set (Giersz & Heggie, 1994a) in most of our runs and (Giersz & Heggie, 1996) in a part of runs for multimass clusters.
3 Results
3.1 Comparison with body models: singlemass clusters
First we compare FP models with the full tidal field models of Baumgardt (2001) and additional body runs performed for this comparison. All the model clusters are composed of equalmass stars and move on circular orbits round a pointmass galaxy. The initial distribution of stars is given by King models (King, 1966).
Results are presented in body units, where the initial total mass and energy of a cluster are equal to 1 and , respectively, and the gravitational constant . The same units are used throughout this paper.
Here we will refer to FP models with the boundary condition of equation (14) as “FPf” models, which aim to model clusters in a selfconsistent full tidal field. FP models with equation (3) will be called “FPd” models, where stars beyond the tidal cutoff radius are removed on the dynamical timescale.
Fig. 1 compares FPf and body models concerning the evolution of the total mass of bound stars. The initial models are King models with the number of stars , 4096, 16384 and 65536. The new treatment of escapers described by equation (14) with the apocentre criterion of equation (18) is employed in the FPf models. The agreement between the FPf and body models is good in all the cases. In fact the value of the parameter in equation (14) has been determined so that good agreement is obtained by performing test runs with different values of as was done by Takahashi & Portegies Zwart (2000). We have finally chosen the value of . All the FPf models shown in Fig. 1 are calculated with this value.
Fig. 2 shows the evolution of the ratio of the mass of potential escapers to the total cluster mass for the runs shown in Fig. 1. The agreement between the FPf and body models is fairly good also in this comparison. Note that here for the FPf models is defined as the mass of stars with , although the apocentre criterion is used in the simulations. The mass of stars satisfying the apocentre criterion is smaller than that of stars with , but shows a similar trend in time variation.
(body)  (FPf)  (FPd)  

128  
256  
512  
1024  
2048  
4096  
8192  
16384  
32768  
65536  
131072  —  
262144  —  
524288  —  
1048576  —  
2097152  — 
Fig. 3 shows the halfmass time , which is the time required for a cluster to lose a half of its initial mass, as a function of the initial halfmass relaxation time . Here the halfmass relaxation time is defined by
(19) 
(Spitzer 1987, Chapt. 2) with (Giersz & Heggie, 1994a). The results are summarised also in Table 1. The FPf and body models show good agreement over the whole range of where the comparison is made. The scaling gives a reasonable fit to the results of these models as Baumgardt (2001) found.
The results of FPd models are also shown in Fig. 3. In these models the parameter is used for equation (3) (Takahashi & Portegies Zwart, 2000). The FPd models show clearly a different scaling from the other models; expect for models with very short (i.e. small ).
In Fig. 4 FPf models with the energy criterion are compared with those with the apocentre criterion as well as the body models. We have set in the energycriterion models so that their mass evolution reasonably agrees with that of the body models for small . There is no significant difference between the energycriterion models and the other models for , but the energycriterion models tend to lose mass much faster as increases. This indicates that the apocentre criterion is a better escape criterion for FPf models.
1.69  
1.57  
1.48  
1.41  
1.35  
1.32  
1.30  
1.28  
1.27 
As stated above, the results of the FPf models shown in Fig. 3 are reasonably well described by the scaling law . However we should not expect this scaling continues to hold in the limit of large . If this scaling continues, the halfmass time measured in the units of the halfmass relaxation time, , would go to zero as . This must be impossible because the mass loss is driven by twobody relaxation. In order to see the scaling of in the limit of large , we have calculated FPf models with very large , to , which are much larger than typical numbers of stars in globular clusters. The results of these models are shown in Table 2 and Fig. 5. In Fig. 5 we see that is nearly proportional to for very large clusters, say, for or . This trend is more qualitatively shown in Fig. 6, where the change in the logarithmic slope,
(20) 
is plotted. The slope approaches one as increases. The ratio for our largest models.
(body)  (FPf)  

1024  
2048  
4096  
8192  
16384  
32768  
65536  —  
131072  —  
262144  —  
524288  —  
1048576  —  
2097152  — 
We have performed simulations also for the initial conditions of King models. The halfmass times of body and FPf models for are summarised in Table 3 and are plotted in Fig. 7. Here we find good agreement again. The same parameter is used for both the and clusters. In Fig. 7 the slope of the – relation seems to be in between and 1. This point is further examined in subsection 3.3.
3.2 Dependence on the escapetime function
Baumgardt (2001) argued that the scaling can be explained by a steady state solution of a simple model for the evolution of potential escapers (see equation (12) of his paper). His model adopts the escape timescale of equation (13). If a different function is assumed for , his model predicts a different scaling law. It is shown that the scaling
(21) 
is obtained for (see Appendix A). It is interesting to see if this prediction is confirmed by the results of our FPf models.
We have performed FP runs using a generalized form of equation (15),
(22) 
with and 3. Fig. 8 plots the halfmass time against the initial halfmass relaxation time for these runs as well as for the standard runs, where King models with are used as initial conditions. The value of has been adjusted so that the nonstandard models should have roughly the same halfmass times with those of the standard ones for lower ; and for and 3, respectively.
The results of the FPf models actually depend on , but the degree of the dependence is weaker than predicted by equation (21). While this equation predicts the slopes 2/3, 3/4 and 4/5 for , 2 and 3, respectively, linear leastsquares fitting of the data in Fig. 8 gives the slopes 0.69, 0.72 and 0.75. When the fitting is done only for , the slopes are 0.75, 0.75 and 0.77. Thus the scaling law is not a bad approximation in all the cases investigated here. This is not consistent with equation (21).
3.3 Dependence on the strength of the tidal field
Tanikawa & Fukushige (2005) found that the dependence of on is affected by the strength of the tidal field and that the logarithmic slope , defined by equation (20), approaches unity as the strength of the tidal field decreases. In order to confirm their findings, we have calculated FPf models for the initial conditions where the initial tidal radius is greater than the King cutoff radius (i.e. the radius at which the density drops to zero) for each value of . On the other hand, all the models presented above are calculated for the initial conditions with .
Table 4 lists the halfmass times for King models with , 2, 4 and 6, and Fig. 9 illustrates these results. In this figure the results for and King models with are also plotted. Note that the ratio . Fig. 10 shows the variation of with .
The results shown in Figs. 9 and 10 confirm the findings of Tanikawa & Fukushige (2005). The dependence of on does depend on the strength of the tidal field. In the limit of and , it is expected that .
Note that the curve for King models with lies very close to that for King models with in each of Figs. 9 and 10. This indicates that the massloss timescale does not depend very much on the initial concentration of the cluster but is mainly determined by the strength of the tidal field, as was found by Tanikawa & Fukushige (2005).
()  ()  ()  ()  

128  
256  
512  
1024  
2048  
4096  
8192  
16384  
32768  
65536  
131072  
262144  
524288  
1048576  
2097152 
3.4 Comparison with body models: multimass clusters
So far we have concentrated on singlemass clusters. Here we consider the evolution of multimass clusters comparing our FP models with the body models of Gieles & Baumgardt (2008). They performed body simulations of clusters on circular orbits around a pointmass galaxy. In their simulations the initial mass function (IMF) is given by with the ratio . Stellar evolution is not considered in their simulations. The clusters initially have the density distribution of King models with . The ratio of the initial tidal radius to the King radius is varied from 1 to 8. The results of the simulations of Gieles & Baumgardt (2008) are summarised in their Table 1. Note that they use different notations from ours: is for the tidal (Jacobi) radius and is for the King radius.
(body)  (FPf)  (FPf)  (FPf)  
(0.11, 7)  (0.02, 7)  (0.02, 40)  
1024  
2048  
4096  
8192  
16384  
32768 
FPf models are calculated for the same initial conditions as those of Gieles & Baumgardt (2008). The results for are summarised in Table 5 and Fig. 11. There the results of the FPf models with three different sets of parameters and are reported. Giersz & Heggie (1994a) estimated the best value of for singlemass clusters by comparing body models with FP and gas models. Similarly Giersz & Heggie (1996) obtained for multimass with the IMF . We have calculated FPf models for multimass clusters using these two values of .
Fig. 11 shows that the parameter set adopted for singlemass clusters gives good fit to the body models also for multimass clusters. On the other hand the parameter set results in a clear deviation from the body models. If we stick to , the value of needs to be increased to about 40 in order to obtain good agreement with the body models. We will discuss in more detail what values of the parameters we should choose in the next section.
()  ()  ()  

1024  
2048  
4096  
8192  
16384  
32768 
()  ()  ()  

1024  
2048  
4096  
8192  
16384  
32768 
The results for the initial conditions with are shown in Tables 6 and 7 and Fig. 12. The results of Gieles & Baumgardt (2008) are not shown in these tables (see their Table 1). Fig. 12 shows that the FPf models with are in good agreement with the body models for and 4. The FPf models with are a little farther to the body models but still follow them rather well. However, for , a noticeable difference is observed between the FPf and body models; in Fig. 12 the curve for the body models is approximately linear but the slopes of the curves for the FPf models apparently change with . Neither parameter set reproduces the results of the body models as well as in the cases of . The reason for this discrepancy is not clear at present, but there is a possibility that very early corecollapse in the models with is, at least partially, responsible for it. The FPf model with and experiences core collapse (bounce) at for , and at for . The Coulomb logarithm may take different values for precollapse and postcollapse stages (see the next section), which affects the timescale of the evolution of FP models.
4 Discussion
We have shown that FP models can well follow the mass evolution of star clusters in a tidal field if a new scheme for treating potential escapers is implemented. This is the first time the effect of rescattering of potential escapers has been taken into account in FP models. Although Takahashi & Portegies Zwart (1998, 2000) showed that anisotropic FP models are in good agreement with body models for the mass evolution of star clusters in a galaxy, the tidal field is treated as a tidal cutoff rather than an actual force field. In the present study we have found that our new FP models are in good agreement with body models calculated with the inclusion of the tidal force field. Thus the new scheme has improved the accuracy of FP models.
Baumgardt (2001) argued that some potential escapers are scattered back to lower energies before they leave the cluster and that this complicates the scaling of the massloss time. The success of our models is consistent with his argument. Actually our equation for potential escapers, equation (14), can be regarded as a generalization of the equation of his toy model, his equation (12), used for explaining the scaling .
The toy model of Baumgardt (2001) is useful for giving us insight into the effect of potential escapers on the cluster evolution. On the other hand, the results presented in subsection 3.2 have revealed the limitation of the model. When the energy dependence of the escape time is artificially changed from the true one, the toy model does not correctly explain the results of our FP models. This failure of the toy model is not a big surprise, because it is only a simplified model based on many assumptions, some of which are not very realistic. For example, our simulations show that an exact steady state is never established, but the toy model assumes a steady state. In addition, the scaling of the cluster lifetime depends on the strength of the tidal field, as found by Tanikawa & Fukushige (2005) and confirmed by the present study, but the toy model does not take account of the strength of the tidal field.
Our FP models show good agreement with body models not only for singlemass clusters but also for multimass clusters. However, we have encountered a difficulty in determining proper values of the two parameters, and , in the FP models. As shown in subsection 3.4, the parameter set brings good agreement for both singlemass and multimass clusters. Since the escape timescale given by equation (15) is expected to be independent of stellar mass, it is natural that the same value of the parameter is applicable to both singlemass and multimass clusters.
On the other hand, the value of is expected to depend on the stellar mass function. Hénon (1975) argued theoretically that the value of is generally smaller in multimass clusters than in singlemass clusters. Based on the results of body simulations, Giersz & Heggie (1994a) obtained a value of for isolated singlemass clusters, and Giersz & Heggie (1996) obtained a much smaller value, , for isolated multimass clusters having an IMF similar to the IMF used in our simulations.
When we adopt the value of for multimass clusters, we have to use a much larger value of , , than the best value of for singlemass clusters, in order to obtain good agreement with body models. Thus we have not found a parameter set satisfying both the independence of on the mass function and the dependence of on it. It needs further investigation to solve this incompatibility, but even the determination of itself is not a simple task. For example, Giersz & Heggie (1994b) obtained the best value of by examining the postcollapse evolution of body models of isolated singlemass clusters. This value is much smaller than the value of obtained for precollapse singlemass clusters. These results suggest that the value of changes along with the evolution of clusters. It may also change with radius within a cluster (Giersz & Heggie, 1994a).
Fukushige & Heggie (2000) theoretically estimated not only the energy dependence of the escape timescale but also its numerical coefficient, which is given in their equation (9). If we ignore the difference between energy and the Jacobi integral , their estimate for a King model leads to a value of . This is about four times larger than our best value of for singlemass clusters. However, Fukushige & Heggie (2000) also did numerical experiments and found that their theoretical estimate of is too small; escape timescales obtained from the numerical experiments are more than a few times larger than the theoretical one. Therefore our value is not inconsistent with the result of Fukushige & Heggie (2000). On the other hand, our value of for multimass clusters with is a little larger than their theoretical estimate.
Another issue not addressed in the present paper is how the mass profile of the parent galaxy affects the results. In all the simulations presented here we assume that the parent galaxy is represented by a point mass. On the other hand, Tanikawa & Fukushige (2010) showed that the massloss timescale depends on the mass profile of the parent galaxy; the timescale increases as the mass profile gets shallower. Therefore we expect that the parameter depends on the mass profile of the parent galaxy. This issue will be examined in a future study.
5 Conclusion
In this paper we have developed new FP models of globular clusters in a steady galactic tidal field. Our FP models are novel in the method of treating escapers: potential escapers are allowed to experience gravitational scattering with other stars before they really leave clusters. The new method has been devised in order to construct more realistic models of star clusters in a tidal field compared to simple tidalcutoff models as in previous studies. The mass evolution of clusters in a tidal field does not simply scale with the relaxation time, and our FP models are in good agreement with body models in this respect.
Our FP models include two parameters and ; is the numerical factor in the Coulomb logarithm and adjusts the speed of the tidal mass loss. We have determined the best values of for given values of by comparing FP results with body results. For singlemass clusters the best parameter set is . This parameter set is applicable to multimass clusters as well, but another set does work equally well as long as multimass clusters are concerned. The parameter is expected to depend on the mass profile of the parent galaxy, though a pointmass galaxy is assumed in all the simulations of the present paper. Further investigation is required for the determination of the best values of the parameters and under various conditions.
While FP models are generally thought to be less faithful models of globular clusters than body models, the present study has significantly improved the accuracy of FP models. An advantage of FP models is that they can be calculated much faster than body models. Therefore FP models are particularly useful when we need to calculate a huge number of models. For example, when we try to specify the initial conditions of individual clusters, we have to perform simulations for many sets of the initial conditions, because the parameter space to be searched is very large. We believe that our FP models is quite useful for such searching.
Acknowledgments
Part of the work was done while the authors visited the Center for Planetary Science (CPS) in Kobe, Japan, during a visit that was funded by the HPCI Strategic Program of MEXT. We are grateful for their hospitality. HB acknowledges support by the Australian Research Council (ARC) through Future Fellowship Grant FT0991052. The numerical calculations of the FokkerPlanck models were carried out on Altix3700 and SR16000 at YITP in Kyoto University.
References
 Baumgardt (2001) Baumgardt H., 2001, MNRAS, 325, 1323
 Chernoff & Weinberg (1990) Chernoff D. F., Weinberg M. D., 1990, ApJ, 351, 121
 Cohn (1979) Cohn H., 1979, ApJ, 234, 1036
 Fukushige & Heggie (2000) Fukushige T., Heggie D. C., 2000, MNRAS, 318, 753
 Gieles & Baumgardt (2008) Gieles M., Baumgardt H., 2008, MNRAS, 389, L28
 Giersz & Heggie (1994a) Giersz M., Heggie D. C., 1994a, MNRAS, 268, 257
 Giersz & Heggie (1994b) Giersz M., Heggie D. C., 1994b, MNRAS, 270, 298
 Giersz & Heggie (1996) Giersz M., Heggie D. C., 1996, MNRAS, 279, 1037
 Heggie (2001) Heggie D. C., 2001, in Steves B. A., Maciejewski A. J., eds., The Restless Universe. Institute of Physics Publishing, Bristol, p. 109
 Hénon (1975) Hénon M., 1975, in Hayli A., ed., Proc. IAU Symp. 69, Dynamics of Stellar Systems. Reidel, Dordrecht, p. 133
 King (1966) King I., 1966, AJ, 71, 64
 Lee & Ostriker (1987) Lee H. M., Ostriker J. P., 1987, ApJ, 322, 123
 Lee et al. (1991) Lee H. M., Fahlman G. G., Richer H. B., 1991, ApJ, 366, 455
 Shin et al. (2008) Shin J., Kim S. S., Takahashi K., 2008, MNRAS, 386, L67
 Spitzer (1987) Spitzer L. Jr., 1987, Dynamical Evolution of Globular Clusters. Princeton University Press, Princeton
 Takahashi (1995) Takahashi K., 1995, PASJ, 47, 561
 Takahashi (1997) Takahashi K., 1997, PASJ, 49, 547
 Takahashi & Portegies Zwart (1998) Takahashi K., Portegies Zwart S. F., 1998, ApJ, 503, L49
 Takahashi & Portegies Zwart (2000) Takahashi K., Portegies Zwart S. F., 2000, ApJ, 535, 759
 Takahashi et al. (1997) Takahashi K., Lee H. M., Inagaki, S., 1997, MNRAS, 292, 331
 Tanikawa & Fukushige (2005) Tanikawa A., Fukushige T., 2005, PASJ, 57, 155
 Tanikawa & Fukushige (2010) Tanikawa A., Fukushige T., 2010, PASJ, 62, 1215
Appendix A Estimation of the scaling of the cluster lifetime
We follow the arguments given by Baumgardt (2001) and Heggie (2001) in order to derive the scaling law of equation (21).
Let and assume that the escape timescale has energydependence such as
(23) 
Then Baumgardt’s toy model is modified as
(24) 
where is the number of stars with energies in the range and is a constant. If we assume that the distribution of escapers is nearly in equilibrium, equation (24) shows that the width of the distribution is approximately given by
(25) 
and the number of escapers . The escape rate is estimated to be
(26) 
Therefore the scaling of the halfmass time is given by
(27) 