text:
# Formation of Massive Black Holes in Dense Star Clusters. I. Mass Segregation and Core Collapse

## 1. Introduction

### 1.1. Astrophysical Motivation

It is now well established that the centers of most galaxies host supermassive black holes (BH) with masses in the range $`M_{}10^610^9M_{}`$ (Kormendy & Gebhardt, 2001; Ferrarese et al., 2001). The evidence is particularly compelling for a BH of mass $`4\times 10^6M_{}`$ at the center of our own Galaxy (Ghez et al., 2000; Eckart et al., 2002; Schödel et al., 2002; Ghez et al., 2003). Dynamical estimates indicate that, across a wide range, the central BH mass is about 0.1% of the spheroidal component of the host galaxy (Ho, 1998). A related correlation may exist with the total gravitational mass of the host galaxy (basically the mass of its dark matter halo) Ferrarese (2002). An even tighter correlation is observed between the central velocity dispersion and the central BH mass (Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002).

Theoretical arguments and recent observations suggest that a central BH may also exist in some globular clusters (van der Marel, 2001, 2003). In particular, recent HST observations and dynamical modeling of M15 by Gerssen et al. (2002, 2003) yielded results that are consistent with the presence of a central massive BH in this cluster. Similarly Gebhardt et al. (2002) have argued for the existence of an even more massive BH at the center of the globular cluster G1 in M31. However, $`N`$-body simulations (Baumgardt et al., 2003a, b) suggest that the observations of M15 and G1 could be explained equally well by the presence of many compact objects near the center without a massive BH (van der Marel, 2001).

When the correlation between the mass of the central BH and the spheroidal component in galaxies is extrapolated to smaller stellar systems like globular clusters, the inferred BH masses are $`10^310^4M_{}`$, much larger than a $`10M_{}`$ stellar-mass BH, but much smaller than the $`10^610^9M_{}`$ of supermassive BH. Hence, these are called intermediate-mass black holes (IMBH). If some globular clusters do host a central IMBH the question arises of how these objects were formed (for recent reviews see van der Marel, 2003; Rasio et al., 2003). One natural path for their formation in any young stellar system with a high enough density is through runaway collisions and mergers of massive stars following core collapse. These runaways could easily occur in a variety of observed young star clusters such as the “young populous clusters” like the Arches and Quintuplet clusters in our Galactic center and the “super star clusters” observed in all starburst and galactic merger environments (see, e.g., Figer et al., 1999a; Gallagher & Smith, 1999). The Pistol Star in the Quintuplet cluster (Figer et al., 1998) may be the product of such a runaway, as demonstrated recently by direct $`N`$-body simulations (Portegies Zwart & McMillan, 2002). A similar process may be responsible for the formation of seed BH in proto-galactic nuclei, which could then grow by gas accretion or by merging with other IMBH formed in young star clusters (Ebisuzaki et al., 2001; Hansen & Milosavljević, 2003). Further observational evidence for IMBH in dense star clusters comes from recent Chandra and XMM-Newton observations of ultra-luminous X-ray sources, which are often (although not always) associated with young star clusters and whose high X-ray luminosities in many cases suggest a compact object mass of at least $`10^2M_{}`$ (Kaaret et al., 2001; Ebisuzaki et al., 2001; Miller et al., 2003), although beamed emission by an accreting stellar-mass BH provides an alternative explanation (King et al., 2001; Zezas & Fabbiano, 2002).

When they are born, star clusters are expected to contain many young stars with a wide range of masses, from $`0.1M_{}`$ to $`100M_{}`$, distributed according to a Salpeter-like initial mass function (IMF) (Clarke, 2003). Inspired by an early version from Rees (1984) for the formation of supermassive BH, we show in Figure 1 a diagram illustrating the two main scenarios leading to the formation of an IMBH at the center of a dense star cluster. The early stages of dynamical evolution are dominated by the stars in the upper part of the mass spectrum. Through dynamical friction, these heavy stars tend to concentrate toward the center and drive the system to core collapse. Successive collisions and mergers of the massive stars during core collapse can then lead to a runaway process and the rapid formation of a very massive object containing the entire mass of the collapsing cluster core. Although the fate of such a massive merger remnant is rather uncertain, direct “monolithic” collapse to a BH with no or little mass loss is a likely outcome, at least for sufficiently small metallicities (Heger et al., 2002). An essential condition for this runaway to occur is that the core collapse must occur before the most massive stars born in the cluster end their lives in supernova explosions (Portegies Zwart & McMillan, 2002; Rasio et al., 2003).

The accumulation at the center of a galaxy of many IMBH produced through this runaway process in nearby young star clusters (like the Arches and Quintuplet clusters in our Galaxy) provides an interesting new way of building up the mass of a central supermassive BH (Portegies Zwart & McMillan, 2002). It is possible that this process of accumulation is still ongoing in our own Galaxy (Hansen & Milosavljević, 2003). These ideas have potentially important implications for the study of supermassive BH by the Laser Interferometer Space Antenna (LISA), since the inspiral of IMBH into a supermassive BH provides the best source of low-frequency gravitational waves for direct study of strong field gravity (Cutler & Thorne, 2002).

In the alternative scenario where massive stars evolve and produce supernovae before the cluster goes into core collapse, a subsystem of stellar-mass BH will be formed (Figure 1). As demonstrated in Section 4.5, the mass loss from supernovae provides significant indirect heating of the cluster core, delaying the onset of core collapse until much later, after the stellar remnants undergo mass segregation. The final fate of a cluster with a component of stellar-mass BH remains highly uncertain. This is because realistic dynamical simulations for such clusters (containing a large number of black holes and ordinary stars with a realistic mass spectrum) have yet to be performed. For old and relatively small systems (such as Galactic globular clusters), complete evaporation is likely (with all the stellar-mass BH ejected from the cluster through 3-body and 4-body interactions in the dense core). This is expected theoretically on the basis of simple qualitative arguments (Kulkarni et al., 1993; Sigurdsson & Hernquist, 1993) and has been demonstrated recently by direct $`N`$-body simulations for very small systems containing only $`10`$ BH (Portegies Zwart & McMillan, 2000). However, for larger systems (more massive globular clusters or proto-galactic nuclei), contraction of the cluster to a highly relativistic state could again lead to successive mergers (driven by gravitational radiation) and the formation of a single massive BH (Quinlan & Shapiro, 1989; Lee, 1995, 2001). Moreover, it has been recently suggested that, if stellar-mass BH are formed with a broad mass spectrum (a likely outcome for stars of very low metallicity; see Heger et al. 2003), the most massive BH could resist ejection, even in a system with low escape velocity such as a globular cluster. These more massive BH could then grow by repeatedly forming binaries (through exchange interactions) with other BH and merging with their companions (Miller & Hamilton, 2002). However, as most interactions will probably result in the ejection of one of the lighter BH, it is unclear whether any object could grow substantially through this mechanism before running out of companions to merge with.

### 1.2. Core Collapse and the Spitzer Instability

The physics of gravothermal contraction and core collapse is by now very well understood for single-component systems (containing all equal-mass stars). In particular, the dynamical evolution of an isolated, single-component Plummer sphere to core collapse has been studied extensively and has become a testbed for all numerical codes used to compute the evolution of dense star clusters (Aarseth et al., 1974; Giersz & Spurzem, 1994). This evolution can be visualized easily using Lagrange radii, enclosing a fixed fraction of the total mass of the system see, e.g., Fig. 3 of Joshi et al., 2000 and Fig. 5 of Freitag & Benz, 2001. As the system evolves, the inner Lagrange radii contract while the outer ones expand. This so-called gravothermal contraction results from the negative heat capacity that is a common property of all gravitationally bound systems (Elson et al., 1987). In the absence of an energy source in the cluster core, and in the point-mass limit (i.e., neglecting physical collisions between stars), the contraction of the cluster core becomes self-similar and continues indefinitely. This phenomenon is known as core collapse and its universality is very well established (see, e.g., Freitag & Benz, 2001, Table 1 and references therein).

When the system contains stars with a variety of masses, the evolution to core collapse is accelerated. This has been demonstrated by the earliest $`N`$-body and Monte Carlo simulations (Aarseth, 1966; Hénon, 1971a; Wielen, 1975). To illustrate this behavior, and as a preview of the results presented later in Sec. 4, we show in Figure 2 the evolution of a cluster described initially by a simple Plummer model containing $`1.25\times 10^6`$ stars with a broad Salpeter IMF between $`m_{\mathrm{min}}=0.2M_{}`$ and $`m_{\mathrm{max}}=120M_{}`$. Core collapse occurs after $`0.1t_{\mathrm{rh}}(0)`$, where $`t_{\mathrm{rh}}(0)`$ is the initial half-mass relaxation time (see eq. (6) below). In sharp contrast, core collapse in a single-component Plummer model occurs after $`10t_{\mathrm{rh}}(0)`$ (a well-known result). Thus the presence of a broad IMF can dramatically accelerate the evolution of the cluster to core collapse.

This acceleration of the evolution to core collapse is due to the changing nature of energy transfer in the presence of a wide mass spectrum. Relaxation processes tend to establish energy equipartition (see, e.g., Binney & Tremaine, 1987, Sec. 8.4). In a cluster where the masses of the stars are nearly equal, this can be (very nearly) achieved. The core collapse is then a result of energy transfer from the inner to the outer parts of the cluster, leading to gravothermal contraction (Lynden-Bell & Wood, 1968; Larson, 1970). A large difference between the masses of the stars allows a more efficient mechanism for energy transfer. In this case, energy equipartition would tend to bring the heavier stars to lower speeds. However, as a result, the heavier stars sink to the center, where they tend to gain kinetic energy, while the lighter stars move to the outer halo. This process is called “mass segregation.” As the mass segregation proceeds, the core contracts and gets denser, leading to a shorter relaxation time, which in turn increases the rate of energy transfer from heavier to lighter stars. In typical cases this evolution eventually makes the heavier stars evolve away from equipartition (Spitzer, 1969).

The fundamental inability of the heavier stars to establish energy equipartition with the lighter stars in a system with a continuous mass spectrum is similar to the Spitzer “mass-segregation instability” in two-component clusters. Spitzer (1969), using analytic methods and a number of simplifying assumptions, determined a simple criterion for a two-component system to achieve energy equipartition in equilibrium. If the mass of the lighter (heavier) stars is $`m_1`$ ($`m_2`$) and the total mass in the light (heavy) component is $`M_1`$ ($`M_2`$), then Spitzer’s criterion can be written

$$S\left(\frac{M_2}{M_1}\right)\left(\frac{m_2}{m_1}\right)^{3/2}<0.16.$$
(1)
Watters et al. (2000), using numerical simulations, obtained a more accurate empirical condition,

$$\mathrm{\Lambda }\left(\frac{M_2}{M_1}\right)\left(\frac{m_2}{m_1}\right)^{2.4}<0.32.$$
(2)
When this stability criterion is not satisfied, energy equipartition cannot be established between heavy and light stars. Spitzer (1969) noted that the equilibrium would not be achieved for realistic mass spectra, because there is always sufficient mass in high-mass stars that, through mass segregation, they can form a subsystem (near the center of the cluster) that decouples dynamically from the lower-mass stars. This was later supported by more detailed theoretical studies (Saslaw & De Young, 1971; Vishniac, 1978; Inagaki & Saslaw, 1985).

Here we carry out a simple calculation to show that a typical system, with a Salpeter IMF (see Sec. 3.2) extending from $`0.2M_{}`$ to $`120M_{}`$, when viewed as a two-component cluster of “light stars” and “heavy stars,” does not satisfy any of the simple stability criteria. Let us separate the stars into these two components according to some arbitrary boundary: we group all stars lighter than $`m_\mathrm{b}`$ in the first component, and all stars heavier than $`m_\mathrm{b}`$ in the second component. We use $`m_1`$ and $`M_1`$ to denote the values of the average and total mass in the first component, and $`m_2`$ and $`M_2`$ similarly for the second component. We obtain, with $`m_\mathrm{b}`$ in solar mass,

$`{\displaystyle \frac{M_1}{M_2}}`$ $`={\displaystyle \frac{0.2^{0.35}m_\mathrm{b}^{0.35}}{m_\mathrm{b}^{0.35}120^{0.35}}},\mathrm{and}`$ (3)
$`{\displaystyle \frac{m_1}{m_2}}`$ $`={\displaystyle \frac{M_1}{M_2}}{\displaystyle \frac{N_2}{N_1}}={\displaystyle \frac{0.2^{0.35}m_\mathrm{b}^{0.35}}{m_\mathrm{b}^{0.35}120^{0.35}}}{\displaystyle \frac{m_\mathrm{b}^{1.35}120^{1.35}}{0.2^{1.35}m_\mathrm{b}^{1.35}}}.`$ (4)
When these values are used in equations (1) and (2) we see (Figure 3) that the stability criteria are almost never satisfied for any value of $`m_\mathrm{b}`$, except when it is extremely close to the maximum mass (as expected, since the model reduces artificially to a single-component system in this limit, with all stars having the average mass).

Vishniac (1978), devised a criterion genuinely adapted to clusters with a continuous mass spectrum, under the ad-hoc assumption that the shape of the density distribution does not depend on the stellar mass. He derives the following necessary condition for stability:

$$\beta \left(\frac{m_\mathrm{b}}{m_1}\right)^{3/2}\frac{M_2}{M_1+M_2}<1,\text{with}\beta 0.5,$$
(5)
for all $`m_\mathrm{b}`$ (his eq. 17). Figure 3 shows that this condition cannot be satisfied either111The applicability of Vishniac’s criterion appears somewhat questionable. Through FP simulations, Inagaki & Saslaw (1985) have shown that an IMF exponent $`\alpha 6.0`$ (see Sec. 3.2) is required for central equipartition to set in before core collapse while Vishniac’s criterion predicts that $`\alpha 3.5`$ is sufficient.. These results suggest that for any Salpeter-like IMF, one can always find a collection of stars that have large enough average and total mass that they will dynamically decouple from the lighter stars in the cluster. Once these decouple they form a subsystem with a much shorter relaxation time and consequently evolve to core collapse very rapidly.

The half-mass relaxation time (relaxation time at the half-mass radius $`r_h`$) for a cluster of $`N`$ stars is given by (Spitzer, 1987, eq. 2.63):

$$t_{\mathrm{rh}}(0)=\frac{0.138N}{\mathrm{ln}(\gamma _\mathrm{c}N)}\left(\frac{r_\mathrm{h}^3}{GM}\right)^{1/2}\frac{N/\mathrm{ln}(\gamma _\mathrm{c}N)}{\rho ^{1/2}},$$
(6)
where $`\rho `$ is the mass density and $`\gamma _\mathrm{c}0.01`$ in the Coulomb logarithm. Let us assume that a collapsing subsystem is formed by stars that constitute 1% of the total mass, and that all of them come from the uppermost part of the mass spectrum. For a Salpeter IMF with $`m_{\mathrm{min}}=0.2M_{}`$ and $`m_{\mathrm{max}}=120M_{}`$, the number of stars in this subsystem is then $`N_{\mathrm{sub}}10^4N`$. At the time of dynamical decoupling the central density of the subsystem must be comparable to the central density of the overall cluster. Therefore we conclude that the relaxation time is around three orders of magnitude smaller for the collapsing subsystem of heavy stars at the onset of instability. So, essentially, the heavy stars will go into core collapse as soon as they start dominating the mass density near the center.

If the cluster starts its evolution with heavy stars distributed throughout, the timescale for core collapse will be determined by their mass segregation timescale. The process of mass segregation for the heaviest stars in the cluster is driven by dynamical friction (Binney & Tremaine, 1987, Section 7.1). The timescale for mass segregation and the onset of core collapse therefore depends primarily on the mass ratio between the dominant heavy stars and the lighter background stars (Spitzer, 1969; Chernoff & Weinberg, 1990). For simple dynamical friction in an effective two-component model, one would then expect the core collapse time to be comparable to $`t_{\mathrm{DF}}(m_\mathrm{b}/m)t_{\mathrm{rh}}`$. Here $`m_\mathrm{b}`$ ($`m_{\mathrm{max}}`$, the upper mass limit of the IMF) should be the mass above which the IMF contains a large enough number of massive stars to define a “collapsing subsystem.” If this number is $`10^3`$ and the total number of stars $`N10^6`$ we get $`m_\mathrm{b}20M_{}`$ for our standard Salpeter IMF and this simple analysis would then predict the correct order of magnitude for the core collapse time ($`t_{\mathrm{cc}}/t_{\mathrm{rh}}(0)0.1`$; see Figure 2). Note, however, that this simple dynamical friction picture corresponds to one or a few massive objects traveling through a uniform background of much lighter particles, a description that provides, at best, only a rough approximation of the real situation under consideration here (see Appendix).

### 1.3. Goals of this Study

In this paper, the first of a series, we consider a wide variety of initial cluster models and we investigate the evolution of all systems until the onset of core collapse. Results from numerical simulations that include stellar collisions and track the growth of a massive object through successive mergers of massive stars during core collapse will be presented in a subsequent paper. This work is in progress and preliminary results have already been reported elsewhere (Rasio et al., 2003). Here we determine how the core collapse time is related to various initial parameters including the IMF (Sec. 3.2 and 4.2) and the central concentration and tidal radius of the cluster (Sec. 3.1 and 4.3), and we also include the possibility of initial mass segregation (Sec. 3.3 and 4.4). We then make a comparison with the stellar evolution time and derive limits on cluster initial conditions to allow the development of runaway collisions and the possible subsequent formation of a central massive BH (Sec. 5). We also derive from our calculations an estimate of the total mass of the stars that participate in the runaway collisions, which provides an upper limit to the final BH mass.

We do not include stellar evolution in our calculations since our aim here is to investigate dynamical processes taking place before even the most massive main-sequence stars in the cluster have evolved. We do, however, study the effects of mass loss from stellar winds and their dependence on metallicity (Sec. 4.5). In a subsequent paper (Gürkan & Rasio, 2003), we will study the evolution of “post-collapse” star clusters in which a central IMBH is assumed to have formed early on (through the runaway collision process). In particular, we will study the possibility that mass loss from the stellar evolution of the remaining massive stars (those that have escaped the central collapse and runaway) could disrupt the cluster (cf. Joshi et al., 2001), thereby producing a “naked” IMBH. This is motivated by observations of ultraluminous X-ray sources in regions of active star formation (e.g., in merging galaxies) containing many young star clusters, but with the X-ray sources found predominantly outside of those clusters (Zezas & Fabbiano, 2002).

An important factor that could affect significantly the dynamics of core collapse in young star clusters is the presence of primordial binaries (Fregeau et al., 2003a) or the dynamical formation of binaries through three-body interactions (Hut, 1985; Giersz, 1998). As pointed out by Inagaki (1985), the formation rate of hard three-body binaries would be accelerated in clusters with a mass spectrum (Heggie, 1975). In this first paper we do not take into account the presence of binaries in the cluster. This is justified because we do not expect binaries to play an important role until after the onset of core collapse, when the central density increases suddenly. However, collisions should become dominant immediately after the onset of core collapse. Our expectation, based on several previous theoretical studies of physical collisions during interactions of hard primordial binaries and for three-body binary formation, is that the presence of binaries in the core will in fact increase collision rates, thereby helping to trigger the runaway (Chernoff & Huang, 1996; Bacon et al., 1996; Fregeau et al., 2003b). This is also supported by the results of direct $`N`$-body simulations showing that, in smaller systems containing $`N10^4`$ single stars, collisions indeed occur predominantly through the interactions of three-body binaries formed at core collapse (Portegies Zwart & McMillan, 2002).

## 2. Numerical Methods and Summary of Previous Work

Numerical methods for investigating the dynamical evolution of star clusters include direct $`N`$-body integration, solutions of the Fokker-Planck equation by direct (finite-difference) or Monte Carlo methods, and gaseous models (for a review, see Heggie & Hut, 2003). Here we refer to direct $`N`$-body integrations simply as “$`N`$-body simulations,” direct integrations of the Fokker-Planck equation as “Fokker-Planck (FP) simulations,” and Fokker-Planck simulations based on Monte Carlo techniques as “Monte Carlo (MC) simulations.” Note, however, that our MC simulations, in which the cluster is modeled on a star-by-star basis, are in fact another type of “$`N`$-body simulations.” Each approach offers different advantages and disadvantages for understanding core collapse and massive black hole formation in star clusters with a realistic mass spectrum.

### 2.1. Summary of Previous Numerical Work

Portegies Zwart & McMillan (2002) carried out $`N`$-body simulations starting from a variety of initial conditions for clusters containing up to $`6\times 10^4`$ stars. They found that runaway collisions driven by the most massive stars can happen in sufficiently dense clusters. Their results apply directly to small star clusters containing $`10^410^5`$ stars. However, in such a small cluster, any realistic IMF typically contains only a very small number of massive stars. For example, a Salpeter IMF with minimum mass $`m_{\mathrm{min}}=0.2M_{}`$ and maximum mass $`m_{\mathrm{max}}=120M_{}`$ contains a fraction $`3\times 10^4`$ of its stars above $`60M_{}`$. The dynamical role played by massive stars can therefore depend strongly on the total number of stars in the cluster. In addition, in small systems, the dynamical evolution might be dominated by the random behavior or the initial conditions of just a few very massive stars. Consequently, to investigate runaway collisions in larger systems such as super star clusters or proto-galactic nuclei, realistically large numbers of stars must be used in numerical simulations. The computational time required for direct integration of an $`N`$-body system over one crossing time $`t_{\mathrm{dyn}}`$ scales as $`N^2`$ ($`N^{1.25}`$ on parallel machines, if one can adjust the number of CPUs to $`N`$ optimally; see Spurzem & Baumgardt 2003). Since the relaxation time is $`(N/\mathrm{ln}N)t_{\mathrm{dyn}}`$ (see Sec. 1.2), the scaling of the total CPU time of $`N`$-body simulations is nearly as steep as $`N^3`$. Currently, using a state-of-the-art GRAPE-6 board to accelerate the computations (Makino, 2001, 2002), the evolution of a cluster containing $`10^5`$ stars with $`m_{\mathrm{max}}/m_{\mathrm{min}}=1000`$ and no primordial binaries can be integrated up to core collapse in about one day (H. Baumgardt, private communication). However, including primordial binaries or increasing the number of stars would still lead to prohibitively long computation times, especially for a parameter study where a large number of integrations are required.

A wide mass spectrum also leads to increased computation times for FP simulations and gaseous models222We have been informed (H. Cohn & B. Murphy 2004, private communication) that latest isotropic FP codes require only a few seconds of computation to reach core collapse, for as many as 20 mass components.. In these methods, a continuous mass spectrum is approximated by discrete mass bins. In the gaseous model, the most time-consuming operation of the algorithm is the inversion of large matrix whose dimension is proportional to the number of equations, itself proportional to $`N_{\mathrm{comp}}`$. Hence the computing time increases like333Note that, in principle, one could split each step into separate “Poisson” and “Fokker-Planck” parts, in a way similar to what is done in FP codes, hence reducing the cost to $`T_{\mathrm{CPU}}N_{\mathrm{comp}}^2`$ (R. Spurzem, private communication). $`T_{\mathrm{CPU}}N_{\mathrm{comp}}^3`$. For FP codes, there is one diffusion term coupling each component to all others, leading to $`T_{\mathrm{CPU}}N_{\mathrm{comp}}^2`$. These steep scalings limit $`N_{\mathrm{comp}}`$ to at most $`20`$ to avoid computation times longer than a few days.

Even with a small number of mass bins, FP simulations have yielded important qualitative results. Inagaki and collaborators investigated two- and multi-component systems (Inagaki & Wiyanto, 1984; Inagaki, 1985; Inagaki & Saslaw, 1985). They found that the energy transfer within the core is an important process for determining the onset of core collapse. FP simulations by Inagaki (1985) and Chernoff & Weinberg (1990) show that a mass function with $`m_{\mathrm{max}}/m_{\mathrm{min}}1015`$ needs to be discretized into about 15 components or more to obtain an accurate value of the core collapse time. A coarser discretization leads to an artificially slow evolution. Chernoff & Weinberg (1990) successfully tracked the energy transfer between different components and demonstrated the importance of this process. Quinlan (1996) used FP simulations to follow the evolution of single- and two-component systems. He used a unique setting for his two-component clusters, with each component having a different initial spatial distribution. His aim was to study the interaction of galactic nucleus containing mainly dark (compact) objects with the bulge of the galaxy. In one model he considered a structure with “inverse initial mass segregation,” where the central nucleus is made of objects much lighter than the normal stars composing the bulge. He showed that this situation also leads to highly accelerated evolution as the more massive stars get trapped by the nucleus through dynamical friction and undergo rapid core collapse. Takahashi (1997) also investigated the evolution of clusters with a mass spectrum using FP simulations. His simulations were two-dimensional (in phase space) and therefore he was able to study the development of velocity anisotropy.

Gaseous models have the advantage of being very fast (for $`N_{\mathrm{comp}}20`$) but they include the greatest number of simplifying assumptions and so require independent checking. The good agreement shown with our MC simulation for a Plummer model with Salpeter IMF in Figure 2 is encouraging. Note, however, that the gaseous model calculation shown in Figure 2 used 50 mass components. This number of components is exceptionally high, leading to a total computation time in excess of three weeks, much longer than most MC runs. It was chosen to be able to follow core collapse to a very advanced stage, until the most massive component dominates the collapsing core. We note that the evolution of the gaseous model to core collapse is faster by some 30 % and exhibits a more gradual contraction of the inner region. The reasons for these small discrepancies will be investigated in future work, where additional comparisons between Monte Carlo and gaseous model calculations will also be presented (Freitag et al., 2003a). For the time being, suffice it to mention that for multi-mass clusters, there are basically two adjustable (dimensionless) parameters in the gaseous model, one setting the effective thermal conductivity and the other the timescale for energy exchange between components ($`\lambda `$ and $`\lambda _{\mathrm{eq}}`$, see Louis & Spurzem 1991; Giersz & Spurzem 1994; Spurzem & Takahashi 1995). For the model plotted in Figure 2, we used the standard values of these parameters ($`\lambda =0.4977`$ and $`\lambda _{\mathrm{eq}}=1`$) which were established for clusters with a different structure and mass spectrum and some adjustment (preferentially through comparisons with $`N`$-body runs) may be required.

A direct comparison between FP and gaseous models was also carried out by Spurzem & Takahashi (1995), for two-component clusters, and also resulted in good agreement. Recently, a hybrid code has been developed by Giersz & Spurzem (2003), combining a gaseous model with MC techniques. In this approach the single stars are represented by the gaseous model, while primordial binaries are followed with a Monte Carlo treatment. A similar hybrid treatment could be applied to the problem we are studying here, but with massive (single and binary) stars included in the Monte Carlo component, and lower-mass stars represented by a gaseous model.

### 2.2. The Monte Carlo Code

The solution of the Fokker-Planck equation with an orbit-averaged Monte Carlo method is an ideal compromise for the problem at hand. It can handle a suitably large number of stars and a wide mass spectrum can be implemented with very little additional difficulty. Most importantly, as in direct $`N`$-body integrations, MC simulations can implement a star-by-star description of the cluster. This allows the inclusion of many important processes such as collisions, binary interactions (including primordial binaries), stellar evolution (and the accompanying mass loss), as well as the effects of a massive central object with relative ease and much higher realism compared to direct FP simulations or gaseous models. We will incorporate the effects of all these additional processes on cluster evolution in the subsequent papers of this series.

The MC code we have used to obtain the main results of this paper is described in detail by Joshi et al. (2000). It is based on the ideas of Hénon (1971a, b, 1973) and in many respects it is very similar to MC codes developed by Stodółkiewicz (1982, 1986), Giersz (1998, 2001), and Freitag & Benz (2001). In the rest of this section, we give a brief summary of the numerical method.

The main simplifying assumption is the Fokker-Planck approximation, in which relaxation processes are assumed to be dominated by small distant encounters rather than strong encounters with large deflections (Spitzer, 1987; Binney & Tremaine, 1987). The dynamical evolution can then be treated as a diffusion process in phase space. Following the individual interactions between the stars, as in direct $`N`$-body simulations, is computationally expensive. However, the average cumulative effect on a star in a given amount of time can be characterized by diffusion coefficients in phase space. To compute the relaxation, the timestep for the numerical evolution has to be chosen smaller than the relaxation time. Since the relaxation time is normally shortest at the center, we choose our timestep to be a fraction of the central relaxation time. This ensures that the relaxation is followed accurately throughout the cluster.

Another important simplification in the MC method is the assumption of spherical symmetry. The position of the particles is represented by a single radial coordinate, $`r`$, and the velocity is represented by radial, $`v_r`$, and tangential, $`v_t`$, components. The specific angular momentum, $`A`$, and the specific energy, $`E`$, of a star with index $`i`$, are given by:

$$A_i=v_{t,i}r_i,E_i=U(r_i)+\frac{1}{2}(v_{r,i}^2+v_{t,i}^2),$$
(7)
where $`U(r)`$ is the gravitational potential at a given point. The assumption of spherical symmetry implies that the potential at all points can be computed in a time proportional to the number of particles $`N`$, rather than $`N^2`$.

In Hénon’s algorithm, the evolution is simulated by reproducing the effect of the cluster on each star by a single effective encounter in each timestep. At every iteration the two integrals of the motion $`A_i`$ and $`E_i`$ characterizing the orbit of each star are perturbed in a way that is consistent with the value of the diffusion coefficients (Hénon, 1973). To conserve energy, this perturbation is realized by a single effective scattering between two neighboring stars. The square of the scattering angle, $`\beta `$, is proportional to the timestep chosen. We choose our timestep such that $`\mathrm{sin}^2(\beta /2)0.05`$ at the center of the cluster. Choosing too large a timestep will lead to a saturation effect and the relaxation will proceed artificially slowly. Using too small a timestep, on the other hand, not only increases the computation time, but also can lead to spurious relaxation (see the Appendix for a discussion of spurious relaxation effects).

After the perturbation of $`A_i`$ and $`E_i`$, the stars are placed at random positions between their apocenters and pericenters using a probability distribution that is proportional to the time spent at a given location on their new orbit. For a star with index $`i`$, the apocenter and pericenter distances are calculated by finding the roots of

$$2E_i2U(r)\frac{A_i^2}{r^2}=0.$$
(8)
This random placement is justified by the assumption of dynamical equilibrium, i.e., the evolution of the system does not take place on the crossing (or dynamical) timescale, but rather on the relaxation timescale. The only important point about assigning a specific position to a particle on its orbit is that its contribution to density, potential, interaction rates, etc., has to be estimated correctly. After all stars have been placed at their new positions, the potential is recalculated and the whole cycle of perturbation is repeated.

This method can be modified so that the timestep is a fraction of the local relaxation time (Hénon, 1973; Freitag & Benz, 2001). Stodółkiewicz (1982, 1986) and Giersz (1998, 2001) divided the system into zones resulting in an approach intermediate between a fixed and smoothly varying timestep. Dividing the system into radial zones also allows the MC method to be parallelized efficiently for use on multi-processor machines (Joshi et al., 2000). Another possible modification uses a scaling of the units such that each particle in the simulation can represent an entire spherical shell of many identical stars rather than a single star (Hénon, 1971a; Freitag & Benz, 2001, 2002). It should also be noted that, although the effects of strong encounters between stars on relaxation are assumed to be negligible compared to weak, more distant encounters, they can be incorporated by estimating their rate of occurrence in a way similar to physical collisions (Freitag & Benz, 2002) or interactions with binary stars (Fregeau et al., 2003a).

Our MC code has been used previously to study many fundamental dynamical processes such as the Spitzer instability (Watters et al., 2000) and mass segregation (Fregeau et al., 2002) for simple two-component systems, as well as the evolution of systems with a continuous but fairly narrow mass spectrum of evolving stars (Joshi et al., 2001). A difficulty introduced by a broad continuous IMF (with a large $`m_{\mathrm{max}}/m_{\mathrm{min}}`$ ratio) is the necessity of adjusting the timestep to treat correctly encounters between stars of very different masses. When pairs of stars are selected to undergo an effective hyperbolic encounter as described above, one has to make sure that the deflection angle remains small for both stars. In situations where the mass ratio of the pair can be extreme, one has to decrease the timestep accordingly (Stodółkiewicz, 1982). In practice, for the simulations described here, we find that the timestep has to be reduced by a factor of up to $`500`$ compared to what would be appropriate for a cluster of equal-mass stars. We discuss further the applicability of orbit-averaged MC methods to systems with a continuous mass spectrum and large $`m_{\mathrm{max}}/m_{\mathrm{min}}`$ ratio in the Appendix.

## 3. Initial Conditions and Units

The characteristics of core collapse and the subsequent runaway collisions depend on the initial conditions for the cluster. These initial conditions include the total number of stars, the IMF, the initial spatial distribution of the stars (density profile and, possibly, initial mass segregation), and the position of the cluster in the galaxy, which determines the tidal boundary. As we shall see, the most important initial parameters are the slope of the IMF, the maximum stellar mass, and the initial degree of central concentration of the cluster density profile. We have used a wide variety of initial conditions, both to test the robustness of our findings and to establish the dependence of our results on these parameters. As a typical reference model we use an isolated Plummer sphere with a Salpeter IMF and stellar masses ranging from $`m_{\mathrm{min}}=0.2M_{}`$ to $`m_{\mathrm{max}}=120M_{}`$. We then explore variations on this model by changing the initial cluster structure, the IMF or the number of stars.

### 3.1. Density Profile

We have examined three families of models: Plummer and King models (Binney & Tremaine, 1987, Sections 2.2 and 4.4), which have a core-halo structure, and $`\gamma `$-models (Dehnen, 1993; Tremaine et al., 1994), which have a cusp near the center. All these models have a characteristic radius given by

$`a_\mathrm{P}`$ $`=(2^{2/3}1)^{1/2}r_\mathrm{h}0.766r_\mathrm{h},`$ (9)
$`a_\mathrm{K}`$ $`=\sqrt{{\displaystyle \frac{9\sigma ^2}{4\pi G\rho _0}}},`$ (10)
$`a_\gamma `$ $`=(2^{1/(3\gamma )}1)r_\mathrm{h},`$ (11)
for Plummer, King, and $`\gamma `$-models respectively. In these formulae, $`r_\mathrm{h}`$ is the half-mass radius, $`\rho _0`$ is the central density, and $`\sigma `$ is a King model parameter. We show the density profiles corresponding to these various models in Figure 4. Here $`W_0`$ is the dimensionless central potential, related to the concentration parameter (Binney & Tremaine, 1987, Fig. 4-10). Other useful quantities characterizing the various density profiles are given in Table 1. The initial model used for each of our MC simulations is listed in the second column of Table 2. A general procedure for producing these models is given by Freitag & Benz (2002); a less general but simpler procedure for the Plummer model is given by Aarseth et al. (1974).

### 3.2. Initial Mass Function

Since we expect the conditions leading to core collapse to depend sensitively on the IMF, we have carried out simulations with a wide range of IMF parameters. However, it is generally established that, at least for the high-mass end of the spectrum, a universal IMF very close to the simple Salpeter-like power-law is obeyed. This is indicated both by observations and by theoretical calculations (Kroupa, 2002; Clarke, 2003, and references therein).

In our simulations we assign the masses of individual stars using a sampling procedure. We first choose a random number, $`X`$, from a uniform distribution between 0 and 1. For a simple power-law IMF with $`dNm^\alpha dm`$ between $`m_{\mathrm{min}}`$ and $`m_{\mathrm{max}}`$ we calculate the corresponding mass using

$$m(X)=m_{\mathrm{min}}\left(1+X\left[\left(\frac{m_{\mathrm{max}}}{m_{\mathrm{min}}}\right)^{1\alpha }1\right]\right)^{1/(1\alpha )},$$
(12)
where the value $`\alpha =2.35`$ would correspond to a Salpeter IMF.

In addition to simple power laws, for two models we have used the Miller-Scalo (1979) and Kroupa (Kroupa et al., 1993) IMFs, which are steeper at the high-mass end of the spectrum and shallower at the low-mass end. Both of these IMFs can be represented by broken power laws. For the Miller-Scalo IMF,

$$\frac{dN}{dm}\{\begin{array}{cc}m^{1.4}\hfill & 0.1m1.0\hfill \\ m^{2.5}\hfill & 1.0m10.0\hfill \\ m^{3.3}\hfill & 10.0m\hfill \end{array},$$
(13)
and for the Kroupa IMF,

$$\frac{dN}{dm}\{\begin{array}{cc}m^{1.3}\hfill & 0.08m0.5\hfill \\ m^{2.2}\hfill & 0.5m1.0\hfill \\ m^{2.7}\hfill & 1.0m\hfill \end{array},$$
(14)
where all numerical values are in solar mass.

However, for the sake of computational convenience, we prefer to implement a somewhat different parametrization of these distributions. To generate the mass spectra corresponding to these IMFs we use the functions (Eggleton et al., 1989; Kroupa et al., 1993)

$$m(X)=\frac{0.19X}{(1X)^{0.75}+0.032(1X)^{0.25}}$$
(15)
and

$$m(X)=0.08+\frac{0.19X^{1.55}+0.05X^{0.6}}{(1X)^{0.58}},$$
(16)
for the Miller-Scalo and Kroupa mass functions, respectively. In both cases if $`m(X)<m_{\mathrm{min}}`$ or $`m(X)>m_{\mathrm{max}}`$ the result is discarded and a new random number is generated. The values in equation (16) correspond to the choice of $`\alpha _1=1.3`$ in Table 10 of Kroupa et al. (1993).

In most of our simulations we have used $`N=1.25\times 10^6`$ stars. For all mass functions, this implies the presence of many massive stars, allowing us to resolve fully the higher end of the mass spectrum. For example, for $`m_{\mathrm{min}}=0.2M_{}`$ and $`m_{\mathrm{max}}=120M_{}`$, with $`N=1.25\times 10^6`$, we have $`60`$ stars with $`m>100M_{}`$ for a Salpeter IMF. The results from our simulations are therefore not much affected by random fluctuations in a small number of very massive stars. Note, however, that in much smaller systems, containing perhaps only $`10^4`$ stars (as in the Arches and Quintuplet clusters near our Galactic center), these small number effects and random fluctuations may indeed play a dominant role in determining the dynamical fate of the few most massive stars in the system.

### 3.3. Initial Mass Segregation

Initial mass segregation in star clusters (i.e., the tendency for more massive stars to be formed preferentially near the cluster center) is expected to result from star formation feedback in dense gas clouds (Murray & Lin, 1996) or from competitive gas accretion onto proto-stars and mergers between them (Bonnell et al., 2001; Bonnell & Bate, 2002). There is also some observational evidence for initial mass segregation in both open and globular clusters (Bonnell & Davies, 1998; Raboud & Mermilliod, 1998; de Grijs et al., 2003).

We have considered the possibility of initial mass segregation in a few of our MC simulations. We adopt a simple prescription whereby we increase the average stellar mass within a certain radius $`r_{\mathrm{ms}}`$. For $`r<r_{\mathrm{ms}}`$, rather than sampling from an IMF with fixed $`m_{\mathrm{min}}`$, we randomly choose between two values of $`m_{\mathrm{min}}`$, one that is used for the outer part of the cluster and another one that is larger. For $`r>r_{\mathrm{ms}}`$, we follow a similar procedure, this time changing $`m_{\mathrm{max}}`$. The average mass within $`r_{\mathrm{ms}}`$ is larger by a factor $`C_{\mathrm{ms}}`$ and outside $`r_{\mathrm{ms}}`$ is smaller by a factor $`C_{\mathrm{ms}}^{}`$, with respect to a cluster without initial mass segregation. We choose $`C_{\mathrm{ms}}`$ and $`C_{\mathrm{ms}}^{}`$ such that the overall average stellar mass in the cluster does not change with these modifications. Changing the average stellar mass within any region of the cluster would of course in general leave the system out of dynamical equilibrium. To maintain virial equilibrium, the mass density profile must also be preserved. We achieve this by modifying the number density of stars appropriately.

Our initial conditions for models with initial mass segregation are summarized in Table 3. Here $`q`$ is the initial cluster mass fraction within $`r_{\mathrm{ms}}`$. This implies $`C_{\mathrm{ms}}^{}=(1q)/(1qC_{\mathrm{ms}})`$. Increasing $`q`$ or $`C_{\mathrm{ms}}`$ represents more extended or more pronounced initial mass segregation.

Our prescription for initial mass segregation allows the formation of massive stars in the outer parts of the cluster, as well as the formation of lighter stars in the inner parts, but the more massive stars are more likely to be found near the center. This makes our approach different from that of Bonnell & Davies (1998), who put all massive stars closer to center.

### 3.4. Units

For all our numerical calculations, we adopt the standard $`N`$-body units (Hénon, 1971a): we set the initial total cluster mass $`M=1`$, the gravitational constant $`G=1`$, and the initial total energy $`E_0=1/4`$. In the tabulation of the results we also use the initial half-mass relaxation time $`t_{\mathrm{rh}}(0)`$, given by equation (6), as the unit of time for comparison with other work in the literature.

The conversion to physical units is done by evaluating the initial half-mass relaxation time in years. For example, for the Plummer model, we can write

$$\begin{array}{c}t_{\mathrm{rh}}(0)330\mathrm{Myr}\hfill \\ \hfill \times \left(\frac{N/10^6}{\mathrm{ln}(\gamma _cN)/\mathrm{ln}(10^4)}\right)\left(\frac{a_\mathrm{P}}{1\mathrm{pc}}\right)^{3/2}\left(\frac{M}{10^6M_{}}\right)^{1/2}.\end{array}$$
(17)
Similar expressions can also be obtained for King models and $`\gamma `$-models by use of the quantities in Table 1. For a single-component model (containing equal-mass stars) the value of $`\gamma _\mathrm{c}`$ in the Coulomb logarithm can be calculated theoretically (Farouki & Salpeter, 1994), but for a system with a wide mass spectrum it must be determined by comparing to direct $`N`$-body integrations. Giersz & Heggie (1996, 1997) carried out such comparisons and found $`\gamma _\mathrm{c}0.015`$. Our own comparison with a recent $`N`$-body result for $`10^5`$ stars with a wide mass spectrum led us to adopt the value $`\gamma _\mathrm{c}=0.01`$ (H. Baumgardt, private communication).

For processes occurring in the central parts of the cluster, the central relaxation time is a more relevant quantity,

$$t_{\mathrm{rc}}(0)\frac{\sigma _{3\mathrm{D}}^3}{4.88\pi G^2\mathrm{ln}(\gamma _\mathrm{c}N)nm^2},$$
(18)
where $`\sigma _{3\mathrm{D}}`$, $`n`$ and $`m`$ are the 3D velocity dispersion, number density and average stellar mass at the cluster center (Spitzer, 1987, eq. 3.37). The $`N`$-body unit system only specifies unambiguously dynamical times, which are independent of $`N`$. In this system, relaxation times are proportional to $`N/\mathrm{ln}(\gamma _\mathrm{c}N)`$. It is therefore useful to define also the so-called “Fokker-Planck time unit,” which is the $`N`$-body time unit ($`t_{\mathrm{dyn}}GM^{5/2}(4E_0)^{3/2}`$) multiplied by $`N/\mathrm{ln}(\gamma _\mathrm{c}N)`$.

The most important physical properties of all our initial cluster models, including $`t_{\mathrm{rh}}(0)`$ and $`t_{\mathrm{rc}}(0)`$, are given in Table 1 (in $`N`$-body and FP units).

## 4. Results

The initial conditions and main results of all our MC simulations are summarized in Tables 2 and 3. All models have initially $`N=1.25\times 10^6`$ stars, except for Models 2r, s, b, c, and d, which have varying $`N`$ between $`3\times 10^5`$ and $`10^7`$. The maximum value of $`N`$ is set in practice by the available computer memory and $`N10^7`$ corresponds to 2 GB of available memory for our code. The run for Model 2d took about two weeks of CPU time to complete on a 2.8 GHz Pentium 4 Linux workstation. More typical runs for $`N=1.25\times 10^6`$ took around $`2030`$ CPU hours. The close agreement between the outcomes of Models 2r – 2d confirms the expectation that our results should be independent of the number of stars in the system, at least for sufficiently large $`N`$ to avoid small-number effects in the cluster core.

### 4.1. Mass Segregation and Core Collapse

As expected, all models with a Salpeter-like IMF and a wide mass spectrum undergo core collapse considerably faster than any single-component cluster (cf. Sec. 1.2). In Figure 5, we show the evolution of various Lagrange radii, as well as the average stellar mass inside these radii, for our reference model (Model 2; same as shown in Figure 2). In contrast to the evolution of a single-component model, the inner Lagrange radii remain almost constant until the very abrupt onset of core collapse at $`t/t_{\mathrm{rh}}(0)0.07`$. A more detailed view of core collapse is shown in Figure 6. Here the time axis has been replaced by the central potential depth, which increases monotonically in time and provides a natural stretch near core collapse. Note that the core collapse time for a particular run can be determined very accurately, to within $`0.01\%`$. However, random fluctuations in the realization of each initial condition for a particular model lead to a much larger physical “error bar” on $`t_{\mathrm{cc}}`$ (see below).

The rapid increase of the average stellar mass inside the innermost Lagrange radii seen in Figures 5 and 6 is an indication of significant mass segregation. In fact, very significant mass segregation takes place throughout the evolution of the system. This can be seen in the upper panel of Figure 7, which shows the evolution of half-mass radii for stars in various mass bins (compare with Fig. 1 of Spitzer & Shull, 1975). It is clear that the rate of mass segregation in each mass bin, measured by the slope $`\alpha _{\mathrm{ms}}`$ of the corresponding half-mass radius $`r(t)`$ in that bin, is very nearly constant from $`t=0`$ all the way to core collapse (the least-square straight-line fits, constrained to $`r/r_\mathrm{h}(0)=1`$ at $`t=0`$, are shown in the figure). This rate can be positive or negative. For Model 2 we find that stars more massive than about $`5M_{}`$ drift inward on average, while less massive stars drift outward. Even for the most massive stars, we do not find that the mass segregation rate is proportional to the average mass in the bin (in contrast to what would be expected for a tracer population of massive stars driven by simple dynamical friction; cf. Sec. 1.2 and Fregeau et al. 2002). Instead, the following simple expression provides a good fit (to within a few percent) to the observed mass-dependence for Model 2:

$$\alpha _{\mathrm{ms}}t_{\mathrm{rh}}(0)=\alpha _0\mathrm{exp}(m/m_\mathrm{f})+\alpha _2,$$
(19)
where $`m`$ is the average stellar mass and the best-fit parameters are $`\alpha _0=9.45`$, $`m_\mathrm{f}=21.9M_{}`$, and $`\alpha _2=8.07`$. Note that the mass segregation rate actually approaches a constant for large $`m`$ (but of course it is unphysical to extrapolate beyond $`m_{\mathrm{max}}`$). In the lower panel of Figure 7 we show the average mass within $`r_\mathrm{h}`$, $`r_\mathrm{h}/2`$ and $`r_\mathrm{h}/4`$. The steady increase of the average mass in each region is further indication that the mass segregation not only starts immediately, but also continues until core collapse. For a smaller number of stars, the average mass within a given radius can reach saturation (cf. Figure 2 of Bonnell & Davies, 1998). Our results for larger $`N`$ do not show this saturation.

Initially and throughout the evolution until core collapse, the cluster maintains a core-halo structure. However, when the heaviest stars start dominating the core and the Spitzer instability occurs, this initial structure is lost. We demonstrate this behavior in Figure 8, where the evolution of the gravitational potential profile is shown. We have checked that the final structure and, in particular, the formation of a cusp are not dominated by small number effects but are instead the result of many massive stars participating in the core collapse. We illustrate this in Figure 8 where close to 200 innermost stars, which constitute 0.2% of the total mass, are shown explicitly.

Another important point, which can be deduced from the evolution of the Lagrange radii, as well as from Figure 8 (see the lower two curves), is that the outer parts of the cluster are not much affected by mass segregation and core collapse, except for the disappearance of the most massive stars.

As expected, the stars in the collapsing core have a mass distribution much richer in heavy stars than the IMF. We compare the mass function at core collapse with the IMF, for $`m1M_{}`$, in Figure 9. The mass function shown in this figure is obtained by choosing the innermost 200 stars at three different times all before and within $`5\times 10^6t_{\mathrm{rh}}(0)`$ of core collapse time.

Mass segregation continues until the core is dominated by massive stars. At this point, the massive stars near the center dynamically decouple from the rest of the cluster and go into collapse as a separate subsystem. The stars outside this subsystem act as an energy sink, and, as a result of their energy gain, expand away from the center of the cluster. In a single-component system, as the collapse proceeds, the subset of stars that participate in the collapse becomes smaller and smaller. Here instead the stars that participate in the collapse at the edge of the core are heavier than the surrounding stars so they can more efficiently give their energy away and remain together. As a result, in systems with a wide mass spectrum and a steep enough IMF, there is a clear separation between stars in the collapsing core and stars outside the collapsing core, almost suggesting a real condensation process. The total mass in the collapsing stars is a crucial property of these systems, clearly representing an upper limit to the mass of any BH that could form eventually through the gravitational collapse of a runaway merger remnant.

To estimate the total mass in the collapsing core, $`M_{\mathrm{cc}}`$, we proceed as follows. At every timestep we calculate and record the Lagrange radii for various mass fractions. For each Lagrange radius, this provides typically a few thousand data points per run. The innermost radii exhibit a great amount of noise for the lowest mass fractions. We remove this noise in two steps. First, we take arithmetic averages over 60 points and reduce the total number of points to $`50100`$. Then for each point and using the closest 6 points we fit a cubic polynomial using least squares and we evaluate this polynomial at the corresponding time. This way we both smooth the Lagrange radii further and estimate their derivatives. When the system is going into core collapse, a decrease in the collapse rate is an indication that the results provided by our code are beginning to be dominated by numerical errors. We stop the simulation when this happens and define the core collapse time as the last point before this behavior is observed. We then find the innermost Lagrange radius that has a positive derivative at the time of core collapse and estimate the mass of the stars that participate in core collapse as the mass enclosed by this radius. An investigation of Figure 6 by eye verifies that this method produces reasonable results. In practice, we track many more Lagrange radii than shown in this figure, at intervals of 0.02% between 0.08% and 0.36%.

Because of small-number fluctuations (and, in particular, the intrinsic noise in the innermost Lagrange radii) there is a statistical uncertainty in all the numbers quoted in Table 2. This uncertainty is not unphysical, as real systems will be affected similarly by fluctuations in the number and specific properties of a relatively small number of massive stars. To test the robustness of our results and estimate this uncertainty we have repeated 10–20 simulations with different random seeds for Model 2 and some of its variants with different numbers of stars (Models 2r, 2s, and 2b). For these models the numbers in Table 2 are obtained by averaging the results of the various (physically equivalent) runs. The standard deviations obtained for the values of $`t_{\mathrm{cc}}/t_{\mathrm{rh}}(0)`$ and $`M_{\mathrm{cc}}/M_{\mathrm{tot}}`$ for our typical Model 2 are about 5% and 20%, respectively. The larger uncertainty in $`M_{\mathrm{cc}}/M_{\mathrm{tot}}`$ is a result of the high sensitivity of this quantity to noise in the innermost Lagrange radii.

We cannot study the evolution of the cluster past core collapse without treating in detail the dynamics of the central stars during collapse, which is beyond the scope of this paper. In future work we will include a detailed treatment of stellar collisions in the core (PaperII, Freitag et al. 2003b). Alternatively, one could use a hybrid method and treat the central part of the cluster with a direct $`N`$-body approach (Lightman & McMillan, 1985). However, the simplest approximation may be to introduce an effective boundary condition at some very small radius. Stodółkiewicz (1982), in his model B2, took this approach, but his implementation was not conserving energy.

### 4.2. Dependence on the IMF

The sharp onset of core collapse for our reference model, shown in Figure 5, is a result of the Spitzer instability. This instability is driven by the segregation of the heaviest stars toward the center and their dynamical decoupling from the rest of the system. As indicated by the results of Section 4.1, the ratio of maximum to average stellar mass in the IMF, $`m_{\mathrm{max}}/m`$, is an important parameter setting the timescale for the onset of instability. There are various ways to study the dependence of our results on this parameter. One can use an IMF different from Salpeter, e.g., Miller-Scalo or Kroupa. In addition, when using a power-law IMF, changing the slope $`\alpha `$ or the maximum mass $`m_{\mathrm{max}}`$ will obviously alter the value of the ratio $`m_{\mathrm{max}}/m`$.

Our results from simulations for a large number of models exploring these various alternative IMF are presented in Figure 10. Note that for $`m_{\mathrm{max}}/m<40`$ our results suggest a relation $`t_{\mathrm{cc}}(m_{\mathrm{max}}/m)^{1.3}`$, also obeyed by all Fokker-Planck models from Inagaki (1985) and Takahashi (1997) with $`\alpha <3.5`$ and a number of mass components sufficient to ensure proper sampling of the IMF. However our computations extend to much higher values of $`m_{\mathrm{max}}/m`$ than those works. Beyond $`m_{\mathrm{max}}/m50`$, a domain reached by any realistic IMF, the core collapse time approaches a constant $`0.15t_{\mathrm{rc}}(0)`$. Therefore our main conclusions appear to be independent of the details of the IMF as long as the number of massive stars in the system is large enough.

For small values of $`m_{\mathrm{max}}/m`$, not only the timescale but also the very nature of the collapse changes. In these systems the evolution timescale, i.e., the relaxation time, of the subsystem that can decouple from the rest of the cluster is no longer small enough that the evolution of lighter stars can be neglected. The Lagrange radii for these models behave similarly to the single-component case, i.e., there is no clear separation between collapsing and expanding Lagrange radii. Consequently, we do not give values of $`M_{\mathrm{cc}}/M_{\mathrm{tot}}`$ for these models in Table 2. Quinlan (1996), using FP simulations, also found that, for a clear decoupling, the relaxation time of the subsystem has to be short compared to that of the other stars.

### 4.3. Dependence on Initial Cluster Concentration

One expects the core collapse time to depend strongly on the initial density profile of the cluster, and, in particular, on the central concentration. A simple and systematic way to examine this dependence is to use different King models with varying concentration parameter or $`W_0`$ (Binney & Tremaine, 1987, Fig. 4-10). For single-component clusters this has been done using a variety of methods (Quinlan, 1996; Einsel, 1996; Joshi et al., 2001; Kim et al., 2002).

We have carried out a number of simulations using King models with $`W_0=19`$ as initial structure. Our simulations include both isolated clusters (with no tidal boundary enforced, even though the initial models are truncated) and clusters with a tidal boundary (assuming a circular orbit in a spherical galactic potential). This tidal boundary is initially chosen to be at the tidal radius of the King model and then adjusted as the cluster loses mass (Chernoff & Weinberg, 1990; Joshi et al., 2001). In Figure 11a we plot the ratio of core collapse time to initial half-mass relaxation time for our simulations, along with results for single-component clusters. The ratio $`t_{\mathrm{cc}}/t_{\mathrm{rh}}(0)`$ is much smaller for systems with mass spectrum, as would be expected. We also see that the core collapse time for these systems is independent of the presence of a tidal boundary, as expected when core collapse is driven by local processes within the core rather than by global energy transfer. This idea is strongly supported by the results in Figure 11b, which shows the comparison of the ratio of core collapse time to central relaxation time, for our simulations and single-component models. With respect to single-component models, the ratio $`t_{\mathrm{cc}}/t_{\mathrm{rc}}(0)`$ for systems with a broad mass spectrum shows very little variation.

We have also run two simulations for $`\gamma `$-models, using $`\gamma =1.0`$ (Hernquist model), and $`\gamma =1.5`$. Initially, models with $`\gamma <2`$ have vanishing central velocity dispersion and, for $`\gamma 2`$, they exhibit a central “temperature inversion” (Tremaine et al., 1994). For single-component clusters, the central region at first undergoes rapid gravothermal expansion until it becomes isothermal and “normal” core collapse can start (Quinlan, 1996). Models with $`\gamma 2`$ also have zero initial central relaxation time444Models with $`2.0<\gamma <3.0`$ have infinite $`t_{\mathrm{rc}}(0)`$ because their velocity dispersion rises like $`1/r`$ near the center (Tremaine et al., 1994).. Our results for King models would then suggest that core collapse induced by central mass segregation should be extremely fast, if it were not for the opposite effect of the temperature inversion. For $`\gamma =2`$, mass segregation seems to have the upper hand in the competition and it proves impossible to resolve core collapse because extreme mass segregation appears nearly instantaneously during the MC runs. Models with $`\gamma =1`$ and $`\gamma =1.5`$ expand at first and then evolve to core collapse very rapidly. The core collapse time and collapsing core mass are very hard to determine numerically for these models. In Table 2, we give only an upper limit on their core collapse times. Obviously, these values are so short that their physical relevance is unclear. It is hard to imagine through which process such a cluster could be created if one wants it to be virialized but without initial mass segregation because these conditions impose constraints on the formation timescale unlikely to be satisfied in real systems. In addition, note that the local dynamical time in a $`\gamma =1`$ model approaches a constant non-zero value near the center, while, formally, the relaxation time goes to zero there. This suggests that such models can only provide an approximate description of real clusters, where finite-number effects will play a key role near the center, so that the question of determining the evolution in the limit of very large $`N`$ is not well posed.

### 4.4. Initial Mass Segregation

Naturally we expect that any initial mass segregation in a cluster should lead to an even shorter core collapse time as the heavier stars are starting their life closer to the center on average, and therefore do not need as much time to concentrate there through mass segregation.

Our results from MC simulations with initial mass segregation (Sec. 3.3) are presented in Table 3 and Figure 12. As expected, we find that $`t_{\mathrm{cc}}/t_{\mathrm{rh}}(0)`$ decreases with increasing $`C_{\mathrm{ms}}`$ (i.e., when stars become more and more massive on average in the inner region), and with increasing $`q`$ (i.e., when the increase in average stellar mass takes place over a larger central region). However, the values of $`M_{\mathrm{cc}}/M_{\mathrm{tot}}`$ do not appear to be affected significantly (within expected random fluctuations from run to run). The implication for BH formation through runaway collisions is that the final BH mass may not be affected by initial mass segregation, while the condition for runaway growth to occur (see Sec. 5) will be relaxed as the core collapse takes place even earlier.

### 4.5. Effects of Stellar Evolution

The simulations described so far consider only the dynamical evolution of the cluster, neglecting the stellar evolution entirely. We have also carried out a number of simulations to understand how stellar evolution and the accompanying mass loss can modify the dynamical evolution. The stellar evolution treatment we have adopted is that of Belczynski et al. (2002), based on the approximations of Hurley et al. (2000).

Stellar evolution introduces a new physical clock in the system, independent of relaxation. Therefore it is necessary to specify the physical scale of an initial cluster model (e.g., the half-mass radius in pc) before starting a simulation so that the relaxation time $`t_{\mathrm{rh}}(0)`$ can be calculated in years and related to the stellar evolution timescale. The only other parameter to be specified is the metallicity $`Z`$, which plays an important role in the calculation of stellar evolution mass loss. In our treatment of stellar evolution, all wind mass loss rates are proportional to $`Z^{1/2}`$. Since stellar evolution introduces two new parameters in our initial models, a full exploration of the initial parameter space is clearly impossible. Fortunately, we will see that the effects of mass loss on the evolution to core collapse are rather unimportant, and so a systematic study is unnecessary at this point.

We first illustrate how post-main-sequence evolution of massive stars can prevent core collapse. In Figure 13 we show the evolution of a system similar to Model 2 with $`t_{\mathrm{rh}}(0)=60\mathrm{Myr}`$ and $`Z=10^4`$. At such low metallicity, there is little mass loss on the main sequence, so one expects $`t_{\mathrm{cc}}0.07\times 60\mathrm{Myr}=4.2\mathrm{Myr}`$ without evolution beyond the main sequence. However, massive stars evolve off the main sequence after about $`3`$Myr, before core collapse has occurred. Once they evolve off the main sequence, they typically lose up to $`2/3`$ of their mass rapidly. At the end of their life, which extends about 10 % beyond their main-sequence lifetime, these stars undergo supernova explosions and their cores collapse to BH. In the upper panel of Figure 13 we plot the number of these BH in the cluster. The significant mass loss during the late stages of stellar evolution causes the cluster core to expand so that collapse is prevented and the system then goes into long-term dynamical evolution.

Varying $`Z`$ from $`10^4`$ to 0.02 (the solar value), we found that this reversal of core contraction occurs for any metallicity when massive stars are allowed to evolve off the main sequence. Higher metallicities of course lead to increased mass loss and an even stronger tendency for the cluster core to expand. We have not considered the case of $`Z=0`$, i.e., a cluster of Pop III stars. These stars may collapse to BH that incorporate essentially the entire initial stellar mass (Heger et al., 2003). The dynamical evolution of the cluster would then be unaffected by mass loss. Whenever a runaway is avoided, a dense cluster of relatively massive primordial BH would then form near the center. This is an intriguing possibility to keep in mind for future consideration, although the IMF of Pop III stars is essentially unknown (but see, e.g., Nakamura & Umemura 2001; Abia et al. 2001) and recent hydrodynamic simulations suggest that the first stars may actually form isolated rather than in clusters (Abel et al., 2002).

We have also examined the possibility that mass loss on the main sequence may already be a source of indirect heating strong enough to reverse core collapse or, at least, delay it significantly. Our results indicate that wind mass loss from main-sequence stars alone can never prevent core collapse and that the delay introduced remains very small, even for high metallicities. It is of course always possible to fine-tune the relaxation time so that the mass loss slows down the evolution just enough to allow a few massive stars to evolve off the main sequence and stop contraction. But such models can only represent a very small domain in the parameter space of initial conditions. For example, in a system like Model 2 with $`t_{\mathrm{rh}}(0)=40\mathrm{Myr}`$, the core collapse time increases only by about 5% for $`Z=0.02`$. For $`t_{\mathrm{rh}}(0)=50\mathrm{Myr}`$, the evolution is just delayed long enough to prevent core collapse for $`Z=0.02`$, but not for $`Z=0.001`$. One has to increase the relaxation time to $`t_{\mathrm{rh}}(0)=60\mathrm{Myr}`$ for stellar evolution to prevent core collapse for all metallicities.

In conclusion, we find that post-main-sequence mass loss can be strong enough to prevent early core collapse, but that it does not significantly tighten the condition (on $`t_{\mathrm{rh}}(0)`$) for core collapse to occur. Mass loss from main-sequence stars always plays a negligible role.

## 5. Summary and Discussion

It has long been realized that a mass spectrum accelerates the dynamical evolution of dense star clusters through the process of mass segregation (Aarseth, 1966; Hénon, 1971a; Wielen, 1975). This is a consequence of the statistical tendency of two-body gravitational encounters to establish energy equipartition between stars of different masses. For any realistic IMF, equipartition can never be achieved because the massive stars quickly form a separate self-gravitating subcluster (a system with negative effective heat capacity) as they segregate near the center by transferring kinetic energy to lighter stars (Spitzer, 1969; Vishniac, 1978; Inagaki & Wiyanto, 1984; Inagaki & Saslaw, 1985). This instability leads to a rapid core collapse: after a finite time $`t_{\mathrm{cc}}`$, the central density of stars would actually become infinite in the absence of finite size effects (both physical collisions between stars and the finite number of stars in the core).

With the exception of a few recent $`N`$-body simulations (Portegies Zwart et al., 1999; Portegies Zwart & McMillan, 2002), all previous investigations of core collapse in clusters with a mass spectrum have considered a relatively narrow range of stellar masses ($`m_{\mathrm{max}}/m_{\mathrm{min}}100`$). This restriction is appropriate for old globular clusters in which stellar evolution probably had time to remove all massive stars before significant relaxation took place. In contrast, in the present study, we ignored stellar evolution to concentrate on systems in which, by assumption, core collapse occurs before even the most massive stars (with $`m100M_{}`$) leave the MS.

Through a set of high-resolution Monte Carlo simulations of clusters with a variety of IMF and structural parameters (concentration, presence or absence of tidal truncation or initial mass segregation), we established two important results. The first concerns the time needed to go deep into core collapse ($`t_{\mathrm{cc}}`$). As the evolution to core collapse is driven by relaxation processes, $`t_{\mathrm{cc}}`$ can always be written as proportional to the initial half-mass or central relaxation time ($`t_{\mathrm{rh}}(0)`$, $`t_{\mathrm{rc}}(0)`$). The proportionality constant does not depend on the size of the cluster or the number of stars, as long as they are numerous enough to avoid small-number effects. For isolated clusters with an IMF of realistic slope ($`\alpha 2.5`$ in the high-mass range), we find

$$t_{\mathrm{cc}}/t_{\mathrm{rh}}(0)\mu ^\delta \text{with}\mu =\frac{m_{\mathrm{max}}}{m}\text{and}\delta 1.3,$$
(20)
as $`\mu `$ increases from 1 (single-mass) to $`50`$. This result extends and agrees nicely with previous work based on FP simulations (Inagaki, 1985; Takahashi, 1997). From simple arguments about dynamical friction, one would naively expect a linear relation with $`\delta =1`$ (Binney & Tremaine, 1987). The steeper dependence found in our numerical results may be related to the fact that the instability causing the core collapse is triggered only after some critical number, $`N_{\mathrm{cr}}`$, of sufficiently massive stars have drifted to the core by dynamical friction. If one has to go out to radius $`R_{\mathrm{cr}}`$ to find this number of heavy stars, the required timescale will be of order $`\mu ^1t_{\mathrm{rel}}(R_{\mathrm{cr}})`$. For systems with increasing $`\mu `$, it is reasonable to expect $`N_{\mathrm{cr}}`$ and, consequently, $`R_{\mathrm{cr}}`$ and $`t_{\mathrm{rel}}(R_{\mathrm{cr}})`$ to decrease. While this provides a plausible explanation, further investigations will be needed to confirm it, in particular studies making use of the gaseous model of cluster dynamics (Freitag et al., 2003a).

For $`\mu >50`$, the core collapse time saturates to

$$t_{\mathrm{cc}}0.15t_{\mathrm{rc}}(0),$$
(21)
a key result given that any realistic IMF is likely to have $`\mu >100`$. This value of $`t_{\mathrm{cc}}`$ is at least two orders of magnitude shorter than the core collapse time for single-component clusters. Moreover, we find this $`t_{\mathrm{cc}}/t_{\mathrm{rc}}(0)`$ ratio to hold for all the models we have considered that have a finite $`t_{\mathrm{rc}}(0)`$. This includes King models (with or without tidal truncation) with $`W_0`$ ranging from 1 to 9, a sequence along which $`t_{\mathrm{cc}}/t_{\mathrm{rh}}(0)`$ decreases by more than $`100`$. Our result that the core collapse time is set fundamentally by $`t_{\mathrm{rc}}(0)`$ rather than $`t_{\mathrm{rh}}(0)`$ reflects an important difference between single-component clusters and systems with a broad mass spectrum: in systems with a mass spectrum, core collapse is driven ultimately by energy transfer occurring locally in the core; in contrast, relaxation in a single-component system is a global process taking place on all scales, resulting in a timescale comparable to $`t_{\mathrm{rh}}`$ (Inagaki, 1985). For Plummer models or King models with moderate concentration ($`W_04`$), we find $`t_{\mathrm{cc}}0.070.15t_{\mathrm{rh}}(0)`$, in good agreement with the results of $`N`$-body simulations by Portegies Zwart & McMillan (2002)555Although the fundamental relation is between $`t_{\mathrm{cc}}`$ and $`t_{\mathrm{rc}}(0)`$, $`t_{\mathrm{rh}}`$ is, by far, easier to estimate observationally..

The main-sequence lifetime $`t_{}`$ of very massive stars approaches a constant value, of about $`34`$Myr, with only weak dependence on metallicity, rotation, or the mass $`m`$ of the star, as long as $`m50M_{}`$ (Schaller et al., 1992; Meynet & Maeder, 2000; Maeder & Meynet, 2001). It exhibits little variation with $`m`$ because such massive stars are nearly Eddington-limited, hence $`Lm`$, where $`L`$ is the luminosity, and $`t_{}f_\mathrm{c}m/L`$constant if the fractional mass of the convective core, $`f_\mathrm{c}`$, does not depend too strongly on the stellar mass666If the mass segregation time scale $`t_\mathrm{s}1/m`$ is compared to the “usual” relation $`t_{}1/m^3`$, one could conclude erroneously that massive stars never play an important role in core collapse (Applegate, 1986). However, the approximate $`1/m^3`$ scaling of the stellar lifetime applies only for $`m10M_{}`$.. Therefore, the necessary condition for core collapse not being stopped by stellar evolution can be written

$$t_{\mathrm{rh}}(0)30\mathrm{Myr}.$$
(22)
Interestingly, we also find that core collapse cannot be prevented by mass loss through stellar winds on the MS. Such mass loss can only delay the evolution to core collapse very slightly. Even for high metallicities and hence strong winds, only clusters for which the core collapse time would otherwise be very close to the critical value (i.e., just below $`t_{}`$) can see their fate changed through wind mass loss if $`t_{\mathrm{cc}}`$ increases to become slightly longer than $`t_{}`$. These clusters must represent a very small domain in the parameter space of initial conditions.

Our second important finding is that, as the collapse proceeds, the mass contained in the ever-shrinking core converges to a non-zero value,

$$M_{\mathrm{cc}}0.0010.002M_{\mathrm{tot}},$$
(23)
which depends only weakly on the properties of the initial cluster. This contrasts again with single-component clusters which, exhibiting self-similar collapse, have $`M_{\mathrm{cc}}0`$ as $`tt_{\mathrm{cc}}`$.

This value of $`M_{\mathrm{cc}}/M_{\mathrm{tot}}`$ is in remarkable agreement with the fractional mass of a dark object, presumably a massive BH, discovered at the center of galaxies (in which case, $`M_{\mathrm{tot}}`$ stands for the mass in the spheroidal component) and, possibly, globular clusters like M15 and G1. This suggests that the collapsing core may well be the progenitor of a central IMBH. However, for a galactic nucleus, the requirement that $`t_{\mathrm{cc}}0.15t_{\mathrm{rc}}(0)<34`$Myr would imply that, either it evolved from a cluster with short (central) relaxation time to its present-day state with $`t_{\mathrm{rh}}10^8`$yr (Lauer et al., 1998), or that the growth of the central BH was predominantly from accretion of smaller star clusters with short relaxation times (each one carrying an IMBH at its center; Ebisuzaki et al. 2001; Hansen & Milosavljević 2003). M15 and G1 also have long half-mass relaxation times ($`1`$Gyr for M15, Dull et al. 1997, and $`3050`$Gyr for G1, Meylan et al. 2001; Baumgardt et al. 2003b) indicating that the process we envision could only have taken place if these clusters were born with much more centrally concentrated density profiles than observed today777It is not clear that the present-day value of $`t_{\mathrm{rh}}`$ for a cluster, which can be estimated from observations, must be close to $`t_{\mathrm{rh}}(0)`$. (for M15, an initial value of $`W_08`$ seems required; see Fig. 11). However, there is no doubt that young star clusters with half-mass relaxation times shorter than the critical value of 30 Myr do exist. The Arches cluster near the center of our Galaxy is the best established example, with $`t_{\mathrm{rh}}`$ probably shorter than 10 Myr (Figer et al., 1999a), while the Quintuplet cluster appears to lie very close to the critical value (Figer et al., 1999b). It is also possible that some super star clusters, the birthplaces of most stars in starburst environments, have sufficiently short relaxation times (e.g. Ho & Filippenko, 1996).

Although the final fate of the stars that participate in core collapse is not completely certain, very high rates of physical collisions are expected in the core. Primordial binaries are probably unable to stop the collapse and prevent these collisions. In smaller systems like globular clusters interactions with hard binaries may in fact increase the collision rate. In more extreme environments like (proto-)galactic nuclei, the velocity dispersion is so high ($`10^210^3\mathrm{km}\mathrm{s}^1`$) that most primordial binaries are soft and will have been disrupted before they can play a role in the collapse. The main uncertainty is then whether these collisions, occurring at relatively high velocity, allow a massive star to gain mass in a runaway fashion, as suggested by analytical models (based on the dependence of the cross section on the mass and neglecting collisional mass-loss; Malyshkin & Goodman 2001, and references therein), or, on the contrary, grind it down progressively. Our preliminary results from MC simulations including a realistic treatment of stellar collisions (Rasio et al., 2003) indicate that the formation of a very massive star through runaway collisions is a likely outcome. More extensive calculations, as well as a detailed discussion of the likely fate of the merger remnant, will be presented in Paper II.

We thank Holger Baumgardt for discussions and for providing data from his $`N`$-body simulations. We are grateful to Rainer Spurzem for making the code for gaseous models (SPEDI) available and for his assistance in using it. Piero Spinnato and Simon Portegies Zwart kindly gave us the initial condition data for the “dynamical friction” run presented in the Appendix. We thank K. Belczynski for providing us the population synthesis code StarTrack and for his help in integrating it into our MC code, and to John Fregeau for useful discussions. This work was supported by NASA ATP Grant NAG5-12044 and NSF Grant AST-0206276 to Northwestern University. The work of MF is funded by the Sonderforschungsbereich (SFB) 439 ‘Galaxies in the Young Universe’ (subproject A5) of the German Science Foundation (DFG) at the University of Heidelberg. MF thanks the Theoretical Astrophysics Group at Northwestern University for financial support and hospitality during several visits which made this collaboration possible.

## Appendix A Applicability of the Monte Carlo Method

The Monte Carlo method has been applied previously to two-component systems (see, e.g., Watters et al., 2000) and to systems with a continuous mass function covering a relatively narrow range (with $`m_{\mathrm{max}}/m_{\mathrm{min}}30`$; see the references in Section 2). However, one may question the applicability of the method to systems with a much broader mass spectrum. In particular, the nature of the energy transfer in these systems is different than in systems with small $`m_{\mathrm{max}}/m_{\mathrm{min}}`$ ratio (See the discussion in Sec. 1.2). In addition, the gravitational potential is assumed to change smoothly within the cluster in the MC method. However, in a system containing a broad IMF, a relatively small number of very massive stars could introduce significant discrete changes in the slope of the potential. As the positions of these steps change randomly from one iteration to the next, one may worry that this could introduce a significant amount of spurious relaxation in the system. In this Appendix, we provide numerical tests and theoretical arguments showing that the MC method is indeed applicable to systems with a very broad mass spectrum, such as those under consideration in this paper.

### A.1. Spurious Relaxation

In MC simulations, the orbits of the stars are calculated using a potential assumed to be spherical and smooth. This potential is constructed numerically from a radial distribution of point masses and therefore includes random fluctuations that can lead to spurious relaxation effects (Hénon, 1971a). To test the importance of these effects, we carried out simulations in which we turned off physical relaxation (by skipping the perturbation of stars). Ideally, in this case, the cluster should maintain a constant state indefinitely and none of its properties should change over time. In practice, however, it is sufficient to make sure that the changes are not significant compared to changes occurring in the presence of physical relaxation.

In one of our tests, we let a MC simulation, for our typical Model 2, run for 3600 iteration steps, which is twice the number of steps normally required for this model to reach core collapse. We plot various relevant quantities showing the evolution of the system in Figure 14. It is clear that the deviations in these quantities are much smaller than the changes due to physical relaxation, so we conclude that the effects of spurious relaxation is not significant.

### A.2. Relaxation and Dynamical Friction

The energy transfer from massive to lighter stars and the segregation of massive stars to the center of the cluster are processes similar to simple dynamical friction for a heavy test particle embedded in a background of lighter stars (Binney & Tremaine, 1987, Sec. 7.1). Some of the processes that govern the dynamics of dense stellar systems cannot be treated simply as part of two-body relaxation. Well-known examples include collisions, binary interactions, and three-body formation of binaries. A natural question to ask is whether dynamical friction is a process that can be modeled correctly by our implementation of two-body relaxation in the MC method. The answer is central to the applicability of the MC method to our problem.

The frictional drag on a massive object is proportional to the mass density of light background stars, and is independent of the masses of individual stars (see Binney & Tremaine, 1987, eq. 7-18). So the effect of dynamical friction would be unaltered, for example, if every light star were replaced by two stars of half the mass. However, the treatment of relaxation in the MC method (Sec. 2.1) does depend on the mass ratio of the two particles considered and would be altered by such a replacement. The ratio of relaxation to crossing time is given by

$$\frac{t_\mathrm{r}}{t_{\mathrm{dyn}}}\frac{N}{\mathrm{ln}N}$$
(A1)
which would evidently change if the stars in a system were replaced by many more stars with smaller masses. This might suggest that the MC method cannot follow dynamical friction correctly. However, note that the relaxation time can be written (Spitzer, 1987, eq. 3-37)

$$t_\mathrm{r}\frac{v^3}{nm^2}$$
(A2)
and the timescale for dynamical friction is proportional to (Binney & Tremaine, 1987, eq. 7-18)

$$t_{\mathrm{DF}}\frac{v_M^3}{nmM}.$$
(A3)
Here, $`v`$ is the average velocity of background stars, of mass $`m`$ and number density $`n`$, and $`v_M`$ is the velocity of the massive test particle, of mass $`M`$. When $`v_M`$ is comparable to $`v`$ (which is true, for example, if the massive particle is on a circular orbit; see below), we conclude that

$$t_{\mathrm{DF}}\frac{m}{M}t_\mathrm{r}\frac{v^3}{\rho M}.$$
(A4)
The first proportionality implies that dynamical friction can be treated as a relaxation process. In the second, we show explicitly the dependence on the total mass density $`\rho `$ of the background stars. This result is also established by several numerical calculations (Spitzer & Shull, 1975; Bonnell & Davies, 1998; Fregeau et al., 2002).

To demonstrate that the MC algorithm is indeed able to handle simple dynamical friction, we have simulated a system consisting of a single massive object (mass $`m_2`$) in a cluster of much lighter particles (mass $`m_1`$, number $`N_1`$). The initial conditions are identical to those used for the $`N_1=\mathrm{400\hspace{0.17em}000}`$ model of Spinnato et al. (2002), with $`m_2/m_1=211`$, hence the total number of particles used for this test is 400 001. Spinnato et al. were interested in the spiral-in of an IMBH in a galactic nucleus. The massive particle is initially on a circular orbit at a relatively large distance from the center, $`r(0)`$. For the power-law density profile used to represent the cluster, $`\rho (r)r^\gamma `$, one can predict analytically the decay $`R(t)`$ using Chandrasekhar’s formula for dynamical friction (Spinnato et al., 2002, eq. 5). In Figure 15, we compare the spiral-in of the massive particle with the analytical solution. In the formula, we have set $`\gamma =2.3`$ to fit (visually) the density profile of the cluster. The initial $`N`$-body set-up was provided by Spinnato et al. and converted directly into MC particles. The agreement between the MC run and the analytical formula is very satisfactory, especially considering that the density profile does not follow the exact power law and in fact evolves slightly during the simulation as a result of relaxation between light particles. Note that this good agreement does not validate the assumptions used in deriving Chandrasekhar’s formula (e.g., neglecting large-angle scatterings, correlations between light particles, and gradients in the properties of the cluster) because the MC method relies basically on the same set of assumptions. In particular, $`\gamma _\mathrm{c}`$ is a free parameter in both approaches. What this result demonstrates is that these assumptions are correctly implemented in our MC code.