A COMPARISON BETWEEN THE STELLAR AND DYNAMICAL MASSES OF SIX GLOBULAR CLUSTERS

Among the stellar system zoo there is a widely recognized discrepancy between photon-emitting mass (stars and gas) and dynamical mass (including the contribution of dark matter). From galaxy clusters to dwarf galaxies, the mass of the baryonic component, estimated through the conversion of light into mass, has been found to be significantly smaller than that estimated via kinematical considerations (see Tollerud et al. 2011 for a recent review). The commonly accepted solution to this problem invokes the presence of a large amount of dark, non-baryonic matter driving the kinematics of these systems without contributing to their luminosity. Within this framework, globular clusters (GCs) are considered a notable exception: although they are the next step down from the smallest stellar systems containing dark matter (the Ultra Compact Dwarf galaxies; Mieske et al. 2008) and brighter than the most dark matter dominated systems (the Ultra Faint Galaxies; Simon et al. 2011), they are not believed to contain dark matter (Baumgardt et al. 2009). This conclusion is mainly justified by evidence that dynamical mass-to-light (M/L) ratios of GCs are ∼25% smaller than those predicted by simple stellar population models that assume a canonical initial mass function (IMF; McLaughlin & van der Marel 2005; Strader et al. 2009, 2011). However, these models do not account for dynamical effects such as the preferential loss of low-mass stars due to energy equipartition which alters the shape of the present-day mass function (PDMF) and consequently the stellar M/L ratio. In fact, all the past determinations of PDMFs in GCs derived slopes shallower than that of a canonical (e.g., Kroupa 2001) IMF (Piotto & Zoccali 1999; Paust et al. 2010; De Marchi et al. 2010). In a recent paper, Kruijssen & Mieske (2009) modeled the evolution of the mass function of 24 GCs including the effects of dynamical dissolution, low-mass star depletion, stellar evolution, and stellar remnants and found substantial agreement between the prediction of their model and the dynamical M/L ratio estimated through cluster kinematics. They conclude that dynamical effects are likely responsible for the observed discrepancy between stellar and dynamical masses. The flood of available data from deep photometric and spectroscopic surveys now make it possible to test this issue observationally at a higher level, sampling the mass function of GCs down to the hydrogen burning limit (Paust et al. 2010) and their velocity dispersion profile up to many half-mass radii (Lane et al. 2010b). This gives us an unprecedented opportunity to derive stellar masses from direct star counts instead of converting light into mass, thus obtaining an estimate which does not require the assumption of an M/L ratio.

In this paper, we compare the results of deep Hubble Space Telescope (HST) photometric observations and wide-field radial velocities available from recent public surveys with a set of dynamical models with the aim of deriving the amount of stellar and dynamical mass in a sample of six Galactic GCs.

For the present analysis, we make use of three different databases.

• 1.  
  The set of publicly available deep ACS@HST photometric catalogs of the "globular cluster treasury project" (Sarajedini et al. 2007);⁴
• 2.  
  The radial velocities obtained by the AAOmega@AAT survey of GCs (Lane et al. 2011);⁵
• 3.  
  The surface brightness profiles derived by Trager et al. (1995).⁶

The photometric database consists of high-resolution HST observations of a sample of 65 Galactic GCs. The database has been constructed using deep images secured with the Advanced Camera for Surveys (ACS) Wide Field Channel through the F606W and F814W filters. The field of view of the camera (202'' × 202'') is centered on the cluster's center with a dithering pattern to cover the gap between the two chips, allowing a full coverage of the core of all the six GCs considered in our analysis. This survey provides deep color–magnitude diagrams (CMDs) reaching the faint main sequence (MS) of the target clusters down to the hydrogen burning limit (at M_V ⩽ 10.7) with a signal-to-noise ratio S/N > 10. The results of artificial star experiments are also available to allow an accurate estimate of the completeness level and photometric errors. A detailed description of the photometric reduction, astrometry, and artificial star experiments can be found in Anderson et al. (2008). Since all the considered GCs have relaxation times shorter than their ages (McLaughlin & van der Marel 2005), the effects of mass segregation are expected to be important. It is therefore useful to constrain the radial variation of the mass function outside the half-mass radius. Unfortunately, the ACS treasury survey covers only the central part of the cluster at distances <2' from the cluster centers. Additional deep data sampling the external part of the cluster are available only for NGC 288. They consist of a set of ACS images observed under the program GO-12193 which are centered at ∼5' away from the cluster center. In particular, 3 × 200 s + 1 × 15 s long exposures have been taken through the F606W filter and 3 × 150 s + 1 × 10 s long exposures through the F814W one. All images were passed through the standard ACS/WFC reduction pipeline. Data reduction has been performed on the individual pre-reduced (.flt) images using the SExtractor photometric package (Bertin & Arnouts 1996). For each star, we measured the flux contained within a radius of 0.1 arcsec (corresponding to 2 pixel ∼FWHM) from the star center. We performed the source detection on the stack of all images while the photometric analysis was performed independently on each image. Only stars detected in two out of three long exposures or in the short ones have been included in the final catalog. We used the most isolated and brightest stars in the field to link the aperture magnitudes at 0.5 arcsec to the instrumental ones, after normalizing for exposure time. Instrumental magnitudes have been transformed into the VEGAMAG system by using the photometric zero points by Sirianni et al. (2005). Finally, each ACS pointing has been corrected for geometric distortion using the prescriptions by Hack & Cox (2001). We preformed artificial star experiments on the science frames following the prescriptions described in Sollima et al. (2007). In summary, a set of artificial stars have been simulated using the Tiny Tim model of the ACS PSF (Krist et al. 2010) and added to all images at random positions within 36 × 36 pixel cells centered on a grid of 100 × 50 knots along the x and y directions of the two ACS chips (a single star for each cell). The F814W magnitude of each star has been extracted from a luminosity function simulated by adopting a Kroupa (2001) mass function which has been converted to F814W magnitude using the mass–luminosity relation of Dotter et al. (2007). The corresponding F606W magnitude was derived by interpolating along the cluster mean ridge line of the F814W–(F606W–F814W) CMD. We performed the photometric reduction on the simulated frames, with the same procedure adopted for the science frames producing a catalog of ∼10⁵ artificial stars.

The radial velocity database has been derived by Lane et al. (2009, 2010a, 2010b, 2011) from a large number of spectra of red giant branch (RGB) stars observed in fields centered on 10 Galactic GCs. Spectra were obtained with the multi-fiber AAOmega spectrograph mounted at the Anglo Australian Telescope with the 1700D and 1500V gratings on the red and blue arms, respectively. With this configuration, spectra covering the Ca ii triplet region (8340–8840 Å) and the swathe of iron and magnesium lines around ∼5200 Å were obtained with a resolution of R = 10, 000 and R = 3700 for the red and blue arms, respectively. A detailed description of the reduction procedure and radial velocity estimates can be found in Lane et al. (2009, 2010a, 2010b, 2011). For the present work, we adopted the radial velocities extracted with the RAVE pipeline since they provide a better estimate of the radial velocity uncertainty when compared with the available high-resolution spectroscopic studies (see Bellazzini et al. 2012). In the analysis, we used only stars satisfying the same membership criterion adopted by these authors and within d < 2r_h from the cluster center (see Section 4); there are 132 in NGC 288, 164 in NGC 5024, 190 in NGC 6218, 357 in NGC 6752, 676 in NGC 6809, and 193 in NGC 7099.

The surface brightness profiles were taken from the Trager et al. (1995) database. They were constructed from generally inhomogeneous data based mainly on the Berkeley Globular Cluster Survey (Djorgovski & King 1984). The surface brightness profile of each cluster has been derived by matching several sets of data obtained with different techniques (aperture photometry on CCD images and photographic plates, photoelectric observations, star counts, etc.). Although these inhomogeneities may represent a drawback of the analysis (see Section 6.1), this database represents the unique collection of wide-field surface brightness profiles covering the entire extent of the target clusters up to their tidal radius.

The target clusters analyzed in this work were selected from the 10 objects in common among the above three databases. We require that (1) the photometric observations reach the magnitude of stars with masses M ∼ 0.15 M_☉ with a completeness level ⩾50%, and (2) the systemic rotation does not exceed 10% of the cluster velocity dispersion. This latter requirement is necessary since rotation affects the determination of the actual velocity dispersion (and consequently the dynamical mass) by providing a source of rotational kinetic energy and introducing an azimuthal variation of star counts and kinematics (Wilson 1975), which will not be included in our models. Six clusters satisfy the above criteria, namely: NGC 288, NGC 5024, NGC 6218, NGC 6752, NGC 6809, and NGC 7099.

Following the prescriptions of Gunn & Griffin (1979), the phase-space distribution of stars is given by the contribution of eight groups of masses with the distribution function:

$$\begin{eqnarray} f(E,L)&=&\sum _{i=1}^{8} k_{i} \exp \left(\! -\frac{A_{i}L^{2}}{2\sigma _{K}^{2}r_{a}^{2}}\right) \left[ \exp \left(\! -\frac{A_{i}E}{2\sigma _{K}^{2}}\right)-1 \right]\nonumber \\ \sum _{i=1}^{8}f_{i}(r,v_{r},v_{t})&=&\sum _{i=1}^{8}k_{i} \exp \left[\! -A_{i}\frac{v_{t}^{2}}{2\sigma _{K}^{2}}\left(\frac{r}{r_{a}}\right)^{2}\right] \nonumber \\ &&\times \left[\exp \left(\! -\frac{A_{i}\big(v_{r}^{2}+v_{t}^{2}+2\psi\big)}{2\sigma _{K}^{2}}\right)-1 \right]\!, \end{eqnarray} \tag{ 1 }$$

where E and L are, respectively, the energy and angular momentum per unit mass, v_r and v_t are the radial and tangential components of the velocity, the effective potential ψ is the difference between the cluster potential ϕ at a given radius r and the potential at the cluster tidal radius ψ ≡ ϕ − ϕ_t, A_i and k_i are scale factors for each mass group, and σ_K is a normalization term. For each cluster, we defined the eight mass bins (all covering equal-mass intervals at different ranges) to sample the entire mass range from the hydrogen burning limit to the mass at the end of the asymptotic giant branch evolution set by comparison between the CMD and the isochrones by Dotter et al. (2007; see Figure 1). The dependence on mass of the coefficients A_i determines the degree of mass segregation of the cluster. Since the half-mass relaxation times of all the six GCs of our sample are significantly smaller than their ages (McLaughlin & van der Marel 2005), we adopted A(i)∝m_i (Michie 1963). The parameter r_a determines the radius at which orbits become more radially biased. As r_a → ∞, models become isotropic. As usual, the above distribution function has been integrated to obtain the number density and the radial and tangential components of the velocity dispersion of each mass group:

$$\begin{eqnarray} n_{i}(r)&=&4\pi \int _{0}^{\sqrt{-2\psi }}\int _{0}^{\sqrt{-2\psi -v_{r}^{2}}} v_{t} f_{i}(r,v_{r},v_{t}) dv_{t}dv_{r},\nonumber \\ \sigma _{r,i}^{2}&=&\frac{4\pi }{\rho _{i}(r)}\int _{0}^{\sqrt{-2\psi }}v_{r}^{2}\int _{0}^{\sqrt{-2\psi -v_{r}^{2}}} v_{t} f_{i}(r,v_{r},v_{t}) dv_{t}dv_{r},\nonumber \\ \sigma _{t,i}^{2}&=&\frac{4\pi }{\rho _{i}(r)}\int _{0}^{\sqrt{-2\psi }}\int _{0}^{\sqrt{-2\psi -v_{r}^{2}}} v_{t}^{3} f_{i}(r,v_{r},v_{t}) dv_{t}dv_{r}.\nonumber \\ \end{eqnarray} \tag{ 2 }$$

The above equations can be written in terms of dimensionless quantities by substituting

$$\begin{eqnarray*} \zeta &=&\frac{v_{r}^2}{2\sigma _{K}^{2}}, \eta =\frac{v_{t}^2}{2\sigma _{K}^{2}},\\ \tilde{\rho }_{i}&=&\frac{m_{i} n_{i}}{\sum _{j} m_{j} n_{0,j}}, \tilde{r}=\frac{r}{r_{c}},\\ W&=& -\frac{\psi }{\sigma _{K}^{2}}, \tilde{r}_{a}=\frac{r_{a}}{r_{c}}, \end{eqnarray*}$$

where ρ₀ is the central cluster density and

$$r_{c}\equiv \left(\frac{9 \sigma _{K}^{2}}{4 \pi G \rho _{0}}\right)^{1/2}$$

is the core radius (King 1966). The potential at each radius is determined by the Poisson equation

$$\begin{equation} \nabla ^{2}\psi =4\pi G \rho. \end{equation} \tag{ 3 }$$

Equations (2) and (3) have been integrated once we have assumed a value of the adimensional potential at the center W₀ outward to the radius r_t at which both density and potential vanish. The total mass of the cluster M is then given by the sum of the masses of all the groups:

$$\begin{eqnarray*} &&M=\sum _{i=1}^{8}m_{i}N_{i} = \frac{9 r_{c} \sigma _{K}^{2}}{G}\sum _{i=1}^{8} \int _{0}^{r_{t}}\tilde{\rho }_{i}\tilde{r}^{2}d\tilde{r}. \end{eqnarray*}$$

As a last step, the above profiles have been projected on the plane of the sky to obtain the surface mass density,

$$\Sigma _{i}(R)=2\int _{R}^{r_{t}}\frac{\rho _{i}r}{\sqrt{r^{2}-R^{2}}}dr,$$

and the line-of-sight velocity dispersion,

$$\sigma _{v,i}^{2}(R)=\frac{1}{\Sigma _{i}}\int _{R}{r_{t}}\frac{\rho _{i}\big[2\sigma _{r,i}^{2}(r^{2}-R^{2})+\sigma _{t,i}^{2}R^{2}\big]}{r\sqrt{r^{2}-R^{2}}}\,dr.$$

In the above models, the shape of the density and velocity dispersion profiles is completely determined by the parameters ($W_{0},\tilde{r}_{a},N_{i}$), while their normalization is set by the pair of parameters (r_c, M).

Zoom In Zoom Out Reset image size

For each choice of the above parameters, a large number (N_tot = 10⁶) of particles have been randomly extracted from a continuous mass function which reproduces the number counts in each mass bin and distributed across the cluster according to their masses by interpolating through the density profiles of the various mass groups. For each particle, the tangential and radial components of its velocity have then been extracted from the distribution function defined by Equation (1) and projected in an arbitrary direction.⁷ The masses of the particles have been converted into F814W and Johnson V magnitude by adopting the mass–luminosity relation of Dotter et al. (2007). Isochrones has been selected adopting the metallicities by Carretta et al. (2009a)⁸ and their magnitudes have been converted from the absolute to the apparent plane adopting the distances and reddening by Paust et al. (2010),⁹ while ages have been chosen to match the turn-off morphology of each cluster. We adopted the reddening coefficients listed by Sirianni et al. (2005). The F606W magnitude has been instead assigned by interpolating on the MS mean ridge line of each cluster in the F814W–(F606W–F814W) CMD. The same procedure has been performed for a fraction f_b of binaries and f_remn of dark remnants whose masses and magnitudes have been assigned following the prescriptions described in Section 3.1. A synthetic horizontal branch (HB) has been simulated for each cluster using the tracks by Dotter et al. (2007), tuning the mean mass and mass dispersion along the HB to reproduce the observed HB morphology. Photometric errors and completeness corrections have been included in the following way: to each simulated particle, we associated an artificial star within 10'' and with input F814W magnitude within 0.25 mag. If the artificial star is recovered, then its corresponding shifts in color and F814W magnitude are then applied to the particle, while the particle is if the associated artificial star is in the output artificial star catalog. To account for Galactic field interlopers, we used the Galaxy model of Robin et al. (2003). A catalog covering an area of 1 deg² around each cluster center has been retrieved. A subsample of stars has been randomly extracted from the entire catalog scaled to the ACS field of view and the V and I Johnson-Cousin magnitudes were converted into the ACS photometric system by means of the transformations of Sirianni et al. (2005). The correction for photometric errors and incompleteness has been applied to this sample as for cluster stars. So, at the end of the above procedure, according to a given choice of input parameters, a mock cluster with masses, projected velocities, and magnitudes accounting for observational effects has been produced.

3.1. Binaries and Dark Remnants

The presence of binaries can influence both the luminosity and the mass estimate of a cluster. Indeed, while the two components of the binary system contribute to the total budget of light and mass in the same way as if they were single stars, their mutual interaction has two important effects: (1) since the two components are not resolved, the binary system can share the same portion of the CMD with single stars of different mass (Malkov & Zinnecker 2001), and (2) from a dynamical point of view, they behave as a single massive object so that, as energy equipartition is established, they present on average lower velocities and a more concentrated distribution with respect to single stars. The above characteristics make them more resistant to the process of mass loss via evaporation and tidal shocks. It is therefore important to account for the presence of a fraction of unresolved binaries to properly model the PDMF of the cluster. To model the population of binaries, we extracted N_bin = 2f_bN_tot stars from the Kroupa (2001) IMF in the mass range 0.1 < M/M_☉ < 7 and paired them randomly. From the resulting sample, we extracted a subsample of pairs according to a "chance of pairing" calculated by imposing that (1) the distribution of the mass ratios f(q) (where q = m₂/m₁ is the ratio between the masses of the secondary and the primary components of each binary system) in the mass range 1 < M/M_☉ < 7 results in the same distribution observed by Fisher et al. (2005) in their sample of spectroscopic binaries in the solar neighborhood,¹⁰ and (2) the distribution of systemic masses must be the same as that of single stars in the common mass range. The accepted pairs were added to a library of binaries while the rejected ones are "broken" and their components used for the next iteration. The procedure was repeated until all stars were included in the library. Then, all binaries whose primary component has a mass exceeding the maximum allowed for an evolving cluster star were removed. The location and velocity of binaries within the cluster were assigned according to their systemic mass (see Section 3). The relative projected velocity of the primary component was added to the projected velocity of each binary according to the formula

$$\begin{eqnarray*} &&v=\frac{2\pi a_{1} \sin i}{P(1-e^{2})^{1/2}}[\cos (\theta +\omega)+e \cos \omega ]\nonumber\\ && {\rm (McConnachie\ \& \ C{\hat{o}}t{\acute{e}}\ \hbox{\protect\hyperlink{apj436802r46}{2010}}),} \end{eqnarray*}$$

where a₁ is the semimajor axis of the primary component, P is the orbital period, e is the eccentricity, i is the inclination angle to the line of sight, θ is the phase from the periastron, and ω is the longitude of the periastron. We followed the prescriptions of Duquennoy & Mayor (1991) for the distribution of periods and eccentricities. The semimajor axis was calculated using Kepler's third law:

$$a_{1}=\frac{1}{1+\frac{m_{1}}{m_{2}}}\left[\frac{P^{2}G(m_{1}+m_{2})}{4\pi ^{2}}\right]^{1/3},$$

where m₁ and m₂ are the masses of the primary and secondary components. We removed all those binaries whose corresponding semimajor axes lie outside the range a_min < a₁ + a₂ < a_max, where a_min is linked to the radius of the secondary component (according to Lee & Nelson 1988) and a_max is set physically by the maximum separation beyond which the binary becomes unbound due to stellar interactions within the cluster (see Equation (1) of Hills 1984). The distribution of the angles (i, θ, ω) was chosen according to the corresponding probability distributions (${\rm Prob}(i)\propto \sin i; {\rm Prob}(\theta)\propto \dot{\theta }^{-1}; {\rm Prob}(\omega)={\rm constant}$). The magnitudes of the two components were then estimated as for single stars and the overall magnitude calculated from the sum of their fluxes in the different passbands. The photometric errors and completeness correction were then applied following the same procedure adopted for single stars (see Section 3).

The number of dark remnants (white dwarfs, neutron stars, and black holes) was estimated, adopting for the precursors a Kroupa (2001) IMF normalized to the number of stars in the most massive stellar bin.¹¹ We then adopted the initial–final mass relation and the retention fractions (for the initial velocity kicks) of Kruijssen (2009). We applied an additional correction to account for the selective removal of low-mass remnants as a result of tidal effects by imposing the same PDMF in the mass range in common with single stars.

The choice of the best model that reproduces the structure and kinematics of each cluster was made by comparing simultaneously four observables measured on the databases and in the simulated mock cluster.

• 1.  
  The surface brightness profile.
• 2.  
  The F814W luminosity function.
• 3.  
  The fraction of binaries.
• 4.  
  The velocity dispersion profile.

The surface brightness profile of the mock cluster was derived by summing the Johnson V fluxes of particles contained within circular concentric annuli of variable width and dividing the resulting flux by the annulus area. The widths of the annuli were chosen to contain at least 20 stars to minimize Poisson noise while maintaining a good resolution in distance. A χ² of the fit has then been calculated.

A comparison of the F814W luminosity function and of the binary fraction was performed simultaneously by comparing the number counts (N^obs_i) in nine regions of the F814W–(F606W–F814W) CMD defined as follows: eight F814W magnitude intervals were defined corresponding to the mass bins adopted in the cluster modeling (see Section 3) and including all stars with colors within three times the photometric error corresponding to their magnitudes. A binary region was defined within four boundaries: the selected stars are those contained in the magnitudes between the loci of binaries with primary star mass m₁ = 0.45 M_☉ (faint boundary) and m₁ = 0.75 M_☉ (bright boundary), and in color between the MS ridge line (blue boundary) and the equal-mass binary sequence (red boundary), both redshifted by three times the photometric error (see Figure 1; see also Sollima et al. 2007). We excluded the innermost 05 of the core-collapsed clusters NGC 6752 and NGC 7099 as in these extremely crowded regions the completeness level drops below 50% even at relatively bright magnitudes.

The velocity dispersion profile of the cluster was compared using the likelihood merit function:

$$\begin{equation} l=\sum _{i=1}^{N}\left(-\frac{(v_{i}-\overline{v})^{2}}{\sigma _{v}(r_{i})^{2}+\delta v_{i}^{2}}-ln\big(\sigma _{v}^{2}+\delta v_{i}^{2}\big)\right), \end{equation} \tag{ 4 }$$

where r_i, v_i, and δv_i are the distance to the cluster center, the radial velocity, and the corresponding uncertainty of the ith observed star, $\overline{v}$ is the mean velocity of the sample, and σ_v(r_i) is the velocity dispersion predicted for RGB stars by the model at r_i. We excluded the outermost region of the cluster at d > 2r_h to minimize the effect of tidal heating that can inflate the velocity dispersion in the outskirts of our clusters (Johnston et al. 1999; Küpper et al. 2010).

We adopted an iterative algorithm that starts from an initial guess of the parameters ($W_{0},r_{c},M,\tilde{r}_{a},N_{i},f_{b}$), generates the mock cluster, and calculates the observables described above. Corrections to the parameters were then calculated using the method described below and the procedure was repeated until convergence. In each iteration, the following steps were performed sequentially.

• 1.  
  • (a)  
    The value of W₀ and r_c providing the best fit to the surface brightness profile is searched by minimizing the χ² statistic. At the end of this step, a new mock cluster is generated leaving the other parameters unchanged.
  • (b)  
    The relative fraction of stars in the eight mass bins N_i is adjusted by multiplying them for corrective terms which are proportional to the ratio between the number of objects in each bin of the observed and of the simulated CMDs:

    $$\protect\begin{eqnarray*} &&N_{i}^{\prime }=N_{i} \left(\frac{N_{i}^{{\rm obs}}}{N_{i}^{{\rm mock}}} \frac{\sum _{i=1}^{8}N_{i}^{{\rm mock}}}{\sum _{i=1}^{8}N_{i}^{{\rm obs}}}\right)^{\eta }, \protect\end{eqnarray*}$$

    where N^obs_i and N^mock_i are the number of stars counted in the ith bin in the observed CMD and in the simulated one, N'_i is the updated value of N_i, and η is a softening parameter, set to 0.5, used to avoid divergence. An updated value for the binary fraction f_b is also set using the same method. Then, a new mock cluster is generated and the procedure is repeated from point (1) until convergence.
• 2.  
  The above steps are repeated for different values of $\tilde{r}_{a}$ and the merit function l is calculated by comparison with the radial velocities to determine the best-fit value of $\tilde{r}_{a}$.

The procedure usually converges after 10–20 iterations independently of the initial guess of parameters. Convergence is set when the various parameters start to fluctuate around the best-fit value with an amplitude (typically <1%) linked to the degree of degeneracy among the various solutions. The uncertainty in the derived parameters was calculated as the standard deviations of such fluctuations quadratically combined with the Poisson noise of star counts (only for the case of the M, N_i, and f_b parameters). The mass of the model can be constrained in two independent ways: (1) by best-fitting the total number counts (so imposing ∑⁸_{i = 1}N_i^mock = ∑⁸_{i = 1}N_i^obs) or (2) by matching the amplitude of the velocity dispersion. The former estimate gives the stellar mass (M_*), the latter the dynamical one (M_dyn).

In Figure 2, the comparison between the considered observables and the prediction of the best-fit model is shown for the case of NGC 288. It is clear that the model reproduces well the surface brightness and velocity dispersion profiles as well as the F814W luminosity functions measured in the two radial samples considered for this cluster. In this cluster, where ancillary observations are available in an external region, the model that best fits the inner sample observables also fits both the shape and the normalization of the luminosity function in the external sample. Although this agreement cannot guarantee the full adequacy of models and observables in the whole sample of clusters, it represents a sanity check supporting the validity of the assumptions made on mass segregation and on the star-count–surface-brightness conversion.

Zoom In Zoom Out Reset image size

The procedure described above allows us to determine a set of parameters for the target clusters, in particular their structural parameters, degree of anisotropy, binary fraction, PDMF, and their stellar and dynamical masses. The best-fit parameters for each cluster are listed in Table 1. For completeness, the M/L_V ratios, calculated from the V magnitudes by van den Bergh et al. (1991), are also listed. The structural parameter derived here is in good agreement with those estimated by previous studies (e.g., McLaughlin & van der Marel 2005). We note that all the clusters of our sample are best fitted by isotropic (or marginally anisotropic) models. However, the relatively large uncertainties in the binned velocity dispersion measures do not allow a meaningful analysis of this issue: the likelihood merit function l (see Section 4) varies by less than 10% over a relatively large range of $\tilde{r}_{a}$ values in all cases.

Table 1. Best-fit Parameters of the Six Target Clusters

Name	W0	rc	rc	$\tilde{r}_{a}$	fb	fcoreb	α	log Mdyn	log M*	Mdyn/LV	M*/LV
		(')	(pc)		(%)	(%)		(M☉)	(M☉)	(M☉/LV, ☉)	(M☉/LV, ☉)
NGC 288	6.2	1.90	5.40	11	3.3 ± 0.2	4.8 ± 0.3	−1.04 ± 0.05	4.83 ± 0.07	4.73 ± 0.1	1.3 ± 0.4	1.0 ± 0.4
NGC 5024	11	0.55	2.99	+∞	4.3 ± 3.0	11.5 ± 7.2	−1.08 ± 0.07	5.60 ± 0.07	5.61 ± 0.1	1.5 ± 0.3	1.5 ± 0.4
NGC 6218	7.7	1.05	1.50	12	3.2 ± 0.2	5.8 ± 0.4	−0.22 ± 0.02	4.90 ± 0.05	4.74 ± 0.1	1.0 ± 0.3	0.7 ± 0.2
NGC 6752	13	0.42	0.51	+∞	0.6 ± 0.1	3.6 ± 0.6	−1.41 ± 0.09	5.45 ± 0.04	5.20 ± 0.1	2.4 ± 0.5	1.3 ± 0.4
NGC 6809	5.5	2.50	3.32	4	2.5 ± 0.2	3.4 ± 0.3	−0.67 ± 0.01	5.19 ± 0.03	5.08 ± 0.1	1.6 ± 0.2	1.3 ± 0.2
NGC 7099	13.5	0.13	0.34	+∞	1.7 ± 0.1	12.0 ± 0.9	−1.14 ± 0.04	5.28 ± 0.05	5.08 ± 0.1	2.1 ± 0.4	1.3 ± 0.3

Download table as:  ASCIITypeset image

The main results on the binary fraction, PDMF, and mass are discussed in the following sections.

5.1. Binary Fractions

The binary fractions listed in Table 1 are remarkably lower than those estimated in previous studies for the clusters in common (Sollima et al. 2007; Milone et al. 2012). This is related to the distribution of mass ratios adopted in this work that is different from those used in previous studies. The above difference can be responsible for the observed discrepancies since in binaries with low-mass ratios the contribution of the secondary component to the total flux is negligible, and they appear in the CMD almost indistinguishable from single stars. Therefore, the adoption of a mass-ratio distribution skewed toward low (high) mass ratios produces a larger (smaller) correction for hidden binaries. To verify the consistency of our analysis with the previous works, we calculated the fraction of binaries in NGC 288 adopting a flat mass-ratio distribution over the entire mass range instead of the assumptions described in Section 3.1. In this case, the overall binary fraction turns out to be 6.5% ± 0.3% corresponding to 9.6% ± 0.4% in the cluster core, in agreement within the errors of the estimates by Sollima et al. (2007) and Milone et al. (2012). It is worth noting that both the PDMF and the mass estimates are almost unaffected by this change: in this case, in fact, the slope of the PDMF varies of less than Δα ≃ −0.07 and both the dynamical and luminous masses remain the same within Δlog M/M_☉ < 0.01. So, although the adopted mass-ratio distribution is somehow arbitrary, it influences only the absolute fraction of binaries having a negligible effect on the relative cluster-to-cluster ranking, on the luminous and dynamical masses, and on the shape of the PDMF.

In Figure 3, the fraction of binaries of our sample of clusters are plotted as a function of the derived cluster dynamical mass. Once we have excluded NGC 5024, which has a large uncertainty in the derived binary fraction, a clear anticorrelation is noticeable, with massive clusters hosting a smaller fraction of binaries. This result has been already observed by Milone et al. (2008, 2012) and Sollima et al. (2010) on a larger sample of clusters.

Zoom In Zoom Out Reset image size

5.2. Present-day Mass Function

In Figure 4, the global PDMFs of the six clusters are shown. The PDMFs can be well represented by single power laws whose slopes are all within the range −0.2 < α < −1.45, significantly less steep than the IMF of Salpeter (1955; α = −2.35) and, in all but the case of NGC 6752, of Kroupa (2001; α = −1.3 at M < 0.5 M_☉ with a mean slope of 〈α〉 = −1.7 between 0.3 < M/M_☉ < 0.8). The extreme value for NGC 6218 (α = −0.22 ± 0.02) confirms the peculiar flatness of the PDMF of this cluster already noted by De Marchi et al. (2006). For the two clusters NGC 288 and NGC 7099, as in Paust et al. (2010), the derived slopes agree within the combined uncertainties of the two works (∼0.2). In the same way, the slopes derived for NGC 6809 and NGC 7099 agree with those estimated for these clusters by Piotto & Zoccali (1999) and Piotto et al. (1997), respectively. We find no signs of a turnover at low masses in all cases, although a convex shape of the PDMF is visible in NGC 5024. It is however worth noting that the mass at which the PDMF is expected to culminate is ∼0.3 M_☉ (Paresce & De Marchi 2000), so this last conclusion depends only on the relative fraction of stars in the lowest-mass bin. However, the estimate in this bin is by far the most uncertain since (1) the mass–luminosity relation in this mass range is largely uncertain, (2) the completeness in this bin is low and uncertain, and (3) some contamination from faint unresolved galaxies falling in the large selection box of this mass bin (see Figure 1) could be present. So, a firm conclusion on the presence of a turnover in the PDMF in these clusters cannot be established.

Zoom In Zoom Out Reset image size

The main result of this analysis is the measure of the dynamical and stellar masses of the target clusters. Previous estimates of the dynamical mass of GCs have been made by Mandushev et al. (1991), Pryor & Meylan (1993), Lane et al. (2010b), and Zocchi et al. (2012), who have, respectively three, five, six, and three clusters in common with the present study. The mean differences and standard deviations between the dynamical masses estimated by these authors and those estimated here are Δlog (M/M_☉) = 0.13(σ = 0.15), −0.03(0.17), 0.10(0.26), and −0.09(0.14),¹² respectively, indicating good agreement within the errors without significant systematics. In the right panel of Figure 5, the derived stellar masses are compared with those reported by McLaughlin & van der Marel (2005) for the four clusters in common using the M/L_V values predicted by the population synthesis model of Bruzual & Charlot (2003) adopting a Kroupa (2001) IMF. The stellar masses derived here are in reasonable agreement within the errors with those derived by these authors when considering individual clusters, although on average our estimates are smaller by ∼ 20%. This is expected since the Kroupa (2001) IMF is steeper than the PDMFs shown in Figure 4, overestimating the number of low-mass stars that contribute significantly to the clusters' mass budget and only marginally to their light.

Zoom In Zoom Out Reset image size

From an inspection of Table 1, it is suddenly evident that there is a discrepancy between the dynamical and stellar masses of five out of six clusters (with the exception of NGC 5024) with the dynamical masses systematically larger than the stellar ones by ∼40%. This is shown in the left panel of Figure 5 where the two estimates are compared.

The observed difference between the dynamical and stellar masses can also be visualized in Figures 6 and 7 where the velocity dispersion profile and the F814W luminosity function of the six target clusters are compared with the prediction of the best-fit models with the appropriate M_* and M_dyn. It is clear that, apart from the case of NGC 5024, models are not able to fit simultaneously both the velocity dispersion profile and the F814W luminosity function. Such a discrepancy is significant from a statistical point of view: a χ² test indicates that the probability that the above difference occurs by chance in the entire sample of clusters is smaller than 3%. Such a discrepancy can be due to either observational biases in the determination of masses or to a physical reason. In the next sections, we will analyze both possibilities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Possible Observational Biases

6.2. Effect of Assumptions

If the discrepancy is not linked to any observational bias, then it can be related to a physical effect that is not (properly) considered in our models. There are many possible mechanisms that may affect the mass determination. In this section, we evaluate the impact of the various inputs, taking the cluster NGC 288 as a tester. In this cluster, the observed mass discrepancy (Δlog M/M_☉ = 0.1) is smaller than the average of the six analyzed GCs (Δlog M/M_☉ = 0.18) so care should be taken in the general interpretation of the derived results.

As discussed in the previous section, the estimate of the dynamical mass is sensible for all the parameters affecting the velocity dispersion. Radial anisotropy can in principle distort the velocity dispersion profile, resulting in a different normalization of mass. As discussed in Section 5, we take this effect into account in our models, but nevertheless the relatively large uncertainties in the velocity estimates do not allow a firm constraint on this quantity. To test this effect, we repeated the fit, imposing the maximum degree of radial anisotropy that still ensures stability (see Nipoti et al. 2002). The largest effect on the derived dynamical mass has been found to be Δlog M_dyn/M_☉ ⩽ 0.03. Binaries can also inflate the velocity dispersion because of the presence of the additional orbital motion of the binary components. Also, in this case, while our code accounts for such an effect, a different assumption of the distribution of mass ratios might produce a difference. However, the test performed in Section 5.1 indicates that while the above assumption can significantly affect the fraction of binaries, the resulting mass is only marginally affected (Δlog M_dyn/M_☉ < 0.01). Another important effect is the heating due to the tidal interaction of the clusters with the Milky Way: the periodic shocks occurring every disk crossing and pericentric passage transfer kinetic energy to the cluster stars, inflating its velocity dispersion, particularly in the outskirts (Gnedin & Ostriker 1997). While this effect is limited in our case since we excluded by fit the outermost portion of the velocity dispersion profile, Küpper et al. (2010) showed that significant heating can also occur within the half-mass radius of loose clusters. To test this hypothesis, we run an N-body simulation to follow the structural and kinematical evolution of a cluster with the orbital and structural characteristics of NGC 288. This is the less dense cluster of our sample and should therefore be among those where tidal shocks are expected to produce the largest effect. We used the last version of the collisionless N-body code gyrfalcON (Dehnen 2000). The cluster has been represented by 61,440 particles distributed according to the best-fit multimass model of this cluster (see Table 1). The mean masses of the particles have been scaled to match the observed cluster mass.¹³ We adopted a leapfrog scheme with a time step of Δt = 1.45 × 10⁴ yr and a softening length of 0.43 pc (following the prescription of Dehnen & Read 2011). Such a relatively large time step and softening length do not affect the accuracy of the simulation because we are interested in the tidal effects that occur on a timescale significantly shorter than the relaxation time (the timescale at which collisions become important). The cluster was launched within the three-component (bulge + disk + halo) static Galactic potential of Johnston et al. (1995) with the energy and angular momentum listed by Dinescu et al. (1999) and its evolution has been followed for 0.852 Gyr (corresponding to four orbital periods). The projected density and velocity dispersion profiles at the end of the simulation are shown in Figure 8. It is clear that the final velocity dispersion deviates from the prediction of the steady-state model at a distance of ∼4' ≃ r_h to the center. By applying the same technique described in Section 4 to the outcome of the simulation, we overestimate the mass of the cluster by Δlog M_dyn/M_☉ ∼ 0.06.

Zoom In Zoom Out Reset image size

The stellar mass can be underestimated because of an improper conversion between the F814W magnitude and mass or because of the presence of non-luminous matter. Considering that the mass discrepancy is of the order of ∼40%, uncertainties on distance, reddening, and age can be excluded as possible causes: to reproduce the above mass difference, one should adopt for these clusters ages <4 Gyr, distance moduli shifted by ∼2.5 mag, or reddening differences of ΔE(B − V) ∼ 0.2, many orders of magnitudes larger than their statistical uncertainties. Uncertainties in the mass–luminosity relation can also be responsible for part of the discrepancy. In this regard, it is worth noting that while different stellar evolution models all agree for masses M > 0.1 M_☉ (Alexander et al. 1997; Baraffe et al. 1997; Dotter et al. 2007; Marigo et al. 2008), significant differences can be present in the least massive bin. Also in this case, however, the effect of such uncertainty is not expected to exceed a few per cent in the estimated mass. A larger effect is instead produced by the assumption on dark remnants. As described in Section 3.1, we followed the prescription by Kruijssen (2009) for the retention fraction of these objects and a Kroupa (2001) IMF for their precursors. To estimate the extent of the impact of this assumption, we considered the extreme case of a 100% retention fraction for all the remnants. In this ideal case, the stellar mass would increase by Δlog M_*/M_☉ ∼ 0.08. The adoption of a different slope for the IMF has an even larger effect: the entire mass discrepancy could be completely removed by assuming a slope of α = −1.5 (instead of −2.3) for masses M > 0.8 M_☉, because the shallower the slope of the IMF is, the larger the relative fraction of massive stars which terminated their evolution and are now in the form of dark remnants.

By comparing the most extensive photometric and spectroscopic surveys with suitable dynamical models, we derived the global binary fraction, the PDMF of six Galactic GCs, as well as their stellar and dynamical masses. The approach adopted here has the advantage of not requiring an a priori assumption of the M/L ratio to estimate the stellar mass (as usually done in previous works; e.g., McLaughlin & van der Marel 2005) which is instead estimated via direct star counts in the CMD.

We confirm the anticorrelation between binary fraction and cluster mass already found in previous studies (Milone et al. 2008, 2012; Sollima et al. 2010) and explained it as the consequence of the increased efficiency of the binary ionization process in massive clusters. Indeed, in massive clusters, both the number of collisions and the mean kinetic energy of stellar encounters are larger, resulting in a higher probability of ionization of multiple systems (Sollima 2008).

The PDMFs of the six clusters are well represented by single power laws, although we cannot exclude a turnover at masses M < 0.2 M_☉, because of the uncertainty in the relative fraction of the least massive stars.

The stellar masses derived for these clusters are on average smaller than those predicted by population synthesis models. This is a consequence of the fact that the IMF adopted by these models is steeper than the actual PDMF of the considered clusters. This has been previously predicted by Kruijssen & Mieske (2009) to be the result of the preferential loss of low-mass stars during the cluster evolution.

We found a discrepancy between the stellar and dynamical masses of five clusters. In particular, dynamical masses are on average ∼40% larger than stellar ones. Such a discrepancy is statistically significant when considering the entire sample of six clusters and could not be detected by previous studies because of the overestimate of the stellar mass mentioned above. Unfortunately, we cannot exclude that such a discrepancy is due to the presence of systematics affecting our estimates. In particular, a possible source of systematics is the inhomogeneity of the surface brightness profile sample of Trager et al. (1995). However, if the detected difference were real, then a number of physical interpretations could be given. After a careful test on the effect of various assumptions and physical processes, we found that a significant spurious increase of the estimated dynamical masses is given by the tidal heating, which can reduce (but not eliminate) the observed discrepancy. On the other hand, the assumption on the distribution and retention of dark remnants has the largest impact on the estimated stellar masses. Part (up to 75%) of the observed mass difference can be indeed explained by assuming that most of the remnants, including the most massive neutron stars and black holes, are retained. This could be the case if the cluster mass at the moment of the supernovae II explosion were large enough to trap the remnants, the kinetic energy of which has been increased by the velocity kick following the explosion. The mass necessary to retain the majority of massive remnants is of ∼10⁷ M_☉ (see Figure 1 of Kruijssen 2009). Such a large mass is comparable with the initial mass predicted by D'Ercole et al. (2008) to ignite the self-enrichment process observed in most GCs (see Carretta et al. 2009b). Moreover, massive remnants are expected to be more resistant to tidal stripping since they tend to sink toward the cluster core as a result of mass segregation. It is however worth noting that a retention fraction larger than 80% would require initial masses and a quite extreme early mass-loss history, which has indeed been proposed in theoretical works to be the result of severe tidal shocking in the primordial environment of GCs (Elmegreen 2010; Kruijssen et al. 2012). An even larger effect is provided by the slope of the high-mass end of the IMF of the precursors which can drastically increase the amount of non-luminous mass. In this regard, the recent analysis by Marks et al. (2012) suggests that the slope of the IMF in the range M > M_☉ could actually be shallower than the canonical Salpeter (1955) and Kroupa (2001) values with a dependence on density and metallicity. In this case, a slope α ≃ −1.5 (instead of α = −2.3 by Kroupa 2001) could account for the entire difference between the estimated stellar and dynamical masses. Such a variation would have strong implications on the evolution of the GC MF: In fact, while the largest mass contribution would be due to the increased number of white dwarfs, the consistent increase of the fraction of massive remnants (black holes and neutron stars) makes these objects the main drivers of the MF evolution. It is not clear whether such an IMF could actually evolve in the PDMF observed in these clusters.

Another possibility is that GCs contain a modest fraction of non-baryonic dark matter. In this case, it is possible that what we see now is the remnant of a larger halo lost during the interaction with the Milky Way (see Saitoh et al. 2006). Some amount of dark matter is also expected if some GCs are the remnants of past accretion events, as previously suggested by Freeman & Bland-Hawthorn (2002). However, given the involved uncertainties, all these hypotheses are merely speculative. Further studies will be needed to verify the above possibilities.

Based on observations made with the NASA/ESA Hubble Space Telescope, which is operated by the association of Universities for Research in Astronomy, Inc., under the NASA contract NAS 5-26555, under programs GO-10775 (PI: Sarajedini) and GO-12193 (PI: Lee). We thank the anonymous referee for helpful comments and suggestions. A.S. acknowledges the support of INAF through the 2010 postdoctoral fellowship grant. M.B. acknowledges the financial support of INAF through the PRIN-INAF 2009 grant assigned to the project Formation and evolution of massive star clusters, PI: R. Gratton. J.-W.L. acknowledges financial support from the Basic Science Research Program (grant No. 2010-0024954) and the Center for Galaxy Evolution Research through the National Research Foundation of Korea.

Facilities: HST (ACS) - Hubble Space Telescope satellite, AAT (AAOmega) - Anglo-Australian Telescope