© ESO, 2013

1. Introduction

The dynamical evolution of globular clusters (GCs) in a tidal field causes them to gradually dissolve (e.g. Ambartsumian 1938; Spitzer 1958; Baumgardt 2001). After GCs have migrated out of their gas-rich birth environment (Shapiro et al. 2010; Elmegreen 2010; Kruijssen et al. 2012), this dissolution is driven by two-body relaxation, that is, the repeated effect of soft encounters between stars, which leads to energy equipartition and populates the high-velocity tail of the Maxwellian velocity distribution (e.g. Portegies Zwart et al. 2010). Because escape is driven by a tendency toward energy equipartition between stars of different masses and a mass-dependent encounter rate, the escape probability is not the same for all stars and generally decreases with mass (Hénon 1969; Vesperini & Heggie 1997; Takahashi & Portegies Zwart 2000; Baumgardt & Makino 2003; Kruijssen 2009). As a result of mass segregation, high-mass stars sink to the center of a star cluster, while low-mass stars are pushed to the outskirts, where they are easily removed by the Galactic tidal field. This changes the slope of the mass function (MF) and the mass-to-light ratio (M/L) of the GC (e.g. Richer et al. 1991; De Marchi et al. 2007; Mandushev et al. 1991; Kruijssen & Mieske 2009; Baumgardt & Makino 2003).

How exactly the content of evaporating clusters evolves depends on their initial composition and global properties. Kruijssen (2009) showed that the most massive objects in a cluster play a key role in determining which bodies escape. It is therefore important at which point the bulk of the dynamical mass loss occurs, that is, early on (when massive stars are still present) or after ~1 Gyr (when all stars with masses m > 2  M_⊙ will have died). In the latter case, the cluster evolution will be dominated by stellar remnants, because at such old ages they generally have masses that exceed the masses of the most massive luminous stars. During their formation process, remnants receive velocity kicks (e.g. Lyne & Lorimer 1994; Pfahl et al. 2002; Moody & Sigurdsson 2009; Fregeau et al. 2009), and hence the fraction of remnants that is retained depends on the depth of the potential well M/r_h (see Kruijssen 2009 for a quantitative analysis), as well as on the magnitude of the kick velocity.

If the central density of a young GC is high enough, runaway collisions between stars may be able to produce an intermediate-mass black hole (IMBH, see Portegies Zwart et al. 2004; Gürkan et al. 2004). However, recent studies have shown (e.g. Glebbeek et al. 2009) that stellar mass-loss complicates the formation of IMBHs through runaway collisions. In the context of the above discussion, it is clear that the presence of an IMBH could potentially have an observable influence on the evolution of the GC composition. Thus far, the search for IMBHs has focused on direct detection through kinematics (e.g. Baumgardt et al. 2003; van den Bosch et al. 2006; Noyola et al. 2008; Lützgendorf et al. 2011), surface brightness profiles (e.g. Noyola & Gebhardt 2006), or the search for Bondi-Hoyle accretion in the centers of GCs (e.g. Strader et al. 2012). An analysis of how the presence of IMBHs changes the stellar composition of GCs would provide an independent tracer, or at the very least facilitate obtaining better criteria for the selection of target GCs in observational work.

Previous numerical work on IMBHs in GCs has had a natural emphasis on the internal evolution of GCs. Baumgardt et al. (2005) and Noyola & Baumgardt (2011) found that the surface brightness profiles of clusters hosting an IMBH exhibit weak central cusps, in contrast to core-collapsed clusters with very steep profiles and pre-core collapsed systems with no cusp at all. However, Trenti et al. (2010) and Vesperini & Trenti (2010) showed that the shallow cusp may also form as a transition state of GCs undergoing core collapse and therefore cannot be a sufficient criterion for a GC hosting an IMBH. Furthermore, it has been shown that IMBHs prevent a cluster from undergoing core collapse and reduce the degree of mass segregation compared with non-IMBH clusters (Baumgardt et al. 2004; Gill et al. 2008). Trenti et al. (2007) predicted that clusters with a high ratio of core radius to half-mass radius are good candidates for hosting an IMBH. This was challenged by Hurley (2007), who showed that the ratios observed for Galactic GCs can be explained without the need for an IMBH when treating model data as if they were observational data. Most of these studies focused on internal observables such as the surface brightness profile and hence did not consider tidal effects. Because of the importance of a tidal boundary for the eventual composition of a GC, our simulations do include a tidal field. A first effort in this direction was made by Trenti et al. (2010), who included a single N-body run with an IMBH in their model grid of tidally dissolving clusters.

In this paper, we address the dynamical evolution of GCs with direct N-body simulations, serving two particular aims. Firstly, we study the evolution of the composition of dissolving clusters itself, covering a substantial part of the parameter space. In doing so, we vary quantities that are important but have not received much prior attention, such as the remnant retention fraction and the initial binary fraction. Secondly, we intend to identify global tracers of the presence of IMBHs or massive remnants. We discuss the evolution of the characteristics of GCs and compare them with observational results. Wherever possible, we identify potential targets for future observational searches for IMBHs.

The structure of the paper is as follows: in Sect. 2, we introduce the N-body code and the initial conditions of the simulations. Section 3 is dedicated to the influence of different IMBH masses, black-hole retention fractions, and binary fractions on the properties of tidally dissolving clusters. In Sect. 4, we compare our numerical results with the observed distribution of GCs in the mass function slope versus concentration plane. The paper concludes in Sect. 5, where we summarize our work and present our conclusions.

2. Models

In this section we describe the survey of N-body simulations that are treated in this paper.

2.1. N-body simulations

We ran N-body simulations based on the GPU (Graphic Processing Unit)-enabled version of the collisional N-body code NBODY6 (Aarseth 1999; Nitadori & Aarseth 2012) on GPU graphic cards at the Headquarters of the European Southern Observatory (ESO) in Garching and the University of Queensland in Brisbane. This code uses a Hermite integration scheme with variable time steps. Furthermore, it treats close encounters between stars by applying KS (Kustaanheimo & Stiefel 1965) and chain regularizations and accounts for stellar evolution (Hurley et al. 2000). The regularization procedures are crucial for following orbits of tightly bound binaries over a cluster lifetime accurately and treat strong binary-single and binary-binary interactions properly. The simulations were carried out with particle numbers of N = 32  768 (32k), 65  536 (64k), and 131  072 (128k) stars.

The initial density profile is given by a King (1966) model with a central concentration W₀ = 7 and an initial half-mass radius of r_h = 1 pc. This is preferable to a tidally limited cluster where the tidal radius from the King model is set equal to the Jacobi radius of the tidal field, which has proven to be unrealistic for modeling GC-like densities, since the resulting initial half-mass radii are on the order of 4−5 pc. This has the effect that the clusters are too extended at the beginning of their evolution and undergo core collapse only at the very end of their evolution. Furthermore, after expansion due to mass loss through stellar evolution, these clusters attain half-mass radii of 8−10 pc, much larger than what is observed in GCs today. For this reason we fixed the initial half-mass radius to 1 pc. This is more time-consuming because the smaller crossing times increase of the numbers of time steps.

All models started with a Kroupa (2001) initial mass function (IMF, ξ(m)   dm ~ m^α   dm) that has a Salpeter-like power law slope of α = −2.3 for stars more massive than 0.5  M_⊙ and a slope of α = −1.3 for stars between 0.08 and 0.5  M_⊙. For the lower and upper mass limit we chose 0.1 and 100  M_⊙, respectively. This configuration leads to an initial mean mass of ⟨ m ⟩  = 0.63  M_⊙. Primordial binaries were included for some models (see Sect. 2.3) by randomly sampling the desired fraction of stars in the cluster to binaries with a thermal eccentricity distribution (f(e) = 2e) and a flat period distribution between 1 and 10⁶ days. The median ratio of the binary binding energy to the mean particle kinetic energy for these particular initial conditions is e_bin/e_kin = 2−4 and decreases with increasing particle number.

Our clusters were set on a circular orbit around an external galaxy that follows the potential of an isothermal sphere (Φ(r) ∝ lnr). The distance to the Galactic center was chosen to be 8.5 kpc and the circular velocity of the galaxy V_G was set to 220 km   s^-1  for all our simulations because we study Galactic GCs only. The particle integration was performed in an accelerated, but non-rotating reference frame where the cluster remained in the center and the galaxy moved around the cluster on a circular orbit. The forces generated by the Galactic potential and the stellar gravity were both applied to each star when it had advanced on its orbit through the cluster.

Table 1 lists the initial parameters for our set of models, which from left to right are the particle number (N), IMBH-mass fraction (M_•/M₀), black-hole retention fraction (f_ret), binary fraction (f_bin), and dissolution time (t_diss, defined as the time the cluster needs to lose 95% of its initial mass, see Sect. 3.1). Due to an error in our mass-randomizing routine, some of the models with the same particle number contain the exact same mass. However, since every model starts with different initial parameters, this does not affect the individuality of each model.

2.2. Reference models

In Fig. 1 we depict the evolution of our reference models with N = 32k, 64k, and 128k stars. Shown are the bound mass fraction (a), the mean slope of the mass function α in the mass interval 0.3  M_⊙ < m < 0.8  M_⊙ (b), the Lagrangian radii r_10% (c) and r_50% (d), which include 10% and 50% of the cluster mass, respectively (the latter is the half-mass radius r_h), the remnant mass fraction (e), and the number of stellar-mass black holes (f) as a function of time.

The evolution of the bound mass fraction in Fig. 1a shows a similar behavior for all models. After a steep drop from 1.0 to 0.7 in the first Gyr, which results mainly from the mass loss due to stellar evolution, the mass fraction drops almost linearly until the last few Gyr of the cluster lifetime. After this, the mass loss slows down. Apart from the increasing lifetime with particle number we did not find significant differences among these models in terms of their mass evolution.

The next quantity we show is the evolution of the mass function slope α. This is achieved by binning the mass intervals logarithmically, going from m = 0.2  M_⊙ to 10  M_⊙ in 20 bins. We determined the slope of the mass function by fitting a power law to the points with 0.3  M_⊙ < m < 0.8  M_⊙, because this is the same fitting area as used for the observational data. The result is shown in Fig. 1b. As expected, the slope of the mass function increases with the age of the cluster because the cluster loses low-mass stars through two-body relaxation. As with the bound mass fraction, there are no significant differences between models with different particle numbers other than the time scaling (when these values are plotted as a function of the fraction of the total lifetime t/t_diss all three models agree).

Fig. 1 Evolution of the reference models. Shown are a) the bound mass fraction, b) the slope of the mass function α (for 0.3  M⊙ < m < 0.8  M⊙), c) the Lagrangian radii r10% and d) r50%, e) the remnant mass fraction, and f) the number of stellar-mass black holes as a function of time. We stress that the absolute time on the x-axis is arbitrary because it depends on the total mass and the galactocentric distance and is therefore only shown for relative comparison. Dashed lines in panel f) mark the logarithmic fits to the evolution of the black hole numbers.

Plots 1c and d depict the evolution of the two Lagrangian radii that include 10% and 50% of the cluster mass, respectively. The core collapse is indicated by a drop in the 10% Lagrangian radius followed by an immediate rise due to binary heating. Surprisingly, the exact occurrence of that event does not seem to depend on the particle number. While the 32k model undergoes core collapse first, the 128k core collapse still occurs before the 64k model reaches this point. We assume the reason for the absence of a clear trend for the occurrence of core collapse with particle number arises from stochastic effects caused by the low number of massive stars, which mainly drive the dynamical evolution of the cluster even for the 128k models.

The last two plots (Fig. 1e and f) show the evolution of the stellar remnants in the clusters. In Fig. 1e the remnant mass fraction increases with time. This is because low-mass stars are ejected more efficiently than high-mass stars. Consequently, the fraction of high-mass stars to low-mass stars constantly increases. The steep rise at the beginning is due to initially rapid stellar evolution and the consequent formation of stellar remnants. The evolution of the stellar-mass black hole numbers in the clusters in Fig. 1f shows a roughly exponential behavior. For a model with 128k stars, the number of black holes already drops to half of the initial value within 1 Gyr. After 5 Gyr all black holes have been ejected from the cluster. To lead the eye, we applied rough logarithmic fits to the evolution of black-hole numbers (dashed lines). This agrees well with previous results (e.g. Sigurdsson & Hernquist 1993; Kulkarni et al. 1993; Miller & Hamilton 2002; O’Leary et al. 2006), which show that stellar mass-black holes eject themselves from the cluster core by building a dynamically decoupled subcore through mass segregation (Spitzer 1969). A recent study by Sippel & Hurley (2013) has shown, however, that it is possible to retain a large number of stellar-mass black holes in the cluster after a Hubble time with even small retention fractions of 10%. The initial half-mass radius in their simulations was a factor of 6 larger than our values. This leads to longer crossing and relaxation times and consequently to a delayed dynamical evolution of the cluster. This is shown by the fact that the models of Sippel & Hurley (2013) do not undergo core collapse after 10 Gyr. Therefore, the ejection rate of the stellar-mass black holes is also low. We note, however, that models with initial half-mass radii of several parsec end up with physical sizes too large compared with those of the majority of galactic GCs after their dynamical evolution. A comparison of panel c) and f) indicates that the core collapse in all models occurs after the cluster has lost its stellar-mass black holes.

2.3. Model families

We ran models with different IMBH masses, black-hole retention fractions, and binary fractions, which we sorted to three different families. To use realistic values for the IMBH masses we considered several scaling relations such as the M_• − M_bulge relation (Magorrian et al. 1998; Häring & Rix 2004) and the M_• − σ relation (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Gültekin et al. 2009). When these empirical correlations are extrapolated, they predict IMBH masses of M_• ~ 0.001−0.03  M_tot in GCs (Goswami et al. 2012). This corresponds to black-hole masses M_• ~ 10²−5 × 10³  M_⊙ in an average GC with a mass of M_tot = 1.5 × 10⁵  M_⊙. These masses were found by Portegies Zwart et al. (2004), who used N-body simulations to model the runaway merging of massive stars in the dense star cluster MGG-11. Their values, however, should be treated as upper limits, since stellar mass loss of the “runaway” star was not included (see Glebbeek et al. 2009). Taking into account these predictions, we used two different black-hole masses: 1% and 3% of the total cluster mass. For models with black-hole masses <1% the actual black-hole mass is <800  M_⊙ even in simulations of N = 128k stars and the central black hole has a high chance of being ejected by encounters with massive stars.

The retention fraction gives the fraction of stellar remnants that remain in the cluster after their formation. It is assumed that very compact objects such as black holes and neutron stars suffer a kick velocity during their formation that is often higher than the escape velocity of the cluster (Lyne & Lorimer 1994). Theoretical estimates for the retention fraction of neutron stars range from 5% to 20% (Drukier 1996). For black holes, however, the range of possible values is not well-constrained since the physics of the formation of black holes is still not fully understood. For this reason, we assumed a constant retention fraction of 10% for neutron stars and varied the retention fractions of the black holes from 0% (no black hole remains in the cluster) to 100% (all black holes that are formed remain in the cluster). For white dwarfs we always assumed a retention fraction of 100% since they do not suffer any kick velocities (Davis et al. 2008).

Binaries are important for the dynamical evolution of GCs and need to be considered in our models. To test how primordial binaries affect the models, we introduced a third family, in which we varied the fraction of binaries in the initial conditions of the model from 10% to 30%. This covers the current typical binary fraction in GCs, which is about 10% (e.g. Rubenstein & Bailyn 1997; Davis et al. 2008; Dalessandro et al. 2011). Since calculating the individual time steps of the binary systems is very time-consuming, we limited this exercise to a maximum particle number of 64k for the models with primordial binary fractions.

3. Observables

The main goal of this study is finding observables in GCs that give information about the presence of a possible IMBH at their centers. In this section we study the effects of IMBHs on their environment to pin down important tracers for central black holes in GCs. We especially concentrate on mimicking observations, that is, measuring cluster properties in the same way as observers would do to compare our results with observational studies.

3.1. Lifetimes

The lifetime, or dissolution time t_diss, gives information about the dynamical evolution of the model as a function of its internal properties. In this work, we define the dissolution time as the time when the cluster has lost 95% of its initial mass.

In Fig. 2 we compare the lifetimes of our model families with the reference models (green bars). As expected, there is a clear trend with particle number and slight differences (of about ~10−20%) between models with different initial conditions. The higher the mass of the central IMBH, the shorter the lifetime of the cluster. The opposite trend is observed for models with different black-hole retention fractions. The model with the longest lifetime is the model with a black-hole retention fraction of 100%. However, it is remarkable that the difference in lifetimes is only about 20% and therefore does not show a strong dependence on whether the cluster hosts an IMBH, a large number of stellar-mass black holes, or primordial binaries.

Fig. 2 Lifetimes of all clusters as a function of their initial particle number.

The primordial binary fraction weakly affects the lifetime of the cluster model. Encounters with binaries can result in flybys, exchanges, ionizations, and sometimes mergers (see Lützgendorf et al. 2012). Some of these encounters result in high-velocity stars that exceed the escape velocity of the cluster. This can lead to a higher ejection rate than in clusters without tight binaries.

We compared our results with the lifetimes found in Baumgardt & Makino (2003), who studied N-body simulations of GCs in tidal fields with different galactocentric distances and eccentricities. With a neutron star retention fraction of 100% and no stellar-mass black holes, the initial conditions assumed by these authors slightly differ from ours. However, since their stellar masses are cut off at 15  M_⊙, the effective retention fraction of neutron stars is lower and the models from their family II are comparable with our f_ret = 0% models (see Table 1). Except for the 32k model, the lifetimes of the models in this study are ~10% longer than those found in of Baumgardt & Makino (2003). These discrepancies are most likely caused by the underfilling initial conditions of our simulations in contrast to the tidally filling clusters of Baumgardt & Makino (2003).

3.2. Remnant fractions

Another quantity we considered is the fraction of the cluster mass constituted by stellar remnants. As remnants we counted all stars within one tidal radius classified as either a white dwarf, a neutron star, or a black hole. The results are displayed in Fig. 3, where we compare the remnant fractions for models with different black-hole masses, retention fractions, and primordial binary fractions as a function of the completed lifetime fraction (i.e. t/t_diss). The remnant fractions increase monotonically. This occurs because of the increasing number of remnants as a result of stellar evolution and also because low-mass stars are ejected first. For all models, this increase continues until the end of the cluster’s lifetime. For models with non-zero IMBH masses and primordial binary fractions, however, the remnant fraction shows a slower increase than that of the reference model. For models with a high fraction of stellar-mass black holes, the remnant fraction does not exhibit a strong variation. An explanation for the lower remnant fraction in models with IMBHs and primordial binaries could be that the lower degree of mass segregation caused by the IMBH distributes massive stars and remnants farther outside in the cluster where they can be more easily ejected. For binaries, the effect might be caused by the higher survival chance of low-mass stars in binaries. Owing to their higher total mass, they remain in the cluster longer than to single low-mass stars, and the remnant fraction increases more slowly.

In summary, remnants are ejected more efficiently from a cluster with a central massive black hole, while no notable effect is found for clusters with different black-hole retention fractions. This is in contrast with the result found by Kruijssen (2009), who used semi-analytic models to study the evolution of the stellar mass function in star clusters. These authors proposed that the remnant fraction could also decrease toward the end of a cluster’s lifetime if a large portion of stellar-mass black holes were retained.

Fig. 3 Remnant mass fractions as a function of the completed lifetime fraction for simulations with different a) IMBH masses, b) black-hole retention fractions, and c) primordial binary fractions. All models contain 128k particles, except for the binary-fraction models, which are limited to 64k stars.

Fig. 4 Mass-function slopes (for 0.3  M⊙ < m < 0.8  M⊙) as a function of the completed lifetime fraction for simulations with different a) IMBH masses, b) black-hole retention fractions, and c) primordial binary fractions. All models contain 128k particles, except for the binary-fraction models, which are limited to 64k stars. The dashed lines indicate the true mass function slope, while the solid lines indicate the mass function slopes as measured by an observer, that is, counting binaries as single stars.

3.3. Mass function slope

The slope of the mass function is a measure for the evolutionary state (or dynamical age) of a GC because it gives information about how many stars and which types of stars were ejected. We used the method described in Sect. 2.2 to derive the mass function slopes in the mass range of 0.3  M_⊙ < m < 0.8  M_⊙ for our models. We stress that our goal is to compare the derived properties of the simulations with observations. For binaries, we therefore took the total luminosity of the system and transferred this into mass using the existing M(L) relation of single stars in the cluster. In this way we ensured that the binaries contribute to the mass function as they would in observations where the binary components cannot be separated. In Fig. 4 we plot the slope of the mass function as a function of the completed lifetime fraction (solid lines). For comparison we overplot the “true” evolution of α (dashed lines) as it would be measured if all binary stars could be resolved. As expected, low-mass stars gain high velocities due to two-body relaxation and escape the cluster. Therefore, the slope of the mass function at the low-mass end changes from negative to positive and increases monotonically. This behavior is observed for all our simulations. However, the details of the change in the slope differ from model to model.

As shown in panel a) of Fig. 4, the IMBH mass has a significant influence on the mass function slope. This most probably arises from the fact that a central black hole ejects high-mass stars more efficiently than the non-IMBH model, as already demonstrated in the previous section. In panel b) of Fig. 4 the evolution of α seems to be very similar irrespective of how many stellar-mass black holes are retained by the cluster. Only for a retention fraction of 100% a slight difference in α becomes noticeable.

Binaries also cause a slower rise of the mass function slope and therefore a lower depletion in low-mass stars. As shown in panel c) of Fig. 4, the simulation that starts with a primordial binary fraction of 30% in the end has the lowest true value of α of all simulations (dashed lines). As described in the previous section, this can be explained by the lower ejection rate of low-mass stars that reside in binaries. However, this is only measurable when all binaries are resolved. A large portion of binaries can bias the observed value of the mass function slope to higher values (solid lines). The differences increase with cluster lifetime and range from Δα = 0.05−0.5 and Δα = 0.2−0.7 for the 10% and the 30% binary fraction model, respectively.

Fig. 5 Characteristic radii ratios (r10%/r50%, i.e., deprojected radii containing 10% and 50% of the stars) as a function of the completed lifetime fraction for simulations with different a) IMBH masses, b) black-hole retention fractions, and c) primordial binary fractions. All models contain 128k particles, except for the binary-fraction models, which are limited to 64k stars. The dashed lines indicate the rL,10%/rL,50% (deprojected radii containing 10% and 50% of the mass) evolution when taking the Lagrangian radii.

3.4. Structural properties

The final quantities of interest are the cluster’s structural properties, such as the characteristic radii containing 10% and 50% of its stars. Again, we aim to compare these quantities with observations and therefore derive the values in the same way as we do for the observed clusters. To obtain these observables in a non-parametric and consistent way, we applied several steps: First, we derived a projected surface brightness profile for stars brighter than M_v = 8.5 mag in 15 logarithmic bins from r_min = 1″ (by placing all our clusters at a distance of 10 kpc, this corresponds to ~0.05 pc) out to the tidal radius of the cluster. In each bin the stellar number density is computed for three orthogonal projected directions. For the final density we took the average of the three projected directions and used the scatter as the error. With this method we derived a projected density profile similar to that observed for a real cluster. As a next step we deprojected the stellar density profile of each cluster using the projected profiles and a multi-Gaussian expansion method (Emsellem et al. 1994) to compare this with observations. This parametrization has the advantage that a smooth deprojected profile can be easily obtained since the deprojection of any Gaussian again is a Gaussian. By integrating the deprojected profiles, we obtain radii that include different percentages of the total number of stars, such as r_50% and r_10%.

In Fig. 5 we show the evolution of the ratio r_10%/r_50% with time. One prominent feature of this plot is the strong dip in the ratio at about t/t_diss = 0.2 (~5 Gyr) that indicates core collapse. The plots demonstrate that for the standard model with a conservative black-hole retention fraction of 30% and no primordial binaries or IMBH, core collapse is a very prominent feature in the evolution of the structural parameters (see Fig. 5a, green line). The presence of an IMBH at the center of the cluster entirely prevents core collapse and r_10%/r_50% is nearly constant over the cluster’s lifetime (Fig. 5a). A similar effect arises if a large number of stellar-mass black holes is present. As shown in Fig. 5b, the model with a retention fraction of 100% undergoes core expansion instead of core collapse, which is delayed until the final stages of its evolution. This underlines the dependence of the time of the core collapse on the black-hole retention fraction. Surprisingly, the binary fraction does not seem to have a strong effect on the concentration of the cluster. Core collapse occurs at the same time as for the model without any primordial binaries, and the overall evolution of r_10%/r_50% is very similar. However, this result should be treated with care because of the low-number statistics of the 64k models.

To estimate the degree of mass segregation in our models, we furthermore considered the evolution of the Lagrangian radii r_L,10% and r_L,50%. These quantities indicate the radii that contain 10% and 50% of the total cluster mass (i.e., also including remnants), respectively. Because of two-body relaxation and mass segregation, remnants sink to the center of the cluster, causing a higher concentration of mass. This causes the ratio of r_L,10% and r_L,50% to be lower than r_10% and r_50% that were derived from the stellar distribution alone. Therefore, studying the difference between the two ratios gives an estimate of the degree of mass segregation. In Fig. 5 the evolution of r_L,10%/r_L,50% is indicated by the dashed lines. It is notable that the ratios of the characteristic radii differ for all our models after some time of cluster evolution. This indicates that mass segregation is present in all the clusters, but its degree differs from model to model. Especially after core collapse, the remnant segregation is much higher than observed in the visible stars for clusters without IMBHs or low stellar-mass black hole retention fractions. The lowest mass segregation is observed for models with non-zero IMBH masses. This result is of great interest for understanding the connection between observed quantities and the actual dynamical stage of the GC.

4. Comparison with observed GCs

This section describes the comparison of GC N-body simulations with observed properties. This may help to explain empirical correlations and identify new cluster candidates for hosting IMBHs.

4.1. α − c plane

De Marchi et al. (2007) and Paust et al. (2010) studied the mass function slopes α as a function of the concentration parameter c for a total sample of 36 Galactic GCs. As mentioned in the previous section, the slope of the mass function is a measure of the dynamical age of the cluster. The higher the slope (i.e. higher value of α), the more low-mass stars were already lost through two-body relaxation in the external tidal field. Therefore, one would expect a direct correlation between this slope and other quantities that are correlated with the dynamical time scale of the cluster (see Fig. 4). One of these quantities is the concentration c, which is a measure for the central density of the cluster. The more advanced the cluster evolution, the closer it is to core collapse, and in turn, the higher is its central concentration. The observations, however, contradict this picture by showing the opposite relation between mass-function slope and concentration¹.

To measure the concentration parameter c that we can compare with the observations, we computed a projected density profile (as described in Sect. 3.4) and fit an isotropic King (1966) model to the simulation data. This allowed us to derive the structural parameters (e.g. c, r_c, and r_h) from the best fit.

Fig. 6 α – c plane (left panels) and α − r10%/r50% plane (right panels) of observed GCs compared with N-body simulations. Shaded areas mark the 3σ limit of the simulations. Clusters falling within these limits and not overlapping with other models are considered to be well-represented by these models. The error bar in panel a) is representative for the uncertainties in α for the observed clusters. The labeled clusters are discussed in Sect. 5.

To compare with the work of De Marchi et al. (2007) and Paust et al. (2010), we derived α (for 0.3  M_⊙ < m < 0.8  M_⊙) and c at points in time spaced by 500 Myr from the cluster models. Figure 6 (left panels) displays the tracks on top of the observed data of De Marchi et al. (2007, filled circles) and Paust et al. (2010, open circles). As before, we compare models with different IMBH masses, black-hole retention fractions and primordial binary fractions. The figure shows that our standard model (without an IMBH and retention fraction of 30%) undergoes core collapse after several Gyr and moves back to lower concentrations during its evolution (green line). Furthermore, the figure shows that the cluster without a black hole behaves in the expected way and c increases with α before it undergoes core collapse. After core collapse, however, the dependence of these two values is reversed because the cluster core slowly expands again due to binary heating. Comparing the green line with the data, it seems plausible that clusters with medium-to-high concentrations and relatively unevolved mass functions can be represented by a no-IMBH model with a medium retention fraction. Clusters with central massive black holes are presented by the yellow and the red lines. These models do not undergo core collapse and represent the data points with low concentrations very well. A similar agreement was found by Trenti et al. (2010) for a single IMBH run.

The suppression of core collapse can also be produced by many stellar-mass black holes, as demonstrated by the model with the 100% black-hole retention fraction (Fig. 6c, blue line). The high concentration of black holes in the center also enlarges the core of the visible stars and prevents core collapse because it functions as an extra heating source. As already observed for the structural properties, binaries do not affect the concentration of the cluster and the tracks are similar for binary and non-binary models.

4.2. α − r_c/r_h plane

Another and preferable way of comparing these data with the models is using the derived non-parametric structural parameters r_10% and r_50% instead of the concentration parameter. The ratio of these radii is a more reliable indicator than the ratio provided by the parametric approach, because a King model does not always provide a good fit to the data after core collapse.

To extract the characteristic radii from the observational data we took all clusters that are included in the sample of Paust et al. (2010) and De Marchi et al. (2007) and applied the same procedure as for the N-body data. We matched their clusters with the sample of Trager et al. (1995) and McLaughlin & van der Marel (2005) to obtain surface-brightness profiles and masses for as many clusters as possible. The overlap is substantial and we obtained a sample of 22 clusters. For each cluster we then deprojected the surface brightness profile using the smoothed profiles from Trager et al. (1995) and the multi-Gaussian expansion method described in Sect. 3.4. We integrated the deprojected surface-brightness profiles to obtain the radii r_10% and r_50% and used the total masses obtained by McLaughlin & van der Marel (2005).

From Fig. 6 (right panels) we conclude that the ratio of the two characteristic radii is a good replacement for the concentration parameter. The observed data follow a similar distribution as in the α − c plot (i.e., there appears to be a lack of clusters with low r_10%/r_50% value and high α). As already concluded from the α − c analysis, models with large cores and low values for the mass function slope are good candidates for hosting a central IMBH. However, the degeneracy with the models with a 100% black-hole retention fraction needs to be treated with caution. The agreement between the observations and the different N-body runs can be exploited to identify candidate GCs for hosting an IMBH. We describe this in the next section.

5. Summary and discussion

We have investigated the effect of intermediate-mass black holes, stellar-mass black hole retention fractions, and primordial binary fractions on the properties of GCs evolving in a tidal field. We ran N-body simulations using the GPU-enabled version of NBODY6 with 32k, 64k, and 128k stars. We studied the effect of the different initial conditions on the cluster lifetime, remnant fraction, mass function, and structural parameters. In addition, we compared the results of the simulations with observational data from the literature and found good agreement. Owing to the specific shape of the King profile, the concentration parameter c is a poor representation of the cluster’s internal properties. Especially after core collapse, a King model is not able to reproduce the central cusp in any of our models and the concentration is systematically underestimated. For that reason we also computed the ratio of the radius containing 10% and 50% of the stars in the cluster r_10% and r_50%, which is a more accurate quantity than the parametric King fit.

When comparing the different models, we found that a cluster with central IMBH has a shorter lifetime than the rest of the models. This is caused by the enhanced ejection of stars due to lower mass segregation. On the other hand, models with high black-hole retention fractions live the longest of all simulations. A possible explanation for the extended lifetimes is that large numbers of stellar-mass black holes would expand the core and prevent core-collapse, thereby slowing down dynamical evolution. We stress that the difference in cluster lifetime is about 20% and the effects of IMBHs, stellar-mass black holes, and binaries are therefore modest.

Our study also showed that a central IMBH in a GC causes the remnant fraction to increase more slowly than in a cluster without a central black hole. However, a similar effect is reached with a high (30%) fraction of primordial binaries. Again, the reason for this behavior might be found in the lower degree of mass-segregation that causes high-mass stars and remnants to be distributed in the outskirts of the cluster where they can be easily ejected. For the evolution of the mass function in the simulations, we found that the central IMBH as well as a primordial binary fraction of 30% lead to a delay of the depletion of low-mass stars because high-mass stars are ejected. Therefore, the mass function slope has a lower value (i.e. is less depleted in low-mass stars) at the end of the cluster’s lifetime than for models without an IMBH or binaries.

When comparing the models with the observations in the plane of the mass function slope versus the cluster’s concentration, we found that the IMBH models clearly explain the observed data with low concentrations. By contrast, the models without an IMBH immediately undergo core collapse. Our more extensive parameter survey thus confirms the initial suggestion by Trenti et al. (2010). The analysis also showed that the same effect can be achieved when considering a black-hole retention fraction of 100%. In this case the black holes in the cluster center survive long enough to prevent core collapse over a long timescale. However, it seems very unrealistic to retain all black holes in the cluster after formation. We therefore conclude that GCs that lie in the upper left side of Fig. 6, that is, those with low concentrations but depleted mass functions, are good candidates for hosting (or having hosted) an IMBH at their center.

To find the best candidates for clusters hosting an IMBH we selected those that agreed within 3σ (derived from the spread in the values of α and c within 500 Myr) with one of the IMBH runs. Figure 6 shows the comparison with the models. The color-shaded areas mark the 3σ limits of the simulations. Only a few clusters can be identified as clear IMBH or non-IMBH candidates. The majority of the clusters lie in the white areas and are therefore not clearly defined. For the α − c plot, the best IMBH candidates are NGC 288, NGC 5466, and NGC 6656 because they directly overlap with the 3σ regions of the IMBH simulations. But NGC 5053, NGC 6809, NGC 6496, and Pal 5 are also promising because their values are closer to those predicted from the IMBH models than the non-IMBH models. Furthermore, NGC 6712 could be considered if the lifetime of the IMBH models were extended to reach higher values for α. NGC 5272, NGC 6341, NGC 5927, NGC 6093, and NGC 3201 agree with the non-IMBH model in their values of α and c. Analogous to the IMBH candidates, the high-concentration models to the right of the no-IMBH models, e.g., NGC 6397 or NGC 7078, can also be considered as clusters that most likely do not host an IMBH at their center. NGC 5286 agrees with all three models, which makes it challenging to compare with previous studies. Feldmeier et al. (2013) found the inner kinematics of NGC 5286 to be consistent with a 1500  M_⊙ black hole at its center. However, as mentioned earlier, the α − c plane for GCs might not be the best parameter to compare with our simulations because the King model is a poor representation of an evolved cluster’s structural parameters.

Comparing the non-parametric radii r_10% and r_50% in the same way for we find that NGC 1261, NGC 3201, NGC 5286, NGC 5272, NGC 6254, and NGC 6656 overlap with the 3σ regions of the IMBH models. The IMBH candidates NGC 5286 and NGC 5139 (ω Centauri) both are more consistent with models that host an IMBH. It is conceivable that clusters with different initial conditions undergo more extreme excursions to low concentrations, and we therefore also considered the clusters with high r_10%/r_50% to the left of the IMBH models as possible candidates. We note that this is almost the entire sample for which surface brightness profiles were available. However, as shown in panel d), some of the clusters are equally well reproduced by a model that contains a high black-hole retention fraction. This degeneracy may be lifted with more detailed N-body simulations for both clusters, so that one can distinguish between the two possible explanations. As mentioned before, a black-hole retention fraction of 100%, that is, when none of the formed stellar-mass black hole are ejected, does not agree with our knowledge of the violent process of black-hole formation. At least some fraction of the black holes will be ejected because of high kick velocities, and the model with the black-hole retention fraction of 100% should therefore be treated as an extreme case.

To test the results of this study, more kinematic observations have to be performed for the clusters that are considered IMBH candidates. A new sample of inner kinematics for Galactic GCs is desired and planned for the future. A larger sample of GCs with available mass function slopes and surface brightness profiles would also help to better constrain the behavior of GCs in the α − c and α − r_10%/r_50% plane. As we have shown, such observations may provide a powerful diagnostic for identifying GCs that possibly host IMBHs. Furthermore, future N-body simulations with an extended parameter space are needed to reproduce the observations that were not covered by the simulations presented in this study.

Acknowledgments

This research was supported by the DFG cluster of excellence Origin and Structure of the Universe (www.universe-cluster.de). We thank the anonymous referee for constructive comments that helped to improve this manuscript.

References

All Tables

All Figures

Fig. 1 Evolution of the reference models. Shown are a) the bound mass fraction, b) the slope of the mass function α (for 0.3  M⊙ < m < 0.8  M⊙), c) the Lagrangian radii r10% and d) r50%, e) the remnant mass fraction, and f) the number of stellar-mass black holes as a function of time. We stress that the absolute time on the x-axis is arbitrary because it depends on the total mass and the galactocentric distance and is therefore only shown for relative comparison. Dashed lines in panel f) mark the logarithmic fits to the evolution of the black hole numbers.

In the text

	Fig. 2 Lifetimes of all clusters as a function of their initial particle number.
In the text

	Fig. 3 Remnant mass fractions as a function of the completed lifetime fraction for simulations with different a) IMBH masses, b) black-hole retention fractions, and c) primordial binary fractions. All models contain 128k particles, except for the binary-fraction models, which are limited to 64k stars.
In the text

Fig. 4 Mass-function slopes (for 0.3  M⊙ < m < 0.8  M⊙) as a function of the completed lifetime fraction for simulations with different a) IMBH masses, b) black-hole retention fractions, and c) primordial binary fractions. All models contain 128k particles, except for the binary-fraction models, which are limited to 64k stars. The dashed lines indicate the true mass function slope, while the solid lines indicate the mass function slopes as measured by an observer, that is, counting binaries as single stars.

In the text

Fig. 5 Characteristic radii ratios (r10%/r50%, i.e., deprojected radii containing 10% and 50% of the stars) as a function of the completed lifetime fraction for simulations with different a) IMBH masses, b) black-hole retention fractions, and c) primordial binary fractions. All models contain 128k particles, except for the binary-fraction models, which are limited to 64k stars. The dashed lines indicate the rL,10%/rL,50% (deprojected radii containing 10% and 50% of the mass) evolution when taking the Lagrangian radii.

In the text

Fig. 6 α – c plane (left panels) and α − r10%/r50% plane (right panels) of observed GCs compared with N-body simulations. Shaded areas mark the 3σ limit of the simulations. Clusters falling within these limits and not overlapping with other models are considered to be well-represented by these models. The error bar in panel a) is representative for the uncertainties in α for the observed clusters. The labeled clusters are discussed in Sect. 5.

In the text