Matching Globular Cluster Models to Observations

The CMC Cluster Catalog comprises 148 models, spanning a realistic and comprehensive range of initial virial radii, tidal radii, metallicities, and masses (Table 1), integrated via CMC. Within CMC, stellar evolution is modeled using the single-star evolution (sse, Hurley et al. 2000) and binary-star evolution (bse, Hurley et al. 2002) algorithms, updated to include the most current prescriptions for compact object formation (see, e.g., Breivik et al. 2020). These prescriptions describe the evolution of stars/stellar objects through various evolutionary stages, distinguished in the code by "stellar type" (see Section 1 of Hurley et al. 2000 for a list of stellar types and discussion). CMC also incorporates the physics of three/four-body encounters, integrated using the fewbody package (Fregeau et al. 2004; Fregeau & Rasio 2007), and updated to include post-Newtonian terms (Rodriguez et al. 2018). Our models assume that GCs experience a constant tidal field throughout their lifetimes. Of course, in general, GCs undergo complicated orbits, characterized by periapse passages which may affect their dynamics in a nonlinear fashion, and there is evidence to suggest that the galactic potential has itself varied over time in the history of the Milky Way (e.g., Kruijssen et al. 2019). Further exploration of the effects of such time-dependent tidal fields falls beyond the scope of this work.

At various snapshots in time, separated by multiples of the estimated dynamical time of the cluster, CMC writes a catalog of stellar and kinematic properties for all stars in the cluster. In this study, we are primarily interested in cluster models similar to the old GCs observed in the Milky Way; we therefore restrict our attention to snapshots for which t > 10 Gyr (all models are run to t ≈ 14 Gyr, or until dissolution), of which there are 7537 throughout the entire CMC Cluster Catalog.

In general, the parameters most germane to dynamical structure are the initial number of stars, N, and the virial radius, r_v . The galactocentric distance, R_g , is most impactful through its influence at the outskirts of the cluster as it defines the tidal radius, while the metallicity, Z, primarily influences stellar evolution. Moreover, both R_g and Z are more easily estimated empirically than N and r_v , which often change dramatically over the course of a cluster's lifetime (e.g., Kremer et al. 2020). We therefore limit our fitting procedure for each cluster to only those models where R_g and Z are closest (in linear and logarithmic scales, respectively) to their observed present-day values, as reported by Baumgardt et al. (2019b), and Harris (1996, 2010 edition), respectively. Hence, for any individual cluster, we only optimize over N and r_v . For simplicity, a constant tidal radius is assumed, although we caution that GC orbits in the Galaxy generally induce time-dependent tidal forces (including possible close pericenter passages). Furthermore, the modern-day distance of a GC to the Galactic center may not be representative of the average tidal force (Baumgardt et al. 2019b).

In order to match an observed GC with a best-fit model on our grid, we identify models whose dynamical properties most closely match the observed cluster features. The most direct dynamical observables of a GC are its surface brightness and velocity dispersion profiles; we therefore extract a simulated SBP and VDP from each model snapshot for comparison to the corresponding observed profiles.

Since GCs are observed only in projection on the sky, we project our simulated stellar positions and velocities onto a two-dimensional plane by assuming spherical symmetry. In particular, a star with a three-dimensional radius, r, has a probability p(a, b; r) of lying within projected radial distances d = a and d = b > a, given by

$$\begin{eqnarray}&&p(a,b;r)=\sqrt{1-\min {\left(a,r\right)}^{2}/{r}^{2}}-\sqrt{1-\min {\left(b,r\right)}^{2}/{r}^{2}}.\end{eqnarray} \tag{ 1 }$$

We calculate the surface brightness, Σ_V (a, b), and one-dimensional velocity dispersion, σ_v (a, b), for 80 logarithmically spaced bins, with an inner bin of 10⁻³ pc, and an outer bin given by the maximum radial position of a star in the catalog. For σ_v (a, b), we only include evolved bright stars (sse/bse with stellar types 2-9), so as to mimic the use of bright stars in real VDP measurements (e.g., Kamann et al. 2017; Ferraro et al. 2018b).

Given the two-dimensional distribution described by Equation (1), the average V-band surface brightness, Σ_V (a, b), in a projected radial bin bounded by d ∈ (a, b] can be calculated as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}_{V}(a,b) & = & \\ & & -2.5{\mathrm{log}}_{10}\left(\displaystyle \frac{100\,{\mathrm{pc}}^{2}\,{\mathrm{arcsec}}^{2}}{\pi ({b}^{2}-{a}^{2})}\displaystyle \sum _{i}\displaystyle \frac{p(a,b;{r}_{i})}{{10}^{{M}_{V,i}/2.5}}\right)+{A}_{V},\end{array}\end{eqnarray} \tag{ 2 }$$

where r_i and M_V,i are the three-dimensional radial distance and absolute V-band magnitude of the ith star in the simulation, respectively, and A_V is the V-band extinction. For simplicity, all stars are assumed to be blackbodies, which should reasonably approximate their actual magnitudes, particularly for more massive stars, where molecular lines are less prominent. Stellar magnitudes thus take the form

$$\begin{eqnarray}&&{M}_{V,i}=-2.5{\mathrm{log}}_{10}\left(\displaystyle \frac{{R}_{i}^{2}}{{(10\mathrm{pc})}^{2}{F}_{\mathrm{ZP},\lambda }}\displaystyle \frac{\int f(\lambda ,{T}_{i}){ \mathcal T }(\lambda )d\lambda }{\int { \mathcal T }(\lambda )d\lambda }\right),\end{eqnarray} \tag{ 3 }$$

where F_ZP,λ = 3.57453 × 10⁻⁹ erg cm⁻² s⁻¹ Å⁻¹ is the zero-point spectral flux, f(λ, T) is the wavelength-space Planck distribution for temperature T, and ${ \mathcal T }(\lambda )$ is the transmission function for the filter (Casagrande & VandenBerg 2014). Photometric magnitudes are derived using the generic Johnson V-band filter function and Vega zero point from the SVO Filter Profile Service, a public repository for astronomical filter parameters (Rodrigo & Solano 2013; Rodrigo et al. 2020). The V-band extinction, A_V , of a GC is computed using the standard Cardelli et al. (1989) extinction law as A_V = 3.1 × E(B − V), where E(B − V) is taken from the Harris (1996, 2010 edition) catalog. While the blackbody approximation deviates significantly for cool M dwarfs, which have significant molecular absorption lines (Allard et al. 1994; Baraffe et al. 1995), the approximation is a reasonable estimate for brighter and hotter stars, whose continuum emission dominates the profile. Equation (3) is also applied in Section 4.2 to select binaries, although the cut is restricted to relatively bright main sequence stars, where the blackbody approximation should be expected to hold.

The velocity dispersion, σ_v (a, b), in the same radial bin is given by

$$\begin{eqnarray}&&{\sigma }_{v}(a,b)=\sqrt{\displaystyle \frac{\displaystyle \sum _{i}| {{\boldsymbol{v}}}_{i}{| }^{2}p(a,b;{r}_{i})}{3\cdot \displaystyle \sum _{i}p(a,b;{r}_{i})}},\end{eqnarray} \tag{ 4 }$$

where v _i is the three-dimensional velocity vector of the ith star. This expression for σ_v (a, b) assumes an isotropic velocity dispersion. Although, in principle, the tangential and radial components of the velocity dispersion may differ, in most cases the ratio of the two is very close to one, particularly for cases near to a cluster's center (Watkins et al. 2015). Furthermore, although some have claimed detection of coherent rotation within many GCs (e.g., Kamann et al. 2017), this subtle behavior is not captured in our models, and we do not consider departures from spherical symmetry in this work. Such rotation is, in any case, usually much smaller than the velocity dispersion, and would not be expected to change it significantly, especially in the central regions of the GC, where CMC is expected to be most accurate. While spherical asymmetry is the subject of very interesting work (see, e.g., Gieles et al. 2021, where measurements of Palomar 5's tidal tails are used in fitting), it lies beyond the scope of this work.

Here, we assess the relative likelihood that a given CMC model fits observed SBPs/VDPs by computing χ² statistics between this data and linearly interpolated model SBPs/VDPs (Heggie & Giersz 2008, 2014; Giersz & Heggie 2009, 2011; Kremer et al. 2018, 2019). The fitness of a model with a given GC is assessed using ${\tilde{\chi }}_{\mathrm{SBP}}^{2}={\chi }_{\mathrm{SBP}}^{2}/{N}_{\mathrm{SBP}}$ and ${\tilde{\chi }}_{\mathrm{VDP}}^{2}\,={\chi }_{\mathrm{VDP}}^{2}/{N}_{\mathrm{VDP}}$, the χ² statistics between the model SBP/VDP, and the observations, normalized by the number of observational data points. For a given GC, we consider well-fitting snapshots to have the fitting heuristic $s\equiv \max \left({\tilde{\chi }}_{\mathrm{SBP}}^{2},{\tilde{\chi }}_{\mathrm{VDP}}^{2}\right)\,\lt 10$. As such, to be a "good fit" to the data, a snapshot must be a reasonably good fit to both the SBP and VDP. For diagnostic reasons, we also report the best-fitting snapshot of each cluster ${\tilde{\beta }}_{\mathrm{SBP}}^{2}$ and ${\tilde{\beta }}_{\mathrm{VDP}}^{2}$, defined as the reduced χ² sums with terms weighted by the sign of the residual. These statistics parameterize the extent to which a given snapshot overestimates or underestimates a cluster's surface brightness or velocity dispersion, and are included to guide the creation of future models in order to better fit particular observed GCs.

The SBPs are taken from ground-based observations by Trager et al. (1995). ⁷ While other data sets, such as the Noyola & Gebhardt (2006) SBPs for 38 GCs, may better probe the surface brightnesses of GC cores, particularly with respect to core-collapsed clusters, we opt to exclude this data, in order to avoid assigning ad hoc relative weights between the Trager and Noyola profiles. The VDPs are taken from the radial velocity measurements of Kamann et al. (2017), Ferraro et al. (2018b), and Baumgardt & Hilker (2018), as well as the proper motion measurements of Watkins et al. (2015), and Baumgardt et al. (2019b). The VDPs measured using radial velocities and proper motions are generally observed to be consistent with one another; cluster–VDP uncertainties across all data sets are therefore taken as reported, without any rescaling or homogenization, and we fit to a combination of these VDPs for maximal constraint.

As Trager et al. (1995) do not calculate formal uncertainties on their SBP measurements, but instead provide a Chebyshev polynomial fit, and coarse quality weights, w_i , for each data point, we follow the procedure adopted in McLaughlin & van der Marel (2005) in order to estimate uncertainties. In particular, we assume that the measurement uncertainties are inversely proportional to w_i , and that the third-order Chebyshev polynomial fits have ${\tilde{\chi }}^{2}=1$ exactly. We then estimate the uncertainty of the ith data point to be δΣ_V = δΣ_V,0/w_i , where δΣ_V,0 is estimated separately for each cluster as $\delta {{\rm{\Sigma }}}_{V,0}\,=\sqrt{{({N}_{\mathrm{SBP}}-4)}^{-1}{\displaystyle \sum }_{i=1}^{N}{w}_{i}^{2}{e}_{{{\rm{\Sigma }}}_{V},i}^{2}}$, where ${e}_{{{\rm{\Sigma }}}_{V},i}$ is the surface brightness residual of the ith data point with respect to the Chebyshev polynomial fit. Since this effectively discards observations for which Trager et al. (1995) assign w_i = 0, we omit these points when normalizing ${\tilde{\chi }}_{\mathrm{SBP}}^{2}$ by N_SBP. As the publicly available SBP for NGC 2419 appears to be multi-valued, we restrict our profile to data sets given in Trager et al. (1995), following the branch shown in the plot in their paper, and we estimate the SBP uncertainties using these points alone.

Figure 1 shows the GCs for which both the SBP and VDP are sampled with at least five points, in both the core radius–brightness (r_c –M_V ) and galactocentric distance–metallicity (R_g –[Fe/H]) planes, with M_V , r_c , and [Fe/H] taken from Harris (1996, 2010 edition), and R_g from Baumgardt et al. (2019b). Points are color coded by the fitting statistic, s, with well-fit clusters circled. Despite the wide range of Milky Way GC properties, we are able to satisfactorily fit a wide range of clusters across the observed parameter space. As expected, well-fit clusters are concentrated at lower brightnesses, specifically at dimmer M_V ≳ − 9.5, indicating a lack of grid coverage at larger masses. Unsurprisingly, the quality of the fit does not obviously correlate with present-day galactocentric distance or metallicity, as the bulk cluster dynamics are less sensitive to these parameters.

Zoom In Zoom Out Reset image size

In the sections below, we examine best-fitting models for seven specific GCs, during which we make comparisons to the observed population of low-mass X-ray binaries, millisecond pulsars, and cataclysmic variables. Although we do not use these stellar exotica as factors in our goodness-of-fit measurements, due to their large uncertainties, we can use the rough numbers as guidelines to further constrain and explore our models.

An X-ray binary (XRB) is a mass-transferring binary where the donor, typically a main sequence star, accretes onto a compact object, either a neutron star or a black hole. While short-lived high-mass X-ray binaries (with OB-type donors) dominate the X-ray sources in young, star-forming populations, low-mass X-ray binaries (LMXBs) are believed to form in the dense cores of GCs (Pooley et al. 2003; Fabbiano 2006). In our models, we consider as XRBs any main sequence star (sse stellar type 0–1) in a mass-transferring binary with a neutron star or black hole. Our models contain characteristically between zero and a few XRBs.

One possible outcome of a disrupted LMXB is a millisecond pulsar. Millisecond pulsars (MSPs) are rapidly rotating pulsars with periods on the order of milliseconds. Unlike standard pulsars, which are young neutron stars rotating fast enough to beam, MSPs are older neutron stars which have been "recycled" by accretion from a companion. Whereas the former have periods P ≃ 0.1–3 s, mass transfer onto the latter spins MSPs up to periods as small as 1.5 ms (Lorimer 2008). Such objects are expected to be formed at an enhanced rate in the high-density environments at the center of GCs, particularly in core-collapsed clusters with very few black holes (Ye et al. 2019). In our models, we identify MSPs as neutron stars with periods P < 10 ms.

Cataclysmic variables (CVs) are usually low-mass main sequence stars undergoing mass transfer onto a white dwarf via Roche-lobe overflow. As their name suggests, they are characterized by large, often rapid flux variability, in some cases undergoing violent eruptions in the form of novae or dwarf novae (Robinson 1976). GCs are believed to harbor significant CV populations, which may help to guide studies of their evolution, as well as their possible role as SNe Ia progenitors (Ivanova et al. 2006; Knigge 2012; Maoz et al. 2014). Similarly to XRBs, we identify CVs as sse stellar type 0 stars undergoing mass transfer onto white dwarfs. Our models contain a wide range of CVs, spanning from a few to ∼100 CVs. Interestingly, given that most CVs originate from primordial binary progenitors, the number of CVs in a GC actually correlates inversely with its central density. This scaling is in contrast to that of other objects (e.g., MSPs), whose formation is predominantly dynamical.

Blue straggler stars are unusually bright/blue main sequence stars which have been rejuvenated, either via accretion from another star (e.g., Chen & Han 2004), or a collision at some point in the GC's history (e.g., Glebbeek et al. 2008). With links to both standard binary evolution and cluster dynamics, the radial distribution of blue stragglers within GCs has gained significant attention as a possible tracer of a cluster's dynamical history (e.g., Ferraro et al. 2012, 2018a). We leave the modeling of these interesting sources to a future work.

Reaching about 0.5 kpc at perigalacticon, NGC 6293 is a low-metallicity, core-collapsed GC which lies very close to the Galactic center—an intense tidal environment which spatially flattens the cluster's shape (Chen & Chen 2010; Baumgardt et al. 2019b). Within the CMC Cluster Catalog, we locate 59 well-fitting snapshots from two core-collapsed models, n8-rv1-rg2-z0.01 (8 snapshots), and n8-rv0.5-rg2-z0.01 (51 snapshots), both with an initial N = 8 × 10⁵. While the SBPs are generally matched quite well within the cores, some snapshots slightly overestimate the core brightness, and some underestimate it. As such, well-fitting snapshots should produce predictions for cluster properties (e.g., total mass) which "surround" their true values. While all well-fitting snapshots underestimate the VDP somewhat, observational uncertainties are comfortably large enough to be consistent with predictions.

NGC 6293 is associated with at least one soft X-ray source (XTE J1709-267, Jonker et al. 2003). Of the 59 well-fitting snapshots, 54 snapshots do not have any XRBs, and five snapshots contain a single XRB (three from the r_v = 1 pc model pc, and two from the r_v = 0.5 pc model). Furthermore, all of these XRBs have low-mass donors with sse stellar type 0 (M ≲ 0.7 M_⊙). Despite possible observational biases, and incompleteness on the total number of XRBs in the cluster (especially for a cluster near the Galactic Center), our models appear consistent with the capacity of NGC 6293 to produce a small number of X-ray sources.

NGC 6397 is a nearby (D = 2.44 ± 0.04 kpc, Baumgardt et al. 2019b) metal-poor, core-collapsed GC, whose close proximity has attracted significant study of its white dwarf and low-mass stellar populations (Paresce et al. 1995; Cool et al. 1996; Taylor et al. 2001; Hansen et al. 2007). In addition, NGC 6397 has 15 known CV candidates, a quiescent LMXB (qLMXB), and 1–2 MSPs (Cool et al. 1995; Grindlay et al. 2001; Cohn et al. 2010; Dieball et al. 2017). We locate 11 well-fitting snapshots from the model, n4-rv1-rg8-z0.01. While all of these snapshots tend to overestimate the surface brightness, and underestimate the velocity dispersion somewhat, they do so within our tolerances. By eye, it appears that the slight overestimation of the SBP occurs primarily at ≈5 pc, where the model is slightly brighter than the data.

The well-fitting snapshots for this cluster each have between 11 and 13 CVs. In reasonable agreement with observations, these snapshots also show between 1 and 2 XRBs each (all of which have M ≲ 0.7 M_⊙ donors). Finally, these snapshots all contain two MSPs. Of course, while these small numbers should not be taken as precise predictions, they provide a measure of reassurance that our models generate these populations in reasonable numbers.

NGC 6681 (M70) is a core-collapsed cluster, which has occasionally been used to derive distortion solutions for instruments operating in the far-ultraviolet (e.g., Sohn 2018). Like NGC 6293, NGC 6681 lies quite close to the Galactic center (R ∼ 0.8 kpc, Baumgardt et al. 2019b), and observations also suggest some degree of tidal deformation (Han et al. 2017). We find 49 well-fitting snapshots from the core-collapsed model, n8-rv0.5-rg2-z0.01. As this particular model also provides snapshots consistent with the SBP and VDP of NGC 6293, one can similarly view NGC 6681 as somewhat similar to NGC 6293 from a dynamical perspective. As with NGC 6293, all snapshots slightly underestimate the VDP, while remaining within tolerance in terms of observational uncertainties.

A relatively metal-rich, non-core-collapsed GC, NGC 288 garnered much attention in the 1990s and early 2000s, as its similar metallicities and distances to those of NGC 362 provided promising avenues for constraining the age difference between the two clusters (Green & Norris 1990; Sarajedini & Demarque 1990; Bellazzini et al. 2001). There is a large number of subsequent studies into its dynamics (e.g., Piatti 2018) and stellar populations (e.g., Roh et al. 2011). Within our model grid, we identify 16 well-fitting snapshots from a large r_v = 4 model, n4-rv4-rg8-z0.1. Both the SBPs and VDPs from these snapshots fit the observations exceptionally well.

Using Chandra, Kong et al. (2006) reported between two and five possible CVs or other chromospherically active binaries within the half-mass–radius of NGC 288. Meanwhile, our models predict between 27 and 33 CVs. This discrepancy may be attributable to the high temporal variability in the activity of CVs, which could make quiescent accreting binaries difficult to detect, or to some other source of observational incompleteness.

While Trager et al. (1995) do measure the V-band SBP for this cluster, they do not report data outside the 175 core radius reported by Harris (1996, 2010 edition). As such, the SBP neither contains the characteristic turnover in brightness, nor constrains particularly well the size of the cluster's core. Therefore, while the CMC Cluster Catalog includes 30 well-fitting snapshots each from the N = 8 × 10⁵ models, n8-rv4-rg8-z0.01, and n8-rv2-rg8-z0.01, the initial virial radius of the cluster is largely uncertain. Although the snapshots come from models with different initial r_v , their similar initial masses and dynamical states produce a relatively small spread in the predicted mass of the cluster. Kacharov et al. (2014) estimated a cluster mass M = (2.0 ± 0.5) × 10⁵ M_⊙, compatible with the estimated mass from snapshots at the 1σ level. They also estimate a mass-to-light ratio of between 1.4 and 2.3 for this cluster, which is broadly compatible with predicted values for all seven of our clusters of interest (see Section 4.1).

Within NGC 4372, Kaluzny & Krzeminski (1993) identify a candidate CV, although atmospheric effects imply a high degree of incompleteness in the data. Using X-ray observations from XMM-Newton, Servillat et al. (2008) were unable to detect specific CVs—although these were limited by X-ray luminosity—although they report unresolved emission, consistent with a population of ∼20 CVs. As our well-fitting snapshots come from two models with different initial r_v , we find a bimodal distribution, with modes approximately around 26 and 84 CVs, and with well-fitting snapshots having between 21 and 98. The lower and higher modes correspond to snapshots from the r_v = 2.0 pc and r_v = 4.0 pc models, respectively. Intuitively, this correlation between the number of CVs and r_v (inversely related to the central density) is due to the origination of most CVs in primordial binaries, which are less likely to be disruptive in a less dense cluster (Kremer et al. 2020). The observed number of CVs is consistent with the predicted number from the r_v = 2 pc model.

NGC 5897 is a non-core-collapsed cluster at a distance D ≈ 12.6 kpc (Baumgardt et al. 2019b). While the VDP has been measured using both radial velocities (Baumgardt & Hilker 2018) and proper motions (Baumgardt et al. 2019b), these measurements have not extended into the core of the cluster, and thus do not constrain the dark mass distribution particularly well. We find 14 well-fitting snapshots in the model, n4-rv4-rg8-z0.01, which differs from the well-fitting model for NGC 288 only in that its metallicity is lower. Thus, between the two clusters, one might expect similarity in their dynamics, but not necessarily in their stellar populations.

One of a handful of clusters with γ-ray emission > 100 MeV (Tam et al. 2011), NGC 6624 is an interesting, high-metallicity ([M/H] ≈ − 0.44; Harris 1996, 2010 edition) GC, known to contain at least four MSPs, in addition to two young pulsars (Biggs et al. 1994; Lynch et al. 2012). In the past, it has been argued that the relatively large spin-period derivatives of an MSP near the cluster center is evidence for an IMBH (Peuten et al. 2014; Perera et al. 2017), although these signals have since been found to be consistent with dynamical interactions alone (Gieles et al. 2018; Baumgardt et al. 2019a). Moreover, the cluster is known to contain at least one well-studied LMXB (4U 1820-30; Biggs et al. 1994). Significantly, the cluster lacks any well-fitting snapshots from the unmodified CMC Cluster Catalog—the best-fitting snapshot belongs to the model n8-rv0.5-rg2-z1.0, for which ${\tilde{\chi }}_{\mathrm{SBP}}^{2}=11.95$. We present NGC 6624 as a test case for introducing new models, in order to fit a known SBP and VDP, not satisfactorily fit by the main CMC Cluster Catalog.

Without the addition of any other models, the CMC Cluster Catalog already provides a satisfactory fit for 26 out of 59 of the GCs for which we attempt to locate an analogous model. Of the remaining GCs, there are broadly three reasons why a GC may not be fit well by the CMC Cluster Catalog. Firstly, the GC's parameters may lie entirely outside of the range of parameter space probed by the CMC Cluster Catalog, although further simulations outside of this parameter space (e.g., new models with larger N than the largest N on the grid) may fit such clusters. Secondly, the GC may lie between points in parameter space sampled by the grid; for example, a GC's initial r_v may lie between 1 and 2 pc, but not especially close to either. In this case, the inclusion of models which increase the resolution of the grid would fit these clusters. Finally, the GC may differ from the models in other parameters besides those which have been varied in the grid, e.g., r_v , R_g , Z, and N. Such parameters include the initial binary fraction and initial mass function, which are not varied in the CMC Cluster Catalog, but can generally be varied in CMC to potentially improve GC fits.

NGC 6624 does not have a well-fitting snapshot, according to s < 10. However, as its closest-fitting snapshots on the main grid do not appear to deviate substantially from the observed data, it is likely to lie within the second aforementioned category: the lack of good fits is probably due to the coarseness of the grid, rather than its limited range. Nevertheless, its closest-fitting snapshots provide guidance as to the manner in which the grid ought to be extended so as to obtain a good fit. Guided by these pre-existing models, we supplement the model grid to better fit NGC 6624, in order to demonstrate targeted cluster fitting with further CMC simulations.

Notably, ${\tilde{\chi }}_{\mathrm{SBP}}^{2}=11.95$ for this cluster, indicating that some model SBPs are only slightly discrepant with the observed SBP. Upon inspection, this disagreement arises because the model's results for the SBPs' core brightnesses are higher than observed, while the outer halo brightnesses are slightly lower than observed. Modification of the initial N can narrow these discrepancies (either directly, by altering visible mass, or indirectly, by changing the relaxation time), as can modifying r_v (which has been shown to influence progress toward core collapse, see Kremer et al. 2018).

To better fit NGC 6624, we accordingly run six new models, and examine the extent to which they improve (or fail to improve) the fit (Table 2). We consider separately a decrease in initial N to 7 × 10⁵ particles (n7-rv0.5-rg2-z1.0), and an increase in the virial radius to r_v = 0.7 pc (n8-rv0.7-rg2-z1.0). We also consider an increase in the virial radius to r_v = 1.0 pc under three distinct initial particle counts, N = (6, 7, 9) × 10⁵, (n6-rv1-rg2-z1.0, n7-rv1-rg2-z1.0, and n9-rv1-rg2-z1.0, respectively). Finally, for N = 7 × 10⁵, we also consider a decrease in metallicity to Z = 0.35 Z_⊙ (n7-rv1-rg2-z0.35), i.e., the observed metallicity reported by Harris (1996, 2010 edition). As the metallicity is coarsely rounded up to Z = 1.0 Z_⊙ in the main CMC Cluster Catalog, we consider this latter model, in order to ascertain whether or not the structure of the GC displays fine sensitivity to Z.

Table 2. New Models for NGC 6624

Model	Nsnap	Ngood	rv	Rg	[M/H]	N	${\tilde{\chi }}_{\mathrm{SBP}}^{2}$	${\tilde{\chi }}_{\mathrm{VDP}}^{2}$	${\tilde{\beta }}_{\mathrm{SBP}}^{2}$	${\tilde{\beta }}_{\mathrm{VDP}}^{2}$
			(pc)	(kpc)		(×105)
n7-rv0.5-rg2-z1.0	570	0	0.5	2	0	7	14.51	2.00	6.23	0.78
n8-rv0.7-rg2-z1.0	437	0	0.7	2	0	8	12.24	3.37	9.11	−2.70
n7-rv1-rg2-z0.35	306	13	1	2	−0.46	7	5.22	2.41	3.04	−2.17
n6-rv1-rg2-z1.0	233	0	1	2	0	6	11.38	1.76	7.83	−0.15
n7-rv1-rg2-z1.0	218	0	1	2	0	7	13.60	4.05	10.94	−3.90
n9-rv1-rg2-z1.0	317	0	1	2	0	9	21.66	13.20	21.35	−13.20

Note. Parameters of additional models generated to better fit NGC 6624, alongside goodness-of-fit measures for closest-fitting snapshots.

Download table as:  ASCIITypeset image

Of these six models, we find that only n7-rv1-rg2-z0.35 provides well-fitting snapshots (13, Figures 6 and 7). Within the other models, the best-fitting snapshots from n7-rv0.5-rg2-z1.0, n8-rv0.7-rg2-z1.0, and n9-rv1.0-rg2-z1.0 all contain the same central overdensity as the best-fitting snapshots from the CMC Cluster Catalog proper. The models n7-rv1-rg2-z1.0 and n6-rv1-rg2-z1.0 appear to decrease the core overdensity in the models, although they appear to be slightly overbright at r_v ≈ 6'' ≈ 0.2 pc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Of the 13 snapshots that fit NGC 6624 well, all lack MSPs. Most of these snapshots also lack any XRBs, with one snapshot containing a single XRB, and two snapshots containing two XRBs. However, 16 additional snapshots pass the slightly relaxed fitting criterion, s < 15, with six from different models: n6-rv1-rg2-z1.0 (with one MSP and two to four XRBs), n7-rv1-rg2-z1.0 (with two MSPs, and one XRB), and n8-rv0.7-rg2-z1.0 (four to five MSPs, and zero to one XRBs). Our models are consistent with NGC 6624's single observed LMXB. Although we find fewer MSPs in our well-fit models, the difference is small, and may be explained by moderate sensitivity of this number to the initial size, compactness, and metallicity of the cluster. Since current observations are not necessarily complete, this could also suggest that our models may somewhat under-produce these types of stellar exotica. However, it is also important to keep in mind that NGC 6624 exists in a more complex high-metallicity regime, and lies behind a moderate amount of extinction. Coupled with uncertainty in terms of the distance from Gaia (D = 7.19 ± 0.37 kpc, Baumgardt et al. 2019b) and the initial mass function (which may play a significant role; see Weatherford et al. 2021), it remains plausible that both observational and modeling uncertainties could account for the deficit in MSPs and XRBs in the best-fitting models.

Although cluster brightness and extent are readily observable quantities, their precise translation to total cluster mass is complex, and generally requires dynamical modeling. Using scaled-up versions of N-body simulations, where N = (1–2) × 10⁵, Baumgardt & Hilker (2018) estimated the total masses and mass-to-light ratios of a number of GCs using radial velocities, including the specific clusters discussed in this paper. They further refined the mass-to-light calculations with additional measurements of L_V , where they found typical values of M/L_V ∼ 1.8 (Baumgardt et al. 2020). Here, we perform analogous estimates for the total cluster mass, M, and mass-to-light ratio, M/L, where L refers to bolometric luminosity (Figure 8).

Zoom In Zoom Out Reset image size

The seven clusters we examine are estimated to have present-day masses ranging between ∼ 9.2 × 10⁴ M_⊙ and ∼ 1.1 × 10⁵ M_⊙, and M/L between ∼1.5 and ∼1.9 (Table 3). Our mass estimates are very consistent with those of Baumgardt et al. (2020), except in the case of NGC 5897, where they estimate a significantly larger mass, M ∼ 2 × 10⁵ M_⊙, and NGC 6624, where they estimate a significantly lower mass M ∼ 7 × 10⁴ M_⊙. Our mass-to-light ratios are very consistent with those of Baumgardt et al. (2020) in all cases, and are narrowly scattered around M/L ∼ 1.8.

Table 3. Masses and Mass-to-light Ratios for Seven GCs

	Mass M (M⊙)		Mass-to-light Ratio M/L (M⊙/L⊙)
Cluster	This work	Baumgardt & Hilker (2018)	This work	Baumgardt & Hilker (2018)	Baumgardt et al. (2020)
NGC 6293	$\left({1.26}_{-0.15}^{+0.61}\right)\times {10}^{5}$	(1.88 ± 0.18) × 105	${1.49}_{-0.53}^{+0.24}$	1.67 ± 0.29	1.75 ± 0.31
NGC 6397	$\left({9.25}_{-0.16}^{+0.65}\right)\times {10}^{4}$	(8.89 ± 0.16) × 104	${1.85}_{-0.30}^{+0.25}$	2.18 ± 0.34	1.58 ± 0.10
NGC 6681	$\left({1.20}_{-0.10}^{+0.47}\right)\times {10}^{5}$	(1.13 ± 0.02) × 105	${1.49}_{-0.53}^{+0.18}$	2.00 ± 0.28	1.84 ± 0.11
NGC 288	$\left({1.13}_{-0.03}^{+0.03}\right)\times {10}^{5}$	(1.16 ± 0.03) × 105	${1.93}_{-0.32}^{+0.44}$	2.39 ± 0.17	2.14 ± 0.15
NGC 4372	$\left({2.39}_{-0.08}^{+0.09}\right)\times {10}^{5}$	(2.49 ± 0.25) × 105	${1.79}_{-0.42}^{+0.31}$	1.89 ± 0.19	2.10 ± 0.18
NGC 5897	$\left({1.12}_{-0.03}^{+0.03}\right)\times {10}^{5}$	(2.03 ± 0.21) × 105	${1.67}_{-0.29}^{+0.23}$	3.05 ± 0.43	2.19 ± 0.29
NGC 6624	$\left({1.53}_{-0.06}^{+0.08}\right)\times {10}^{5}$	(7.31 ± 0.20) × 104	${1.89}_{-0.22}^{+0.23}$	1.02 ± 0.13	1.50 ± 0.16

Note. Cluster mass and mass-to-light ratio for NGC 6293, NGC 6397, NGC 6681, NGC 288, NGC 4372, NGC 5897, and NGC 6624. The uncertainty bars reported here are taken to span the entire range of values for M and M/L appearing in well-fitting snapshots for a given cluster. Values of M from Baumgardt & Hilker (2018), as well as M/L from Baumgardt & Hilker (2018), and Baumgardt et al. (2020) for these clusters are also reproduced above.

Download table as:  ASCIITypeset image

Within GCs, the binary fraction is a photometrically observable property, sensitive to cluster dynamics, particularly in GC cores. The dense environments provided by GCs frequently scatter and eject binary systems, and the dynamical formation of binaries is generally thought to be the halting mechanism for collapse in core-collapsed clusters after the expulsion of their black holes (e.g., Chatterjee et al. 2013). Moreover, as binary systems almost always have total fluxes equal to the sum of their component fluxes, main sequence binary systems can be found on the color–magnitude diagram in predictably brighter sequences above the main sequence defined by their mass ratios. In a GC, binaries can either persist from the initial formation of the cluster, or be dynamically generated over time—in the CMC Cluster Catalog, it is assumed in all cases that the primordial binary fraction is 5%, with a flat mass ratio distribution of between q = 0.1 and 1. Using the ACS Globular Cluster Survey, Milone et al. (2012) photometrically measured the binary fraction for 59 GCs, where the mass ratio q > 0.5, 0.6, and 0.7 individually within three radial regions, r < r_c , r_c < r < r_h , and r > r_h , where r_c and r_h are as given by Harris (1996, 2010 edition; Figure 9). Of these, five overlap with our seven clusters of interest (they did not report binary fractions for NGC 6293 or NGC 4372).

Zoom In Zoom Out Reset image size

We calculate the binary fraction, subject to the same minimum mass ratios and radial ranges, for all best-fitting snapshots for each of our seven clusters of interest (Figure 10). To mimic the magnitude cuts applied by Milone et al. (2012), we restrict our sample to a locus on the color–magnitude diagram consistent with binaries whose primary has an F814W magnitude between 0.75 and 3.75 mag below the main sequence turn-off. On the blue edge, we enforce that included sources must have an F606W–F814W color lying no bluer than 0.1 mag of the main sequence turn-off. To calculate a binary fraction for stars above a mass ratio, q, we use the sse main sequence prescription to define a locus on the color–magnitude diagram corresponding to the sum of fluxes due to two main sequence stars with mass ratio q. Sources on the red side of the locus are then considered to be binaries. The binary fraction is then calculated by dividing the weight of sources identified as binaries by the total weight of all sources in the magnitude and radial range, with weights defined by Equation (1). Notably, while Milone et al. (2012) additionally applied a general inner radius cut for a number of clusters, of our seven clusters of interest, only NGC 6681 and NGC 6624are affected by such a cut, and, in particular, only the r_c < r < r_h annulus (they do not report binary fractions for r < r_c ). However, in order to keep applied cuts relatively consistent between the clusters, we omit this cut for these particular GCs.

Zoom In Zoom Out Reset image size

We find reasonable consistency between the model binary fractions and the data in most cases. One notable exception is the binary fraction within the inner spatial bins of NGC 6397, where our models appear to predict binary fractions of up to ∼3 times the observed value—even here, the data and models come back into agreement in the outermost bin. Another is in the r_c < r < r_h bin of NGC 6624, where Milone et al. (2012) reported very low binary fractions, curiously lower than the binary fraction of the outskirts of the cluster (r > r_h ). However, in general, measurement of the binary fraction within core-collapsed clusters is relatively difficult for a number of reasons. Firstly, formally, core-collapsed clusters do not have a well-defined core radius, and observational definitions of the core radius tend to be quite small (e.g., r_c ∼ 0.03 pc for NGC 6397), with conversions from angular units to physical distances being very sensitive to heliocentric distance. As binaries are dynamically generated and then burned in large numbers in the cores of core-collapsed clusters, the binary fraction is expected to vary substantially with respect to distance from the GC center. Moreover, in the central regions of such clusters, stellar density is extremely large, implying relatively low completeness.

Overall, we note acceptable agreement between our predictions and observations for non-core-collapsed clusters, as well as the outer regions of core-collapsed clusters. This provides a degree of reassurance that CMC can sensibly replicate the binary populations of realistic GCs, although the fixed assumed binary fraction of 5% in a flat mass ratio distribution, combined with a complex initial to final binary fraction mapping, complicates the picture somewhat.

Over thermodynamically long timescales, massive stars are expected to sink to the center of the cluster through dynamical friction. It can be shown that a population of stars with mass m will segregate within a cluster on a timescale ${t}_{m,\mathrm{MS}}\sim \left(\left\langle m\right\rangle /m\right){t}_{\mathrm{rlx}}$, where $\left\langle m\right\rangle$ is the mean mass of the cluster (Portegies Zwart et al. 2010). As t_rlx ∼ few × 10⁹ yr for typical GCs, mass segregation in GCs is readily identified as a preferential clustering of massive stars closer to the center of the cluster. Moreover, given that "dark" objects, such as black holes, and other stellar remnants also participate in mass segregation, they may influence the observed metrics of mass segregation in nontrivial ways.

Here, we examine the capacity of the SBP and VDP alone to predict the degree of mass segregation within a cluster. In particular, in accordance with Weatherford et al. (2020), we define Population II stars as the "high-mass" population, where L_MSTO/5 < L < L_MSTO, and Population IV stars as the "low-mass" population, such that L_MSTO/125 < L < L_MSTO/25, and where L_MSTO is the luminosity of the main sequence turn-off. We then parameterize the mass segregation for each cluster, Δ₂₄, defined as the difference between the median projected radial distance of Population II stars and Population IV stars, normalized by the half-light radius of the cluster.

Note that Δ₂₄ is straightforwardly calculated both for simulations and observed clusters. In particular, using the ACS Globular Cluster Survey, we reproduce this calculation for four of our seven clusters of interest with suitable observations, following the procedure of Weatherford et al. (2018, 2020). We take into account incompleteness in the observed catalogs, estimated from artificial star tests. We then compare these to the distribution of Δ₂₄ within well-fitting snapshots for our seven clusters (Figure 11). Within the simulated data, we mimic the limited field of view of the data by restricting it to stars within a projected radius of 61% of the half-light radius (which is the highest value of the four clusters that can be accommodated by the data). In cases where comparison is possible, we find very strong agreement between the simulated and observed Δ₂₄ except in NGC 6397, where the simulations slightly overestimate the degree of mass segregation.

Zoom In Zoom Out Reset image size

While black holes (BHs) in GCs are very difficult to detect directly, their presence and number can be indirectly inferred by examining their effect on a GC's dynamical state. For example, by considering a grid of CMC models finely gridded over initial virial radius, Kremer et al. (2019) demonstrated the importance of BHs in the halting of core collapse in GCs NGC 3201, NGC 6656 (M22), and NGC 6254 (M10), and noted their likely absence in NGC 6752, which is core-collapsed. These BH populations are, in turn, intimately related to the cluster's initial size. Along a similar vein, trends in the core radius with age for massive clusters in Milky Way satellite galaxies in the seminal Mackey & Gilmore catalogs (Mackey & Gilmore 2003a, 2003b, 2003c) have been interpreted as evidence for the role of BHs in the clusters' bulk evolution (see, e.g., Mackey et al. 2008, which reproduces these observed trends in N ∼ 10⁵-body simulations).

Evidence of the effect of BH populations on the structure of a GC is also apparent in the spatial distribution of stars in different mass ranges. In particular, using CMC models, Weatherford et al. (2020) constrained the number of BHs in GCs by taking advantage of an anticorrelation between the extent of mass segregation in a cluster and its BH population, a trend quantifiable in CMC models. Intuitively, this anticorrelation arises from the rapid segregation of a GC's black hole population, followed by dynamical heating of the massive star population, and their typical distances from the cluster center.

For a similar reason, the presence of a BH population in a GC halts core collapse—large BH populations transfer significant energy to their host cluster's stellar population through binary-mediated dynamics, preventing core collapse of the bulk stellar population (e.g., Kremer et al. 2020). This manifests in terms of a flattened core surface brightness, as well as a heightened dynamical temperature in the core of the cluster. This motivates the use of observed SBPs and VDPs to constrain the BH populations of GCs, which can in turn be accomplished by analyzing the simulated stellar populations of their best-fitting CMC models.

For each of the seven GCs of interest, we calculate the median, 18th (−1σ) and 84th (+1σ) percentiles, and minimum and maximum number of BH in well-fitting snapshots (Table 4). As expected, all three of the core-collapse clusters examined have fully single-digit BH counts. Of their sample of 50 GCs, four of the GCs for which Weatherford et al. (2020) have estimated BH counts coincide with our seven: NGC 6397, NGC 6681, NGC 288, and NGC 6624. Reassuringly, our BH predictions are consistent with theirs when the metric parameterizing mass segregation is consistent with the definition in Section 4.3. Moreover, both Weatherford et al. (2020) and our work broadly reflect the tendency of core-collapsed clusters to have fewer BHs, reinforcing the scenario whereby black holes provide the dominant mechanism for halted core collapse in the majority of GCs today. Although both this work and that of Weatherford et al. (2020) calibrate N_BH to the same grid of models, we obtain constraints from two distinct observables (the SBP and VDP versus the degree of mass segregation in the cluster). This indicates at least concordance with the idea that both the (suppressed) degree of mass segregation and cluster dynamics are broadly driven by a single BH population at the center of a GC (or the lack thereof). Nevertheless, given a lack of direct observations of N_BH, the actual size of this population remains highly uncertain.