MONTE CARLO SIMULATIONS OF GLOBULAR CLUSTER EVOLUTION. VI. THE INFLUENCE OF AN INTERMEDIATE-MASS BLACK HOLE

As recently as 10 years ago, it was generally believed that black holes (BHs) occur in two broad mass ranges: stellar (M_BH⋍3–20 M_☉), produced by the core collapse of massive stars, and supermassive (M_BH ∼ 10⁶–10¹⁰ M_☉), believed to have formed in the centers of galaxies at high redshift and grown in mass as a result of galaxy mergers (see, e.g., Volonteri et al. 2003). However, the existence of BHs with masses intermediate between these two ranges could not be established by observations until recently, although intermediate-mass BHs (IMBHs) were discussed by theorists more than 30 years ago (see, e.g., Wyller 1970). Indirect evidence for IMBHs has accumulated over time from observations of so-called ultraluminous X-ray sources (ULXs), objects with fluxes that exceed the angle-averaged flux of a stellar mass BH accreting at the Eddington limit. An interesting result from observations of ULXs is that many, if not most, of them are associated with star clusters. It has long been speculated (e.g., Frank & Rees 1976) that the centers of globular clusters (GCs) may harbor BHs with masses ∼10³ M_☉. If so, these BHs affect the distribution function of the stars, producing velocity and density cusps (Bahcall & Wolf 1976). While the detection of ULXs can only give indirect evidence of the presence of IMBHs, observations of cuspy velocity profiles would make it possible to directly determine the BH mass. However, the radius of influence of an IMBH, defined as the radius where the orbital velocity around the BH equals the velocity dispersion of the cluster, is very small. For example, at a distance of 10 kpc, a 10³ M_☉ BH would influence orbits within ≈1'', making observations very challenging.

Studies of the surface density profile of GCs offer a complementary method of constraining the effects of an IMBH on the host GC stars. A recent study by Noyola & Gebhardt (2006) obtained central surface brightness profiles (SBPs) for 38 Galactic GCs from Hubble Space Telescope (HST) WFPC2 images. They showed that half of the GCs in their sample have slopes for the inner SBPs that are inconsistent with simple isothermal cores, which may be indicative of an IMBH. However, it is challenging to explain the full range of slopes with current models. While analytical models can only explain the steepest slopes in their sample, recent N-body models of GCs containing IMBHs (Baumgardt et al. 2005), might explain some of the intermediate surface brightness slopes.

However, the disadvantage of current N-body simulations is that for realistic cluster models, which take into account stellar evolution and a realistic mass spectrum, the number of stars is restricted to typically less than ∼10⁵ as these simulations require a large amount of computing time. However, many GCs are known to be very massive, with masses reaching up to 2 × 10⁶ M_☉, resulting in a much larger number of stars one has to deal with when modeling these objects. In previous N-body simulations, such large-N clusters have been scaled down to low-N systems. Scaling down can be achieved in two ways (e.g., Baumgardt et al. 2005): either the mass of the central IMBH M_BH is kept constant and N is decreased, effectively decreasing the total cluster mass M_C, or the ratio M_BH/M_C is kept constant, while lowering both M_BH and M_C. As both M_BH/M_C and the ratio of M_BH to stellar mass are important parameters that influence the structure and dynamics of a cluster, but cannot be held constant simultaneously when lowering N, it is clear that only with the real N can a fully self-consistent simulation be achieved. Scaling becomes even more difficult once other physical processes are included in the simulations such as stellar evolution or stellar collisions. Using the correct number of stars in a dynamical simulation ensures that the relative rates of different dynamical processes (which all scale differently with N) are correct.

It is clear that in order to study the evolution of old GCs that might harbor IMBHs, more approximate methods have to be employed. They fall roughly into two categories, methods that treat the cluster as a continuum distribution function (Amaro-Seoane et al. 2004; Takahashi 1997; Giersz & Spurzem 1994; Murphy et al. 1991; Cohn & Kulsrud 1978; Lightman & Shapiro 1977; Bahcall & Wolf 1977) and Monte Carlo (MC) methods that use a particle-based approach (see, e.g., Fregeau & Rasio 2007; Freitag & Benz 2002, 2001; Shapiro & Marchant 1978; Hénon 1971 and references therein). While the former methods have been successfully used to study the effect of a massive BH on the cluster structure (Amaro-Seoane et al. 2004; Murphy et al. 1991; Cohn & Kulsrud 1978; Lightman & Shapiro 1977), the model clusters were highly idealized, consisting only of equal-mass stars, and did not incorporate stellar evolution. Although including additional physical processes is not impossible, it remains nevertheless highly non-trivial for these methods. The MC method, on the other hand, relies on a star-by-star description of the cluster and has, therefore, the great advantage that additional processes are easily incorporated.

Our group has been developing over many years a state-of-the-art MC code, which treats many relevant processes in sufficient detail, making direct comparison with GC observations feasible (Fregeau & Rasio 2007, and references therein). This paper is the sixth in a series studying the fundamental aspects of cluster dynamics using this code. Here, we will describe the changes we made to our code in order to incorporate the effect of a massive central IMBH and carry out comparison runs with idealized models as well as more realistic cluster simulations published previously in the literature. In Section 3, we briefly describe our method and the changes we made to the MC code. We validate our code by comparing our results to previously published results in the literature, using idealized cluster models (Section 4), as well as more realistic ones that include stellar evolution (Section 5). In Section 5.2, we present SBPs from our large-N runs and compare them to observations of Noyola & Gebhardt (2006). We conclude in Section 6.

The dynamical effect of an IMBH on the surrounding stellar system was first described by Peebles (1972), who argued that the bound stars in the cusp around the BH must obey a shallow power-law density distribution to account for stellar consumption near the cluster center. Analyzing the Fokker–Planck equation in energy space for an isotropic stellar distribution, Bahcall & Wolf (1976) obtained a density profile with n(r)∝r^−7/4, which is now commonly referred to as the Bahcall–Wolf cusp. The extension of the cusp solution is given by the radius of influence of the BH, r_i, which is defined as

$$\begin{equation} r_{i}\equiv \frac{GM_{{\rm BH}}}{2\,\sigma ^{2}}\,, \end{equation} \tag{ 1 }$$

where G is the gravitational constant and σ is the velocity dispersion of the core. As shown later by Shapiro & Lightman (1976), such a solution can be readily obtained using simple scaling arguments. The key is to realize that in the region delimited by the tidal radius, r_t, within which stars get tidally disrupted, and r_i within which the stars are bound to the BH, the net energy diffusion timescale, t_U, is proportional to the local relaxation time, t_r, which is the shortest timescale on which any physical quantity can be transported. Furthermore, the quasi-steady state of the cusp region is characterized by a dynamic equilibrium, with a constant net energy flow into the core region that should scale as n(r) r³E(r)/t_U, where E(r) is the mean specific energy at radius r, so E(r) ∼ GM_BH r⁻¹ in the cusp region. Setting then t_U ∼ t_r ∼ σ^3/2/[G²m²n(r)], $\sigma\, {\sim}\, \ssty\sqrt{GM_{{\rm BH}}\, r^{-1}}$, and simple substitutions immediately lead to n(r)∝r^7/4.

The formation of such a cusp has been confirmed subsequently by numerical studies. Most of them employed methods based on the Fokker–Planck equation, solving it directly (Bahcall & Wolf 1976; Lightman & Shapiro 1977; Cohn & Kulsrud 1978), or indirectly, either using the statistical MC approach (Shapiro & Marchant 1978; Freitag & Benz 2002) or a fluid-dynamical approach based on velocity moments (Amaro-Seoane et al. 2004). Being derived from the Fokker–Planck equation, all these methods share essentially the same set of underlying assumptions: (1) the cluster potential has spherical symmetry, (2) the cluster is in dynamical equilibrium at all times, and (3) the evolution is driven by diffusive two-body relaxation. Through direct N-body simulations that do not rely on any a priori assumptions, Baumgardt et al. (2004a) confirmed the cusp solution of Bahcall & Wolf (1976), therefore also providing important justification to the Fokker–Planck approach.

Based on the cusp solution, Frank & Rees (1976) calculated the stellar disruption rate taking into account that stars inside a critical radius, r_crit, are efficiently accreted by the BH as they diffuse quickly into low angular momentum orbits with periastron distances, r_peri, smaller than r_t. As for a given radial position, the velocity vectors that lead to orbits with r_peri < r_t form a cone; these orbits are also called loss-cone orbits. Outside r_crit, stars are always able to leave the loss cone during one orbital period due to two-body relaxation while inside r_crit the orbital changes are smaller so that stars on loss-cone orbits are more likely to reach the tidal radius and get disrupted before they have a chance to get scattered out. Consequently, inside r_crit the loss cone should be nearly empty while outside r_crit it always remains full. In addition, Frank & Rees (1976) argue that the disruption rate is mainly given by the cluster conditions at r_crit. This is because, on the one hand, for r > r_crit the fraction of stars populating the loss cone decreases as the loss-cone angle decreases with increasing radius, while for r < r_crit the total net flux of stars, and therefore the flux of stars that can diffuse into the loss cone, decreases rapidly (Lightman & Shapiro 1977). Calculations by Amaro-Seoane et al. (2004) confirm the loss-cone picture, showing that, for a cluster in dynamical equilibrium, the disruption rate is strongly peaked at r_crit, and the fraction of stars on loss-cone orbits is rapidly approaching zero for r < r_crit, while for r > r_crit it is always close to one. Similarly, Baumgardt et al. (2004a) find generally good agreement between the disruption rates in their simulations and the disruptions rates based on the approximate expression of Frank & Rees (1976; their Equation (22)) and using the cluster conditions at r_crit from their N-body models.

While earlier work on the dynamics of clusters with IMBHs mainly focused on the equilibrium state of the cusp surrounded by a static isothermal core, Shapiro (1977) considered the effect of an IMBH on the global evolution of the cluster. Using a homological model for the dynamical evolution, he calculated the core size in response to evaporation of high-velocity stars and tidal disruption of stars tightly bound to the IMBH. While stars that leave the cluster by evaporation carry away very little energy-driving core contraction (Spitzer & Saslaw 1966), which would ultimately lead to core collapse in the absence of an IMBH, the tidal disruptions close to the central IMBH that remove stars with highly negative specific energies provide an energy source that causes the core to expand. Shapiro (1977) shows that for low initial IMBH masses and large initial core radii, stellar evaporation first dominates and drives core contraction until, due to the increasing core density, the tidal disruption rate becomes large enough to reverse core collapse. The time of this reversal roughly coincides with the time of core collapse for the cluster without IMBH. Tidal disruptions then drive the re-expansion of the core, and the core size increases asymptotically to infinity. This expansion is a generic feature of a stellar system where energy is generated within a very small central volume and the mass contained within this region is very small compared to the cluster mass (Hénon 1965; Shapiro 1977).

The qualitative behavior of the core size evolution for lower mass IMBHs was later confirmed by numerical studies (Marchant & Shapiro 1980; Murphy et al. 1991; Freitag & Benz 2002; Amaro-Seoane et al. 2004). Amaro-Seoane et al. (2004) in particular showed that the core size increases asymptotically as ∝t^2/3, which was also predicted by Shapiro (1977). This expansion also causes the disruption rate to decrease with time as approximately ∝t^−6/5 (Amaro-Seoane et al. 2004) and will ultimately lead to the complete dissolution of the cluster as the outer stars are removed by tidal forces (Wielen 1971). For larger IMBH masses and small initial core sizes, tidal disruptions will prevent any initial core contraction and the core expands from the beginning (Shapiro 1977). This case was calculated by Baumgardt et al. (2004a) using direct N-body simulations, confirming that core expansion starts almost immediately and follows a t^2/3 power law.

Most of the studies mentioned so far considered the evolution of clusters containing a central massive BH and comprised of stars of equal mass. A stellar mass spectrum was first considered by Bahcall & Wolf (1977) extending their previous work in Bahcall & Wolf (1976). They find that, due to mass segregation, lower mass stars have shallower density profiles than more massive ones. For old GCs this means that the observable surface brightness cusp must be much shallower as the more massive dark stellar remnants are concentrated toward the center while the lower-mass main-sequence stars that contribute most of the light are much less centrally concentrated. It follows that, although a cusp in the velocity and density profile provides strong evidence for the presence of an IMBH in a cluster, such cusps might not be easily detectable in real star clusters. Using direct N-body simulations including an initial stellar mass spectrum and stellar evolution, Baumgardt et al. (2004b) find, indeed, flat luminosity density profiles almost indistinguishable from a standard King profile. Carrying out similar simulations but with M_BH/M_c < 1%, Baumgardt et al. (2005) find surface brightness cusps with power-law slopes ranging from α = −0.1 to α = −0.3. Based on these results they identified nine candidate clusters from the sample of galactic GCs of Noyola & Gebhardt (2006) that might contain IMBHs.

MC simulations of realistic clusters with a central IMBH were mainly done in the context of galactic nuclei. A recently developed and well-tested code is that of Freitag & Benz (2002). Similar to our code, it is based on the method of Hénon (1971) but is modified to evolve each star individually on a fixed fraction of the local relaxation time, as opposed to the original shared time-step scheme. In addition to the implementation of loss-cone physics, which we will describe in detail in Section 3, it also incorporates stellar collisions interpolating between results from detailed hydrodynamical simulations. Collisions between stars is an important physical process in dense galactic nuclei. As already pointed out by Frank & Rees (1976), the radius r_coll outside which stellar encounters responsible for relaxation can be treated as elastic encounter can be larger than r_crit for large M_BH and typical sizes and densities for galactic nuclei. Inside r_coll stars cannot deflect each other significantly without colliding. As the relative velocity within r_coll is larger than the escape velocity from the stellar surface the collision can be very disruptive and might under certain conditions provide a significant source to fuel an active galactic nucleus (see Freitag & Benz 2002 and references therein). However, for typical conditions inside GCs stellar collisions are unlikely to play a significant role and are therefore not further considered for the present study.

In addition to the formation of cuspy profiles, an IMBH influences the SBP by producing rather large cluster cores as measured by the core-to-half-light radius, such that larger cores are produced by more massive IMBHs (Heggie et al. 2007). The large core sizes are simply a result of the energy flow from the central cusp region to the core which causes the cluster to expand. Constructing generalized King models (King 1966) including the effect of an IMBH and a stellar mass spectrum, Miocchi (2007) finds that the core size and the cusp slope are related such that clusters with larger slopes, s, have lower concentrations, c = log (r_c/r_t), where r_c is the core radius of the cluster. More specifically, they find that s and c are related by

$$\begin{equation} 11.6s-4.85\lesssim \log \left(\frac{M_{{\rm BH}}}{M_{c}}\right)\lesssim -1.14c-0.694\,, \end{equation} \tag{ 2 }$$

where M_c is the cluster mass. Based on this criterion and data from Noyola & Gebhardt (2006) for s as well as the Harris catalog for c, they identified seven candidate clusters that might contain IMBHs with mass >100 M_☉, with four of them also identified by Baumgardt et al. (2005).

In contrast to the rather shallow surface brightness cusps, the stronger cusp in the stellar velocity dispersion appears to be a much better diagnostic to infer the presence of IMBHs in GCs. However, for GCs this signature turns out to be difficult to detect as any velocity dispersion measurement inside such a cusp has to rely on only a few bright stars for expected IMBH masses ≲ 1000 M_☉ and typical GC masses (Baumgardt et al. 2005). These IMBH mass estimates are based on extrapolating the well-known M_BH–σ relation (Magorrian et al. 1998; Gebhardt et al. 2000; Ferrarese & Merritt 2000) between central BH mass and velocity dispersion in the central bulges of galaxies down to velocity dispersions typical for GCs (∼10 km s⁻¹). In that case, M_BH should be ∼10³–10⁴ M_☉. Furthermore, simulations of collisional runaways by Gürkan et al. (2004) show that the mass of the Spitzer unstable subcluster, which provides the mass reservoir for forming the BH progenitor, is ∼10⁻³ M_c, implying IMBH masses of at most a few 10³ M_☉ for typical cluster masses. However, it is important to point out that this does not mean that M_BH/M_c for an old GC has to be always significantly less than 1%, as a cluster can later lose a substantial amount of mass due to tidal stripping in the Galactic field.

3.1. Monte Carlo Method with IMBH

Our MC code shares some important properties with direct N-body methods, which is why it is also regarded as a randomized N-body scheme (see, e.g., Freitag & Benz 2001). Just as in direct N-body codes, it relies on a star-by-star description of the GC, which makes it particularly straightforward to include additional physical processes such as stellar evolution. Contrary to direct N-body methods, however, the stellar orbits are resolved on a relaxation timescale t_r, which is much larger than the crossing time t_cr, the timescale on which direct N-body methods resolve those orbits. The specific implementation we use for our study is the MC code initially developed by Joshi et al. (2000) and further enhanced and improved by Fregeau et al. (2003) and Fregeau & Rasio (2007). The code is based on Hénon's algorithm for solving the Fokker–Planck equation. It incorporates treatments of mass spectra, stellar evolution, primordial binaries, and the influence of a galactic tidal field.

The effect of an IMBH on the stellar distribution is implemented in a manner similar to that of Freitag & Benz (2002). In this method the IMBH is treated as a fixed, central point mass while stars are tidally disrupted and accreted onto the IMBH whenever their periastron distances lie within the tidal radius, r_t, of the IMBH, which is given by

$$\begin{equation} r_{t}=\left(2\frac{M_{{\rm BH}}}{M_{*}}\right)^{1/3}R_{*}\,, \end{equation} \tag{ 3 }$$

where R_* and M_* are the stellar radius and mass, respectively. Stars are removed from the system and their masses are added to the BH as soon as their velocity vector, v, enter the loss cone, θ_LC, approximately given by

$$\theta _{{\rm LC}}^{2}\simeq 2\frac{GM_{{\rm BH}}r_{t}}{v^{2}r^{2}}\,.$$

However, as the star's removal happens on an orbital timescale, one would need to use time steps as short as the orbital period of the star in order to treat the loss-cone effects in the most accurate fashion. This would, however, slow down the whole calculation considerably. Instead, during one MC time step a star's orbital evolution is followed by simulating the random walk of its velocity vector, which approximates the effect of relaxation on the much shorter orbital timescale. The random-walk procedure is as follows.

• 1.  
  After a gravitational encounter between two stars is calculated in the standard MC fashion, resulting in a deflection angle of δθ_step, the orbital period, P_orb, is calculated using Gauss–Chebychev quadrature, and a "representative" diffusion angle during a single orbit, δθ_orb, is estimated as

  $$\delta \theta _{{\rm orb}}=\left(\frac{\delta t}{P_{{\rm orb}}}\right)^{-\frac{1}{2}}\delta \theta _{{\rm step}},$$

  where δt is the time step.
• 2.  
  The star's velocity vector with respect to the encounter reference frame is calculated and a variable L₂, which represents the remaining quadratic deflection, is set to δθ²_step.
• 3.  
  The star is tested for entry into the loss cone, and, if this is the case, is removed and its mass added to M_BH, whereupon the random walk is terminated. Otherwise, we proceed with the next step.
• 4.  
  If L₂ ⩽ 0, the random walk is terminated. The star's position and velocity are reset to its values before the random-walk procedure.
• 5.  
  A new random-walk step is carried out with amplitude $\Delta =\max (\delta \theta _{{\rm orb}},\min (0.1\pi,\Delta _{{\rm safe}},\sqrt{L_{2}}))$, and a random direction on the velocity sphere, where Δ_safe is set to roughly half the angular distance to the loss cone. This way, Δ becomes progressively smaller down to δθ_orb when approaching the loss cone in order to keep the risk of missing a disruption minimal. According to the new step size, the direction of the star's velocity is changed and L₂ is updated: L₂ := L₂ − Δ. The random walk continues at step 3.

Although many of our results turn out to agree very well with previously published data, there are discrepancies when comparing the disruption rates. One possible reason for these differences might be related to the fact that our code uses a shared time-step scheme. In an individual time-step scheme, as in Freitag & Benz (2001), the time step is some constant fraction of the local relaxation time, i.e., dt_i = ft_r(r_i), where f is a constant and the subscript i refers to the individual star. In a shared time-step scheme the smallest of these dt_i is chosen for all stars. This results in much shorter time steps for stars farther out in the cluster compared to an individual time-step scheme. As has been noted by Freitag et al. (2006), the time-step size must be chosen to be small enough in order to achieve good agreement with N-body simulations, with f ≲ 0.01. While choosing such a small time step was still feasible in the code of Freitag & Benz (2002), to enforce such a criterion for all stars in our code would lead to a dramatic slowdown and notable spurious relaxation as the time steps for the stars in the outer cluster regions relative to the local relaxation time become extremely small (see Gürkan et al. 2004). In order to reduce the effect of spurious relaxation we are forced to choose a larger dt resulting in a larger f for the inner regions, up to f ≈ 0.1. Figure 1 shows an example of dt/t_r(r) as a function of radius. One can see that in this case only for r ≳ 0.1 r_h we have dt/t_r(r) ≲ 0.01, while for the inner region it rises quickly to 0.1. As a result, any expansion of the inner cluster region, which limits the growth of the IMBH in our comparison calculations, will be slower and the disruption rates, therefore, larger, due to the higher core densities. Such an expansion occurs either in response to the growing IMBH mass or as an initial expansion of the cluster. Previously, we applied a procedure that tries to compensate for the larger time steps in the inner region (Umbreit et al. 2008) by evolving the inner stars individually on smaller time steps while keeping the cluster potential constant during one shared time step, a procedure that bears some resemblance to the method of Marchant & Shapiro (1980). However, it turned out that in order to achieve good agreement with direct N-body simulations one has to choose the time-step sizes for each cluster configuration separately. This is not only undesirable but might also imply that there are additional processes at work that significantly influence the disruption rate and are not included in the MC scheme. In Section 4.2.2 we discuss several possibilities, among them the wandering of the IMBH. As there is no obvious way to compensate for such a processes in a uniform and consistent way through adjustments in the two-body relaxation timescale, we do not use the sub-time-step scheme for the present paper.

Zoom In Zoom Out Reset image size

3.2. Initial Conditions

Table 1 summarizes the initial conditions and main results for all our single-mass runs and Table 2 for our multi-mass runs that include stellar evolution (implemented by Fregeau et al. 2009 using the SSE code of Hurley et al. 2002).

Table 1. Details of the Performed Monte Carlo Runs for Single-mass Clusters with Initial Conditions as in Baumgardt et al. (2004a; BME) and Amaro-Seoane et al. (2004) and Freitag & Benz (2002) (ASFB)

Name	N	W0	rt	MBH, i	MBH, f	TEND	NDisr
			(Rvir)	(M☉)	(M☉)	(tcr)
BME-1	80, 000	10	1 × 10−7	266	1429 ± 35 (827)	3000	1163 ± 35 (561)
BME-2	80, 000	10	1 × 10−7	800	1690 ± 20 (1388)	2000	890 ± 20 (588)
BME-3	80, 000	10	1 × 10−7	2660	3434 ± 20 (3285)	2000	774 ± 20 (625)
BME-16	178, 800	10	1 × 10−7	461	2171 ± 30 (1368)	2000	1710 ± 30 (907)
ASFB-50	100, 000	Plummer	2.26 × 10−8	50	10500 ± 100 (7450, 13050)	5.4 × 106	10, 450 ± 100 (7.4 × 103, 1.3 × 104)
ASFB-500	100, 000	Plummer	2.26 × 10−8	500	11500 ± 100 (8900, 14500)	5.4 × 106	11, 000 ± 100 (8.4 × 103, 1.4 × 104)

Note. Values in parentheses are results from the corresponding literature, where two values are given: the first one refers to Freitag & Benz (2002) and the second to Amaro-Seoane et al. (2004).

Download table as:  ASCIITypeset image

Table 2. Details of the Performed Monte Carlo Runs for Multi-mass Clusters with Initial Conditions as in Baumgardt et al. (2005; BMH)

Name	N	W0	MBH, i	rh, i	MBH, f	NDisr	Mc, f	TEND	rh, f	log trh
			(M☉)	(pc)	(M☉)		(M☉)	(Gyr)	(pc)	(yr)
BMH-1	131, 072	7	125	4.91	133 ± 3 (137)	4 ± 2	45, 671 ± 6 (45, 534)	12	12.02 ± 0.05 (12.31)	9.86 (9.82)
BMH-2	131, 072	7	250	4.91	260 ± 4 (280)	6 ± 2	45, 677 ± 14 (45, 311)	12	11.98 ± 0.04 (12.60)	9.86 (9.84)
BMH-3	131, 072	7	500	4.91	513 ± 5 (531)	9 ± 3	45, 400 ± 30 (44, 741)	12	12.36 ± 0.04 (13.70)	9.88 (9.89)

Notes. The star masses were chosen according to a Kroupa mass function ranging from 0.1to30 M_☉. Stellar evolution was modeled using the SSE code of Hurley et al. (2002).

Download table as:  ASCIITypeset image

The single-mass clusters consist of N stars all of mass 1 M_☉ with positions and velocities chosen according to a Plummer model or a King model with dimensionless potential W₀ = 10. Radii and times are given in terms of the virial radius, R_vir, and the crossing time, t_cr, respectively, which are defined by

$$\begin{equation} R_{{\rm vir}}=\frac{GM_{c}^{2}}{-4\, E_{0}}\,, \end{equation} \tag{ 4 }$$

and

$$\begin{equation} t_{{\rm cr}}=\frac{GM_{c}^{5/2}}{\left(-4E_{0}\right)^{3/2}}\,, \end{equation} \tag{ 5 }$$

where E₀ is the total gravitational energy of the cluster. The cluster was evolved up to a time T_END, with an initial IMBH mass M_{BH, i} and r_t a constant for all stars. Also shown are the resulting final IMBH mass, M_{BH, f}, and the number of tidally disrupted stars, N_disr.

The initial conditions for multi-mass clusters are chosen as in BMH, with stellar masses drawn from a Kroupa mass function (Kroupa 2001) in the range of 0.1–30 M_☉, stars evolved at a metallicity Z = 0.001, and initial positions and velocities chosen according to a King model with W₀ = 7. The tidal disruption radius for each star is given by Equation (3) with stellar mass and radius provided by the stellar evolution code. In addition to M_{BH, i}, M_{BH, f}, and N_disr the half-mass radius, r_h, cluster mass, M_c, and the final relaxation time at the half-mass radius, t_rh, are shown, where the added subscripts i and f indicate initial and final values, respectively.

For each set of initial conditions nine MC runs were performed with varying random seed, and the values given in the table are the averages and standard deviations of these runs.

4.1. Imprints of IMBHs

In Figure 2 density and velocity profiles from our simulations of single-mass clusters are shown, together with the expected r_i calculated from Equation (1) and using σ ≈ 0.55 in N-body units, appropriate for a W₀ = 10 King model. As can be clearly seen, the density profile of the inner region of the evolved clusters follow closely the expected n(r)∝r^−7/4 power law and the extent of the cusp matches that of the region where the velocity dispersion is Keplerian, which in turn matches the expected r_i. However, contrary to what is seen in direct N-body simulations, the cusp extends down to much smaller radii, especially for BH masses below 1% of the cluster mass. This is mainly because, in our simulations, the central IMBH has a fixed position, while in direct N-body simulations it is allowed to move freely. As a consequence, the density profile flattens inside its wandering radius compared to a pure cusp profile, resulting in fewer stars in the central region. As will be shown later, this might have an influence on the rate at which stars are tidally disrupted.

Zoom In Zoom Out Reset image size

4.2. Disruption Rates

4.2.1. Comparison with Amaro-Seoane et al. (2004) and Freitag & Benz (2002)

Figure 3 compares the growth of the IMBH in our simulations and those of Amaro-Seoane et al. (2004), using a gas code, and Freitag & Benz (2002), using an individual time-step MC code. Here, the evolution of a cluster with 10⁵ stars all of 1 M_☉ and a fixed massive BH in the center with a mass of 50 M_☉ or 500 M_☉ was calculated. The stars were initially distributed according to a Plummer density profile with r_h = 0.707 pc. As can be seen, similar to the results of the MC code of Freitag & Benz (2002), our results match qualitatively the results of the gaseous method of Amaro-Seoane et al. (2004) inasmuch as there is a steep rise of M_BH at the time when the cluster densities are largest, which coincides with the time the cluster would formally go into core collapse if there were no IMBH. After that time the BH growth levels off due to the expansion of the core, described in Section 2, leading to the convergence of the IMBH mass to an asymptotic value. Quantitatively, however, there are differences in the onset of the rapid growth phase as well as in the value of the final IMBH mass. The IMBH masses in Amaro-Seoane et al. (2004) are generally larger at late times, while the rapid growth phase is somewhat delayed. Our results follow more closely, and unsurprisingly, the ones by Freitag & Benz (2002) up until shortly after the onset of the rapid growth phase, while they converge to larger M_BH at late times. Part of the reason for this discrepancy with Freitag & Benz (2002) must be related to the larger time-step size compared to the local relaxation time in the inner region of the cluster and, thus, their slower expansion (discussed in Section 3). This is further demonstrated in Figure 4, where we plot the IMBH growth for two different values of the time-step parameter θ_max, which is the maximum deflection angle for all stars (see Freitag & Benz 2001, their Equation (7)). Here we clearly see that increasing the time-step size increases the IMBH mass, while the onset of the rapid IMBH mass growth is delayed. Both effects can be ascribed to the system becoming more and more under-relaxed for larger time steps. First, the core-contraction phase becomes longer, causing the delay of the onset of the rapid growth phase, and, second, the cluster expands slower, increasing the period of high core densities, and, thus, accretion rates, as discussed in Section 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Despite these differences, the asymptotic IMBH masses differ by less than a factor of two and are, therefore, in reasonable agreement with each other, given the very long integration time.

4.2.2. Comparison with BME

We now compare the growth rate of IMBHs in our simulations with the direct N-body simulations of BME. For this comparison we restrict ourselves to runs with a larger number of particles to ensure that the central cusp around the IMBH is sufficiently populated with stars for the MC method to be applicable (see Table 1).

Figure 5 shows the disruption rate as a function of time for runs BME-2 and BME-16 for our MC and the direct N-body code. The qualitative behavior, i.e., the decrease of the disruption rate due to the expansion of the cluster, is well reproduced. However, our rates are systematically larger, always leading to IMBH masses that are larger than in direct N-body simulations (see Table 1). We find the largest discrepancies for the lowest M_BH/M_c ratios, where the total number of disruptions differ by a factor of 2, whereas the disruption rates differ by approximately a factor of 1.5. The difference in the disruption rates becomes quickly smaller for larger IMBH masses, so that for M_BH/M_c = 1% (Figure 5 (left panel)) there is agreement within the error bars. However, our rates are still systematically larger for that case, leading to approximately 50% more disruptions at the end of the simulation compared to N-body simulations. This difference is again smaller for M_BH/M_c ≈ 3% going down to only 25%.

Zoom In Zoom Out Reset image size

There are several possible reasons that may explain larger disruption rates in our simulations compared to direct N-body simulations. First, one has to note that our results in Figure 5 are the averages over nine runs with different initial seeds to generate the cluster, while the disruption rates of Baumgardt et al. (2004a) come from only one realization of a cluster and are time averages. As the run-to-run variations seem to be larger than the time variations in our MC runs, it might be possible that this also applies to the N-body results, in which case the differences in disruption rates could be statistically less significant or even disappear for larger BH masses.

Furthermore, both calculations start from non-equilibrium cluster configurations, consisting of a central massive IMBH and a flat constant density core (Baumgardt et al. 2004a). As the MC method assumes that the cluster is always in dynamic equilibrium, it cannot adequately model the initial phase until such an equilibrium is reached. It might, therefore, be possible that the difference in the number of tidal disruptions is at least in part due to differences in modeling the initial, violent relaxation process. On the other hand, the number of disruptions during the first 100 t_cr account for at most 5% of the total in all runs, indicating that these differences have a rather minor influence on our results.

Another effect is the wandering of the IMBH due to close passages of stars that are not bound to the BH. Although, one would intuitively think that a wandering IMBH would increase the cross section for stars to enter the tidal radius as it covers a larger volume and thus provides a larger cross section (Chatterjee et al. 2002), the wandering also tends to flatten the density profile (Baumgardt et al. 2004a; Amaro-Seoane et al. 2004). As the disruption rate is proportional to the density near the IMBH, either at r_crit if the loss cone is empty, or r_t if the loss cone is full, such a flattening of the density profile could, therefore, cause a lower disruption rate. Similarly, Magorrian & Tremaine (1999) discuss the influence of BH wandering on the disruption rate, arguing that, depending on cluster parameters and BH mass, the disruption rate can be increased or decreased relative to the disruption rate of a fixed central BH.

In order to estimate this effect more quantitatively we assume for simplicity that, in the case of a wandering IMBH, the IMBH remains inside the cusp and the stellar density, n(r), outside of the wandering radius, r_w, is given by n_cusp ∼ r^−7/4, while inside it has a constant value of n_cusp(r_w). We furthermore make the assumption that, in this case, we are in the full-loss-cone regime, i.e., the loss cone is constantly replenished with stars on a dynamical timescale, contrary to the empty-loss-cone regime, where this happens on a relaxation time. The assumption of a full loss cone for a wandering IMBH is justified as long as its motion is fast enough so that the influx of stars is always sufficiently high to completely replenish the loss cone. The situation is of course different for a fixed BH where the loss cone is not replenished but rather depleted on such a short timescale. In effect, we are comparing the disruption rate of fixed BH in the empty-loss-cone regime with that of a wandering BH in the full-loss-cone regime.

The disruption rate for a wandering IMBH in the full-loss-cone regime is simply given by the product of the disruption cross section of the IMBH, Σ_BH, the stellar number density, n_*, and the relative velocity dispersion between IMBH and stars, v_rel, evaluated at the wandering radius, r_w, as (see Binney & Tremaine 2008, their Equations (8)–(123))

$$\begin{equation} \dot{n}_{w}=n_{*}(r_{w})\,\Sigma _{{\rm BH}}(r_{t},r_{w})\, v_{{\rm rel}}(r_{w}), \end{equation} \tag{ 6 }$$

where Σ_BH is given by

$$\begin{equation} \Sigma _{{\rm BH}}(r_{t},r_{w})=4\sqrt{\pi }r_{t}^{2}\left(1+\frac{GM_{{\rm BH}}}{r_{t}\,2\, v_{{\rm rel}}^{2}(r_{w})}\right)\,, \end{equation} \tag{ 7 }$$

and $v_{{\rm rel}}\,{=}\,\ssty\sqrt{(v_{*}^{2}+v_{{\rm BH}}^{2})/2}$ with v_* and v_BH denoting the stellar and the IMBH velocity dispersion, respectively. Using v_BH ≪ v_*, $v_{*}(r)=\sqrt{{GM_{{\rm BH}}}/{r}}$ in the cusp region, r_w ≫ r_t, we obtain

$$\begin{equation} \dot{n}_{w}=4\sqrt{\pi }\, n_{c}\, v_{c}\, r_{i}\left(\frac{r_{i}}{r_{w}}\right)^{5/4}r_{t}, \end{equation} \tag{ 8 }$$

where the subscript c refers to the core, n_c = n(r_i), and v_c = v(r_i).

For a fixed IMBH the disruption rate is given by the influx of stars into the region inside r_crit and can be similarly estimated using (Frank & Rees 1976)

$$\begin{equation} \dot{n}_{{\rm lc}}\simeq n(r_{{\rm crit}})\,\Sigma _{{\rm BH}}(r_{{\rm crit}})\, v_{*}(r_{{\rm crit}})\,\frac{\theta _{{\rm LC}}^{2}(r_{{\rm crit}})}{2}, \end{equation} \tag{ 9 }$$

where θ²_LC = 2r_t/3r is the loss-cone angle, r_crit is the critical radius, and r_t and r_w are replaced with r_crit in Equation (7). The term θ²_lc/2 represents the fraction of stars with velocity vectors pointing into the loss cone for an isotropic velocity distribution. Using the same substitutions as above, Equation (9) can be written as

$$\begin{equation} \dot{n}_{{\rm lc}}\simeq 8\sqrt{\pi }\, n_{c}\, v_{c}\, r_{i}\,\left(\frac{r_{i}}{r_{{\rm crit}}}\right)^{5/4}r_{t}. \end{equation} \tag{ 10 }$$

From the ratio

$$\begin{equation} \frac{{\dot{n}_{w}}}{\dot{n}_{{\rm lc}}}\simeq \frac{12}{8}\left(\frac{r_{{\rm crit}}}{r_{w}}\right)^{5/4} \end{equation} \tag{ 11 }$$

as well as from Equation (8) we see, as we expected, that, in general, the disruption rate is decreasing with increasing wandering radius. Furthermore, when we calculate this ratio using r_w ≈ 4 × 10⁻³ and r_crit ≈ 2 × 10⁻³ from Figure 1 in Baumgardt et al. (2004a) for a cluster with M_BH ≈ 700 M_☉ and N = 8 × 10⁴ (run BME-1), we find a value of ≈0.6. This is very similar to what we obtain for the ratio of disruption rates in BME-16 (see Figure 5). However, Equation (11) does not seem to hold for larger IMBH masses. If we calculate the ratio for, e.g., the cluster with M_BH ≈ 1395 and N = 8 × 10⁴ in the same figure (run BME-2), Equation (11) predicts that the disruption rate for a fixed IMBH should be lower than for a wandering IMBH by a factor of ≈1.6 while from Figure 5 we find that it is mostly larger. The reason for that is probably related to the assumption of a full loss cone when deriving Equation (8), as it is clear that in the limit of very small r_w, and consequently larger IMBH masses, the empty-loss-cone regime, and thus, Equation (10), must be approached. Indeed, in Figure 5 (left panel) we see that the disruption rates seem to converge at late times. Similarly, we find very good agreement in the total number of disruptions between N-body and our MC runs for IMBHs with even larger masses. As $r_{w}\sim \sqrt{M_{*}/M_{{\rm BH}}}$ (Chatterjee et al. 2002), it appears that the effect of a wandering IMBH has a negligible influence on the disruption rate for M_BH/M_* ≳ 1000. Given that in the core of a multi-mass cluster with a central IMBH M_* ≈ 0.6 M_☉ (Baumgardt et al. 2004b), we would expect that an IMBH can be treated as being fixed at the cluster center as long as M_BH ≳ 600 M_☉.

A third effect that could decrease the disruption rate and that is not included in the current MC scheme is the formation of, and the strong interaction with, an IMBH binary that is able to scatter other stars into the outer cluster regions. This mechanism has been recently suggested by Gill et al. (2008) to suppress mass segregation in a multi-mass cluster with IMBH. Gill et al. (2008) also found that the efficiency of this scattering process is relatively independent of the relative mass of the IMBH. Given the trend we see for the difference in the number of tidal disruptions between our simulations and N-body results as a function of M_BH/M_c, it appears that this mechanism is unlikely to be the main reason.

Finally, there is the possibility that, as already argued in Section 4.2.1, the larger disruption rates are a result of the slower expansion of the inner regions caused by the larger time-step size relative to the local relaxation time. However, while this slower expansion might increase the disruption rate in general, it is not at all clear why this should make a larger difference for clusters with lower-mass IMBHs. For instance, when we compare the ratios of the final masses from our runs to the results of the individual time-step code of Freitag & Benz (2002) for the different M_{BH, i} we see from Figure 3 that they are rather similar (0.70 for M_{BH, i} = 50 M_☉ and 0.71 for M_{BH, i} = 500 M_☉). On the other hand, the final mass ratios for runs BME-1 to BME-3 vary dramatically in comparison (from 0.6 to 0.95) for a similar range in initial IMBH masses. Another reason why it is less likely that the differences to the direct N-body results are caused by differences in the expansion rates is given in Figure 6, where we compare the Lagrange radii for run BME-16. As can be seen, apart from very minor deviations, the expansion of the cluster in our simulation agrees well with the N-body results, even for the regions within 0.1 r_h. However, one should also note that, here, only the radii that contain more than 1% of the total cluster mass are shown, as only for those were N-body data available. Larger differences would be expected for smaller radii, especially because for this particular run the 1% Lagrange radius is just outside the cusp region which mostly determines the disruption rate.

Zoom In Zoom Out Reset image size

In general we can say that, despite the significant differences in the total number of disrupted stars, the disruption rates do not differ greatly from the ones from direct N-body simulations and are even in very good agreement for runs with M_BH/M_c ≳ 0.01 given the error bars. Furthermore, from Figure 6 we see that our MC code can reproduce the evolution of the cluster structure, at least down to the cusp region, quite well. We expect that this also applies to other clusters as long as the wandering radius of the IMBH is inside the cusp.

5.1. Comparison to BMH

Here we re-examine a subset of the N-body models of BMH using our MC code to see if we are able to reproduce their results and, in particular, the inner surface brightness slopes mentioned in Section 2. For this comparison we only include runs with large N in order to populate the high-mass end of the initial mass function (IMF) sufficiently well, which is important for the applicability of the MC method. Furthermore, as in BMH, we only include bright stars, defined as main-sequence stars and giants with masses larger than 90% of the turnoff mass, for the calculation of the SBPs to take into account that those stars contribute most (>80%) of the observed light.

Table 2 gives a summary of the resulting cluster structural parameter after the clusters have been evolved to an age of 12 Gyr together with the corresponding results of BMH. In contrast to our single-mass runs we find lower IMBH masses in our runs compared to the direct N-body results. However, this difference is much less pronounced, if not negligible, as it is at most 20 M_☉, resulting in final IMBH masses that deviate at most by 7% from the N-body results. The reason for that is not only because the tidal disruption radii of the stars are generally lower compared to our single-mass cases but, more importantly, the cluster had a much lower density in the core initially and subsequently expanded because of mass loss due to stellar evolution. After the core started contracting again, after about 1 Gyr, it became quickly dominated by massive dark remnants, mostly BHs and massive white dwarfs, that drove out lower-mass main-sequence stars. This decreases the overall disruption rate, because the main-sequence stars which are more easily disrupted given their much larger radii compared to compact remnants, are, on average, much farther out, while the compact remnants, though much closer to the IMBH, are unlikely to get disrupted.

The fact that for these simulations there are fewer disruptions in our MC runs than in the N-body simulations is most probably again related to the wandering of the IMBH. While for a fixed IMBH there is a well-populated cusp, for a wandering IMBH no such clear cusp can be identified. In the N-body simulations the massive dark remnants in the cusp are ejected through strong gravitational interactions with the IMBH, which allows the main-sequence stars and giants to diffuse inward and to come closer to the tidal radius. In addition, due to its motion, the IMBH is also able to come closer to stars that are just outside of the cusp region, which also increases the number of disruptions.

Comparing the cluster structural parameters we find again very good agreement between the MC and direct N-body runs. There are only minor but systematic differences inasmuch as in our MC runs the clusters have systematically larger masses, though only by less than 1.5%, and are more compact, with their half-mass radii differing at most by 10%. The most likely reason is that, in our simulations, we do not model close encounters with stars tightly bound to the IMBH. Those interactions might frequently lead to the ejection of massive stars or remnants, which also contributes to the cluster expansion. In addition, the somewhat larger disruption rates might have also added to the stronger expansion. However, given the minor differences between the MC and N-body results, close encounters with stars tightly bound to the IMBH do not seem to have a significant influence on the cluster structure as a whole.

We also achieve good agreement for the surface density profiles between the two methods. Figure 7 shows the two-dimensional density profiles of bright stars for two clusters at the end of the simulation. These profiles are obtained by averaging five snapshots obtained at 50 Myr intervals. The profiles show only very shallow cusps with power-law slopes α between −0.2 and −0.3, consistent with the N-body results. In order to reduce the amount of noise and obtain a more reliable fit of the inner slope of the SBP, we re-ran run BMH-2 but with twice as many particles. The resulting profile is shown in Figure 7. As one can see, the power-law fit has a slope of α = −0.23, very close to the average of α = −0.25 found in Baumgardt et al. (2005) and agreeing very well with the relation (2) found by Miocchi (2007).

Zoom In Zoom Out Reset image size

5.2. Comparison to Real Star Clusters

The ultimate goal of any cluster simulation is to reproduce observations of real star clusters as closely as possible. For GCs these observations are mostly in the form of photometric data and SBPs, while there are only a few clusters for which well-measured velocity dispersion profiles are available. In the previous section we found that our MC simulations with IMBH as well as the corresponding N-body simulations of Baumgardt et al. (2005) are able to reproduce the intermediate inner surface brightness slopes seen in some Galactic GCs. As discussed in Section 2, an IMBH also influences the SBP in that it produces rather large cluster cores as measured by r_c/r_h. It is, therefore, interesting to see if the combination of slope and concentration we find in our models matches any observations. However, instead of comparing slopes and concentrations quoted in the literature it is more suitable to directly compare SBPs in this case, given that the tidal radii of observed clusters are very uncertain (see, e.g., discussion in Baumgardt et al. 2010).

For this comparison we choose the cluster NGC 5694 for which Noyola & Gebhardt (2006) report an inner slope of 0.19 ± 0.11, close to our fit in Figure 7 (right panel). We modeled the cluster by carrying out a large parameter survey consisting of approximately 600 model calculations, varying the initial number of stars, concentration, virial radii, and IMBH mass. All models had the same IMF as in BMH and were evolved for 12 Gyr. As the orbit of NGC 5694 in the Galaxy is not known we assumed for simplicity that it moves on a circular orbit at its current distance of ≈29 kpc from the Galactic center. Given this distance, it is a rather isolated cluster, and has been speculated to be of extragalactic origin (Lee et al. 2006). We included in our search also clusters that contain no IMBH but, instead, 10% hard binaries, as it has been shown that a cluster with binaries can possess similar shallow surface brightness slopes as a cluster harboring an IMBH (Vesperini & Trenti 2010). We calculated the SBP by converting the stellar radius and bolometric luminosity for each star obtained from the Cluster Monte Carlo (CMC) stellar evolution module binary star evolution (BSE) to V-band luminosity using the standard stellar library in Lejeune et al. (1998). The luminosities were then radially binned similar to Noyola & Gebhardt (2006) and converted to apparent magnitudes using a distance of 35 kpc from the sun. In order to minimize the large fluctuations caused by the brightest giants, we only consider objects with an absolute V-band magnitude fainter than 3. This value is also low enough to ensure that the shape of the SBP does not become significantly biased (see also Giersz & Heggie 2009). Table 3 summarizes the parameter ranges explored.

Table 3. Parameter Ranges Explored

Parameter	Range
N	0.6–1.7 × 106
rvir	2–4 pc
W0	0.8–7
MBH	500–4500 M☉
RG	29.4 kpc
Z	0.0004

Note. The galactic distance, R_G, and metallicity, Z, are from the Harris catalog (Harris 1996).

Download table as:  ASCIITypeset image

5.2.1. Surface Brightness Profiles

Figure 8 shows the resulting SBPs of our best-fit models along with the data from Noyola & Gebhardt (2006), and in Tables 4 and 5 are their initial and final cluster parameters.

Zoom In Zoom Out Reset image size

Table 4. Evolution of the Characteristics of Our Best-fit Model without IMBH

	t = 0	t = 12 Gyr
N	1.30 × 106	1.28 × 106
Mc	8.08 × 105 M☉	4.83 × 105 M☉
rh	2.1 pc	4.4 pc
trh	0.96 Gyr	3.75 Gyr
rtide	308 pc	260 pc
fbin	10%	8.9%

Notes. Here, N is the number of stars and binaries, r_tide is the tidal, or Jacobi, radius, and f_bin is the binary fraction. The initial density distribution is a King profile with W₀ = 3.

Download table as:  ASCIITypeset image

Table 5. Evolution of the Characteristics of Our Best-fit Model with a 500 M_☉ IMBH

	t = 0	t = 12 Gyr
N	1.30 × 106	1.28 × 106
Mc	7.63 × 105 M☉	4.65 × 105 M☉
rh	1.98 pc	3.98 pc
trh	0.90 Gyr	3.29 Gyr
rtide	308 pc	261 pc

Note. The initial density distribution is a King profile with W₀ = 0.8.

Download table as:  ASCIITypeset image

In general, the model SBPs match well within 60 arcsec, while outside, the observed profile seems to flatten. This excess light is usually attributed to both, the confusion with background stars and stars that are escaping but remain still close to the cluster for a considerable time (Fukushige & Heggie 2000), which we cannot model using a simple tidal cutoff prescription. In this case, however, confusion with background stars is a more likely explanation as the estimated tidal radius from the literature (Harris 1996) and the Jacobi radius we obtain at the end of our models is larger by a factor of four compared to the radius where the observed profile begins to flatten.

Comparing the inner profiles, we find that, indeed, within the photometric errors, the data can be well reproduced even without IMBH, while models with IMBH show good agreement if the IMBH mass is less than ≈1000 M_☉. As the latter mass corresponds to a BH-to-cluster mass ratio of 0.2%, with a total cluster mass of 4.7 × 10⁵ M_☉ at the end, our best-fit IMBH masses are in good agreement with the extrapolated relation for SMBH-harboring galaxy bulges.

The good agreement of the SBPs for the cluster model with IMBH is somewhat surprising as, according to Miocchi (2007), NGC 5694 is not expected to harbor an IMBH based on the fact that its concentration is too large for the relatively steep inner surface brightness slope (see relation (2)). However, this statement is based on values for the concentration quoted from the literature (Harris 1996) which, as we mentioned before, are in general rather uncertain due to difficulties determining the tidal radius of a cluster. An error of 40% in the measured tidal radius would bring the value of the slope and the concentration of the cluster in agreement with relation (2). Furthermore, the observed slope has a rather large uncertainty as well, which, even considered alone, could make the observed parameters consistent with this relation. Therefore, it appears that the close match we obtain for the shape of the SBP in our simulation with the observed one of NGC 5694 does not necessarily contradict the validity of the slope–concentration relation found by Miocchi (2007) but seems to be broadly consistent with it.

5.2.2. Time Variability

Since Giersz & Heggie (2009) and Heggie & Giersz (2009), the importance of fluctuations when comparing cluster models with observations has been better understood. More recently Vesperini & Trenti (2010) investigated the variation of the inner surface brightness slopes for clusters with and without IMBH based on N-body models with up to 65, 536 stars, and found that intermediate slopes in the range expected for clusters harboring an IMBH are also ubiquitous among clusters without IMBH, in particular when there is a non-negligible fraction of dynamically hard binaries. However, the models in Vesperini & Trenti (2010) were, owing to the computational expense of direct N-body simulations including binaries, rather idealized in the sense that they only contained a much lower number of particles, by a factor of 20, than are actually present in NGC 5694. In addition, the surface brightness slopes were derived by considering all main-sequence stars while most of the light is actually contributed by only a small subset of main-sequence stars and giants with masses larger than ≈0.7 M_☉ for clusters with ages of around 12 Gyr (see also BMH and Section 5.1). It is likely that these limitations will influence the derived surface brightness slopes.

Noyola & Baumgardt (2011) raised similar concerns. They analyzed a series of low-N, N-body models with and without IMBHs that were stacked on top of each other to contain in the end up to 6 × 10⁶ stars for constructing corresponding HST-like images and then deriving surface brightness slopes the same way as has been done for observed GCs in Noyola & Gebhardt (2006). Contrary to Vesperini & Trenti (2010), their results show that the inner surface brightness slope is a good diagnostic to discern IMBH-less clusters from likely candidates of IMBH harboring clusters, with intermediate values, ranging from −0.1 to −0.4, indicating the presence of an IMBH.

Given these discrepancies and limitations in the literature, it is worthwhile to investigate the slope variations for our full-scale models, where the cluster as well as single and binary evolution is self-consistently taken into account. Although we have not checked whether our MC code is able to reproduce the density fluctuations of an N-body model, Heggie & Giersz (2009) show that the time variations of the core radius in their MC model, which they identify as to be due to variations of the inner density slope, are of the same magnitude as in a corresponding N-body model, though less coherent due to the random sampling of the stellar orbits. It thus appears that the density fluctuations in a Henon-type MC code reflect the degree of variations present in N-body simulations and can at least be used to get some first insight into the temporal behavior of the SBP. Ideally, we would like to follow the inner SBP with a direct N-body simulation for a brief time period around the current age of NGC 5694, but even for such a relatively short time a direct N-body calculation of the entire cluster with N ≈ 1.3 × 10⁶ stars would take a prohibitive amount of time.

For this analysis, we used a simple least-square fit for the determination of the SBP slope over a range that covers the three innermost points of the Noyola & Gebhardt (2006) profile. This range corresponds to 11%–32% of the core radius, and is, thus, similar to the one Vesperini & Trenti (2010) used for their determination of slopes. In order to reduce the uncertainty of the individual fits we used a much finer binning than in Figure 8, with 50 logarithmically spaced bins in the fitted range.

Figure 9 shows the evolution of the inner surface brightness slopes from five snapshots covering a period of 1 Gyr at around 12 Gyr. As can be seen, the slopes of the cluster models without IMBH do not steepen much beyond −0.1, while for the cluster with IMBH they cover a much wider range, from > − 0.1 to −0.4. Comparing with the slopes derived from observations it appears, at face value, somewhat unlikely that the cluster model without IMBH can match the observed slope, and a cluster with IMBH seems more likely to fit. On the other hand, the uncertainties of the individual slopes are, with 0.1, rather large, so the disagreement is only marginally statistically significant. In addition, it might pretty well be that due to our coarse sampling we have missed possible maxima. For this reason we also analyzed four similar models in our grid, with the same initial particle number and radius but different, slightly lower concentrations, going down to W₀ = 1.8 in steps of 0.2 (Figure 10). As one can see, with the exception of one point, all slopes remain shallower than −0.17, the value from our least-square fit to the observations. Although about a third of the values are within the error bars, it is nevertheless remarkable that the average slope is close to zero while for models with IMBH it is clearly steeper than −0.1 and closer to the observed one. Thus, based on our results it remains, at least on a qualitative level, somewhat more difficult to reconcile a cluster model without IMBH with observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Clearly, given the large uncertainties of not only our slopes, but also of the slope in Noyola & Gebhardt (2006) for NGC 5694, a more detailed statistical analysis, taking the SBP directly into account, is certainly desirable in order to quantify how likely or unlikely it is that NGC 5694 harbors an IMBH. As Giersz & Heggie (2009) already pointed out, the agreement between model and observations has to be determined in terms of the probability that an observational profile be rejected as a member of the ensemble of profiles provided by the model. Such an analysis, however, is beyond the scope of the present paper.

5.2.3. Dynamical Age

The advantage of having modeled the evolution of NGC 5694 self-consistently with a realistic number of stars is that it allows us to obtain a more reliable estimate of its dynamical age. The dynamical age has consequences for the expected inner surface brightness slope for an observed cluster (e.g., Noyola & Gebhardt 2006; Trenti et al. 2010), and is essential when assessing possible mass-segregation signatures for IMBHs (Gill et al. 2008).

Hurley (2007) already emphasized that the relaxation time of a cluster is time-dependent and the actual dynamical cluster state, whether in the core-contraction or past core-collapse phase, is better described in terms of $N_{t_{{\rm rh}}}$, the number of elapsed half-mass relaxation times, defined as

$$\begin{equation} N_{t_{{\rm rh}}}=\int \nolimits _{0}^{\tau }\frac{dt}{t_{{\rm rh}}(t)}\,, \end{equation} \tag{ 12 }$$

where τ is the cluster age. He showed that in his models t_rh decreases by a factor of two over the cluster lifetime and compared $N_{t_{{\rm rh}}}$ with simple estimates, dividing the cluster age by either the initial or current t_rh, finding that the former estimate was much closer to $N_{t_{{\rm rh}}}$ than the latter.

Apparently assuming the same behavior for the evolution of NGC 5694 and using an estimate for the current half-light relaxation time of t_rh = 1.9 Gyr from the updated Harris catalog (Harris 1996), Vesperini & Trenti (2010) estimated its dynamical age as $N_{t_{{\rm rh}}}\approx 3$ (see their Figure 3). However, they also point out that such estimated ages have rather large uncertainties, and detailed models are necessary in order to get more accurate results.

In fact, as inspection of Tables 4 and 5 reveals, the simple extrapolation of the Hurley (2007) results to our full-scale models is not justified. This is mainly because NGC 5694 is not tidally limited and can expand freely which increases its relaxation time rather than decreasing it over time. Calculating the dynamical age of our models for NGC 5694 as in Equation (12), we find $N_{t_{{\rm rh}}}= 4.2$ for the best-fit model without and $N_{t_{{\rm rh}}}= 4.7$ for the one with IMBH. Therefore, NGC 5694 appears to be not far from core collapse and to have evolved considerably. Straightforward estimates based on the initial and final relaxation time, on the other hand, give τ/t_rh(0) ≈ 13 and τ/t_rh(τ) ≈ 3, respectively, which means that, contrary to the case in Hurley (2007), using the current instead of the initial half-mass relaxation time gives a better approximation for the dynamical age of this cluster. Considering, however, that even this estimate differs by as much as 40% and the difference still worsens when the cluster gets older and expands further, neither approximation seems really suitable for a reliable estimate.

Coincidentally, if the value quoted from the Harris catalog is corrected for the fact that the projected half-mass radius is smaller by ≈25% than the unprojected one (see, e.g., Hurley 2007), the resulting current t_rh ≈ 2.8 Gyr would imply $N_{t_{{\rm rh}}}=4.4$, which is very close to the integrated value from Equation (12) for our models. Clearly, one cannot attach any meaning to this as the corrected relaxation time still differs significantly from the one in our model, owing to the underlying simplifying assumptions (see Harris 1996), and the close match is rather accidental.

As this comparison shows, there seems to be no simple way to reliably estimate dynamical ages of clusters without considering the previous cluster evolution in more detail.

In this paper, we studied the influence of IMBHs on the evolution of GCs using our MC code, which has been extended to include the effects of a central IMBH on the stellar distribution. The IMBH is treated as a fixed point mass in the cluster center, and at each time step stars are tested for entry into their loss cone. In order to test our implementation we carried out a large number of idealized, as well as more realistic dynamical simulations of GCs and compared our results with published results from direct N-body and Fokker–Planck calculations.

In general we found that our results agree reasonably well with results from N-body and Fokker–Planck-type codes. Significant differences only exist for idealized models while for more realistic clusters these differences become negligible.

We found that our code is able to reproduce the expected density cusp and Keplerian velocity profile around the IMBH for a single-mass cluster very well. In addition, we also found that the evolution of the cluster density profile outside of the 1% Lagrange radius is in very good agreement with direct N-body simulation. The only notable difference is that the cusp extends farther down toward the IMBH than in the direct N-body simulations of Baumgardt et al. (2004a). The reason is most likely that the IMBH is fixed at the cluster center while in the N-body simulations it is allowed to move freely. Although this has a minor effect on the density profile, causing it mainly to be flatter inside r_w, we argue that it produces significant differences in the disruption rates for single-mass clusters. In particular, as we demonstrated by simple estimates of the disruption rate, the effect of IMBH wandering might be able to explain why we see up to a factor of two larger disruption rates for low-mass IMBHs than for higher-mass ones, compared to N-body results, which is difficult to explain otherwise. For a sufficiently fast-moving IMBH, the loss cone can always be assumed full in which case the disruption rate is a strong function of the density at the disruption radius r_t. Therefore, for larger r_w, and consequently for lower M_BH, the density at r_t becomes lower inside the very steep power-law cusp which, in turn, decreases the disruption rate. In our calculations the disruption rate is determined by the diffusion of stars into the loss cone. This rate turns out to be larger than for the case of a wandering IMBH for our cluster parameters if r_w is significantly larger than r_crit, which agrees with the differences we see for run BME-16. For lower r_w, and, thus, larger M_BH, loss-cone depletion effects appear to become important as we find from our simulations that the rates for these two cases converge for IMBHs with M_BH/M_* ≳ 1000, while Equation (11) predicts that the rate for a wandering IMBH should become larger than for a fixed one.

The agreement for larger M_BH also implies that, in a realistic GC, an IMBH with M_BH > 600 M_☉ can be safely treated as a fixed central point mass, given the scaling of the wandering radius and the fact that the average stellar mass in the core converges quickly to M_* ≈ 0.6 M_☉ (Baumgardt et al. 2004b) in multi-mass clusters with a central IMBH. These larger IMBH masses are expected for clusters with M_c ≳ 2 × 10⁵ M_☉ initially, based on a straightforward extrapolation of the well-known M–σ relation (Magorrian et al. 1998; Gebhardt et al. 2000; Ferrarese & Merritt 2000). Since most GCs are more massive than this, we conclude that our code is able to calculate the evolution of GCs with a central IMBH for a wide range of relevant cluster parameters.

Comparing the IMBH mass growth in the empty-loss-cone regime with different Fokker–Planck-type codes, such as the gas code by Amaro-Seoane et al. (2004) and the MC code by Freitag & Benz (2002), we found better agreement with our results. For all codes the IMBH growth shows the same qualitative behavior, that is, the masses converge to a constant value. Only at late times do the masses differ significantly by factors of 1.3–2. The larger IMBH masses in our simulations compared to the results of Freitag & Benz (2002) are most likely due to the slower expansion of the inner regions in our simulations, which are caused by the larger time steps relative to the local relaxation time, a consequence of the shared time-step scheme used in our code. Contrary to the effect of IMBH wandering, this effect does not introduce a dependence on M_BH and our final masses always differ by 30% regardless of the initial M_BH.

The situation for realistic cluster models is, however, rather different. In this case, the disruption rates are generally lower because the cluster core gets quickly dominated by dark remnants, which have extremely small r_t. Therefore, the influence of the IMBH wandering is, in general, much reduced and does not produce significantly different results compared to a fixed IMBH. There are only minor deviations in the final cluster mass and size, with our clusters being more compact by ≈10% and slightly more massive by ≈1%. The reason here might be related to the fact that we do not model close encounters with stars in the cusp, which could be an important ejection mechanism for the stars (Baumgardt et al. 2004a). Such encounters are also responsible to efficiently remove dark remnants from the cusp (Baumgardt et al. 2004b) and, thus, allow for lower-mass main-sequence stars to get close to the IMBH and disrupted. The somewhat larger number of disruptions in the N-body simulations is also caused by the wandering of the IMBH as it can, this way, get closer to the main-sequence stars that are, due to mass segregation, farther out. Although the total mass the IMBHs gain in the N-body simulations is up to three times larger than in our simulations, the total difference amounts to just 20 M_☉ at most or less than 7% of the total cluster mass and is, thus, very minor, given that this difference comes from not much more than 10 disrupted main-sequence stars.

Using our MC code we are also able to reproduce the intermediate power-law slopes in the SBP seen in all N-body simulations with IMBH in Baumgardt et al. (2005). Although there is considerable scatter in the profile which makes the determination of the individual slopes rather uncertain, we show that they, nevertheless, are within the same region as in the N-body simulations. Using twice as many stars in our simulations than were used in Baumgardt et al. (2005) we have a high enough resolution to obtain a reliable fit, which turns out to be close to the average slope found by Baumgardt et al. (2005). This, once again, confirms that, despite differences in the details of the dynamics in the innermost cusp region, the overall cluster structure and evolution is rather well reproduced.

Finally, by carrying out a large parameter survey, we were able to model, for the first time, the SBP of NGC 5694 in Noyola & Gebhardt (2006) and constrain the maximum mass of a possible IMBH in the cluster center to <1000 M_☉, which, in this case, corresponds to M_BH/M_c ≲ 0.2%. Given the photometric errors in the observed profile, both models with and without an IMBH can be constructed to match observations. However, considering the surface brightness slopes, there is a slight preference for models with an IMBH, as clusters without IMBH rarely had slopes as large as observed. This was in spite of the presence of a significant population of binaries (≈8% in the end), which has produced larger variations in the direct N-body simulations of Vesperini & Trenti (2010) containing a much lower number of stars. Clearly, taking into account the rather large uncertainties in the individual fitted slopes, a more thorough statistical comparison of the SBPs is required to better assess the likelihood of the presence of an IMBH in NGC 5694. So far the evidence for an IMBH is only marginal. Moreover, the existence of an IMBH in this particular cluster seems inconsistent with the static models of Miocchi (2007), with NGC 5694 having a concentration that is too large for the measured inner surface brightness slope. However, considering possible errors in the value for the concentration quoted in the Harris catalog (Harris 1996) and the rather large uncertainty in the inner surface brightness slope, we find that this is not yet conclusive either.

In summary, we find that our MC code, which now includes the effects of a central IMBH on the stellar distribution, is able to reproduce the overall structure and evolution of single-mass as well as of more realistic multi-mass cluster models reasonably well compared to direct N-body simulations. There are differences in the disruption rates and IMBH masses which are most likely related to the wandering of the IMBH, which is not included in the MC code. However, these differences become almost negligible for realistic cluster models as the disruption rates are generally very low for these cases. In addition, since the IMBH wandering radius scales with M_BH/M_* as $r_{w}\sim \sqrt{M_{{\rm BH}}/M_{*}}$ (Chatterjee et al. 2002) the influence of the IMBH motion is reduced for more massive clusters with a similar M_BH/M_c. One process that still plays a role for more massive clusters and which is not included in our MC simulations are the close gravitational interactions with cusp stars which we plan to implement in the near future.

We thank Holger Baumgardt for providing us with data as well as initial conditions of the N-body simulations and helpful discussions. We also thank Rainer Spurzem for very helpful comments and discussions that finally led to the derivation of $\dot{n}_{w}/\dot{n}_{lc}$. This work was supported by NSF Grant AST-0607498 at Northwestern University. J.M.F. acknowledges support from Chandra Postdoctoral Fellowship Award PF7-80047. Support for program HST-AR-11779 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This research was supported in part by the National Science Foundation under grant No. PHY05-51164. This work was completed while all authors were participants of the program "Formation and Evolution of Globular Clusters" at KITP in spring 2009.