Understanding the Dynamical State of Globular Clusters: Core-Collapsed vs Non Core-Collapsed
We study the dynamical evolution of globular clusters using our Hénon-type Monte Carlo code for stellar dynamics including all relevant physics such as two-body relaxation, single and binary stellar evolution, Galactic tidal stripping, and strong interactions such as physical collisions and binary mediated scattering. We compute a large database of several hundred models starting from broad ranges of initial conditions guided by observations of young and massive star clusters. We show that these initial conditions very naturally lead to present day clusters with properties including the central density, core radius, half-light radius, half-mass relaxation time, and cluster mass, that match well with those of the old Galactic globular clusters. In particular, we can naturally reproduce the bimodal distribution in observed core radii separating the “core-collapsed” vs the “non core-collapsed” clusters. We see that the core-collapsed clusters are those that have reached or are about to reach the equilibrium “binary burning” phase. The non core-collapsed clusters are still undergoing gravo-thermal contraction.
keywords:Galaxy: kinematics and dynamics – Galaxies: star clusters: general – globular clusters: general – Methods: numerical
Studying the evolution of dense star clusters, such as old globular clusters (GCs), is of great interest for a variety of branches of astronomy and astrophysics. The high central densities and high masses of GCs make them hotbeds for strong dynamical interactions facilitating formation of many exotic sources (e.g., X-ray binaries, millisecond radio pulsars, type Ia supernovae, and blue straggler stars). GCs are important targets for extragalactic Astronomy. Detailed observations of their spatial distribution in a galaxy can constrain, for example, the potential of the dark matter halo, and give clues to the assembly history of the galaxy. The old ages of GCs provide a direct window into the major star-formation episodes in the early universe.
One long-standing question regarding GCs concerns the nature of their progenitors. The observed young massive star clusters (e.g., Holtzman et al., 1992; Whitmore & Schweizer, 1995; Miller et al., 1997; Scheepmaker et al., 2007, 2009) seem to be potential progenitors of the current GCs. The masses, typical sizes, and inferred dissolution timescales for the so-called “super star clusters”, as observed, for example, in M51 (Scheepmaker et al., 2007), make them especially good candidates. Interestingly, similar to the Galactic GCs (GGC), the sizes of observed super star clusters are typically a few parsecs independent of the cluster mass (e.g., Scheepmaker et al., 2007, 2009). However, the link between these super star clusters and the old evolved GGCs as well as GCs in other galaxies remains speculative.
The main difficulty in establishing a link between the two populations is that almost a Hubble time of evolution separates them. Although the individual qualitative effects of each physical process in a GC has been known from decades of numerical studies (e.g., Heggie & Hut, 2003), it is impossible to estimate the collective effect of these interdependent processes unless a self-consistent simulation is done including them all. For example, two-body relaxation leads to a slow energy diffusion from the core to the halo in a GC leading to a slow and steady contraction of the core until the gravo-thermal instability occurs and the core collapses under gravity. Since, due to equipartition of energy, the low-mass stars are scattered to higher velocities relative to their high-mass counterparts, two-body relaxation also leads to mass segregation in the cluster. The high-mass stars sink to the core and the low-mass stars escape to the halo and preferentially get stripped through the tidal boundary of the cluster. Because of this preferential loss of low-mass stars, the stellar MF changes with time. The stellar MF, in turn, determines the fraction of low and high-mass stars in a cluster as well as the average stellar mass, affecting the mass-segregation timescale (Gürkan et al., 2004). The contraction of the core via two-body relaxation increases its density. The core density determines the interaction rate in the core. These rates affect the binary-single (BS) and binary-binary (BB) interaction probabilities at a given time. The BS and BB interactions in turn generate energy in the core and can support the core stopping further contraction, thus, affecting the central density. The BS and BB interactions also change the orbital properties of the binaries taking part in these scattering interactions. These changes in the binary orbits can alter the evolution pathways that would be taken by a given binary which in turn changes, for example, the rate of formation of exotic stellar populations such as X-ray binaries and blue straggles stars (BSS). Due to this complexity in the evolution of dense star clusters, the only way to learn more about the possible initial properties of the observed GCs is to perform numerical simulations including all of the above physical processes in tandem with reasonable accuracy for a large enough .
The study of GCs has a somewhat long and varied history during which numerical simulations and observations have complemented each other (Heggie & Hut, 2003). In particular, understanding the physical processes in the cores of GGCs has been of prime interest since the evolution of these dense clusters are driven mainly by the energy generation in the core and the transport of this energy from the core. The GGCs are observed to show a clear bimodal distribution of the core sizes that separates the so-called “core-collapsed” clusters from the non core-collapsed clusters (Harris, 1996; McLaughlin & van der Marel, 2005). The core-collapsed clusters in general show a power-law slope in their density profiles near the center. In contrast, the density profiles of the non core-collapsed clusters are well described by a King profile (King, 1966) and show a clear flat part near the center. Theoretical analysis based on the current estimated relaxation times for the GGCs indicate that the majority of the GGCs should have had a deep collapse (Djorgovski, 1993; Harris, 1996). Thus, before it was found that all GGCs contain dynamically significant numbers of binaries, theoretical studies focused on understanding the process of core-collapse via the balance of outward diffusion of energy from the core due to two-body relaxation and post-collapse evolution due to dynamical formation of binaries (e.g., Heggie & Hut, 2003). After it was observed in the early 1990s that all GGCs contain sufficient number of binaries such that they must have been born with substantial primordial binary populations, theoretical studies focused on properties of clusters in the “binary-burning” phase in which the outward diffusion of energy via two-body relaxation is balanced by production of energy via dynamical hardening of binaries (e.g., Vesperini & Chernoff, 1994; Fregeau & Rasio, 2007). It was also realized that even a small primordial binary fraction can support the core from deep collapse for more than a Hubble time (Fregeau & Rasio, 2007). However, comparison between theoretical predictions and observations show that the theoretically predicted core radii during the binary-burning phase for these clusters are at least an order of magnitude smaller than the observed core radii for the bulk of the GGCs (Vesperini & Chernoff, 1994). Additional energy sources were proposed to explain these large core sizes (e.g., Mackey et al., 2007, 2008; Chatterjee et al., 2008; Trenti et al., 2007), however, these scenarios need rather special conditions and seem unlikely to be satisfied by most of the GGCs. Very basic questions remain. Does the bimodal distribution for the GGC core sizes separating the core-collapsed and non core-collapsed clusters indicate a physical difference between these clusters? What dynamical stage are the core-collapsed and non core-collapsed clusters in today? Can the large core sizes of the non core-collapsed (majority of the Galactic population) be explained without any special conditions simply as a result of of evolution from realistic initial conditions?
In this study we present the results from a large number of computer simulations (224) with initial conditions drawn from a multidimensional grid, spanning all relevant parameter ranges as suggested by observations of young star clusters, with the goal of reproducing a population of old clusters similar to the GGCs. We use CMC, a Hénon-type Monte Carlo code including all physical processes such as single and binary stellar evolution, two-body relaxation, strong encounters comprising physical collisions and binary-mediated scattering, and tidal stripping due to the Galactic tidal field. CMC has been extensively tested and the results from CMC show excellent agreement with those from direct -body simulations (e.g., Fregeau & Rasio, 2007; Chatterjee et al., 2010; Umbreit et al., 2012). Our goal is to find whether starting from realistic initial conditions (including , stellar mass function, central density, compactness parameter, binary fraction, and cluster size) typical of the observed young massive star clusters (e.g., Scheepmaker et al., 2007, 2009) and without any special treatment clusters similar in properties to the observed old GGCs are naturally obtained after of evolution. In particular, we focus on understanding the bimodal distribution of the observed GGC core radii to identify the dynamical stages of the so called core-collapsed and non core-collapsed clusters.
In Section 2 we briefly explain our code and introduce working definitions for key structural parameters of a cluster. In Section 3 we describe the multidimensional grid of initial parameters explored in this study. In Section 4 we present our key results. In Section 5 we show our results identifying the dynamical evolutionary state for the observed core-collapsed GGCs. Finally we conclude in Section 6.
2 Numerical Method
We use our Hénon-type Cluster Monte Carlo code (CMC) to numerically model star clusters with single and binary stars including all physical processes relevant in globular clusters such as two-body relaxation, single and binary stellar evolution, strong encounters including physical collisions and binary mediated strong scattering encounters. This code has been developed and rigorously tested over the past decade (Joshi et al., 2000, 2001; Fregeau et al., 2003; Fregeau & Rasio, 2007; Chatterjee et al., 2010; Umbreit et al., 2012; Pattabiraman et al., 2012). Using a large grid of initial conditions over the range of values typical for observed young massive clusters (e.g., Scheepmaker et al., 2007, 2009) we create over detailed star-by-star models and evolve them for . A handful (about 6) of these models reach a very deep collapse phase. For these models the CMC time steps become minuscule and the code grinds to a halt. We stop our simulations at that point for these clusters. In reality, the deep-collapse phase is halted via formation of the so-called “three-body” binaries and the cluster enters into the gravo-thermal oscillation phase. Since in CMC we do not yet include the possibility of creating new binaries via three body encounters, we do not address this phase at this stage. This however, is not a serious limitation. All clusters that reach a very deep collapse stage before the integration is stopped at , had zero primordial binaries which is not realistic (e.g., Davis et al., 2008). Here, we only include these clusters in our analysis as limiting cases. Even a small non-zero primordial binary fraction (; lowest is used in our simulations) can stop the cluster core from collapsing through the dynamical hardening of binaries preventing very deep core collapse. At this stage the core size remains more or less constant, which is also commonly referred to as “binary-burning” stage. A more detailed description of the code, and the various qualitatively different dynamical stages for a cluster’s evolution is presented in Chatterjee et al. (2010).
Since one key goal for this study is to compare the properties of our simulated models at about with the properties of observed GGCs, we have to make sure the same parameter definitions are used. In particular there are different definitions for the core radius and the cluster size commonly used by theoreticians and observers. The three dimensional core radius widely used in -body simulations is a density-weighted measure related to the virial radius in the core (Casertano & Hut, 1985). In contrast, the core radius for observed clusters, , is often defined as the distance from the center where the projected surface brightness profile drops to half its central value. Alternatively, when observations with sufficiently high resolution are available, the core radius can also be based on star-counts representing the radius where the projected surface number density is half the corresponding central value. If the core radius is calculated using the stellar number density profile, we will call it , while the core radius based on the brightness, or luminosity, density profile is denoted by . In real clusters the density profiles can have large scatter, especially close to the centre due to Poisson noise and the presence of bright giants that dominate the light but are only few in number. Hence, it is common practice to exclude the light from the brightest giants when calculating the luminosity density profile for real clusters to reduce noise (e.g., Noyola & Gebhardt, 2006). For the purpose of this study we define as the core radius obtained from a luminosity density profile excluding giants with . Since the estimated value of is strongly dependent on the estimate of the peak luminosity density at the centre of a cluster, values of and can differ significantly. For real clusters even when star counts are available, only stars above a certain brightness can be counted due to completeness issues. Hence, we define calculated using a number density profile including only stars that are on the main-sequence or on the giant branch with masses , representing the currently achievable completeness limit (e.g., Leigh et al., 2012). We find that the values and the values are not much different for our simulated clusters.
The three dimensional half-mass radius is defined as the radius that includes half of the total mass of the cluster. However, for real clusters this quantity is not directly accessible, rather, the two dimensional, projected distance enclosing half of the total light of the cluster is calculated and is defined as the half-light radius . We further define for the half-light radius for our simulated models by calculating the same quantity but excluding bright giants with . We find that the and values show only minor differences.
3 Initial Conditions
Our choice of initial conditions is based on observations of super star clusters. All our simulated clusters have initial virial radii between – independent of other cluster parameters. This range corresponds to initial three-dimensional half mass radii, , ranging from – . These values are in agreement with observations of young, massive clusters that indicate that the effective cluster sizes are rather insensitive to the cluster mass, as well as metallicity (e.g. Ashman & Zepf, 2001; Scheepmaker et al., 2007, 2009) and have a median value of . In addition, observations of old massive LMC clusters, old GCs in NGC 5128, old clusters in M51, as well as the GGCs indicate that the effective cluster radii show only a weak correlation with the distance from the galactic center (Hodge, 1962; Harris et al., 1984; Hesser et al., 1984; Mateo, 1987; van den Bergh et al., 1991; Scheepmaker et al., 2007; Hwang & Lee, 2008).
To restrict the huge parameter space to a certain extent we place all our simulated clusters in a circular orbit at a moderate Galactocentric distance of . We avoid modeling GCs very close to the Galactic center, where the Galactic field is so strong that the tidal stellar loss dominates the cluster’s evolution. Due to the assumption of spherical symmetry, Monte Carlo codes cannot directly model the tear-drop shaped tidal boundary of a star cluster. Instead these codes use some prescription based method (see a detailed discussion and calibration in Chatterjee et al., 2010). If tidal dissolution is the dominant driver of the GC’s evolution, the approximate method may not be accurate enough. Choosing a circular orbit for the simulated clusters is a simplification; however, the results should still be valid for eccentric orbits with some effective Galactocentric distance () (e.g., Baumgardt & Makino, 2003). The Galactic tidal field, and consequently the initial tidal radius, (), for the clusters is calculated following Equation 1 (Baumgardt & Makino, 2003),
using a Galactic rotation speed . Here the cluster mass is denoted by .
For our set of runs we vary the initial number of stars, , between , encompassing the bulk of the GGCs (McLaughlin & van der Marel, 2005). The initial positions and velocities are sampled from a King model distribution function with dimensionless potential, , in the range . We vary the initial binary fraction, , between . The stellar masses of the stars, or primaries in case of a binary, are chosen from the IMF presented in Kroupa (2001, their Equations and ) in the stellar mass range . Secondary binary companion masses are sampled from a uniform distribution of mass ratios in the range , where is the mass of the primary. The semi-major axes, , for stellar binaries are chosen from a distribution flat in log within physical limits, namely between the physical contact of the components and the local hard-soft boundary. Although initially all binaries in our models are hard at their respective positions, some of these hard binaries can become soft during the evolution of the cluster. The cluster contracts under two-body relaxation and the velocity dispersion increases making initially hard binaries soft. Moreover, binaries sink to the core due to mass segregation where the velocity dispersion is higher than that at the initial binary positions. We include these soft binaries in our simulations until they are naturally disrupted via strong encounters in the cluster.
4 Basic Structural Parameters of Our Model Clusters
Here we present the evolution of some structural properties of the simulated clusters and compare them with the same properties of the observed GGCs. For each of these comparison plots, the evolution of a certain cluster property is shown together with the distribution of the corresponding observationally derived values for the GGC population. Since all our simulated models are at a Galactocentric distance of , we also show the observed distributions for a subset of GGCs satisfying . Our goal here is to simply test whether, starting from observationally motivated initial conditions typical for massive, young clusters, the final properties of the cluster ensemble naturally attain ranges of values as observed in the GGCs. We do not intend to reproduce the present day distribution for these properties since this would require to introduce the distribution of cluster initial conditions as another parameter, and, consequently, significantly more simulations which is beyond the scope of this study.
Figure 1 shows the evolution of the total mass of our simulated clusters. The initial sharp decrease in the cluster mass () is because of the high mass loss rate via winds and compact object formation of the massive stars in the GC. Later on, decreases at a slower, nearly constant rate caused by a steady stellar mass loss through the tidal boundary of the GC (Baumgardt & Makino, 2003; Davis et al., 2010). The histograms on the right show the distribution of observed GGC masses. The observed GGC masses are estimated from the absolute visual magnitudes () given in Harris (1996) using Equation 2 assuming a uniform mass to light ratio for all clusters.
This is an approximation. The can vary from cluster to cluster and can also depend on the evolutionary stage of the cluster in question (Anders et al., 2009). However, since here we are only interested in reproducing the cluster mass ranges for most of the GGCs and do not try to model a particular cluster, this estimate should be appropriate. Figure 1 shows that our model clusters have final masses typical for the bulk of the observed GGCs. Since these clusters are modeled with the total and masses typical of the observed GGCs the model properties can be directly compared with the overall GGC population without any need for scaling.
Figure 2 shows the evolution of the central density in our simulated models. As a dynamically important GC property, the central density () determines the interaction cross-sections for strong scattering inside the core, for example, for BB and BS interactions, and stellar collisions. These strong interactions in turn modify the properties of the core, through, e.g., binary burning. The values sharply decrease during the first of evolution as the high-mass stars lose mass via stellar winds, and compact object formation. Followed by the sharp decrease increases almost linearly over time during the gravo-thermal core contraction stage. The histograms show the central densities of the observed GGCs. Here we convert the bolometric luminosity densities presented in Harris (1996) and assume . The range of central densities for our collection of simulated models is compatible with the range for observed GGCs.
Figure 3 shows the evolution of the density weighted three dimensional core radius . The for all clusters sharply increases initially for up to about , because of mass loss due to stellar evolution. This mass loss happens mainly at the deepest part of the cluster potential since the high-mass stars reside near the cluster center due to mass segregation and are affected by mass loss the most. The resulting loss of gravitational binding energy expands the cluster core. Once the rate of mass loss goes down, the core contracts via diffusive energy transport from the core to the outside through two-body relaxation. The gravo-thermal contraction ceases when this energy flow is balanced by the production of energy via dynamical hardening of binaries, the binary-burning phase. In our collection of models about 20 clusters show clear binary-burning end stages exhibited by a near constant core radius at late stages (). The bulk of our models are still contracting at . A few of our models with do not show the binary-burning stage and enter deep-collapse after gravo-thermal contraction. This is because we do not account for dynamical creation of binaries. However, as mentioned before, is a limiting case and not representative for observed GCs. Even with , we find that the clusters that complete the slow contraction phase do not suffer deep collapse, rather reach a steady binary-burning energy equilibrium phase.
Another reliably measured and frequently used parameter reflecting the dynamical state of the evolution of GCs is the ratio between the core radius and the half-mass radius (e.g., Fregeau & Rasio, 2007). Figure 4 shows the evolution of for all our clusters. The range of final values of the simulated clusters agree well with the observed ones in the GGC population, producing values at close to the peak of the observed distribution.
We should remind the readers, however, that these and values are not exactly the quantities observed directly. A more careful comparison between these structural parameters in our models and the observed values in the GGCs requires “observing” the models and defining the model quantities such as and as the observers define them for a real cluster (Section 2). For example, at our model run12 has the -body defined . If the same model is observed and the core radius is calculated using the surface number density distribution using all main-sequence stars satisfying and low-luminosity () giant stars, the core radius (Table 1). Figure 5 shows the distribution of values for the ratio using values obtained in 4 different ways. The values for the different definitions of and the resulting ratios follow different distributions. Thus, if one wants to estimate what the three dimensional -body defined value is from the observed value of and vice versa, one should be careful about how for that cluster was calculated. In general, calculated using the number density profile of the cluster show lower levels of errors. Moreover, the distributions do not change significantly based on which stars were included in the sample to calculate the surface number density profile. The code defined is typically between – times the for all cluster models in our collection. The luminosity density profiles are subject to a lot more noise compared to the number density profiles. As a result the calculated using the luminosity density profiles show a larger spread. In addition, which stars were included in calculating the luminosity density profile matters in this calculation significantly. Typically, the code-defined values are between – times and – times the and , respectively.
Similarly, only the half-light radius projected on the sky () is observed in reality and not the three dimensional half mass radius (). For example, for the simulated cluster run12 the sky projected half light radius including all stars is . If a luminosity cut-off is used for the same cluster the half light radius becomes . The theoretically calculated three dimensional half-mass radius for the same cluster at the same age is (Table 1). This difference is not simply due to the projection effect. This difference depends intricately on the positions of the bright giant stars that dominate the light but are few in number. Figure 6 shows the distribution of the ratio calculated using two separate samples of stars, i.e., with and without applying a luminosity cut-off. Typically, the three dimensional half mass radius is between – times the observed if the luminosity cut described in Section 2 is used. Note that this range includes the expected geometric factor of due to projection effect. The spread in values is moderately larger for . This is due to the increased statistical fluctuations created by the bright and low number of giants with .
We are not the first to point out this discrepancy in definitions. For instance, Hurley (2007) and Trenti et al. (2010) already found that the -body definition of can differ from an observed by a factor of a few, and our results confirm that. However, Figures 5 and 6 are expected to be useful for observers and theorists modeling GCs since the three dimensional or can be be easily estimated using the presented distributions if the observed values for these are known.
We now compare the values of the structural parameters for our simulated models calculated in a similar way as they are calculated for observed GGCs. Figure 7 shows the as a function of . The observed values for all GGCs and GGCs in the Solar neighborhood are also shown in the histograms along the respective axes. The two clusters of points show vs values for clusters with initial and , respectively. Our collection of models shows very similar half-light radii and central densities compared to those properties for the bulk of the GGCs. Note that the final values are strongly dependent on the initial for a given Galactocentric distance () as evidenced in the clustering of the points in Figure 7. Interestingly, for GCs at the half-light radii remain close to the initial according to our models. Thus, the present observed half-light radii can be used as an effective indicator for what the initial was for these clusters. The other parts of the observed histograms for can be easily populated if more simulations are performed using smaller initial . In addition, using smaller will also lead to smaller final values for these clusters as they cannot expand further once they fill their tidal radii.
Figure 8 compares the final values of and for all our simulated models with the observed core radius distribution for the GGCs. As expected, the smaller the , the higher the . Our simulated models populate a very similar range in values as is observed in the bulk of the GGCs. This is also true when other definitions of are used for the comparison.
Figure 9 shows the scatter plot for the and for our models together with the corresponding histograms for the observed GGCs. The values for this dynamically important and dimensionless measure for the compactness of the cluster show excellent agreement with the values typically found in all observed GGCs.
Figure 10 shows the distributions of relaxation times at () for our models, all GGCs, and a subsample of the GGCs in the solar neighborhood (). The relaxation time is a very important dynamical quantity because it is the time scale on which the global cluster properties evolve. However, it is difficult to derive this quantity accurately for observed clusters since it depends on several dynamical cluster properties that are not directly observable. Traditionally, multiple assumptions are made to calculate for observed GGCs (Harris, 1996, see the bibliography in the online database). To remain consistent with the derived values in Harris (1996) we adopt the same assumptions to calculate the values for our models. At the end of the simulation () we calculate the total luminosity () of the model cluster, and compute its total mass using (actually can have values between – in our models). The total number of stars is then calculated assuming an average stellar mass (actually can have values between – in our models). We use as a proxy for . We estimate the values for our models using Equation 11 of Djorgovski (1993) with the corrected coefficient as mentioned in Harris (1996). The values for our models calculated this way agree well for the GGCs in the Solar neighborhood. The bin showing an unusually low value just above contains a single cluster E3. A quick search indicates that E3 is a sparse unusual cluster. Again we remind the readers that here we are only interested in comparing the ranges of from our models and the observed GGCs. Note that the values for based on the actual cluster parameters can be a few times higher depending on the particular cluster properties. The differences between the actual dynamical values and those derived for the observed clusters come essentially from the many assumptions listed above.
A long standing puzzle has been the apparent discrepancy between the theoretically predicted values and the values observed for the GGCs. Early numerical simulations as well as analytical studies expected that most GGCs are in the binary-burning stage. However, the simulated values resulting from binary burning have been found to be about an order of magnitude smaller than that for the bulk of the observed population (e.g., Vesperini & Chernoff, 1994; Fregeau & Rasio, 2007) and it was already known that this amount of discrepancy cannot be entirely coming from differences in definitions (e.g., Hurley, 2007; Trenti et al., 2010). Consequently, additional energy sources in the core to expand have been investigated. Several studies proposed different additional energy generation mechanisms to explain the large observed values (e.g. Trenti et al., 2007; Chatterjee et al., 2008; Fregeau, 2008; Mackey et al., 2008). However, these additional energy sources require rather special conditions. For example, high central densities are required to create a population of high-mass stars via physical collisions that can then suffer expedited mass loss via compact object formation (Chatterjee et al., 2008). On the other hand, presence of an intermediate mass black hole also cannot be common for the GGCs (Trenti et al., 2007). Our models, in contrast, are generated without any special assumptions using observationally motivated initial conditions, and they naturally create a population of model clusters with properties in excellent agreement with the bulk of the observed GGC properties. These results indicate that the progenitors of today’s GGCs were very similar in properties to the young massive clusters observed, for example, in M 51 (Scheepmaker et al., 2007, 2009). Our results also indicate that the majority of today’s GGCs have not yet reached the binary-burning, energy-equilibrium stage, and are still contracting under energy transport via two-body relaxation.
5 What is a “core-collapsed” cluster?
A lot of early theoretical work was devoted to understand the gravo-thermal collapse and subsequent evolution of the core and the cluster as a whole due to hardening of primordial binaries. However, analytical results as well as numerical simulations showed that the values in the binary-burning phase are about an order of magnitude too small compared to the bulk of the GGC cores.
Our simulation results presented in Section 4 show that this is simply because the bulk of the observed GGCs are not in binary-burning equilibrium stage. Rather they are still contracting. Now we focus on understanding the clearly bimodal distribution of the core radii of the GGCs. Depending on the shape of the observed cluster density profiles, all GGCs are divided into two categories, namely core-collapsed and non core-collapsed clusters (e.g., Harris, 1996; McLaughlin & van der Marel, 2005). The so called core-collapsed clusters exhibit a power-law increase in the density profile until the limit of resolution of the observation. The non core-collapsed clusters show a very clear flattening of the density profile and the profile for these clusters are well fitted with a King profile (King, 1966). Is there a distinct difference dynamically between the two categories of clusters?
Figure 11 shows examples of two representative clusters chosen randomly from our large collection of models, one is in a clearly binary-burning stage and the other is still contracting. The evolution of is shown for the two clusters as well as their respective surface density profiles. The binary-burning cluster can be clearly identified by the near constant value of starting at about . The core contracting cluster shows a constant rate of contraction until . The surface density profiles for the two clusters are very different. The surface density profile for the binary-burning cluster looks very much like a so called core-collapsed cluster and exhibits a clearly power-law slope to very small radius below which the profile is noisy. In contrast, the surface density profile for the core contracting cluster shows a clear King profile with a distinct flat part near the center below a few parsecs.
We divide the full sample of our simulations in two subgroups: 1) binary-burning: models showing a clear binary-burning stage before integration is stopped, and 2) core-contracting: models still contracting due to two-body relaxation until integration stopping time. Figure 12 shows the distributions for for the two subgroups of our theoretical models. The distributions peak at different values, with the binary-burning clusters having much lower core radii. The distributions for for the observed core-collapsed and non core-collapsed GGCs are also shown. Note that the values for the binary-burning clusters from our models and the core-collapsed clusters from the observed GGCs show a very similar range of values. Similarly, for the contracting clusters in our models show very similar values compared to those of the observed non core-collapsed GGCs. In addition, most of the core-contracting models with small values are in fact about to start binary-burning, but was classified by contracting since the constant stage of evolution was not clearly seen. The observed GGCs show larger ranges for the core radii values compared for both categories of clusters compared to our models. This is simply because we are forced to limit the ranges of the grid of initial conditions to constrain the number of required cluster calculations to a tractable amount. We remind the readers that the heights of the histograms between the model clusters and the observed GGCs are not compared here, since that depends directly on the distributions of initial cluster parameters, determination of which is beyond the scope of this study, and are chosen arbitrarily. Only ranges in values are compared.
We have presented a large () collection of cluster models created with the Northwestern group’s Hénon-type Monte Carlo code in star-by-star detail. We start our models with initial parameters including the stellar mass spectrum, cluster size, concentration, and primordial binary fraction over large ranges (Section 3) guided by observed young massive clusters (e.g., Scheepmaker et al., 2007, 2009). We find that these initial clusters very naturally produce a population of model clusters with structural properties including cluster mass, , , , and in excellent agreement with the bulk of the GGC properties after about of evolution without any special considerations or fine tuning (e.g., very high density to aid collisional expedited stellar mass loss via compact object formation; Chatterjee et al. 2008; or intermediate mass black holes; Trenti et al. 2007). We pay attention to the various different commonly used definitions of the structural parameters and and calculate these quantities from our models as an observer would for real clusters. These parameters are then compared and found to agree well with the ranges from observed GGCs. Using our large collection of models we also show the distribution of the ratio of the three dimensional code-defined and to the corresponding “observed” values (Figures 5, 6). We expect that these distributions of ratios for the and values will be valuable for observers and theorists alike to convert the values of these parameters from one set of definitions to another.
From the evolution of the code-defined three-dimensional structural parameters of all our models, we find that all qualitatively different evolutionary stages are observed, in particular, the initial expansion due to stellar evolution driven mass loss, core contraction driven by two-body relaxation, and the binary-burning equilibrium stage (for clusters with ) driven by a balance between energy production via dynamical hardening of binaries in the core and outward diffusion of energy from the core due to two-body relaxation. Our results indicate that the progenitors of today’s GGCs were very similar in properties to the present day young massive clusters (observed, for example, in M 51 Scheepmaker et al., 2007, 2009). Of course, the metallicities of these progenitors must have been much lower compared to today’s massive young clusters.
After establishing that our collection of cluster models are representative of the observed GGCs we investigate the apparently bimodal distribution of the observed core radii of the GGCs created by the core collapsed clusters and non core-collapsed clusters. In particular, we answer the question if the core collapsed clusters are dynamically different from the non core collapsed clusters. We find that the surface brightness profile for the binary-burning cluster shows a prominent power-law slope near the center (Figure 11). The core collapsed clusters are observationally defined as the clusters that show this distinct feature in their surface brightness profiles (e.g., McLaughlin & van der Marel, 2005). In contrast, a model cluster that is still contracting at shows a surface brightness profile that has a clear flat central part and is fitted well by a King profile (King, 1966). We further divide our models into two subsets, one containing clusters in the binary-burning stage, and the other containing clusters that are still contracting at . We compare the ratio between the core and half light radii of our binary-burning and contracting clusters with those for the core-collapsed and non core-collapsed clusters in the observed GGCs, respectively. We find that the binary-burning values in our models are in agreement with those of the core-collapsed GGCs. Similarly, the contracting values in our models are in agreement with those of the non core-collapsed GGCs (Figure 12). Thus our results clearly indicate that the so called core-collapsed GGCs are in fact at the binary-burning stage whereas, the non core-collapsed GGCs are still contracting under two-body relaxation. This also indicates that the majority of the GGCs (since most GGCs are non core-collapsed) are not in energy equilibrium as was expected by some earlier theoretical models (e.g., Fregeau & Rasio, 2007). One key implication for this finding is that analytical estimates of interaction rates in a GGC must take into account the fact that the present day observed structural parameters including has not been constant and is still evolving. Hence, to calculate a correct estimate one must integrate the time dependent cross-sections (based on the changing values of these parameters) over an appropriate length of time as has been done in, e.g., Fregeau (2008).
Acknowledgement: We thank the anonymous referee for his help rectifying some mistakes and useful suggestions. This work was supported by NASA ATP Grant NNX09AO36G at Northwestern University. SC acknowledges support from the Theory Postdoctoral Fellowship from UF Department of Astronomy and College of Liberal Arts and Sciences.
- Anders et al. (2009) Anders P., Lamers H. J. G. L. M., Baumgardt H., 2009, A&A, 502, 817
- Ashman & Zepf (2001) Ashman K. M., Zepf S. E., 2001, AJ, 122, 1888
- Baumgardt & Makino (2003) Baumgardt H., Makino J., 2003, MNRAS, 340, 227
- Casertano & Hut (1985) Casertano S., Hut P., 1985, ApJ, 298, 80
- Chatterjee et al. (2008) Chatterjee S., Fregeau J. M., Rasio F. A., 2008, in IAU Symposium Vol. 246 of IAU Symposium, Effects of Stellar Collisions on Star Cluster Evolution and Core Collapse. pp 151–155
- Chatterjee et al. (2010) Chatterjee S., Fregeau J. M., Umbreit S., Rasio F. A., 2010, ApJ, 719, 915
- Davis et al. (2008) Davis D. S., Richer H. B., Anderson J., Brewer J., Hurley J., Kalirai J. S., Rich R. M., Stetson P. B., 2008, AJ, 135, 2155
- Davis et al. (2010) Davis O., Clarke C. J., Freitag M., 2010, MNRAS, p. 844
- Djorgovski (1993) Djorgovski S., 1993, in S. G. Djorgovski & G. Meylan ed., Structure and Dynamics of Globular Clusters Vol. 50 of Astronomical Society of the Pacific Conference Series, Physical Parameters of Galactic Globular Clusters. p. 373
- Fregeau (2008) Fregeau J. M., 2008, ApJ, 673, L25
- Fregeau et al. (2003) Fregeau J. M., Gürkan M. A., Joshi K. J., Rasio F. A., 2003, ApJ, 593, 772
- Fregeau & Rasio (2007) Fregeau J. M., Rasio F. A., 2007, ApJ, 658, 1047
- Gürkan et al. (2004) Gürkan M. A., Freitag M., Rasio F. A., 2004, ApJ, 604, 632
- Harris et al. (1984) Harris H. C., Harris G. L. H., Hesser J. E., MacGillivray H. T., 1984, ApJ, 287, 185
- Harris (1996) Harris W. E., 1996, AJ, 112, 1487
- Heggie & Hut (2003) Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press, 2003
- Hesser et al. (1984) Hesser J. E., Harris H. C., van den Bergh S., Harris G. L. H., 1984, ApJ, 276, 491
- Hodge (1962) Hodge P. W., 1962, PASP, 74, 248
- Holtzman et al. (1992) Holtzman J. A., Faber S. M., Shaya E. J., Lauer T. R., Groth J., Hunter D. A., Baum W. A., Ewald S. P., Hester J. J., Light R. M., Lynds C. R., O’Neil Jr. E. J., Westphal J. A., 1992, AJ, 103, 691
- Hurley (2007) Hurley J. R., 2007, MNRAS, 379, 93
- Hwang & Lee (2008) Hwang N., Lee M. G., 2008, AJ, 135, 1567
- Joshi et al. (2001) Joshi K. J., Nave C. P., Rasio F. A., 2001, ApJ, 550, 691
- Joshi et al. (2000) Joshi K. J., Rasio F. A., Portegies Zwart S., 2000, ApJ, 540, 969
- King (1966) King I. R., 1966, AJ, 71, 64
- Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
- Leigh et al. (2012) Leigh N., Umbreit S., Sills A., Knigge C., de Marchi G., Glebbeek E., Sarajedini A., 2012, MNRAS, 422, 1592
- Mackey et al. (2007) Mackey A. D., Wilkinson M. I., Davies M. B., Gilmore G. F., 2007, MNRAS, 379, L40
- Mackey et al. (2008) Mackey A. D., Wilkinson M. I., Davies M. B., Gilmore G. F., 2008, MNRAS, p. 374
- Mateo (1987) Mateo M., 1987, ApJ, 323, L41
- McLaughlin & van der Marel (2005) McLaughlin D. E., van der Marel R. P., 2005, ApJS, 161, 304
- Miller et al. (1997) Miller B. W., Whitmore B. C., Schweizer F., Fall S. M., 1997, AJ, 114, 2381
- Noyola & Gebhardt (2006) Noyola E., Gebhardt K., 2006, AJ, 132, 447
- Pattabiraman et al. (2012) Pattabiraman B., Umbreit S., Liao W.-K., Choudhary A., Kalogera V., Memik G., Rasio F. A., 2012, ArXiv e-prints
- Scheepmaker et al. (2009) Scheepmaker R. A., Gieles M., Haas M. R., Bastian N., Larsen S. S., 2009, The Radii of Thousands of Star Clusters in M51 with HST/ACS. p. 103
- Scheepmaker et al. (2007) Scheepmaker R. A., Haas M. R., Gieles M., Bastian N., Larsen S. S., Lamers H. J. G. L. M., 2007, A&A, 469, 925
- Trenti et al. (2007) Trenti M., Ardi E., Mineshige S., Hut P., 2007, MNRAS, 374, 857
- Trenti et al. (2007) Trenti M., Heggie D. C., Hut P., 2007, MNRAS, 374, 344
- Trenti et al. (2010) Trenti M., Vesperini E., Pasquato M., 2010, ApJ, 708, 1598
- Umbreit et al. (2012) Umbreit S., Fregeau J. M., Chatterjee S., Rasio F. A., 2012, ApJ, 750, 31
- van den Bergh et al. (1991) van den Bergh S., Morbey C., Pazder J., 1991, ApJ, 375, 594
- Vesperini & Chernoff (1994) Vesperini E., Chernoff D. F., 1994, ApJ, 431, 231
- Whitmore & Schweizer (1995) Whitmore B. C., Schweizer F., 1995, AJ, 109, 960