TIDAL DISRUPTION, GLOBAL MASS FUNCTION, AND STRUCTURAL PARAMETER EVOLUTION IN STAR CLUSTERS

Numerical investigations of the dynamics of globular clusters are computationally very challenging because of the need to resolve a large dynamic range in timescales, from a few hours typical of a tight binary to the billion of years over which the system evolves due to two-body relaxation. In addition, the computational complexity of a run carried out for a constant number of relaxation times scales as N³/log(N), because the relaxation time t_rh increases with the number of particles N (e.g., see Spitzer 1987). The computational resources required are further increased by the presence of a significant number of binaries, which can be up to 50% of the mass of the core (e.g., see Rubenstein & Bailyn 1997; Albrow et al. 2001; Pulone et al. 2003), although some clusters, like NGC 6397, might also be almost binary-free (Davis et al. 2008). Hard binaries have an important dynamic effect on the evolution of the core size (McMillan et al. 1990, 1991; Vesperini & Chernoff 1994; Heggie et al. 2006). In the past, numerical simulations have been performed either using approximate algorithms such as Monte Carlo methods (e.g., see Gao et al. 1991; Fregeau & Rasio 2007) or direct N-body simulations with a modest number of particles (N ≈ 10³ in McMillan et al. 1990 and Heggie & Aarseth 1992) until a recent improvement of one order of magnitude (N = 16384—Heggie et al. 2006; Trenti et al. 2007a, 2007b, 2008; Hurley 2007). These investigations highlighted two key results: (1) the presence of a significant population of primordial binaries drives the evolution of the core radius toward a value that is 5%–10% of the half-mass radius; and (2) the ratio of the core to half-mass radius (r_c/r_h) does not depend either on its initial value or on the primordial binary fraction f, once a few relaxation times have passed and provided that f ≳ 0.1. In addition, runs with a galactic tidal field show the central concentration parameter c decreases toward the end of the simulation, due to the combined effect of binaries, that keep a large core, and of the evaporation, that progressively reduces the tidal radius (Trenti et al. 2007b). However, the simulations of Trenti et al. (2007b) were obtained using simplified initial conditions with equal mass particles and therefore cannot be directly applied to the study of the mass function evolution or analyzed consistently with observations. In this paper, we extend our earlier work to include a realistic mass spectrum, as discussed below.

The simulation framework that we use here has been extensively described in Trenti et al. (2007b); Gill et al. (2008); and Pasquato et al. (2009). In short, we follow the evolution of star clusters using the direct summation code NBODY6 (Aarseth 2003), which guarantees an exact treatment of the multiple interactions between stars by employing special regularization techniques, without resorting to the introduction of softening.

Our N-body simulations are carried out in natural (dimensionless) units (Heggie & Mathieu 1986), in which:

$$\begin{equation} G=M_T=-4E_T=1, \end{equation} \tag{ 1 }$$

where M_T is the total mass of the system and E_T the total energy (potential plus kinetic). The corresponding unit of time is

$$\begin{equation} t_d = \frac{(GM_T)^{5/2}}{(4E_T)^{3/2}} = 1, \end{equation} \tag{ 2 }$$

which approximately corresponds to the orbital period at the half-mass radius. In these units, the relaxation time can be written as

$$\begin{equation} t_{rh} = \frac{0.138 N r_h^{3/2}}{\log (0.11 N)}, \end{equation} \tag{ 3 }$$

where r_h ≈ 1 is the radius at which the relaxation time is defined, that is the half-mass radius. Equation (3) can be translated in physical units as follows (Djorgovski 1993):

$$\begin{equation} t_{rh} = \frac{8.9\times 10^5 \mathrm{yr}}{\log {(0.11N)}} \times \frac{1 M_{\odot }}{\langle m_* \rangle }\times \left(\frac{M_T}{M_{\odot }}\right)^{1/2} \times \left(\frac{r_h}{1 }\right)^{3/2}, \end{equation} \tag{ 4 }$$

where 〈m_*〉 is the average mass of a star in the system (including dark remnants such as neutron stars).

The individual particle mass is drawn from an initial mass function appropriate to study the late evolutionary stages of star clusters, when stars are ≳10 Gyr old. As in Gill et al. (2008), we consider either a Salpeter (1955) or a Miller & Scalo (1979) initial mass function, that is:

$$\begin{equation} \xi (m) \propto m^{\alpha }, \end{equation} \tag{ 5 }$$

with α = −2.35 and m ∈ [0.2:100]M_☉ for the Salpeter initial mass function (IMF), while for the Miller & Scalo IMF the power-law slope is the following: α = −1.25 for m ∈ [0.2:1]M_☉, α = −2.0 for m ∈ [1:2]M_☉, α = −2.3 for m ∈ [2:10]M_☉, and α = −3.3 for m ∈ [10:100]M_☉. Before starting the run, an instantaneous step of stellar evolution is taken to evolve the mass function to a 0.8 M_☉ turnoff using the Hurley et al. (2000) evolutionary tracks. In our standard models, we have a 100% retention fraction of dark remnants. We also run two simulations with a 30% retention fraction of neutron stars and stellar mass black holes (see Table 1). In all simulations, no kick velocities are given to the remnants. Our idealized treatment of stellar evolution is based on the approximation that most of the relevant stellar evolution occurs on a timescale shorter than a relaxation time. It is thus appropriate to study the evolution of old globular clusters (age of ∼10 Gyr). In fact, most of the impact of stellar evolution on the dynamics of a star cluster happens within the first few hundred Myr, when the most massive stars lose a significant fraction of mass and consequently contribute to a global expansion of the system (e.g., see Hurley 2007; Mackey et al. 2008). Later in the life of a star cluster, two-body relaxation tends to erase the memory of the initial density profile and concentration (Trenti et al. 2007b, 2008).

The initial distribution in the position–velocity phase space is that of a King model with scaled central potential W₀ = 3, 5, 7. Our standard models have a primordial binary ratio of f = N_b/(N_s + N_b) between 0 and 0.1, with N_s and N_b being the number of singles and binaries, respectively. In addition, we also consider a run with a central IMBH (discussed in more detail in Gill et al. 2008). Following Trenti et al. (2007b), our initial distribution of the binaries' binding energies is flat in log scale in the range from _min to 133_min, with _min = 〈m(0)〉σ_c(0)². Here, σ_c(0) is the initial central velocity dispersion, and 〈m(0)〉 is the average stellar mass at t = 0. For a typical velocity dispersion σ_c(0) = 10 km s⁻¹, the semi-axes of binaries in the initial conditions are smaller than 10 AU. This choice is motivated by the fact that wider and therefore softer binaries, which are likely created in young star clusters, are destroyed within a few dynamical times and thus would have no effect on the late-time evolution of a star cluster (Heggie 1975). The stellar mass and type of binary members are drawn at random from our initial mass function and thus binaries containing either one or two dark remnants are possible.

The models considered in this paper are tidally limited. The Galactic tidal field is treated as that of a point mass, and the tidal force acting on each particle is computed using a linear approximation of the field. The tidal cutoff radius of the simulation (which we indicate as {r_t}_s, where the curl-bracket with the subscript "s" is used to denote the definition of a simulation based quantity; see Section 3) is defined as (Spitzer 1987):

$$\begin{equation} \left\lbrace r_t\right\rbrace _s^3 = \frac{M_T}{3M_{{\rm Gal}}}R_{{\rm Gal}}^3, \end{equation} \tag{ 6 }$$

where M_Gal is the galaxy mass, and R_Gal is the distance of the center of the star cluster from the center of the galaxy. We assume R_Gal is constant, i.e., the cluster describes a circular orbit. The initial value of {r_t}_s is fixed to be self-consistent with the King profile used (for example, the W₀ = 7 king model has {r_t}_s = 7.02 in N-body units). Particles are removed from the system when they reach a distance 2{r_t}_s from the center of the cluster (see Trenti et al. 2007b for a full description of the tidal field treatment).

We explore a large parameter space in the initial conditions varying the initial King model concentration, the Roche lobe filling factor, the initial mass function, and the primordial binary fraction, as summarized in Table 1. The NBODY6 code has been run with its Graphic Processing Unit (GPU) extension on the NCSA Lincoln Cluster. Each simulation was assigned to a computing node with eight cores and two Tesla C1060 GPUs. We measured a sustained computational performance above 0.5 Tflops for a single node, which allowed us to simulate systems with up to N = 65, 536 (or N = 64K in a compact notation where 1K ≡ 1024) particles until either t = 8000 was reached (corresponding to several initial relaxation times—t>16t_rh(0)) or at least ∼80% of the initial mass is lost due to evaporation of stars.

The positions and velocities of all particles in our simulations are saved every 10 code time units (one time unit is about one dynamic time; see Heggie & Mathieu 1986). We then proceed to construct a synthetic observation by first ignoring all particles except main-sequence stars. As our simulations only have an instantaneous step of stellar evolution, we do not consider stars along the giant branch. This choice is appropriate for comparison with star counts maps at high angular resolution, such as HST observations, because giants are typically masked out of the analysis (e.g., see De Marchi et al. 2007). Focusing on main-sequence stars allows us to significantly reduce the surface brightness profile fluctuations that arise from the small number of luminous giants that would otherwise dominate the light profile.

We then select a random direction and project the main-sequence stars positions creating a two-dimensional map. Within this map, we search for the center of the system using the Casertano & Hut (1985) density-center method.

A circular surface brightness profile is then constructed using ${\rm Int}[\sqrt{N_{MS}}]$ annuli (where N_MS is the number of main-sequence stars) each containing ${\sim} \sqrt{N_{{\rm MS}}}$ sources. From this profile, we define the half-light radius in projection r_hP as the radius containing half of the total luminosity. The King concentration parameter c = log₁₀(r_t/r_c) is determined by fitting the surface brightness profile with a single-mass King (1966) model. We use a precomputed King profile table with uniform spacing in W₀ (ΔW₀ = 0.1). The fit is carried out imposing that the total luminosity and half-light radius of the fitting model match those of the simulation snapshot within a ⩽1% deviation. From W₀ and r_hP, we derive r_t and r_c. As an alternative definition for the observed core radius, we consider the radius at which the surface brightness profile reaches half of its central value and we indicate this quantity with r_μ.

To construct the global mass function of the system, we proceed to identify close pairs of stars whose light would be blended together, similarly to the procedure described in Gill et al. (2008) and Pasquato et al. (2009). Our main aim is to treat binaries in the simulations consistently with observations. We also capture line-of-sight superpositions and therefore crowding. For reference, we set an HST-like resolution of 005 and we consider the system, of assumed half-light radius 3 pc, at a distance of 5 Kpc from the Sun. This implies that we consider blended those stars separated by less than 250 AU, which includes all our primordial hard binaries. In addition, we get a small number of apparent binaries, typically up to a few tens per snapshot depending on the central concentration of the system. Given a set of blended stars, we only consider the brightest member of the ensemble as observable for the purpose of re-constructing the stellar mass function. As the luminosity of a main-sequence star scales approximately with M^7/2, this is an adequate approximation. With this procedure we define a catalog of observable main-sequence stars where we include stars down to 0.2 M_☉ (for reference, the observed mass function of NCG 2298 has greater than 50% completeness at the center at this mass limit; see Pasquato et al. 2009).

The stellar mass function in main-sequence stars is modeled as a power law. Its index α is obtained by maximizing the likelihood L:

$$\begin{equation} \log {L({\alpha })} = -N_{{\rm MSR}} \log {\left(\frac{m_u^{{\alpha }+1}-m_d^{{\alpha }+1}}{{\alpha }+1} \right)}+{\alpha } \sum _{i=1,N_{{\rm MSR}}} \log (m_i), \end{equation} \tag{ 7 }$$

where N_MSR is the number of resolved main-sequence stars and m_d and m_U are their minimum and maximum mass (here, 0.2 and 0.8 M_☉).

Based on the catalog of resolved main-sequence stars used above to retrieve α, we define the total observed main-sequence mass of the system M as

$$\begin{equation} M = \sum _{i=1,N_{{\rm MSR}}} m_i. \end{equation} \tag{ 8 }$$

In a numerical simulation code such as NBODY6, the tidal radius {r_t}_s is defined based on Equation (6). We define as {M}_s ≡ M_T is the total mass of the cluster (including dark remnants such as black holes, neutron stars and white dwarfs), that is the mass of all the particles within 2{r_t}_s from the center of the cluster.

Several different definitions have been used in numerical simulations for the core radius, so we direct the reader to Trenti et al. (2007b) for a comprehensive discussion. Here, we adopt the standard definition from Casertano & Hut (1985):

$$\begin{equation} \left\lbrace r_c\right\rbrace _s= \frac{\sum _{i=1,N_{{\rm tot}}}{m_i r_i \rho _i}}{\sum _{i=1,N_{{\rm tot}}}{m_i \rho _i}}, \end{equation} \tag{ 9 }$$

where r_i is the distance from the center of the system and the sum is made over all the N_tot particles of the system (including dark remnants). The density ρ_i around each particle is computed from the distance to the fifth nearest neighbor (Casertano & Hut 1985). In passing, we note that this is not the NBODY6 definition of {r_c}_s, which is based on a ρ²_i weight (see Trenti et al. 2007b). The half-mass radius of a simulated system ({r_h}_s) is defined as the three-dimensional radius that contains half of the total mass of the system ({M}_s), thus including again any compact remnant that does not contribute the surface brightness profile of the system.

The evolution of the mass function index α is shown in Figure 1 as a function of the observed remaining mass of the system. As time progresses (and thus the system loses stars and mass), the mass function is preferentially depleted of low-mass stars and becomes flatter (larger α). Models in a stronger external tidal field lose star faster (Trenti et al. 2007b) and consequently have larger variation in α. However, when looked at in terms of the remaining fractional mass, the evolution of α becomes independent of the tidal field strength and almost independent of the binary fraction considered in the runs. A larger binary fraction tends to moderately slow down the flattening of the mass function because low-mass stars can be retained, and observed, if they have a dark remnant as a companion. All the simulations in Figure 1 start with the same IMF below the turnoff (but with different retention fractions of neutron stars and black holes in the right panel). In addition, their initial conditions are such that the tidal field cutoff is self-consistent with the initial density profile of the system. To investigate if a universal evolution of the mass function index is present in a broader range of conditions, we consider in Figure 2 additional runs that either start from a different IMF (Salpeter) or that are initially underfilling their Roche lobe. For the latter, we start from a W₀ = 3 King model but we set {r_t}_s = 6.28, which is twice as large as the self-consistent {r_t}_s(W₀ = 3) = 3.14. To compare the evolution of different mass functions, we consider the change in the power-law fit of the mass function: Δα = α(t) − α(t = 0). From Figure 2, it is immediately clear that the two simulations starting with a Salpeter IMF behave very similarly to the Miller & Scalo runs. The run initially starting with an underfilled Roche lobe does instead evolve marginally faster. Δα is in fact larger at a given fraction of mass loss, especially in the early stages of the evolution of the system. At later times, the Roche lobe is eventually filled, consistent with the results of the simulations by Gieles & Baumgardt (2008). At that point, the evolution of Δα is similar to our standard models. The close correlation between the slope of the stellar mass function and the fraction of the initial mass lost by a cluster was also found in previous studies (Vesperini & Heggie 1997; Baumgardt & Makino 2003). Our simulations confirm this relation and show its validity for a broader range of different initial conditions. Recently, Kruijssen (2009) investigated the evolution of the mass function with an analytic model, that reproduces the universality of Δα for different IMF. The Kruijssen (2009) conclusion that the retention fraction of black holes influences the evolution of the mass function (see Figure 15 and Section 5.2 in that paper) is however not seen in our simulations. This is likely because the analytic model of Kruijssen (2009) does not take into account

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

multiple encounters (three body or higher) of stellar black holes that contribute to their mutual ejection on a relaxation timescale (Merritt et al. 2004). The models of Kruijssen (2009) with different retention fractions start to show differences in the evolution of α when more than half of the mass of the system is lost. At that point, only a few stellar mass BH would have survived in the system, independently of their initial number (e.g., see Merritt et al. 2004; Gill et al. 2008).

The quasi-universal evolution of α as a function of M/M(t = 0) is very interesting because, if the initial mass function of Galactic globular clusters is universal, as for example suggested by Kroupa (2002), then the slope of mass function is a direct tracer of the mass loss of the system, and thus of its dynamical age and of the tidal field it experienced during its life. This diagnostic is most effective to detect a large mass loss as |dα/dM| steepens when M decreases. Our conclusions on the quasi-universality of the α evolution can be applied to old globular clusters, with a turnoff mass below 1 M_☉. For younger systems, both the turnoff mass and the tidal disruption timescale may introduce deviations from the universal behavior observed here (Kruijssen 2009). Finally, it is interesting to note that our conclusions are consistent with Galactic globular clusters observations of α (De Marchi et al. 2007): the clusters more massive than ∼2 × 10⁵ M_☉ have essentially the same α, while lower mass clusters usually show a broader range of α values, likely indicating that mass loss has taken place.

The evolution of the observed structural parameters also tends to attain a universal configuration, as shown in Figures 3 and 4 for the core-to-half-light radius ratio and in Figure 5 for the concentration c. The left panel of Figure 3 shows the r_c/r_hP evolution of N = 64K models with different tidal field configurations and initial density profiles as well as a binary fraction from 0% to 5%. Eventually all these models reach a common value for the observed r_c/r_hP ratio as derived from King profile fitting. However, the observed r_c/r_hP ratio behaves differently from the theoretical and numerical expectation for {r_c/r_h}_s. In fact, core collapse for systems made only of single stars is detectable only for a very short time based on r_c/r_hP. After core collapse, runs with 10% primordial binaries, (left panel of Figure 4) behave essentially as those with single stars only. We interpret this result as a consequence of mass segregation of dark remnants at the center of the cluster. The dark remnants are on average heavier than the main-sequence stars and thus they transfer kinetic energy to them. This results in a heating of the visible component of the cluster, which therefore does not exhibit a typical core-collapsed structure. An independent confirmation of this trend can be inferred from Figure 5 in Hurley (2007), where runs with and without primordial binaries also have a similar core (although Hurley 2007 runs do not reach past core collapse). The core collapse for single-star systems is instead well discernible from the evolution of {r_c/r_h}_s (see Figure 6). Until core collapse, theoretical and observed King-fitted core radii evolve similarly, but immediately after the core bounce the two definitions separate sharply from each other. Post-core-collapsed systems develop a shallow central cusp (μ ∼ R^−ν with ν ∼ 0.4 − 0.7) that is missed by the best-fitting King model (see Figure 7), whose structure is constrained by the global surface brightness profile. In general, a King model becomes a poorer description of the system after core collapse: in fact, the quality of our fits is degraded by about a factor 2, as measured by the ratio of typical χ² before and after core collapse. As a consequence, the fit results in the post-core-collapsed phase also become somewhat dependent on the modeling of photometric errors. We assumed a constant fractional error on the surface brightness profile, which corresponds to a constant error in magnitude, except for the outer points, where we assigned a constant error in flux (to model the effect of the sky background). Assuming a constant error in flux for all the points in the surface brightness profile yields to no differences in the fit before core collapse. After the collapse, slightly more concentrated models are preferred, because in this case the points with the highest surface brightness have an increased weight relative to the outer points.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results in Figure 3 have been obtained assuming a Miller & Scalo (1979) IMF and a 100% retention of neutron stars and stellar mass black holes. A Salpeter (1955) IMF has a higher fraction of dark remnants per unit mass (because of it shallower slope, α = −2.3, compared to the MS IMF, α = −3.3, in the range [10:100] M_☉), thus it sustains a larger core, especially in the earlier stages of the simulation (see the right panel of Figure 4). This effect is primarily due to stellar mass black holes (Merritt et al. 2004; Mackey et al. 2008), so r_c/r_hP returns in line with that of the Miller & Scalo (1979) simulations once the BHs segregate at the center of the cluster and kick each other out of the system via three-body interactions. Reducing the retention fraction of neutron stars and BHs to 30% leads to a decrease of the core radius down to r_c/r_h ∼ 0.05, albeit with large fluctuations. The large fluctuations that appear after the core collapse are also related to the degraded fit quality, which increases the uncertainty on c (that is on W₀).

Interestingly, the main-sequence surface brightness profile has a significantly larger core in the run with a central IMBH compared to all other runs once they reach their long-term equilibrium configuration (see the right panel of Figure 4). The prediction of Trenti et al. (2007a) and Gill et al. (2008) that collisionally relaxed systems with a central IMBH should exhibit a large r_c/r_hP is thus confirmed (see also Baumgardt et al. 2005). However, systems with a large core radius might also be in the pre-core-collapse stage. Unless the cluster relaxation time is short enough to suggest the cluster should have already undergone core collapse, no IMBH (or IMBH-like heating) is required to explain large values of r_c/r_hP. A small value of r_c/r_hP, on the other hand, can be safely used to exclude likely candidates to harbor an IMBH.

Fit of the surface brightness profile with a King model does not appear to provide a very solid determination of the intrinsic tidal radius of the system. The comparison between r_t and {r_t}_s in Figure 8 shows that, when the entire surface brightness profile is used for the fit, then r_t>{r_t}_s. This is however not surprising, because unbound particle escaping from the system are removed from the calculation only when they have r>2{r_t}_s, thus the surface brightness profile at radii beyond {r_t}_s may still be positive. Similarly in actual star clusters, unbound stars remains initially in the proximity of the system, while escaping. The measure of r_t becomes more uncertain, and dependent upon the details of the fit, if only a fraction of the surface brightness profile points is used (see blue dotted line in Figure 8). This is especially true after core collapse (that is at M/M(0) ≲ 0.65 in Figure 8). This is primarily because we carried out the fit over the full surface brightness profile by fixing the total luminosity and the half-light radius within 1% of the values in the profile. If only a fraction of the points is available, these two scales cannot be measured directly and must be left as free parameters of the fit. This leads to increased uncertainty, especially after core collapse, when the surface brightness profiles might present some departures from a King profile (for example, see Figure 9 at R ∼ 0.15). When we do not use the outer points in the fit, the tidal radius tends to be moderately underestimated after core collapse compared to the determination using the full profile (Figure 8). This is consistent with what has been found in the determination of the tidal radius for some Galactic globular clusters (McLaughlin & van der Marel 2005). A different conclusion may however be reached for those globular clusters that do not fill their Roche lobe (Baumgardt et al. 2009).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the difficulty in reliably identifying core-collapsed systems from King model fitting, we explore here alternate definitions for the observed core radius. One possibility is to use the classical core radius r_μ, that is the radius at which the surface brightness falls to half of its central value. A comparison between r_c and r_μ is shown in Figure 10 and suggests that this definition is better suited for identifying core-collapsed systems. r_μ remains in fact close to {r_c}_s (within ∼15% accuracy), confirming for complex dynamic configurations that include mass segregation the results published by Casertano & Hut (1985) for King profiles with a constant mass-to-light ratio. However, the definition of r_μ depends on the choice of bins, especially on the central one. Too small a bin size will lead to large Poisson fluctuations in the central surface brightness value, while an excessively large bin will overestimate the core size by lowering the central value of the surface brightness. A large central bin is not necessarily a choice, but may be imposed in actual observations by a limited angular resolution. To overcome some limitation related to binning, Casertano & Hut (1985) proposed to extend the density-based core radius estimator to the surface brightness profile μ(R) and to define a surface-brightness density radius s_μ. In our notation this reads:

$$\begin{equation} s_{\mu }= \frac{\sum _{i=1,N_{{\rm MSR}}} \mu (R_i) R_i }{\sum _{i=1,N_{{\rm MSR}}} \mu (R_i) }, \end{equation} \tag{ 10 }$$

where the sum is carried out over all the detected main-sequence stars. s_μ is shown for our simulations in Figure 10. As expected from Figure 3 in Casertano & Hut (1985), s_μ deviates from r_μ and {r_c}_s for high-concentration King models, that is mainly in the post-core-collapsed phase. Interestingly, the value of r_c obtained from King fitting for post-core-collapsed systems is approximately bound between r_μ and s_μ.

Fluctuations in the structural parameters beyond those expected from Poisson noise are clearly apparent from our analysis. Heggie & Giersz (2009) reach the same conclusion for their simulation of NGC 6397. Our run with a Miller & Scalo (1979) mass function and low retention of dark remnants has particularly strong fluctuations, as can be seen for all different core radius estimators in Figure 10. It is very interesting to note that the system appears to go through some breathing phases, where subtle changes are made to the surface brightness profiles, affecting primarily the King model fit. The origin of these fluctuations is likely dynamical, especially because we are considering only the main-sequence stars and thus we are relatively unaffected by low-number statistics.

Zoom In Zoom Out Reset image size

Our results on the determination of core radii have profound implications for our understanding of the dynamical state of observed globular clusters: most intrinsically core-collapsed globular clusters in our simulations appears to have a surface brightness profile in main-sequence stars essentially identical to non-core collapsed systems when a King model is used to retrieve r_c. Furthermore, our simulations do not predict very small r_c if the retention fraction of neutron stars and black holes is high. Alternative definitions of the core radius, such as r_μ, are more efficient at identifying core-collapsed clusters, but they may still suffer from large fluctuations and/or biases due to finite resolution at the center of the system. This picture clearly illustrates the potential problems in the observational identification of core collapsed clusters. In fact, the classification of an observed globular as core collapsed or not might be a result of the combination of significant fluctuations in the total surface brightness profiles (enhanced with respect to our main sequence only analysis because the profile is defined primarily by a small number of giant stars; see also Heggie & Giersz 2009), combined possibly with the limited angular resolution of ground-based data (Trager et al. 1995) used to construct the structural parameters presented in the Harris (1996) catalog. The core size in the post-core-collapse phase depends mildly (∝1/log(0.11N)) on the number of particles in the system (Vesperini & Chernoff 1994; Heggie et al. 2006), thus it is expected that smaller core radii can be measured in systems with a larger number of particles. Further exploration of this issue is deferred to a follow-up paper that will include live stellar evolution to follow stars off the main sequence and thus to correctly model the total luminosity of the cluster. This will allow us to carry out a more direct comparison with the ground-based data that define r_c in the Harris (1996) catalog.

An interesting correlation has been recently reported by De Marchi et al. (2007) between the slope α of the global stellar mass function of globular clusters and their central concentration c. Less concentrated clusters are more depleted in low-mass stars than those with higher values of c.

The evolution of a star cluster in the (c; α) plane is driven by collisional (two-body) processes, thus the De Marchi et al. (2007) observations test our understanding of the dynamics of globular clusters. As discussed in Section 1, the observed correlation appears to be in conflict with the standard scenario globular cluster evolution which predicts that relaxation drives an increase in c as the cluster evolves toward core collapse while at the same time preferentially depleting the system of low-mass stars via evaporation (Spitzer 1987). The evolution of our simulated star clusters in the concentration slope of the mass function (c–α) plane is shown in Figure 11, where we also report the values observed by De Marchi et al. (2007) for those clusters that are well relaxed (t_rh ⩽ 1 Gyr) according to the Harris (1996) catalog. In plotting the concentration parameter c for the clusters shown in Figure 11, we have used the values obtained by McLaughlin & van der Marel (2005), which are based on the King (1966) models. Differences from the concentration values reported in De Marchi et al. (2007) are however very small (Δc < 0.1, except for NGC 6712 which has c = 1.05 in McLaughlin & van der Marel 2005 and c = 0.9 in De Marchi et al. 2007). While our simulations never reach a deep core collapse (c ≳ 2.5), their concentration tends to evolve toward c ∼ 1.7, still larger than what observed for the systems with the most depleted stellar mass function. One possibility to explain the system with lower concentrations remains that of a very efficient source of heating at the center of the system, comparable with the heating generated by an IMBH (see cyan line in Figure 11). An IMBH itself appears to be excluded as the preferred explanation, because NGC 2298, where a central IMBH is rejected at high confidence level (Pasquato et al. 2009), is one of the low c points in Figure 11. Stellar mass black holes are also very unlikely to be viable, because the heating source needs to remain efficient in the latest stages of the evolution of the system. Stellar mass BHs lose instead efficiency at later times (see red line in the lower panels of Figure 6; see also Merritt et al. 2004). It has been recently suggested that white dwarf kicks can sustain a large core (Fregeau et al. 2009). However, the effect of the kick has been shown only at the level of {r_c/r_h}_s; and thus its full impact on the observed concentration c is an open question.

Zoom In Zoom Out Reset image size

An alternative explanation for the discrepancy in the (c; α) plane evolution might rest in observational biases. While we carried out our analysis restricting to main-sequence stars, the tidal radius of observed globular clusters is mainly derived from ground-based data (Trager et al. 1995), that thus are mostly dominated by giant-branch stars. Therefore, mass segregation might lead to a systematic underestimate of r_t from these observations compared to our analysis.

Overall, Figure 11 shows that the correlation between α and c observed by De Marchi et al. (2007) is still in mild disagreement with the expectations from numerical simulations. Half of the eight observed clusters are consistent with our models (including NGC 2298, which has c = 1.31; α = 0.5 and could be explained by a core fluctuation). Two of them appear instead to be core-collapsed clusters. They have a conventional concentration c = 2.5, higher than what we obtain from our simulations, as discussed in Section 4.2 (but they could be explained by assuming a low retention fraction of neutron stars and black holes). Two clusters have instead concentrations below our predictions, but here the consistency of simulations and observations might improve if data and models are analyzed with a fully self-consistent approach.