NEW LIMITS ON AN INTERMEDIATE-MASS BLACK HOLE IN OMEGA CENTAURI. II. DYNAMICAL MODELS^*

Imaging observations of a cluster yield an estimate of μ_obs,X(R), the radial profile of the surface brightness in a photometric band X, averaged over concentric circles of radius R around the cluster center. This quantity is related to the projected intensity I_obs(R) according to

$$\begin{equation} \mu _X [{\rm mag} \; {\rm arcsec}^{-2}] = -2.5 \log I_{\rm obs} [L_{\odot } \; {\rm pc}^{-2}] + 21.572 + M_{\odot, {X}}, \end{equation} \tag{ 1 }$$

where M_☉,X is the absolute magnitude of the Sun in the given photometric band. Throughout this paper, units in which quantities are expressed are given in square brackets where relevant. In the presence of A_X magnitudes of foreground extinction, the extinction-corrected profile is given by

$$\begin{equation} I(R) = 10^{0.4 A_{X}} \: I_{\rm obs}(R). \end{equation} \tag{ 2 }$$

For a spherical cluster, I(R) is the LOS integration over the three-dimensional luminosity density j(r). This integration can be inverted through an Abel transform (e.g., Binney & Tremaine 1987) to yield

$$\begin{equation} j(r) = -{1\over \pi } \int _{r}^{\infty } \left[ { {dI}\over {dR} } \right] (R) { {dR} \over {\sqrt{R^2-r^2}} }. \end{equation} \tag{ 3 }$$

To express this quantity in units of L_☉ pc⁻³, one needs to specify the cluster distance D. Sizes in arcsec and pc are related according to

$$\begin{equation} R[{\rm pc}] = R[{\rm arcsec}] \: \pi \: D[{\rm kpc}] / 648. \end{equation} \tag{ 4 }$$

The mass density of stellar objects and their remnants is

$$\begin{equation} \rho (r) = [M/L](r) \: j(r), \end{equation} \tag{ 5 }$$

where [M/L](r) is the average mass-to-light ratio at each radius. Gradients in [M/L](r) can be thought of as arising from two separate ingredients: (1) mass segregation among luminous stars and (2) mass segregation of the heavy compact stellar remnants (heavy white dwarfs, neutron stars, and stellar-mass black holes) with respect to the luminous stars. The first ingredient generally causes a shallow gradient in [M/L](r) over the scale of the entire cluster (e.g., Gill et al. 2008). This has little effect on the model predictions in the area surrounding the very center, where we wish to search for the signature of an IMBH. By contrast, the second ingredient can in some clusters yield a strong increase in [M/L](r) near the very center (e.g., Dull et al. 1997; Baumgardt et al. 2003a) that mimics the presence of a dark sub-cluster.

Our modeling approach allows for exploration of models with arbitrary [M/L](r). However, we have found it useful and sufficient in the present context to restrict attention to a parameterized subset. Based on the preceding considerations, we write the mass density as

$$\begin{equation} \rho (r) \equiv \tilde{\rho }(r) + \rho _{\rm dark}, \qquad \tilde{\rho }(r) \equiv \Upsilon j(r). \end{equation} \tag{ 6 }$$

Here $\tilde{\rho }(r)$ is the mass density obtained under the assumption of a constant mass-to-light ratio ϒ. The excess ρ_dark is treated mathematically as a separate component. For ρ_dark, we have found it convenient to use the parameterization

$$\begin{equation} \rho _{\rm dark} = \big(3 M_{\rm dark} a_{\rm dark}^2 / 4 \pi\big) / \big(r^2 + a_{\rm dark}^2\big)^{5/2}, \end{equation} \tag{ 7 }$$

corresponding to a Plummer density distribution. Here M_dark is the total excess mass in dark objects and a_dark is the characteristic scale length of their spatial distribution. The quantities ϒ, M_dark, and a_dark provide a convenient way for exploring a range of plausible [M/L](r) functions using only three parameters.

By assuming that ϒ is independent of r, we neglect mass segregation among the luminous stars. For ω Centauri, there are in fact several pieces of evidence that it has not yet undergone complete mass segregation, consistent with its long half-mass relaxation time of ∼10^9.96±0.03 yr (McLaughlin & van der Marel 2005). First, Merritt et al. (1997) constructed non-parametric dynamical models for ω Centauri. They found that at the large radii sampled by their data, ρ(r) ∝ j(r) to within the error bars. Second, vdV06 constructed more general dynamical models for a wider range of data, and they too found that at all radii sampled by their study [M/L](r) is consistent with being constant. Third, Anderson (2002) studied luminosity functions (LFs) at different radii in ω Centauri using data obtained with HST. He found that variations in the LFs are less than would be expected if mass segregation had proceeded to the point of energy equipartition, and fourth, in Paper I we studied the proper motion dispersion of stars along the main sequence. This showed that lower mass stars do move faster, but by a smaller amount than would be expected for energy equipartition.

The total gravitational potential of the system is

$$\begin{equation} \Phi _{\rm tot} = \tilde{\Phi }(r) + \Phi _{\rm dark}(r) - G M_{\rm BH}/ r, \end{equation} \tag{ 8 }$$

where G is the gravitational constant and M_BH is the mass of any central IMBH. The gravitational potential corresponding to $\tilde{\rho }(r)$ is

$$\begin{equation} \tilde{\Phi }(r) = -4 \pi G \left[ {1\over r} \int _{0}^{r} \tilde{\rho }(r^{\prime }) r^{\prime 2} \: dr^{\prime } + \int _{r}^{\infty } \tilde{\rho }(r^{\prime }) r^{\prime } \: dr^{\prime } \right]. \end{equation} \tag{ 9 }$$

The Plummer potential corresponding to ρ_dark is

$$\begin{equation} \Phi _{\rm dark}(r) = - G M_{\rm dark} / (r^2 + a_{\rm dark}^2)^{1/2}. \end{equation} \tag{ 10 }$$

The potential Φ_tot could in principle be expanded to also include an extended halo of particulate dark matter, but that is more relevant for galaxies. There is no evidence that star clusters are embedded in dark halos. For the specific case of ω Centauri, vdV06 showed that there is no need for dark matter out to at least two half-light radii.

The Jeans equation for hydrostatic equilibrium in a spherical stellar system is (e.g., Binney & Tremaine 1987)

$$\begin{equation} {{d\big(\nu \overline{v_r^2}\big)}\over {dr}} + 2 {{\beta \nu \overline{v_r^2}} \over {r}} = -\nu {{d\Phi _{\rm tot}}\over {dr}}, \end{equation} \tag{ 11 }$$

where ν(r) is the number density of the stars under consideration. Because of the spherical symmetry, $\overline{v_{\theta }^2}= \overline{v_{\phi }^2}$. The function $\beta (r) \equiv 1 - \overline{v_{\rm tan}^2}/\overline{v_r^2}$ measures the anisotropy in the second velocity moments, where $\overline{v_{\rm tan}^2}\equiv \overline{v_{\theta }^2}= \overline{v_{\phi }^2}$. Models with β = 0 are isotropic, models with β < 0 are tangentially anisotropic, and models with 0 < β ⩽ 1 are radially anisotropic. The models may in principle be rotating if $\overline{v_{\phi }} \not= 0$. However, for the purposes of the present paper, we do not need to address how the second azimuthal velocity moment $\overline{v_{\phi }^2}\equiv \sigma _{\phi }^2 + \overline{v_{\phi }}^2$ splits into separate contributions from random motions and mean streaming.

The measured kinematics in a cluster are generally representative for only a subset of the stars (e.g., stars in that magnitude range for which kinematics can be measured). The number density of this subset is

$$\begin{equation} \nu (r) = j(r) f(r) / \lambda (r), \end{equation} \tag{ 12 }$$

where λ(r) is the average luminosity of a star in the subset, and f(r) is the fraction of the total light that the subset contributes at each radius. Compact stellar remnants play no role in this equation, since they contribute neither to ν(r) (their kinematics are not probed) nor to j(r) (they emit negligible amounts of light). However, f(r) and λ(r) can vary slowly as a function of radius due to mass segregation among the luminous stars. The quantity λ(r) may also be affected by observational selection effects that vary as a function of radius. Nonetheless, variations in λ(r) will generally be smaller than those in f(r), as long as kinematics are measured for stars of similar brightness at each radius. Both f(r) and λ(r) are only expected to have a shallow gradient over the scale of the entire cluster. These gradients are much smaller than those in j(r), which can vary by orders of magnitude. Moreover, in the present paper, we are mostly concerned with the small area surrounding the very center. We therefore assume that

$$\begin{equation} f(r) / \lambda (r) = {\rm constant}, \end{equation} \tag{ 13 }$$

which implies that, as before, we neglect mass segregation among the luminous stars.

It follows from the preceding that one can substitute j(r) for ν(r) in the Jeans equation (11), without actually needing to know the value of the constant in Equation (13). The Jeans equation is then a linear, first-order differential equation for $j \overline{v_r^2}$ with solution

$$\begin{equation} j \overline{v_r^2}(r) = \int _{r}^{\infty } j\left(r^{\prime }\right) \left [ {{d\Phi _{\rm tot}}\over {dr}} \right ](r^{\prime }) \: \exp \left [ 2 \int _{r}^{r^{\prime }} {{\beta (r^{\prime \prime }) \: dr^{\prime \prime }}\over {r^{\prime \prime }}} \right ] \: dr^{\prime }. \end{equation} \tag{ 14 }$$

This is a special case of a more general solution for certain classes of axisymmetric systems (Bacon et al. 1983). We consider anisotropy functions of the form

$$\begin{equation} \beta (r) = (\beta _{\infty } r^2 + \beta _0 r_a^2) / (r^2 + r_a^2), \end{equation} \tag{ 15 }$$

where r_a > 0 is a characteristic anisotropy radius, β₀ ⩽ 1 is the central value of β(r), and β_∞ ⩽ 1 is the asymptotic value at large radii. If β₀ = β_∞ then β(r) is constant. With this β(r), Equation (14) reduces to

$$\begin{eqnarray} j \overline{v_r^2} (r) &=& \int _{r}^{\infty } j(r^{\prime }) \left[ {{d\Phi _{\rm tot}}\over {dr}} \right](r^{\prime }) \: \nonumber \\ & &\times \left[ {{r^{\prime 2} + r_a^2} \over {r^2 + r_a^2}} \right]^{\beta _{\infty }} \left[ \left ({{r^{\prime 2}}\over {r^2}} \right) \Big / \left ({{r^{\prime 2} + r_a^2} \over {r^2 + r_a^2}} \right) \right]^{\beta _0} \!\!dr^{\prime },\quad \end{eqnarray} \tag{ 16 }$$

which is always positive.

The integral in Equation (16) can be used with the equations given for j(r) and Φ_tot(r) in Sections 2.1 and 2.2 to allow numerical evaluation of $\overline{v_r^2}(r)$. This requires knowledge of the observed surface brightness profile μ_obs,X(R) and the known extinction A_X. It also requires choices for the scalar model quantities D, ϒ, M_BH, M_dark, a_dark, β₀, β_∞, and r_a, which can all be chosen to optimize the fit to the available kinematical data.

To calculate quantities that are projected along the LOS, consider an arbitrary point P in the cluster. We adopt a Cartesian coordinate system (pmr, los, pmt) with three orthogonal axis and its origin at the cluster center C. The pmr-axis (proper motion radial) lies in the plane of the sky and points from the projected position of C to the projected position of P. The los-axis (LOS) points from the observer to the point C. The pmt-axis (proper motion tangential) lies in the plane of the sky, and is orthogonal to the pmr- and los-axes in a right-handed sense. The cylindrical polar coordinate system associated with (pmr, los, pmt) is (R, ϕ, z), and the spherical polar coordinate system is (r, θ, ϕ), with θ ∈ [ − π/2, π/2], and θ = π/2 corresponding to the pmt axis. The velocities in the (pmr, los, pmt) system satisfy

$$\begin{eqnarray} v_{\rm pmr}& = & v_r \cos \phi - v_{\phi } \sin \phi, \nonumber \\ v_{\rm los}& = & v_r \sin \phi + v_{\phi } \cos \phi, \nonumber \\ v_{\rm pmt}& = & v_{\theta }. \end{eqnarray} \tag{ 17 }$$

Since $\overline{v_r v_{\phi }} = 0$ in a spherical system, the second moments are

$$\begin{eqnarray} \overline{v_{\rm pmr}^2} (r,\phi) & = & \overline{v_r^2}(r) [1 - \beta (r) \sin ^2 \phi ], \nonumber \\ \overline{v_{\rm los}^2} (r,\phi) & = & \overline{v_r^2}(r) [1 - \beta (r) \cos ^2 \phi ], \nonumber \\ \overline{v_{\rm pmt}^2} (r,\phi) & = & \overline{v_r^2}(r) [1 - \beta (r)], \end{eqnarray} \tag{ 18 }$$

where we have used that $\overline{v_{\phi }^2}= \overline{v_{\theta }^2}= \overline{v_r^2}[1 - \beta (r)]$.

The projected second moments satisfy

$$\begin{eqnarray} I \big\langle v_{\rm pmr}^2 \big\rangle (R) & = & \int _R^{\infty } {{2 r j \overline{v_r^2}(r) [1 - \beta (r) + \beta (r) (R/r)^2) ] \: dr}\over {\sqrt{r^2-R^2}}}, \nonumber \\ I \big\langle v_{\rm los}^2 \big\rangle (R) & = & \int _R^{\infty } {{2 r j \overline{v_r^2}(r) [1 - \beta (r) (R/r)^2 ] \: dr}\over {\sqrt{r^2-R^2}}}, \nonumber \\ I \big\langle v_{\rm pmt}^2 \big\rangle (R) & = & \int _R^{\infty } {{2 r j \overline{v_r^2}(r) [1 - \beta (r)] \: dr}\over {\sqrt{r^2-R^2}}}. \end{eqnarray} \tag{ 19 }$$

The derivations of these equations use the facts that cos²ϕ = (R/r)² and $d\;{\rm los} = d \sqrt{r^2-R^2} = r \: dr /\sqrt{r^2-R^2}$. The factor of 2 arises from the identical contributions from points with los <0 and los >0. The angle brackets around the squared velocities on the left-hand side indicate the combined effect of averaging over velocity space and density-weighted averaging along the LOS. Formally speaking, the weighting is with number density if one is modeling discrete measurements, and with luminosity density if one is modeling integrated-light measurements. However, in the present paper, we are assuming that ν(r) ∝ j(r), so the two cases are equivalent. The integrals in Equation (19) can be evaluated numerically once the Jeans equation has been solved.

The normalized projected velocity distributions ${\cal L}(v,R)$ at projected radius R are generally different for the los, pmr, and pmt directions. The Jeans equation suffices to calculate the second velocity moments of these distributions. To calculate the distributions themselves, it is necessary to know the full phase-space distribution function. This is not generally straightforward (or unique) for anisotropic models. However, for an isotropic model, the distribution function depends only on the energy, and it can be uniquely calculated for given j(r) (being proportional to both ν(r) and $\tilde{\rho }(r)$ under the present assumptions) and Φ_tot(r), with the help of Eddington's equation (Binney & Tremaine 1987). Because of the isotropy, the distribution ${\cal L}_{\rm iso}(v,R)$ is independent of the choice of the velocity direction (los, pmr, or pmt) and satisfies, in analogy with Equation (19),

$$\begin{equation} I(R) {\cal L}_{\rm iso}(v,R) = \int _R^{r_{\rm max}(v)} {{2 r j(r) {\cal L}_{\rm loc}(v,r) \: dr}\over {\sqrt{r^2-R^2}}}. \end{equation} \tag{ 20 }$$

Here r_max(v) is defined to satisfy $\Psi (r_{\rm max}) \equiv {1\over 2}v^2$, with Ψ ≡ −Φ_tot. The quantity ${\cal L}_{\rm loc}(v,r)$ is the local velocity distribution (before projection along the LOS) at radius r, which can be written with the help of Eddington's equation as (van der Marel 1994a)

$$\begin{equation} j(r) {\cal L}_{\rm loc}(v,r) = {1\over {\pi \sqrt{2}}} \int _0^{\Psi (r)-{1\over 2}v^2} \left [ {{d j}\over {d\Psi }} \right ]_{\Psi ^{\prime }}\!\! {{d\Psi ^{\prime }}\over {\sqrt{\Psi (r)-{1\over 2}v^2-\Psi ^{\prime }}}}. \end{equation} \tag{ 21 }$$

The integrals in Equations (20) and (21) can be used with the equations given for j(r) and Φ_tot(r) in Sections 2.1 and 2.2 to allow numerical evaluation of ${\cal L}_{\rm iso}(v,R)$. The result is a symmetric function of v (i.e., ${\cal L}_{\rm iso}(v,R) = {\cal L}_{\rm iso}(-v,R)$). This remains true if the distribution function is given a non-zero odd part in the angular momentum L_z around some axis (corresponding to a rotating model), provided that ${\cal L}_{\rm iso}(v,R)$ is then interpreted as the distribution averaged over a circle of radius R.

Equations (19) and (20) define velocity moments and velocity distributions at a fixed position on the projected plane of the sky. However, to determine these quantities observationally, it is necessary to gather observations over a finite area. In integrated-light measurements this happens naturally through the use of, e.g., apertures, fibers, or pixels along a long slit of finite width. In discrete measurements of individual stars this happens by binning the data during analysis. In either case, the corresponding model predictions are obtained as a density-weighted integral over the corresponding area in the projected (x, y)-plane of the sky. For example, the second LOS velocity moment over an aperture is

$$\begin{equation} \big\langle v_{\rm los}^2 \big\rangle _{\rm ap} = \int _{\rm ap} I\big\langle v_{\rm los}^2 \big\rangle (R) \: dx \: dy {\Big /} \int _{\rm ap} I(R) \: dx \: dy. \end{equation} \tag{ 22 }$$

Similar equations apply to the other velocity moments and to the velocity distribution ${\cal L}_{\rm iso}(v,R)$. In integrated-light measurements, one can also take the observational point-spread function (PSF) into account by convolving both I (x, y) and I〈v²_los〉(x, y) with the PSF before evaluation of the integrals in Equation (22).

In the following, we write

$$\begin{equation} {\bar{\sigma }}_{\rm los}\equiv \big\langle v_{\rm los}^2 \big\rangle _{\rm ap}^{1/2}, \quad {\bar{\sigma }}_{\rm pmr}\equiv \big\langle v_{\rm pmr}^2 \big\rangle _{\rm ap}^{1/2}, \quad {\bar{\sigma }}_{\rm pmt}\equiv \big\langle v_{\rm pmt}^2 \big\rangle _{\rm ap}^{1/2}. \end{equation} \tag{ 23 }$$

It is also useful to consider the quantity

$$\begin{equation} {\bar{\sigma }}_{\rm pm}\equiv \big\langle \big(v_{\rm pmr}^2 + v_{\rm pmt}^2\big) / 2 \big\rangle _{\rm ap}^{1/2} = \sqrt{\big({\bar{\sigma }}_{\rm pmr}^2 + {\bar{\sigma }}_{\rm pmt}^2\big)/2}, \end{equation} \tag{ 24 }$$

which corresponds to the one-dimensional rms velocity averaged over all directions in the plane of the sky. Note that velocities in the plane of the sky are not directly accessible observationally. Instead of ${\bar{\sigma }}_{\rm pmr}$ and ${\bar{\sigma }}_{\rm pmt}$, we measure proper motions ${\bar{\Sigma }}_{\rm pmr}$ and ${\bar{\Sigma }}_{\rm pmt}$ that satisfy

$$\begin{equation} {\bar{\Sigma }}_{\ldots } [{\rm mas}\;{\rm yr}^{-1}] = {\bar{\sigma }}_{\ldots } [{\rm km}\;{\rm s}^{-1}]/(4.7404 \: D[{\rm kpc}]). \end{equation} \tag{ 25 }$$

Transformation to velocities requires an assumption about the cluster distance D.

The overlying bar in the preceding equations is used to differentiate the rms velocity ${\bar{\sigma }}$ from the velocity dispersion σ. These quantities are equal only if the mean streaming is zero (which is always true in the pmr direction).

In general, we will want to model LOS data points ${\bar{\sigma }}_{{\rm los},i} \pm \Delta {\bar{\sigma }}_{{\rm los},i}$ available for i = 1, ..., N, as well as proper motion data points ${\bar{\Sigma }}_{{\rm pm\ldots },j} \pm \Delta {\bar{\Sigma }}_{{\rm pm\ldots },j}$ available for j = 1, ..., M (with each proper motion data point referring to either the pmr or the pmt direction). To calculate model predictions, we use the known surface brightness profile μ_obs,X(R) and the extinction A_X. The questions are then how to assess the fit to the data for given model parameters, and how to determine the parameters that provide adequate fits.

2.7.1. Distance and Mass-to-Light Ratio

Suppose that we have obtained model predictions ${\hat{\sigma }}_{{\rm los},i}$ and ${\hat{\Sigma }}_{{\rm pm\ldots },j}$ for the data points, using trial values ${\hat{D}}$ and ${\hat{\Upsilon }}$ for the distance and mass-to-light ratio, respectively. The fit can then be improved by scaling of the model predictions, which corresponds to rescaling of the mass and distance scales of the model. Let F_los and F_pm be factors that scale the model predictions in the LOS and proper motion directions, respectively. The best-fitting values of these scale factors and their error bars can be calculated using the statistics

$$\begin{eqnarray} \chi ^2_{\rm los} & \equiv & \sum _i [ ({\bar{\sigma }}_{{\rm los},i} - F_{\rm los} {\hat{\sigma }}_{{\rm los},i}) / \Delta {\bar{\sigma }}_{{\rm los},i} ]^2, \nonumber \\ \chi ^2_{\rm pm} & \equiv & \sum _j [ ({\bar{\Sigma }}_{{\rm pm\ldots },j} - F_{\rm pm} {\hat{\Sigma }}_{{\rm pm\ldots },j}) / \Delta {\bar{\Sigma }}_{{\rm pm\ldots },j} ]^2. \end{eqnarray} \tag{ 26 }$$

The best fit in each case is the value with the minimum χ² (i.e., dχ²/dF_... = 0), and the 68.3% confidence region corresponds to Δχ² = χ² − χ²_min ⩽ 1. The inferred scalings imply that the distance and mass-to-light ratio are optimized for $D_{\rm fit} = {\hat{D}} F_{\rm los} / F_{\rm pm}$ and $\Upsilon _{\rm fit} = {\hat{\Upsilon }} F_{\rm los} F_{\rm pm}$. These fits apply at the given values of the other model parameters, some of which are also affected by the scaling. The model parameters M_BH and M_dark are proportional to F³_los/F_pm, and the model parameter a_dark is proportional to F_los/F_pm.

The error bars in D_fit and ϒ_fit (at fixed values of the other parameters) follow from propagation of the uncertainties in F_los and F_pm. The errors in F_los and F_pm are uncorrelated, but the errors in D_fit and ${\hat{\Upsilon }}$ can be highly correlated or anti-correlated. We will not generally quote the correlation coefficient in this paper, but nonetheless, this should be kept in mind where necessary.

The use of scalings as described here provides a means of reducing the parameter space of the models by two. Also, it directly yields the best-fit distance and mass-to-light ratio if the remaining parameters can be determined or eliminated independently. The parameters M_BH, M_dark, and a_dark can be determined as described in Section 2.7.2. They primarily affect the data points near the cluster center. So they can also be effectively eliminated by restricting the definition of the χ² quantities in Equation (26) to the data points at large radius. The anisotropy profile β(r) can be determined as described in Section 2.7.3.

2.7.2. Dark Mass

2.7.3. Anisotropy Profile

Binney & Mamon (1982) showed that there are many different combinations of β(r) and the enclosed mass M(<r) in a spherical system that can produce the same observed projected rms LOS velocity profile ${\bar{\sigma }}_{\rm los}(R)$. The key to constraining the mass distribution, and hence the possible presence of an IMBH, is therefore to constrain the anisotropy profile β(r). This is not possible when only rms LOS velocities are available. However, β(r) can be constrained when information is available on deviations of the LOS velocity distributions (LOSVDs) from a Gaussian, as a function of projected position on the sky (van der Marel & Franx 1993; Gerhard 1993). This information is fairly indirect, and to exploit it one must construct complicated numerical models that retrieve the full phase-space distribution function of the stars (van der Marel et al. 1998; Cretton et al. 1999; Gebhardt et al. 2000; Valluri et al. 2004; Thomas et al. 2004; van den Bosch et al. 2008). This has been explored with considerable success in the context of studies of early-type galaxies (e.g., Cappellari et al. 2007).

In the present context, high-quality proper motion data are available. It is then considerably more straightforward to constrain β(r). This was demonstrated by Leonard & Merritt (1989), who were the first to use proper motion data to examine in detail the mass distribution in a star cluster. The function $[{\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}] (R)$ is a direct measure of anisotropy as seen in the plane of the sky. This ratio is independent of either the distance or the mass-to-light ratio scaling of the model. There is only a weak dependence on the presence of a dark mass (and this only affects the very central region). Therefore, the one-dimensional observed function $[{\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}] (R)$ tightly constrains the one-dimensional intrinsic function β(r). In fact, for a spherical system, there is a unique relation between the two (Leonard & Merritt 1989). This is the primary reason why in the present context it is possible to obtain robust results based on something as simple as the spherical Jeans equation, while much more sophisticated techniques have become standard practice for galaxies.

To build some intuitive understanding of the meaning of $[{\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}] (R)$, consider the case in which β(r) is independent of radius and $j \overline{v_r^2}(r) \propto r^{-\xi }$. The relevant integrals in Equation (19) can then be expressed analytically in terms of beta functions. Their ratio is independent of radius R and simplifies to

$$\begin{equation} [{\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}] = \xi (1-\beta) / (\xi -\beta). \end{equation} \tag{ 27 }$$

For small β (i.e., a system not too far from isotropy), we can use a first-order Taylor approximation, which yields

$$\begin{equation} \beta _{\rm pm} (R) = \beta \big(1 - {1\over \xi }\big), \end{equation} \tag{ 28 }$$

where we have defined the function

$$\begin{equation} \beta _{\rm pm} (R) \equiv 1 - \lbrace [{\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}](R) \rbrace ^2, \end{equation} \tag{ 29 }$$

in analogy with the definition of β(r). This implies that the deviation of $[{\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}]$ from unity is a factor $(1 - {1\over \xi })$ of the deviation of $[\overline{v_{\rm tan}^2}/\overline{v_r^2}]^{1/2}$ from unity. (This uses the Taylor approximation (1 − β)^1/2 ≈ 1 − (β/2) and similarly for β_pm.) Hence, $[{\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}(R)]$ is a "diluted" measure of the anisotropy that is present intrinsically, with the dilution factor $(1 - {1\over \xi })$.

The approximations made above cannot be used in the core of a cluster, where ξ = 0. The assumption that the system remains scale free all along the LOS is then not accurate. However, in the outer parts of a cluster the approximation is a reasonably accurate guide, with the value of ξ depending on the exact cluster structure. For example, an isothermal sphere has ξ = 2 at large radii, while a Plummer model has ξ = 6. The corresponding dilution factors $(1 - {1\over \xi })$ are 0.50 and 0.83, respectively. Real clusters typically fall between these extremes. More extensive analytical calculations of β_pm(R) are presented in, e.g., Leonard & Merritt (1989) and Wybo & Dejonghe (1995). Either way, analytical calculations are mostly useful here for illustration. In practice, we calculate the integrals in Equation (19) numerically.

The extent to which the velocity distribution of the stars (in either the los, pmr, or pmt direction) differs from a Gaussian provides secondary information on the amount of velocity anisotropy in the system. The modeling technique used here is insufficient to exploit this information, since we are not calculating actual phase-space distribution functions for anisotropic models (although this would in principle be possible, e.g., Dejonghe 1987; Gerhard 1993; van der Marel et al. 2000). Nonetheless, we can observationally characterize the Gauss–Hermite moments as a function of projected distance from the cluster center. Also, we can calculate the predicted velocity distributions ${\cal L}_{\rm iso}(v,R)$ for an isotropic model, and from this the corresponding Gauss–Hermite moments. Comparison of the observed and predicted quantities provides an independent assessment of possible deviations from isotropy.

To apply the formalism described above, we augmented a software implementation used previously in, e.g., van der Marel (1994a). The code uses a logarithmically spaced grid, from very small to very large radii (so that a negligible fraction of the cluster mass resides outside the grid). Relevant quantities are tabulated on this grid. Integrals are evaluated through Romberg quadrature (e.g., Press et al. 1992) to high numerical precision. Tabulated quantities are interpolated using cubic splines for use in subsequent integrations. This avoids the need for calculation of higher-dimensional integrals. The derivatives dI/dR and dj/dΨ in Equations (3) and (21), respectively, are calculated using finite differences. The derivative dΦ_tot/dr in Equation (16) can be manipulated analytically so that numerical differentiation is not needed. The accuracy of the code was verified by reproducing analytically known results for special potential-density pairs (e.g., Plummer models).

Trager et al. compiled literature data from various sources. We rejected from their compilation the data points to which they did not assign the maximum weight (unity in their notation). This leaves only the most reliable literature sources for ω Centauri, namely the integrated-light photo-electric measurements of Gascoigne & Burr (1956) and Da Costa (1979), and the number density measurements of King et al. (1968). All these authors presented averages for concentric circular annuli. Da Costa also presented separate measurements using a drift scan technique, which were subsequently calibrated to estimate the profile along concentric circular annuli. King et al. presented results for two separate plates from Boyden Observatory. We follow Trager et al. by referring to the different data sets with the abbreviations GBANN, DANN, DSCAN, ADH-1792, and ADH-1846; the latter two abbreviations refer to the IDs of the photographic plates analyzed by King et al. For our analysis, we used the values provided by Trager et al., who calibrated the published measurements to a V-band surface brightness scale with a common zero point. We assigned to each data point an error bar of 0.142 mag, based on the procedure and analysis of McLaughlin & van der Marel (2005, see their Table 6). This corresponds to the scatter in the data with respect to a smooth curve.

Paper I presented number density measurements for concentric circular annuli based on HST data. Considerable effort was spent to accurately determine the cluster center, and to correct the results for the effects of incompleteness. Error bars were determined from the Poisson statistics of the counted stars. Paper I discussed number density profiles for various magnitude bins and showed that there is no strong dependence on magnitude. For the analysis presented here we have used the profile for all stars in the HST images with B-band instrumental magnitude <−10 (this is defined in Paper I and corresponds to stars measured with signal-to-noise ratio (S/N) ≳100). The choice of magnitude range was motivated by the desire to maximize the number of stars (i.e., to minimize the error bars on the density profile) while maintaining a similar mean magnitude (∼−12) as for the sample of stars for which high-quality proper motions are available from Paper I.

To extend the profile from the Trager et al. compilation to the smaller radii accessible with the HST data, it is necessary to calibrate the projected number density ν_proj(R) to a projected surface brightness μ_obs,V(R). It follows from Equations (1), (4), (12), and (13) that

$$\begin{equation} \mu _{{\rm obs,}V} [{\rm mag} \; {\rm arcsec}^{-2}] = -2.5 \log (\nu _{\rm proj} [{\rm arcsec}^{-2}]) + Z, \end{equation} \tag{ 30 }$$

where the zero point Z is equal to

$$\begin{equation} Z = 21.572 + M_{\odot, V} + 2.5 \log \lbrace (f/\lambda [L_{\odot }]) (\pi D[{\rm kpc}] / 648)^{2} \rbrace, \end{equation} \tag{ 31 }$$

where M_☉,V = 4.83 (Binney & Merrifield 1998). Although (f/λ[L_☉]) can in principle be estimated, we have found it more reliable to calibrate Z so that the surface brightness profiles from Trager et al. and Paper I agree in their region of overlap. We have done this by including Z as a free parameter in the fits described below, which yields Z = 18.21. For a canonical distance D = 4.8 kpc (vdV06) this implies that (f/λ[L_☉]) = 1.0. With this Z there is excellent agreement between the different data sets, as shown in Figure 1 (top panel).

Zoom In Zoom Out Reset image size

NGB08 studied the surface brightness profile of unresolved light in ω Centauri using HST images. In principle, the advantage of using unresolved light over star counts is that one can trace the cumulative effect of all the unresolved faint stars, rather than counting only the less numerous resolved ones. However, we showed in Paper I that even at HST resolution, the unresolved light is dominated by scattered light from the PSF wings of only the brightest stars. So in reality, the inherent statistics of unresolved light measurements are poorer than those of star count measurements. Moreover, we showed in Paper I that the combined photometric and kinematic center of ω Centauri is actually 12'' away from the position where NGB08 identified it to be. There is no simple way to transform their measurements over concentric annuli to a different center position. For these reasons, we have not included the NGB08 surface brightness data points in the data set used for the dynamical modeling.

5.1.1. Line-of-sight Kinematics

vdV06 compiled LOS velocity data of individual stars in ω Centauri. The original data were obtained by Suntzeff & Kraft (1996), Mayor et al. (1997), Reijns et al. (2005), and K. Gebhardt et al. (2010, in preparation). vdV06 performed various operations and cuts to homogenize and correct the data as necessary for dynamical modeling. This included: improved removal of cluster non-members; removal of stars with low enough brightness or large enough velocity error bars to make their reliability suspect; application of additive velocity shifts to correct for offsets in velocity zero points between data sets; multiplicative rescaling of the velocity error bars where necessary, as dictated by actual measurements of scatter between repeat measurements; subtraction of the systemic LOS velocity; and subtraction of the apparent solid-body rotation introduced into the velocity field through the combination of systemic transverse motion and large angular extent ("perspective rotation"). This produced a final homogenized sample of LOS velocities for 2163 stars with uncertainties below 2.0 km s⁻¹.

Glenn van de Ven kindly made the homogenized sample from vdV06 available to us for use in our data–model comparisons. To determine the profile of ${\bar{\sigma }}_{\rm los}(R)$, we binned the stars in radius and then used the statistical methodology described in Appendix A. Radii were calculated with respect to the center determined in Paper I, which differs by 155 from the center used by vdV06 (which was determined by van Leeuwen et al. 2000). We chose a central bin of 20'' radius (with 25 stars) for maximum resolution near the center, and then chose the remaining bins to contain 200 stars each.

Sollima et al. (2009) obtained LOS velocities of 318 cluster members at large radii, which they merged with observations closer to the center by Pancino et al. (2007). This yielded a homogeneous sample of 946 cluster members, from which they derived the profile ${\bar{\sigma }}_{\rm los}(R)$. Seitzer (1983) presented the ${\bar{\sigma }}_{\rm los}(R)$ profile derived from 118 cluster members with known velocities at a range of radii. We included the published profiles from both studies in our data–model comparisons. For each study, we averaged the first two bins in the profile together, since the first bin had a large error bar and did not provide spatial resolution that competes with the vdV06 compilation.

Scarpa et al. (2003) obtained LOS velocities of 75 cluster members at large radii. From this, they determined the velocity dispersion σ_los in four radial bins. The stars were limited to the west side of the cluster, so the results do not represent averages over a complete annulus on the sky. Nonetheless, we included these observations in our data–model comparisons. We subtracted in quadrature from the published σ_los values the average measurement error of 1 km s⁻¹. To obtain ${\bar{\sigma }}_{\rm los}$, we added in quadrature V²/2, where V is the rotation velocity amplitude quoted in Scarpa et al. (2003); the factor 1/2 is the average of a squared sinusoidal variation over a full range of angles.

The profiles obtained from the Sollima et al. (2009), Seitzer (1983), and Scarpa et al. (2003) observations do not pertain to the same center as that inferred in Paper I. However, the smallest radial bin that we use from these studies is the 4 arcmin diameter combined central bin of Sollima et al. A center offset of 155 (as applicable to the Sollima et al. and Scarpa et al. results, who used the van Leeuwen et al. 2000 center) should not significantly affect statistics calculated over bins that extend to such large distances from the center.

At small radii we used as additional data points the integrated-light measurements from NGB08. They inferred velocity dispersions over square fields of 5'' × 5''. Because of the significant size of these fields, we did not include the impact of PSF convolution in our model predictions for these data. We placed the measurements at the distance from the cluster center inferred in Paper I, which is 12'' for both fields. We cannot use the data to calculate a true average over an annulus. However, we do have two measurements at the same distance from the center at a 60° angle. Including both measurements in our model fits is in low-order approximation similar to fitting the average over an annulus.

NGB08 did not measure the mean velocity over their fields as compared to the systemic velocity. Hence, their measurements are dispersions σ_los rather than RMS velocities ${\bar{\sigma }}_{\rm los}$. However, it is known from previous observations (e.g., vdV06) that V_los/σ_los in ω Centauri decreases toward the center. This ratio is small enough at the position of the NGB08 fields that ${\bar{\sigma }}_{\rm los}\approx \sigma _{\rm los}$ to high accuracy (see Section 5.2.3 for a quantitative estimate).

Figure 3 (top panel) shows the collected ${\bar{\sigma }}_{\rm los}$ data from all the sources, as a function of the projected distance R from the cluster center. The different data sets match well in their regions of overlap.

Zoom In Zoom Out Reset image size

5.1.2. Proper Motion Kinematics

5.1.3. Dependence of Kinematics on Stellar Mass

5.2.1. Sample Selection

5.2.2. rms Proper Motions

5.2.3. Proper Motion Rotation

Our data calibration of Paper I (see discussion in Section 3.6.4 of that paper) removed by necessity any possible solid-body rotation component from the data. This is because we did not have access to stationary background sources in the HST fields, as would be required to measure absolute motions. This calibration has the advantage that it automatically removes any perspective rotation present in the actual proper motions (however, this is negligible for the small radii R ⩽ 717 sampled by our data anyway; ⩽0.0035 mas yr⁻¹ = 0.08 km s⁻¹ at the canonical D = 4.8 kpc). But of course, it has the disadvantage that it also removes a solid-body rotation fit to the actual rotation field of the cluster. This is a fundamental limitation of the data. However, we discuss below that this has only negligible importance for our models.

In Paper I, we also used a "local correction" that removed all other mean systematic motions (as opposed to random motions) from the proper motion data. However, that was by choice, and not by necessity. The motivation for this is to calibrate out small time variations in the higher order geometric distortion terms. This improves the measurement of the random motions in the cluster. However, one would of course not want this to remove an important signal that is present in the data, and in particular, any possible differential rotation in the plane of the sky (i.e., the part of the rotation field that remains after subtraction of a solid-body fit). We tested for such differential rotation in Paper I, and found that it is negligible in the central HST field. Proper motion rotation curves (the mean tangential motion as function of radius) derived from data with and without the local correction agree to within 1 km s⁻¹ at all radii. Therefore, solid-body rotation is the only unrepresented rotation component that could affect the modeling.

In principle, the loss of information on solid-body rotation can be partly recovered as in vdV06. They use the fact that for an equilibrium axisymmetric system the equation

$$\begin{equation} V_{\rm los} = 4.7404 \: D[{\rm kpc}] \tan i \: {\cal V}_{\rm y} \end{equation} \tag{ 33 }$$

must hold everywhere on the projected plane of the sky. Here i is the inclination, y is the projected minor axis direction, and the symbol ${\cal V}$ is used to distinguish a mean proper motion from a mean velocity V. However, to use this equation with accuracy for our HST proper motions would require many more stars with LOS velocities than are actually available at R ⩽ 717.

The pmr direction is perpendicular to the direction of absolute rotation, and our estimates of ${\bar{\Sigma }}_{\rm pmr}$ are therefore not affected by our inability to measure absolute rotation. Also, in an equilibrium system there cannot be mean streaming in the pmr direction, when averaged along a circle. However, our calculated values of ${\bar{\Sigma }}_{\rm pmt}= ({\cal V}_{\rm pmt}^2 + \Sigma _{\rm pmt}^2)^{1/2}$ will be systematically low by a multiplicative factor g ⩽ 1, because we are unable to include any actual solid-body contribution from ${\cal V}_{\rm pmt}^2$. The proper motion data presented by vdV06 did include measurements of absolute rotations in the plane of the sky. Figure 4 shows for their data the ratios ${\cal V}_{\rm pmt}/\Sigma _{\rm pmt}$ and $g \equiv \Sigma _{\rm pmt} / {\bar{\Sigma }}_{\rm pmt}$, where the quantities in the nominator and denominator are all averages over an annulus. These measurements allow us to estimate the value of g for the HST data. The central data point shows explicitly that there is very little rotation in ω Centauri at the radii of our central HST field. More generally, in the central 10 arcmin, the rotation rate is approximately linear with radius, ${\cal V}_{\rm pmt}/\Sigma _{\rm pmt} \approx R/750^{\prime \prime }$. Therefore, g ≈ [1 + (R/750'')²]^−1/2. These approximations are shown as solid lines in the figure. This approximation implies that 0.995 ⩽ g ⩽ 1 for the radii R ⩽ 717 pertaining to the HST sample.

Zoom In Zoom Out Reset image size

Given what is already known about proper motion rotation from the vdV06 data, our estimates of ${\bar{\Sigma }}_{\rm pmt}$ should be quite accurate. This is despite the fact that any solid-body part of the ${\cal V}_{\rm pmt}$ field cannot be measured directly in the HST data. In particular, the value of g is insufficient to explain the observed anisotropy ${\bar{\Sigma }}_{\rm pmt}/{\bar{\Sigma }}_{\rm pmr}= 0.983 \pm 0.006$. Also, the data of vdV06 show that V_los/σ_los is about a factor of 2 lower than ${\cal V}_{\rm pmt}/\Sigma _{\rm pmt}$. Hence, the influence of rotation is even more negligible for the LOS measurements at R = 12'' of NGB08 (see the discussion in Section 5.1.1).

5.2.4. Gauss–Hermite Moments

To analyze the shapes of the HST proper motion distributions, we binned the data in concentric circular annuli as in Section 5.2.2. For each bin, we created separate histograms of the proper motions in the pmt and pmr directions. We used bins that correspond to 1 km s⁻¹ for the canonical distance of ω Centauri (D = 4.8 kpc; vdV06). For each histogram, we calculated the best-fitting Gaussian and the corresponding Gauss–Hermite moments, as in van der Marel et al. (2000). The odd moments are zero for distributions that are averaged over circles. Figure 5 shows the lowest order non-trivial even Gauss–Hermite moments h₄ and h₆. These were calculated for annuli that are 8'' wide. The average values over all radii R ⩽ 717 are h_4,pmr = −0.022 ± 0.006, h_4,pmt = −0.024 ± 0.006, h_6,pmr = −0.003 ± 0.006, and h_6,pmt = 0.005 ± 0.006. The 8th- and 10th-order moments were also calculated, but are not shown here. They were found to be consistent with zero, like h₆. For all moments, the values for the pmr and pmt directions were found to be consistent within the errors. The only moment that is non-zero is h₄, with an average value over the two proper motion directions of h₄ = −0.023 ± 0.004. This is indicative of proper motion distributions that are slightly more flat-topped than a Gaussian, and have slightly less extended wings. We will compare this to model predictions in Section 6.7.

Zoom In Zoom Out Reset image size

The observed profile of ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}$ shows statistically significant deviations from unity, both at small and at large radii (see Figure 3 (bottom panel)). It is therefore important to allow for anisotropy in our models for ω Centauri.

To illustrate the relation between the model anisotropy function β(r) and the observed quantity ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}(R)$, we show in Figure 6 (top) the predicted ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}(R)$ for models with constant β as a function of radius r. Models were calculated with $(\overline{v_{\rm tan}^2}/\overline{v_r^2})^{1/2}$ ranging from 2/3 to 3/2, in equal logarithmic steps. At large radii, the proper motion anisotropy is a good tracer of the intrinsic anisotropy, consistent with the arguments presented in Section 2.7.3. For example, at R = 10 arcmin, $({{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}) - 1 = \tau [(\overline{v_{\rm tan}^2}/\overline{v_r^2})^{1/2} - 1]$, with τ = 0.7–0.9 for the models calculated here, and τ ≈ 0.8 for β ≈ 0. At small radii, the behavior of ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}$ is a more complicated function of β. Tangentially anisotropic models (β < 0) all converge to ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}} = 1$ near the center, while radially anisotropic models (0 < β ⩽ 1) predict ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}} < 1$.

Zoom In Zoom Out Reset image size

It is evident from Figure 6 (top) that the observed ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}(R)$ must tightly constrain the intrinsic β(r). We fitted the combined kinematical data from Figure 3 with models that spanned a grid in β₀, β_∞ and r_a. The parameter β₀ determines the behavior of ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}(R)$ near the center, β_∞ determines the behavior at large radii, and r_a is the transition radius. The minimum χ² was obtained for β₀ = 0.13 ±  0.02, β_∞ = −0.53 ± 0.22 and log r_a[arcsec] = 2.86 ± 0.12. The error bars were obtained from the χ² contours of the fit.

The predicted ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}(R)$ for the best-fit model is shown in Figure 6 (bottom). It provides an excellent fit. The predictions shown in the figure were calculated for core models with no dark mass. However, the predicted ratio ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}}(R)$ depends mostly on β(r), and is virtually indistinguishable when considering instead models with a cusp or dark mass.

The error bars of the fit show that velocity anisotropy is detected at high statistical significance in ω Centauri. The parameters β₀ and β_∞ imply that the ratio $\overline{v_{\rm tan}^2}/\overline{v_r^2}$ transitions from 0.935 near the center to 1.235 at large radii. These are of course not necessarily the values at the very center and at infinity, since our model is constrained by proper motion data only over the range 1–1000 arcsec. Either way, there is a transition from radial anisotropy at small radii to tangential anisotropy at large radii. The parameter r_a indicates that the transition occurs at approximately 12 ± 3 arcmin (5 ± 1 core radii). This is broadly consistent with the modeling results of vdV06. They found (their Section 9.2) that ω Centauri is slightly radially anisotropic for r ≲ 10 arcmin and that it becomes increasingly tangentially anisotropic outside this region.

The centers of globular clusters are generally believed to have isotropic velocity distributions because two-body relaxation tends to isotropize the orbits. However, ω Centauri has a long half-mass relaxation time of ∼10^9.96±0.03 yr (McLaughlin & van der Marel 2005). We showed in Paper I that mass segregation in the central region has not yet progressed to the point of energy equipartition. Since full energy equipartition has not yet been achieved, it is not surprising that complete velocity isotropy has not yet been achieved either.

A sufficient number of stars are necessary to measure the proper motion anisotropy with accuracy. Due to the finite number of stars at small radii, the anisotropy is not as well constrained at radii R ≲ 10'' as it is further out. However, because Omega Cen has a very large core, most of the stars that are observed near the projected center are not actually close to the center in three dimensions (see, e.g., Section 6.1 of Paper I). One consequence of this is that the projected kinematics predicted near the center are influenced primarily by the anisotropy at larger radii. Changing the model anisotropy in the central 10'' even quite substantially (e.g., ±30% in σ_r/σ_t) does not change the model predictions for the HST proper motion kinematics by more than a fraction of the error bars, even for the innermost data point. Therefore, uncertainties in anisotropy near the center do not significantly impact our predictions or conclusions.

Figures 7 (top) and (middle) show the predicted and observed RMS projected velocities for the best-fit β(r). Each panel of this figure shows the data for ${\bar{\sigma }}_{\rm los}(R)$ and ${\bar{\sigma }}_{\rm pm}(R)$ together in the same panel, with the observed proper motions ${\bar{\Sigma }}_{\rm pm}(R)$ transformed to km s⁻¹ using the best-fit distance implied by the model. For each model, there are separate curves predicted for ${\bar{\sigma }}_{\rm los}(R)$ and ${\bar{\sigma }}_{\rm pm}(R)$ (dashed and solid curves, respectively), since these quantities are not the same in an anisotropic model. However, the difference is generally smaller than the observational error bars over the radial range where both types of data are available. There is no LOS data available for R ≲ 10'', so for visual clarity we do not show the predicted ${\bar{\sigma }}_{\rm los}(R)$ at those radii. At small radii ${\bar{\sigma }}_{\rm los}(R)$ exceeds ${\bar{\sigma }}_{\rm pm}(R)$ because of the model anisotropy, with ${\bar{\sigma }}_{\rm los}-{\bar{\sigma }}_{\rm pm}\approx 0.9 \;{\rm km}\;{\rm s}^{-1}$ at R = 10'' and ${\bar{\sigma }}_{\rm los}-{\bar{\sigma }}_{\rm pm}\approx 1.1 \;{\rm km}\;{\rm s}^{-1}$ at R = 1''.

Zoom In Zoom Out Reset image size

Figure 7 (top) shows the predictions for the core model fit to the projected intensity and Figure 7 (middle) for the cusp model fit. The curves that predict the lowest ${\bar{\sigma }}$ at small radii (blue) are models with no dark mass. The core and cusp models generally predict the same dynamics at large radii, but differ slightly at small radii.

It is true in many analytical models with a cusped density profile and no central dark mass that the projected RMS velocity decreases toward the center (e.g., Tremaine et al. 1994). The same can be seen in Figure 7 (middle) for our model with γ = 0.05. Even though the cusp is quite shallow in projection, this corresponds to a much steeper three-dimensional luminosity density j ∝ r^−1−γ at asymptotically small radii. Similarly, the decrease in RMS velocities toward the center is much stronger intrinsically in three dimensions than it is in projection. Either way, the observed RMS velocities do not show a significant dip near the center. In the absence of a dark mass, the core models therefore provide a better fit to the data than the cusp models.

The best fits without a dark mass have χ² = 235.3 and 249.2 for the core and cusp models, respectively. There are 202 data points, yielding N_DF = 197 degrees of freedom (there are five free parameters in the models when there is no dark mass). One expects $\chi \approx N_{\rm DF} \pm \sqrt{2 N_{\rm DF}}$. The anisotropic core model is therefore consistent with the data at the 1.9σ level, while the cusp model is consistent at the 2.6σ level. These confidence levels assume that all data points in the fit are statistically independent, which is somewhat of a simplification. It is likely that samples from different LOS velocity studies have some stars in common. However, this cannot be fully explored here since Scarpa et al. (2003) and Sollima et al. (2009) did not list individually which stars they observed.

Much of the excess E(χ²) ≡ χ² − N_DF > 0 for the fits is due to just two data points for ${\bar{\sigma }}_{\rm los}$ that appear anomalously high. The first "discrepant" data point is the measurement ${\bar{\sigma }}_{\rm los}= 23.0 \pm 2.0 \;{\rm km}\;{\rm s}^{-1}$ at R = 124 by NGB08. The second discrepant data point is the average ${\bar{\sigma }}_{\rm los}= 19.9 \pm 1.0 \;{\rm km}\;{\rm s}^{-1}$ for an annulus centered at R = 460, obtained from the compiled data of vdV06. The former contributes 6.8 to the χ² of the fit for the core model, while the latter contributes 11.3. This adds up to 18.1, which is almost half of E(χ²) = 38.3.

The discrepant NGB08 data point pertains to similar radius as the other field that they observed, which is at R = 121 from the cluster center (see Paper I). The measurement ${\bar{\sigma }}_{\rm los}= 18.6 \pm 1.6 \;{\rm km}\;{\rm s}^{-1}$ for that latter field is more consistent with the other measurements in our combined data set from Section 5. The two NGB08 measurements are inconsistent with each other at the 1.7σ level, despite the fact that that they pertain to the same radius, and despite the fact that the proper motions in these fields are consistent with each other (Paper I). We believe that this may be due to underestimates in the NGB08 error bars, since they did not explicitly include the contribution of shot noise (from the finite number of stars) to their error bars. Similarly, the measurement in the vdV06 annulus at R = 460 appears inconsistent with other measurements in that same data set. The weighted average of the two adjacent annuli is lower and inconsistent with it at the 2.1σ level. This suggests that both discrepant measurements that contribute disproportionately to χ² may be spuriously high for unknown reasons.

Evidently, there are some internal discrepancies in the data sets themselves. Hence, the observed values of E(χ²) should not be taken as a sign of shortcomings in the models. In fact, the values of E(χ²) are remarkably low, given that we are fitting kinematical data compiled from 10 different original sources. Potential uncorrected systematics between different studies could easily have induced much larger data–model discrepancies.

One clear result from Figures 7 (top) and (middle) is that the data are well fit by the models out to the largest available radii. This contradicts arguments put forward by Scarpa et al. (2003) that the large radii kinematics of ω Centauri may be inconsistent with traditional theories of gravity. Our conclusion on this issue agrees with that of McLaughlin & Meylan (2003) and Sollima et al. (2009).

To test for the presence of an IMBH, we have also constructed models with M_BH > 0 (and the same β(r) profile). Models with no IMBH have one free parameter less than models with an IMBH. One would therefore on average expect them to fit worse by Δχ² = 1. Moreover, random variations can add another Δχ² = 1 (at 1σ confidence), even when the no-IMBH model is correct. Therefore, an IMBH is statistically significant in the model only if it improves the fit over a no-IMBH model by Δχ² ≳ 2.

Figure 8 shows the curves of χ²(M_BH) that quantify the quality of the fit. For the core model, the best-fit IMBH has mass M_BH = 4.1 × 10³ M_☉. This improves the fit over the no-IMBH model by only Δχ² = 1.8. This is not statistically significant. For the core model, we can therefore determine only an upper limit to the mass of a possible IMBH, namely M_BH ⩽ 7.4 × 10³ M_☉ at 1σ confidence. Figure 7 (top) compares the best-fit core models with and without an IMBH. The inclusion of the IMBH makes very little difference to the predictions. The difference for the innermost data point at R = 11 is only 0.8 km s⁻¹, which is less than a quarter of the observational error bar at that radius.

Zoom In Zoom Out Reset image size

For the cusp model, the inclusion of an IMBH improves the fit by Δχ² = 9.8, which is significant at the 3.1σ level. The best-fit IMBH mass is M_BH = (8.7 ± 2.9) × 10³ M_☉ (with the error bar determined by the requirement Δχ² ⩽ 1). Figure 7 (middle) compares the best-fit cusp models with and without an IMBH. The inclusion of the IMBH provides a small but visible improvement in the fit near the center. The reason that the IMBH is more massive and significant in the cusp model is due to the dip in RMS velocity predicted by these models when there is no dark mass.

In assessing the meaning of formal error bars on M_BH, it is important to keep in mind what one may call "black hole bias." Statistical error bar estimates take into account only the random Gaussian uncertainties in the data. They do not take into account the systematic residuals that are inevitably present in any astronomical data–model comparison. In the present context, these may include oversimplifications in the model assumptions (e.g., sphericity, equilibrium, the adopted density parameterization, no mass segregation among the luminous stars), the data reduction process (e.g., inability to measure absolute proper motion rotation), or the meaning of the data themselves (e.g., inability to separate orbital motion from potential binary motion). By adding M_BH as an extra free parameter to the fit, the model gains in ability to fit these systematic residuals that are not actually due to an IMBH. Because masses are positive definite (M_BH ⩾ 0), this generally tends to bias in the direction of invoking more mass then there actually is, as compared to the opposite. Hence, M_BH tends to be biased high. Marginally significant non-zero M_BH values with no telltale signs of an IMBH (e.g., fast-moving stars or an R^−1/2 increase in ${\bar{\sigma }}$) may simply be indications of small systematic issues in the data–model comparison unrelated to an IMBH.

The core model with no IMBH fits the kinematical data a little better than the cusp model with an IMBH (see Figure 8). The difference is Δχ² = 4.0, which is nominally significant at the 2.0σ level. However, as discussed in Section 4, the cusp model nominally provides a better fit to the photometric data. The difference for a generalized nuker fit is Δχ² = 6.8, which is nominally significant at the 2.6σ level. One could in principle construct a grand total Δχ² that sums the photometric and kinematic values. However, we have not done this here. The arguments about whether a core or cusp model better fits the photometry are complex (see Section 4). They are not easily captured in a single number without having to make arbitrary choices about, e.g., what radii to include in the comparison. Either way, when comparing both the photometric and kinematic information for the core and cusp models we feel that there is no sound statistical basis to prefer one model over the other. In view of this, we conclude that the existing data do not imply that an IMBH is present in ω Centauri. However, IMBHs with masses M_BH ≲ 1.2 × 10⁴ M_☉ cannot be ruled out at 1σ confidence (or ≲1.8 × 10⁴ M_☉ at 3σ confidence).

The best-fit distance and mass-to-light ratio depend only slightly on the particular model we adopt. The core model without a dark mass has D = 4.70 ± 0.06 kpc and ϒ_V = 2.64 ± 0.03 (here and henceforth, mass-to-light ratios are given in V-band solar units). The cusp model with an IMBH has D = 4.75 ± 0.06 kpc and ϒ_V = 2.59 ± 0.03.

Our results for D and ϒ_V are entirely consistent with those inferred by vdV06 from fitting only ground-based projected intensity and kinematical data. They found that D = 4.8 ± 0.3 kpc and ϒ_V = 2.5 ± 0.1. Their modeling technique was considerably more sophisticated than ours, allowing for axisymmetry, inclination, ellipticity variations with radius, and arbitrary profiles of rotation and velocity-dispersion anisotropy. The reason that our error bars are smaller than theirs is likely a combination of two effects. First, we have considerably more data at our disposal than did vdV06, and in particular we have a factor ∼11 more stars with measured proper motions. This produces a significant decrease in the random uncertainties. But second, in our simpler models we cannot explore the full variety of phase-space distributions that might fit the data with non-standard distances or mass-to-light ratios. Therefore, we are not quantifying systematic errors as rigorously as vdV06, and as a result our errors are artificially lowered compared to theirs.

The fact that our D and ϒ_V agree with those of vdV06 indicates that there is no reason to mistrust that our spherical models can address the central mass distribution of ω Centauri with credible accuracy. This is further supported by the facts that both the elongation (Geyer et al. 1983) and the rotation rate (Section 5.2.3) decrease toward the cluster center. To get good results, we did have to make sure to model RMS velocities (which sum rotation velocities and velocity dispersions in quadrature) and not merely velocity dispersions. Also, one should fit only quantities that are properly averaged over annuli on the projected plane of the sky. If these things are not done, then biased estimates are obtained for D and ϒ_V, as demonstrated in Appendix C of vdV06.

Isotropic models predict that ${{\bar{\Sigma }}_{\rm pmt}}/{{\bar{\Sigma }}_{\rm pmr}} = 1$ at all radii, so they are not generally appropriate for ω Centauri. However, isotropic models have the advantage that the full velocity distributions ${\cal L}_{\rm iso}(v,R)$ can be calculated as described in Section 2.5. Also, isotropic models are often used for globular clusters, so it is useful to understand how the predictions of such models deviate from more general anisotropic models.

Figure 8 shows the curve of χ²(M_BH) for isotropic models with an IMBH. The models fit significantly worse than the anisotropic models, as shown by the higher χ². Within the realm of isotropic models, the best fit is provided by M_BH = (1.8 ± 0.3) × 10⁴ M_☉. This model has D = 4.81 ± 0.06 kpc and ϒ_V = 2.61 ± 0.03. Figure 7 (bottom) shows the data–model comparison for the isotropic model with this IMBH. We also show the predictions for M_BH = 0, as well as for the value M_BH = 4.0 × 10⁴ M_☉ advocated by NGB08.

The NGB08 IMBH mass was based on their measurement of ${\bar{\sigma }}_{\rm los}= 23.0 \pm 2.0 \;{\rm km}\;{\rm s}^{-1}$ at a position they believed to be the cluster center. This high IMBH mass is now clearly ruled out, even for an isotropic model, for two reasons. First, Paper I showed that this measurement actually applies to a larger distance R = 124 from the cluster center. The arrow in Figure 7 (bottom) demonstrates how this moved the corresponding LOS velocity measurement. And second, the high measured RMS velocity was not confirmed by proper motion measurements, either at the same position or on the actual cluster center. Note in Figure 7 (bottom) that the isotropic model with the NGB08 IMBH does appear to match well the data at R = 5''–7''. However, it overpredicts the data at both R ⩽ 4'' and R = 8''–10''. Since all data points were derived in the same manner and are equally valid, there is no reason attribute more weight to the data at R = 5''–7''. In particular, the data at R ⩽ 4'' contain most of the information on the very central mass distribution, and they clearly rule out an IMBH mass as high as advocated by NGB08 (independent of the model anisotropy near the center).

It is evident from Figure 8 that isotropic models require a higher IMBH mass than the best-fit anisotropic models. This is not generally true for isotropic models, but is a specific consequence of the fact that ω Centauri is radially anisotropic near the center and tangentially anisotropic at large radii. This affects the general gradient in RMS velocity from small to large radii (see e.g., Figures 10 and 11 of van der Marel 1994a). The anisotropic model without a dark mass has a steeper gradient toward the center than the corresponding isotropic model. The anisotropic model therefore requires less dark mass to fit the observed gradient. This indicates that one should be careful in application of isotropic models to studies of IMBHs. Depending on the exact anisotropy of the cluster, the assumption of isotropy can lead to either an overestimated or an underestimated IMBH mass, or it can require an IMBH when none is present.

When a dark mass is clearly required by the data, as in the case of isotropic models with a cusp, one can ask whether an extended dark cluster may provide a better fit than an IMBH. Figure 9 shows a two-dimensional contour plot of Δχ² in the parameter space of (a_dark, M_dark) for the isotropic cusp model. The best fit is obtained for M_dark = 2.0 × 10⁴ M_☉ and a_dark = 3''. The predictions of this model are shown in Figure 7 (bottom). It has a somewhat shallower increase in RMS velocity toward the center than the IMBH model. The dark-cluster model makes an improvement of only Δχ² = 1.7 over the IMBH model, while increasing the number of free parameters by one. This is not statistically significant. Therefore, one can interpret Figure 9 to mean that if there is a dark mass in ω Centauri, then its extent must be ≲7'' at 1σ confidence.

Zoom In Zoom Out Reset image size

The central density of the best-fitting dark cluster is ρ_0,dark ≡ M_dark/(4πa³_dark/3), according to Equation (7). Using Equation (4) to transform angular units to physical units, this yields for the core model ρ_0,dark = 1.4 × 10⁷ M_☉ pc⁻³. The central density of luminous matter in this model is only ρ_0,lum = 3.3 × 10³ M_☉ pc⁻³. This implies a central M/L = ϒ_V[1 + (ρ_0,dark/ρ_0,lum)] = 1.1 × 10⁴. This applies to an isotropic model with M_dark = 2.0 × 10⁴ M_☉. For an anisotropic model with a lower M_dark, the required M/L scales approximately linearly with M_dark. The (maximum) size a_dark of the dark mass fit is relatively insensitive to the exact anisotropy. Either way, if ω Centauri does indeed have a dark central cluster, this would be indicative of quite an extreme amount of mass segregation or a very top-heavy initial mass function. If a dark mass is invoked to explain the data, then an IMBH may be easier to explain than such a concentrated dark cluster.

For our isotropic models, we calculated the projected velocity distributions ${\cal L}_{\rm iso}(v,R)$ as described in Section 2.5, and from this the Gauss–Hermite moments h_i,iso(R). We did not take the effect of observational proper motion uncertainties into account in the models (essentially a convolution of the model predictions with the error distribution). This has negligible effect on the Gauss–Hermite moments because the observational errors are much less than the intrinsic cluster dispersion. The curves overplotted in Figure 5 show the model predictions for three models: the isotropic core model with no dark mass, the isotropic cusp model with its best-fit M_BH = 1.8 × 10⁴ M_☉, and the isotropic cusp model with no dark mass.

The core model with no dark mass and the cusp model with an IMBH both fit the observed Gauss–Hermite moments within the error bars. Over the range 5'' ≲ R ⩽ 717, the predicted Gauss–Hermite moments for these models are approximately constant. For the core model without dark mass they are h_4,iso(R) ≈ −0.013 and h_6,iso(R) ≈ −0.002, and for the cusp model with an IMBH they are h_4,iso(R) ≈ −0.034 and h_6,iso(R) ≈ 0.009. These averages do deviate in statistically significant manner from the observed averages, h₄ = −0.023 ± 0.004 and h₆ = 0.001 ± 0.004 (see Section 5.2.4). This is not surprising, given that we have already established that ω Centauri does not have an isotropic velocity distribution, not even near its center. However, the differences between the predicted and the average observed Gauss–Hermite moments are quite small |Δh_i| ≲ 0.01. For comparison, models that range in anisotropy from fully tangential to fully radial can easily span the range Δh_i ≈ ±0.2 (e.g., van der Marel & Franx 1993). The measured Gauss–Hermite moments are therefore consistent with the fact that the velocity distribution is not far from isotropic near the center of ω Centauri.

The cusp model without dark mass does not fit the observed Gauss–Hermite moments. As discussed in Section 6.2, this model has a strong decrease in its intrinsic (unprojected) RMS velocity toward the center. The decrease in the projected RMS velocity is relatively mild, as shown in Figure 7 (bottom). However, there is a strong gradient in the RMS velocity along the LOS. This leads to strongly non-Gaussian velocity distributions. In a central aperture of R = 8'' radius, as used for the innermost observational data point in Figure 5, the predictions are h_4,iso = 0.063 and h_6,iso = −0.063. This deviates for each order by ∼3σ from the observed values h₄ = −0.026 ± 0.029 and h₆ = 0.031 ± 0.029 (averaged over the pmr and pmt directions). While this result was derived for isotropic models, the large predicted deviations from Gaussian velocity distributions are likely generic, and would exist in anisotropic models as well. So if ω Centauri has a central cusp in its projected intensity profile, then the observed velocity distributions rule out models without a dark mass. Of course, Figure 7 (middle) and (bottom) show that such models also do not provide good fits to the projected profiles of RMS velocity.

Figure 5 shows that an isotropic model with an IMBH predicts increased h₄ and h₆ in the central R ≲ 5'', corresponding to profiles with broader wings than a Gaussian. To do a more detailed data–model comparison for the wings, we show in Figure 10 the histograms of the observed proper motions for apertures of radii R = 3'' and R = 10''. These apertures contain N = 43 and N = 585 stars with well-measured proper motions, respectively. The proper motion coordinates in the pmr and pmt direction were both included in each histogram, as appropriate for an isotropic velocity distribution. The histograms were symmetrized to decrease the uncertainties. The predictions for the isotropic cusp model with M_BH = 1.8 × 10⁴ M_☉ are overplotted as curves. The predicted velocity distribution for the R = 3'' aperture has more extended wings than the observed histogram. The wings are much reduced in the predicted velocity distribution for the R = 10'' aperture. This is because none of the added stars at 3'' < R ⩽ 10'' are close to the black hole, so they do not move at unusually high velocities.

Zoom In Zoom Out Reset image size

The observed histograms do not extend beyond 60 km s⁻¹. The maximum observed value of either |v_pmt| or |v_pmr| inside R = 10'' is 59.4 km s⁻¹. No stars were rejected as non cluster members inside this radius because of their fast motion (see Section 5.2.1). To assess whether these observations are statistically consistent with the wings of the predicted distributions, let f be the fraction of the normalized model distribution at |v|>60 km s⁻¹. According to the models, this is f = 3.8 × 10⁻³ for the R = 3'' aperture and f = 3.9 × 10⁻⁴ for the R = 10'' aperture. The expectation value for the number of stars with a velocity component at |v|>60 km s⁻¹ is E = 2Nf. The probability that all of the 2N measured velocities have |v| ⩽ 60 km s⁻¹ is P = (1 − f)^2N. This yields for the R = 3'' aperture that E = 0.32 and P = 0.72, while for the R = 10'' aperture E = 0.46 and P = 0.63.

Since less than a single fast-moving star was predicted, the fact none were observed is not statistically inconsistent with the models. This does not change when we use larger apertures, since the models do not predict any fast-moving stars outside R = 10''. It also does not change by using smaller apertures, since in that case the sample of observed stars becomes very small. Therefore, the absence of fast-moving stars in the data cannot be used to independently rule out the presence of an IMBH with mass M_BH = 1.8 × 10⁴ M_☉. To increase the predicted number of fast-moving stars one would need to increase the number of stars with measured proper motions, i.e., extend the kinematical observations to fainter magnitudes. We have experimented with looser cuts on our proper motion catalog. This allows us to increase the sample size by a factor ∼2 (see Paper I). In this extended catalog there is no star with a proper motion coordinate in excess of 65 km s⁻¹ within R = 10''. This provides somewhat tighter constraints on the presence of an IMBH, with the quoted values of P roughly decreasing to their squared values. Even then, it is still not possible to rule out an IMBH of mass M_BH = 1.8 × 10⁴ M_☉ with reasonable confidence.

Alternatively, E increases and P decreases when the BH mass is increased. For the IMBH of mass M_BH = 4.0 × 10⁴ M_☉ suggested by NGB08, and sticking as before with our high-quality proper motion catalog of Section 5.2.1, we obtain E = 2.9 and P = 0.05 for the aperture with R = 10''. This implies that the absence of fast-moving stars rules out all IMBHs with M_BH ≳ 4.0 × 10⁴ M_☉ at 95% confidence. This result is not very competitive for ω Centauri, since we already know from modeling of the RMS velocities that any IMBH must have a mass at least 3.3 times lower than this anyway. However, this method has the advantage that it does not require any accurate measurement of RMS velocities. It may therefore be useful for other clusters, in cases when the observational errors may be too large for accurate measurement of the RMS velocity profile.

We have assumed in the modeling that any IMBH must reside at the cluster center. In reality, an IMBH would exhibit Brownian motion as a result of its energy equipartition with surrounding stars. Chatterjee et al. (2002) derived for the one-dimensional RMS offset from the cluster center that

$$\begin{equation} x_{\rm RMS} = (2 {\tilde{m}} / 9 M_{\rm BH})^{1/2} \: r_0. \end{equation} \tag{ 34 }$$

This was derived for an analytical Plummer model with core radius r₀, which provides a reasonable description of (non core-collapsed) globular clusters. Chatterjee et al. considered single-mass models with stars of mass ${\tilde{m}}$. However, their results should be valid for multi-mass models as well, provided that the quantity ${\tilde{m}} \equiv \langle m^2 \rangle / \langle m \rangle$ is defined appropriately in terms of averages over the stellar mass function (Merritt et al. 2007).

We take for r₀ the King core radius of 141'' (McLaughlin & van der Marel 2005) and for M_BH the upper limit of 1.2 × 10⁴ M_☉ derived in Section 6.3. For a Salpeter mass function from 0.1–10 M_☉, ${\tilde{m}} = 1.27 \;{M_{\odot }}$. We then obtain x_RMS = 069. This is significantly less than the 155 radius of the smallest aperture used in the kinematical analysis (see Section 5.2.2). Therefore, the expected Brownian motion of an IMBH should not affect significantly the analysis that we have presented here.

IMBHs with larger masses than those considered here, placed at considerable offsets from the cluster center (consistent with their expected Brownian motion), might have gone unnoticed in our analysis. However, there is no independent motivation to consider such a situation. In Paper I, we showed that the kinematical center of stars in ω Centauri agrees with the star count center to within its ∼2'' uncertainty.