# New Limits on an Intermediate Mass Black Hole in Omega
Centauri:
II. Dynamical Models^{1}

^{1}

## Abstract

We present a detailed dynamical analysis of the projected density and kinematical data available for the globular cluster Centauri. We solve the spherical anisotropic Jeans equation for a given density profile to predict the projected profiles of the RMS velocity , in each of the three orthogonal coordinate directions (line of sight, proper motion radial, and proper motion tangential). The models allow for the presence of a central dark mass, such as a possible intermediate-mass black hole (IMBH). We fit the models to new HST star count and proper motion data near the cluster center presented in Paper I, combined with existing ground-based measurements at larger radii. The projected density profile is consistent with being flat near the center, with an upper limit on the central logarithmic slope. The RMS proper motion profile is also consistent with being flat near the center. The velocity anisotropy profile, distance, and stellar mass-to-light ratio are all tightly constrained by the data, and found to be in good agreement with previous determinations by van de Ven et al. To fit the kinematics we consider anisotropic models with either a flat core () or a shallow cusp (). Core models provide a good fit to the data with ; cusp models require a dark mass. If the dark mass in cusp models is an intermediate-mass black hole (IMBH), then ; if it is a dark cluster, then its extent must be . Isotropic models do not fit the observed proper motion anisotropy and yield spuriously high values for any central dark mass. These models do provide a good fit to the Gauss-Hermite moments of the observed proper motion distributions (, ). There are no unusually fast-moving stars observed in the wings of the proper motion distribution, but we show that this does not strongly constrain the mass of any possible IMBH. The overall end-result of the modeling is an upper limit to the mass of any possible IMBH in Centauri: at confidence (or at confidence). The limit corresponds to %. We combine this with results for other clusters and discuss the implications for globular cluster IMBH demographics. Tighter limits will be needed to rule out or establish whether globular clusters follow the same black hole demographics correlations as galaxies. The arguments put forward by Noyola et al. to suspect an IMBH in Centauri are not confirmed by our study; the value of they suggested is firmly ruled out.

###### Subject headings:

globular clusters: individual ( Centauri) — stars: kinematics.## 1. Introduction

Intermediate-mass black holes (IMBHs) are of interest in a variety of astrophysical contexts (see, e.g., the reviews by van der Marel 2004; Miller & Colbert 2004). Especially the possible presence of IMBHs in the centers of globular clusters has intrigued astronomers for decades (e.g., Bahcall & Wolf 1976). This subject area underwent a considerable resurgence in recent years due to a combination of two advances. First, theoretical modeling suggested plausible formation scenarios for IMBHs in the centers of globular clusters from realistic initial conditions (Portegies Zwart & McMillan 2002). Second, the quality of observational data sets, especially those obtained with the Hubble Space Telescope (HST), advanced to the point where it is possible to measure velocity dispersions at a spatial resolution comparable to the sphere of influence size for plausible IMBH masses. This led to the identification of three globular clusters with candidate IMBHs, namely M15, G1, and Centauri.

Gerssen et al. (2002) used line-of-sight velocity data for individual stars in the core-collapsed cluster M15 (NGC 7078) to argue for the presence of in dark material near the center, possibly in the form of an IMBH. Dull et al. (1997) and Baumgardt et al. (2003a) showed that this could be a plausible result of mass-segregation, and need not be an IMBH. Moreover, Baumgardt et al. (2005) argued that the observed steep brightness profile of M15, typical of the class of core-collapsed clusters, is in fact not consistent with expectation for a relaxed distribution of stars around an IMBH. Ho, Terashima, & Okajima (2003) and Bash et al. (2008) failed to detect X-ray or radio emission coincident with the cluster center. This is not necessarily inconsistent with an IMBH, since the accretion rate and/or the accretion efficiency are likely to be low. Nonetheless, M15 does not appear to be as strong an IMBH candidate as once thought.

Gebhardt et al. (2002, 2005) used the velocity dispersion measured from integrated light near the center of the M31 cluster G1 to argue for the presence of a dark mass at the cluster center. The relaxation time of G1 is too long to have produced such a concentration of dark mass due to mass segregation. Therefore, this result was taken as evidence of for the presence of an IMBH. Baumgardt et al. (2003b) disagreed with this interpretation, and argued that models with no dark mass (IMBH or otherwise) provided an equally acceptable fit. The crux of this disagreement lies in the extremely small size of the sphere of influence () of the putative IMBH, which is less than the spatial resolution of the observations. Phrased differently, the IMBH induces only a small signature in the data. As a result, the interpretation of the data is sensitive to the exact details of the data-model comparison. The statistical significance of the IMBH detection is not high, especially if one takes into account the possibility of “black hole bias” (discussed later in Section 6.3). On the other hand, the possible presence of an IMBH in G1 recently gained further credence with the detection of weak X-ray and radio emission from near the cluster center (Pooley & Rappaport 2006; Kong 2007; Ulvestad, Greene & Ho 2007). The emission properties are consistent with accretion onto an IMBH, but other explanations cannot be ruled out.

Noyola et al. (2008; hereafter NGB08) argued for the presence of an IMBH in the center of the cluster Centauri (NGC 5139), which is well known for its many enigmatic properties (e.g., Meylan 2002). The dynamics of this cluster were studied previously by van de Ven et al. (2006; hereafter vdV06) using sophisticated axisymmetric models. However, there was insufficient kinematical data at small radii to meaningfully constrain the possible presence of any dark mass. NGB08 performed integrated-light measurements of the line-of-sight velocity dispersion for two fields, one on the cluster center that they had determined, and one at from that center. They obtained and for these fields, respectively. The increase in towards the cluster center was interpreted to imply the presence of an IMBH of mass . They also determined the surface brightness profile of Centauri from HST images, and inferred a shallow central cusp of logarithmic slope . This is consistent with theoretical predictions for the cusp induced by an IMBH (Baumgardt et al. 2005).

Both G1 and Centauri are somewhat atypical globular clusters. They are the most massive clusters of their parent galaxies, M31 and the Milky Way, respectively. Both clusters have been suggested to be the stripped nuclei of larger galaxies (e.g., Freeman 1993; Meylan et al. 2001). Indeed, the nuclei of galaxies have very similar structural properties to globular clusters (these nuclei are therefore also referred to as nuclear star clusters; e.g., Walcher et al. 2005). Accretion activity is known to exist in at least some of these nuclei (e.g., Seth et al. 2008). So if G1 and Centauri contain IMBHs, then these may simply be the low-mass tail of the distribution of super-massive black holes known to exist in galactic nuclei. As such, they would not necessarily be telling us anything fundamentally new about the formation and evolution of globular clusters in general. Nonetheless, G1 and Centauri are probably the best current candidates for globular clusters with IMBHs. Validating the existence of their suggested IMBHs is therefore of fundamental importance.

The (initial) IMBH evidence for M15, G1, and Centauri was all based on studies of line-of-sight velocities. Such studies are complicated because they require spectroscopy in regions of high stellar density. Integrated-light measurements work well for a cluster as distant as G1 (Gebhardt et al. 2002, 2005), but the downside is that the sphere of influence in such clusters is poorly resolved. In closer Milky Way clusters, the limited number of stars (especially brighter ones) introduces significant shot noise in integrated-light measurements. For Centauri, this required special treatment of the data to counteract this (NGB08). For M15, HST/STIS was used to obtain spectroscopy of individual stars (van der Marel et al. 2002). However, the signal-to-noise requirements for a velocity measurement limit the data set to the brightest stars. In turn, this sets the smallest radius inside which a velocity dispersion can be measured.

A significant improvement in data quality is possible with proper motion measurements. Such studies do not require much observing time and allow measurement even of the more numerous faint stars. Moreover, the ratio of the velocity dispersions in the tangential and radial proper motion directions directly constrains the velocity-dispersion anisotropy of the system. This is important because of the well-known mass-anisotropy degeneracy for stellar systems (Binney & Mamon 1982). To determine accurate proper motions close to the center requires an observing platform with high spatial resolution and long-term stability. The HST provides such a platform. McNamara, Harrison, & Anderson (2003) measured proper motions for M15, which were modeled by van den Bosch et al. (2006) and Chakrabarty (2006). This confirmed the presence of dark mass near the cluster center, but did not shed new light on the issue of whether this is due to mass segregation or an IMBH. McLaughlin et al. (2006) obtained and modeled proper motions for 47 Tuc, and used these to set an upper limit of 1000– on the mass of any IMBH in that cluster.

In Anderson & van der Marel (2009; hereafter Paper I) we presented the results of a new observational HST study of Centauri. For the central few arcmin of the cluster we determined the projected number density distribution of some stars, as well as accurate proper motions of some stars. The exact position of the cluster center was determined in three independent ways. First, we determined the center of the number density distribution using various methods, taking care to accurately account for incompleteness. Second, we determined the kinematical center, being the symmetry point of the RMS proper motion velocity field. And third, we determined the center of unresolved light seen in 2MASS data. All determinations agree to within 1–, yielding a final estimate with an uncertainty of .

The NGB08 center position, and hence the position of their central spectroscopic field, is at from the center determined in Paper I (and coincidentally, the second off-center spectroscopic field they studied is also at ). The differences between the centers determined in these two studies is likely due to a combination of subtle effects that biased the NGB08 analysis, as discussed in Paper I. As a result of this, the density cusp found by NGB08 was not measured around the actual cluster center. Also, we showed in Paper I that their use of unresolved light (as opposed to star counts) increased the impact of shot noise. In Paper I we determined the number density profile around the actual center of Centauri. It was found to be well matched by a standard King model, even though a shallow cusp may not be ruled out.

In Paper I we also determined the RMS proper motions of stars detected in the same fields studied spectroscopically by NGB08. For the canonical distance (vdV06), the results translate to one-dimensional values and , for their “central” and “off-center” fields respectively. So we did not confirm the NGB08 finding that one of their fields has a higher dispersion. The observation that there is no significant difference in proper motion dispersion between the two fields is consistent with our finding that they are both at the same radius from the actual cluster center.

Our results of Paper I did not confirm the arguments put forward by NGB08 to suspect the presence of an IMBH in Centauri. However, this does not mean that such an IMBH may not be present after all. The size and quality of our proper motion data set are superb. This provides the opportunity to study the central dynamics of Centauri at a level of detail that is unmatched by other clusters. We therefore present here a detailed study of the dynamics of Centauri, with the primary goal of exploiting the new data from Paper I to constrain the mass of any possible IMBH. We also include existing ground-based data in our data-model comparisons, to obtain the best possible results.

We restrict our modeling efforts in the present paper to spherical models. This is a simplification, given that Centauri is elongated and rotating (in fact, unusually so for a globular cluster). However, both the elongation (Geyer, Nelles, & Hopp 1983) and the rotation rate (vdV06) decrease towards the cluster center. Therefore, the assumption of sphericity is likely to be reasonable in the central region most relevant for constraining an IMBH. We find support for this throughout our paper by showing that spherical models yield an anisotropy profile, distance, and mass-to-light ratio that are in excellent agreement with the values obtained by vdV06. More sophisticated models, such as those of the type described by vdV06, van den Bosch et al. (2006) and Chaname et al. (2008), have the potential to constrain other aspects of the structure of Centauri, such as its inclination, rotation properties, and the presence of any kinematic subcomponents. However, spherical models are sufficient to obtain robust constraints on the central mass distribution.

The structure of this paper is as follows. Section 2 discusses the general modeling approach. Section 3 presents the data on star count and surface brightness profiles used for the analysis. Section 4 discusses parameterized fits to these data. Section 5 presents the kinematical data used for the analysis. Section 6 presents the kinematical data-model comparison, with a determination of the best-fitting model parameters. Section 7 discusses how the new IMBH results for Centauri fit in with our general understanding of IMBH demographics in globular clusters. Section 8 presents and discusses the conclusions.

## 2. Modeling Approach

### 2.1. Mass Density

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

(1) |

where 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 magnitudes of foreground extinction, the extinction-corrected profile is given by

(2) |

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

(3) |

To express this quantity in units of one needs to specify the cluster distance . Sizes in arcsec and pc are related according to

(4) |

The mass density of stellar objects and their remnants is

(5) |

where is the average mass-to-light ratio at each radius. Gradients in can be thought of as arising from two separate ingredients: (i) mass segregation among luminous stars; and (ii) 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 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 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 . 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

(6) |

Here is the mass density obtained under the assumption of a constant mass-to-light ratio . The excess is treated mathematically as a separate component. For we have found it convenient to use the parameterization

(7) |

corresponding to a Plummer density distribution. Here is the total excess mass in dark objects and is the characteristic scale length of their spatial distribution. The quantities , and provide a convenient way for exploring a range of plausible functions using only three parameters.

By assuming that is independent of , 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 years (McLaughlin & van der Marel 2005). First, Merritt, Meylan & Mayor (1997) constructed non-parametric dynamical models for Centauri. They found that at the large radii sampled by their data, 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 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.

### 2.2. Gravitational Potential

The total gravitational potential of the system is

(8) |

where is the gravitational constant and is the mass of any central IMBH. The gravitational potential corresponding to is

(9) |

The Plummer potential corresponding to is

(10) |

The potential 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.

### 2.3. Intrinsic Velocity Moments

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

(11) |

where is the number density of the stars under consideration. Because of the spherical symmetry, . The function measures the anisotropy in the second velocity moments, where . Models with are isotropic, models with are tangentially anisotropic and models with are radially anisotropic. The models may in principle be rotating if . However, for the purposes of the present paper we do not need to address how the second azimuthal velocity moment 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

(12) |

where is the average luminosity of a star in the subset and 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 (their kinematics are not probed) nor to (they emit negligible amounts of light). However, and can vary slowly as a function of radius due to mass segregation among the luminous stars. The quantity may also be affected by observational selection effects that vary as a function of radius. Nonetheless, variations in will generally be smaller than those in , as long as kinematics are measured for stars of similar brightness at each radius. Both and are only expected to have a shallow gradient over the scale of the entire cluster. These gradients are much smaller than those in , 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

(13) |

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

It follows from the preceding that one can substitute for 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 with solution

(14) |

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

(15) |

where is a characteristic anisotropy radius, is the central value of , and is the asymptotic value at large radii. If then is constant. With this , equation (14) reduces to

(16) | |||||

which is always positive.

The integral in equation (16) can be used with the equations given for and in Sections 2.1 and 2.2 to allow numerical evaluation of . This requires knowledge of the observed surface brightness profile and the known extinction . It also requires choices for the scalar model quantities , , , , , , and , which can all be chosen to optimize the fit to the available kinematical data.

### 2.4. Projected Velocity Moments

To calculate quantities that are projected along the line of sight, 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 (line of sight) 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 and the spherical polar coordinate system is , with , and corresponding to the pmt axis. The velocities in the (pmr, los, pmt) system satisfy

(17) |

Since in a spherical system, the second moments are

(18) |

where we have used that .

The projected second moments satisfy

(19) |

The derivations of these equations use the facts that and . The factor of two arises from the identical contributions from points with los and los . 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 line of sight. 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 , so the two cases are equivalent. The integrals in equation (2.4) can be evaluated numerically once the Jeans equation has been solved.

### 2.5. Projected Velocity Distributions

The normalized projected velocity distributions at projected radius 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 (being proportional to both and under the present assumptions) and , with the help of Eddington’s equation (Binney & Tremaine 1987). Because of the isotropy, the distribution is independent of the choice of the velocity direction (los, pmr, or pmt) and satisfies, in analogy with equation (2.4),

(20) |

Here is defined to satisfy , with . The quantity is the local velocity distribution (before projection along the line of sight) at radius , which can be written with the help of Eddington’s equation as (van der Marel 1994a)

(21) |

The integrals in equations (20) and (21) can be used with the equations given for and in Sections 2.1 and 2.2 to allow numerical evaluation of . The result is a symmetric function of [i.e., ]. This remains true if the distribution function is given a non-zero odd part in the angular momentum around some axis (corresponding to a rotating model), provided that is then interpreted as the distribution averaged over a circle of radius .

### 2.6. Observed quantities

Equations (2.4) 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 plane of the sky. For example, the second line-of-sight velocity moment over an aperture is

(22) |

Similar equations apply to the other velocity moments and to the velocity distribution . In integrated light measurements one can also take the observational point-spread-function (PSF) into account by convolving both and with the PSF before evaluation of the integrals in equation (22).

In the following we write

(23) |

It is also useful to consider the quantity

(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 and we measure proper motions and that satisfy

(25) |

Transformation to velocities requires an assumption about the cluster distance .

The overlying bar in the preceding equations is used to differentiate the root-mean-square (RMS) velocity from the velocity dispersion . These quantities are equal only if the mean streaming is zero (which is always true in the pmr direction).

### 2.7. Optimizing and Assessing the Fit to the Data

In general, we will want to model line-of-sight datapoints available for , as well as proper motion datapoints available for (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 and the extinction . 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.

#### Distance and Mass-to-Light Ratio

Suppose that we have obtained model predictions and for the datapoints, using trial values and 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 and be factors that scale the model predictions in the line-of-sight and proper motion directions, respectively. The best-fitting values of these scale factors and their error bars can be calculated using the statistics

(26) |

The best fit in each case is the value with the minimum (i.e., ), and the 68.3% confidence region corresponds to . The inferred scalings imply that the distance and mass-to-light ratio are optimized for and . These fits apply at the given values of the other model parameters, some of which are also affected by the scaling. The model parameters and are proportional to , and the model parameter is proportional to .

The error bars in and (at fixed values of the other parameters) follow from propagation of the uncertainties in and . The errors in and are uncorrelated, but the errors in and 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 , and can be determined as described in Section 2.7.2. They primarily affect the datapoints near the cluster center. So they can also be effectively eliminated by restricting the definition of the quantities in equation (2.7.1) to the datapoints at large radius. The anisotropy profile can be determined as described in Section 2.7.3.

#### Dark Mass

Any dark component in our models can be either a cluster with parameters or an IMBH of mass . The latter is mathematically (but not astrophysically) equivalent to a cluster with parameters . It is therefore not possible to discriminate between an IMBH and a dark cluster of sufficiently small size (within a resolution determined by the properties of the data). Conversely, we only need to explore the parameter space, and this will automatically tell us (at ) what the limits on an IMBH are. In the following we will use the term “dark mass” to refer generically to a dark component that could be either a dark cluster as in equation (7) or an IMBH. We have not explored models that have both a dark cluster and an IMBH, although that would not have been difficult.

At given , we start by calculating initial guesses for and by fitting only the data at large radii. We then calculate for a range of trial values the model predictions and for all the datapoints. For each pair of trial values we calculate and as in Section 2.7.1 by using the quantities and , including observations at all distances from the cluster center. We then map the contours of in the two-dimensional parameter space. The best fit is obtained as the point of minimum . The uncertainties of the fit follow from standard statistics, with the 68.3% confidence region corresponding to (Press et al. 1992). If there are combinations with within the confidence region, then a dark mass is not required to fit the data, for the assumed . If there are combinations with and within the confidence region, then the presence of an IMBH is consistent with the data, for the assumed . If we assume that the dark mass is an IMBH, then its 68.3% confidence region is given by at .

The presence of an IMBH not only increases the RMS velocity towards the center, but it also changes the shape of the velocity distribution. The distribution measured through an aperture around the center will have power-law wings that are significantly broader than those of a Gaussian distribution (van der Marel 1994b). This leads to the possibility of detecting very high velocity stars (Drukier & Bailyn 2003). It is therefore useful to determine and analyze the observed velocity distribution of the stars around the center, to further assess the plausibility of the presence of an IMBH. The observed distribution can be compared to the predicted distribution that can be calculated for an isotropic model. It is also useful to calculate the Gauss-Hermite moments (van der Marel & Franx 1993) of the observed velocity distribution as a function of distance from the center. A sufficiently massive IMBH causes an increase in towards the center (van der Marel 1994b). The absence of broad wings or a central increase in can be used to obtain an upper limit to the mass of any IMBH.

#### Anisotropy Profile

Binney & Mamon (1982) showed that there are many different combinations of and the enclosed mass in a spherical system that can produce the same observed projected RMS line-of-sight velocity profile . The key to constraining the mass distribution, and hence the possible presence of an IMBH, is therefore to constrain the anisotropy profile . This is not possible when only RMS line-of-sight velocities are available. However, can be constrained when information is available on deviations of the line-of-sight velocity distributions (LOSVDs) from a Gaussian, as 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 . 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 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 tightly constrains the one-dimensional intrinsic function . 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 , consider the case in which is independent of radius and . The relevant integrals in equation (2.4) can then be expressed analytically in terms of beta-functions. Their ratio is independent of radius and simplifies to

(27) |

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

(28) |

where we have defined the function

(29) |

in analogy with the definition of . This implies that the deviation of from unity is a factor of the deviation of from unity. (This uses the Taylor approximation , and similarly for .) Hence, is a “diluted” measure of the anisotropy that is present intrinsically, with dilution factor .

The approximations made above cannot be used in the core of a cluster, where . The assumption that the system remains scale free all along the line of sight 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 at large radii, while a Plummer model has . The corresponding dilution factors are and , respectively. Real clusters typically fall between these extremes. More extensive analytical calculations of 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 (2.4) 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 function of projected distance from the cluster center. Also, we can calculate the predicted velocity distributions 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.

### 2.8. Numerical implementation

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 and in equations (3) and (21), respectively, are calculated using finite differences. The derivative 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).

## 3. Projected Density Profile Data

The primary input for our models is the projected radial density profile. To characterize this profile over the full range of radii spanned by Centauri, we have combined the ground-based measurements at intermediate and large radii collected in Trager, King & Djorgovski (1995), with the new HST measurements at small radii presented in Paper I. Much of the measurements come in the form of projected star count profiles . Our equations of Section 2 are mostly expressed in terms of the observed (or implied) surface brightness profile . The profiles and are uniquely related to each other (see eq. [30] below), given our assumption to neglect mass segregation among the luminous stars.

### 3.1. Ground-based Data

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 -band surface brightness scale with a common zeropoint. We assigned to each data point an error bar of 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.

### 3.2. HST Data

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 -band instrumental magnitude (this is defined in Paper I and corresponds to stars measured with signal-to-noise ratio ). 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 () 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 to a projected surface brightness . It follows from equations (1), (4), (12), and (13) that

(30) |

where the zeropoint is equal to

(31) |

where (Binney & Merrifield 1998). Although can in principle be estimated, we have found it more reliable to calibrate 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 as a free parameter in the fits described below, which yields . For a canonical distance (vdV06) this implies that . With this there is excellent agreement between the different datasets, as shown in Figure 1a.

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 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 datapoints in the dataset used for the dynamical modeling.

## 4. Projected Density Profile Models

It is important for the modeling to have a smooth representation of the projected intensity profile , since the derivative is needed in equation (3) for the calculation of . For the dynamical modeling we have used a parametric approach, although this can in principle be done non-parametrically (e.g., Gebhardt & Fischer 1995). We have fit the data using a smooth function of the form

(32) | |||||

This functional form has no particular physical significance, but it is useful because it can fit a wide range of observed profiles. The profile has a power-law cusp at radii . It then breaks to a logarithmic slope at , and further to a logarithmic slope at . The parameters and measure the sharpnesses of the breaks. The intensity , corresponding to a magnitude through equation (1), sets the overall intensity scale. The profile is similar to a so-called “nuker” profile (Lauer et al. 1995), but with two breaks instead of just one. We will refer to it as a generalized nuker profile.

The best fit to the dataset presented in Section 3 was determined using a Levenberg-Marquardt iteration scheme (Press et al. 1992). It has parameters , , , , , (the maximum value allowed during the iteration), , and . The fit and its residuals are shown in Figure 1. The quality of the fit is excellent. For the subset of the data from Paper I, the for 125 datapoints with an RMS residual of 0.04 mag. The total for all 159 datapoints, corresponding to 151 degrees of freedom. This is acceptable at the level for a distribution. The residuals appear generally consistent with random scatter in line with the error bars, with little evidence for systematic trends.

The best fit has a central cusp slope of . To determine the statistical error on this value we also performed a range of fits in which was fixed a priori, while all other parameters were varied to optimize the fit. The central structure of our best fit for is shown in Figure 2 for comparison to our overall best fit. The curve of as function of was used to derive the error bar on , which yields .

When using a generalized nuker fit, a constant surface brightness core is inconsistent with the data at confidence. However, this parameterized approach may not yield an unbiased estimate of the asymptotic slope at the smallest radii. The data at all radii are given equal weight in the of the fit, but it is obviously only the data near the very center that contain most of the relevant information on . Also, the two breaks at radii and in the generalized nuker profile fit correspond to the breaks associated with the traditional core and tidal radii of Centauri. The parameterization does not have sufficient flexibility to reproduce any possible change in profile slope well within the core radius ( arcsec; McLaughlin & van der Marel 2005).

An alternative and simpler approach to determine is to fit a straight line to the data at . For , this yields . The values of inferred for other values of in the range – are in statistical agreement with this. This analysis indicates that the central data points are consistent with a flat core.

An additional assessment of models with a flat core can be obtained from comparison of the new data with previously published models. McLaughlin & van der Marel (2005) presented detailed surface brightness fits for a large sample of galactic and extragalactic globular clusters. They considered three types of models with decreasing amounts of diffuse outer structure, the power-law model, the Wilson (1975) model, and the King (1966) model. The power-law models are similar to the parameterizations that we have used here, but with a more restricted subset of the parameters (namely, , , and were kept fixed). The models from Wilson and King are both based on a parameterized isotropic phase-space distribution function. The Wilson prescription differs from the better known King models in that it uses a slightly different lowering of the exponential energy distribution, so as to produce a model that has a more extended but still finite outer structure. McLaughlin & van der Marel discussed Centauri as a special example (their section 4.2.1) and found that it is well fit by a Wilson model, but less so by a King model. The central structure of both models is shown in Figure 2. The Wilson model is very similar to our best fit model with fixed . The King model is also similar at radii , but it under-predicts the brightness at intermediate radii (before improving again at larger radii, not shown in Figure 2).

The Wilson and King models shown in Figure 2 are the ones that were previously fit by McLaughlin & van der Marel to the Trager et al. data. At small radii these models represent an extrapolation of the ground-based data under the assumption of a constant surface brightness core. Although the models were not fit to the HST data, they do provide a perfect fit to the two innermost HST star count data points. By contrast, our own best generalized nuker model fit with over-predicts both of the innermost datapoints by about . The fact that is preferred in our fits over is therefore because it fits better the data at radii , and not because it provides the best fit closest to the center.

NGB08 inferred from their HST study of integrated light in Centauri that it has a central cusp slope . At first glance, this result appears consistent with that from our generalized nuker fit. However, since they did not measure the brightness profile around the same cluster center, the agreement between our inferred cusp slopes is likely coincidental.

In summary, we conclude that the central cusp slope of Centauri is at confidence. That is, we interpret the from our generalized nuker fit as an upper limit. We do this because it over-predicts the central two data points, and because other statistics do not show evidence for a non-zero slope. In the dynamical analysis that follows we specifically consider two generalized nuker fits that span the range of possibilities: models with a shallow cusp () and models with a constant density core (). We will refer to these as cusp models and core models, respectively. In all calculations we have transformed the observed intensity to an extinction-corrected intensity using equation (2) with for Centauri (McLaughlin & van der Marel 2005).

## 5. Kinematical Data

To constrain the parameters of the dynamical models we need to use both line-of-sight and proper motion data at a range of radii. For this we have combined ground-based line-of-sight and proper motion data at intermediate and large radii with the new HST proper motion data at small radii presented in Paper I.

We construct spherical models to interpret the data. These models yield only a single radial profile for each projected quantity, with no azimuthal variations from the observed major axis of Centauri to its observed minor axis. To get accurate results for the model parameters, it is important that these models be used to reproduce the average radial variations seen in the data. This is achieved by using for each projected quantity that is fit the average of the measurements along concentric annuli. This is indeed what we extracted for the projected density in Section 3, and we therefore do the same here for the projected velocity distributions. When averages are taken over annuli, all the odd moments of the velocity distribution are zero. Of primary interest here are the second moments, which yield the RMS velocity and the RMS proper motions and .

In reality, variations do exist between observed quantities along the major and minor axis. vdV06 discussed this in detail for the existing ground-based data. In Appendix B we address this for the new HST proper motion data. In particular, the RMS proper motions are not the same along the major and minor axes, and the RMS proper motions in the major and minor axis directions on the sky are not the same. Also, the mean proper motion in the transverse (pmt) direction is not zero (vdV06; Section 5.2.3; Appendix C). By modeling these issues it is possible to constrain, e.g., the cluster inclination and its rotation rate. However, it is not essential to model these issues to achieve a good understanding of the central mass distribution, which is the main goal here.

### 5.1. Ground-based Data

#### Line-of-Sight Kinematics

vdV06 compiled line-of-sight 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 Gebhardt et al. (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 line-of-light 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 line-of-sight velocities for 2163 stars with uncertainties below .

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 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 from the center used by vdV06 (which was determined by van Leeuwen et al. 2000). We chose a central bin of 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 line-of-sight 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 . Seitzer (1983) presented the profile derived from 118 clusters 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 line-of-sight velocities of 75 cluster members at large radii. From this they determined the velocity dispersion 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 values the average measurement error of . To obtain we added in quadrature , where is the rotation velocity amplitude quoted in Scarpa et al. (2003); the factor 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 (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 . 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 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 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 rather than RMS velocities . However, it is known from previous observations (e.g., vdV06) that in Centauri decreases towards the center. This ratio is small enough at the position of the NGB08 fields that to high accuracy (see Section 5.2.3 for a quantitative estimate).

Figure 3a shows the collected data from all the sources, as function of projected distance from the cluster center. The different data sets match well in their regions of overlap.

#### Proper Motion Kinematics

Ground-based proper motion data are available from van Leeuwen et al. (2000). vdV06 performed various operations and cuts to optimize and correct the data as necessary for dynamical modeling. This included: improved removal of cluster non-members; removal of stars with low enough brightness, large enough velocity error bars, or neighbors that are nearby enough to make their reliability suspect; and subtraction of the apparent solid-body rotation introduced by perspective rotation, and other spurious solid body rotation components. This produced a final sample of proper motions for 2295 stars with uncertainties below (i.e., at the canonical distance ).

Glenn van de Ven kindly made the homogenized sample from vdV06 available to us for use in our data-model comparisons. Proper motions were given along the major and minor axis directions. For each star we transformed this into proper motions in the radial (pmr) and tangential (pmt) directions. We then determined the profiles of and by binning the stars in radius and using the statistical methodology described in Appendix A. Radii and proper motion coordinate transformations were calculated with respect to the center determined in Paper I. We chose a central bin of radius (with 30 stars) for maximum resolution near the center, and then chose the remaining bins to contain 200 stars each. We followed vdV06 by excluding the 65 outermost stars with radii –30 arcmin.

For fitting models to the data we always use the quantities and . However, for visualizing the results we have found it more convenient to use the quantities and . The RMS proper motion averaged over all directions in the plane of the sky can be calculated directly through its definition in equations (24) and (25). In practice, we found it more convenient to calculate it by applying the statistical methodology described in Appendix A on the combined array of the pmr and pmt data points. For visualization of the ratio we adopted a radial binning scheme with 400 stars in each bin, to obtain smaller and more useful error bars.

Figure 3b shows the data as function of projected distance from the cluster center, and Figure 3c shows the data. The velocity distribution is close to isotropic near the center, but with some radial anisotropy. There is a significant increase in tangential anisotropy at a few core radii. The latter is mostly due to the contribution of rotation to the second moment . The anisotropy in the proper motion dispersions is much less (King & Anderson 2002).

#### Dependence of Kinematics on Stellar Mass

Luminous stars in a cluster undergo two distinct effects as a result of two-body relaxation: mass segregation and energy equipartition. These effects are intimately connected through basic dynamical theory, but from an observational perspective it is useful to think of them as separate effects. Mass segregation causes stars of different masses to have slightly different spatial distributions. In Section 2 we assumed that this could be neglected when parameterizing the run of mass-to-light ratio with radius, and when connecting star count data at small radii to surface brightness data at large radii. This assumption was motivated by the absence of high-quality data to constrain the amount of mass segregation over the full scale of the cluster, and by the desire to minimize the number of a priori model assumptions. Energy equipartition causes stars of low mass to move faster than those of high mass. We measured in Paper I that low-mass stars in Centauri do indeed move faster than high-mass stars, albeit by less than would be predicted for complete energy equipartition (). Since we have a direct measurement of this, there is no reason to ignore this effect in the modeling. We therefore correct for it when comparing the observed kinematics from different data sets.

The available discrete ground-based line-of-sight and proper motion data all pertain to relatively bright stars. The apparent magnitudes of these stars place them on the giant or sub-giant branch. Since stellar evolution proceeds relatively rapidly after the main-sequence turn-off, these stars all have essentially the same mass as the main-sequence turn-off. By contrast, the stars in our new HST proper motion sample from Paper I (discussed in detail in Section 5.2 below) extend from about the main-sequence turn-off to magnitudes below that. From the analysis of Paper I it follows that the RMS proper motion of stars at the average magnitude of the HST sample is % higher than that at the main-sequence turn-off.

For a fair comparison to the new HST proper motion data, we took the kinematics from the discrete velocity studies in Sections 5.1.1 and 5.1.2 and multiplied them by a factor . This implies that we use the observed (sub-)giant star kinematics to estimate the kinematics of lower-mass stars in the same general population. We assume implicitly that the variation of RMS kinematics with mass is the same everywhere in the cluster (although it was only measured near the cluster center), and that it applies also in the line-of-sight direction (although it was only measured in the proper motion direction). In reality, one would expect the relation between mass and velocity dispersion to change with radius, since relaxation proceeds faster near the cluster center than further out. This was neglected in the present context by applying the same small correction throughout the cluster. This was motivated in part by the fact that our study focuses mostly on the central part of the cluster anyway, but also by the fact that we didn’t actually detect a difference in equipartition in Paper I between the different (central and major axis) fields for which we obtained proper motions. We did not apply the multiplicative correction to the NGB08 data. Their measurements are based on integrated light, excluding the brightest regions in their fields. It is therefore not obvious that their measurements would be dominated by (sub-)giant stars. Either way, a % correction would be much less than the error bars on their datapoints.

In summary, we are using our Jeans models to fit the kinematics of stars with masses typical of our HST proper motion sample. We extend our actual HST kinematics to larger radii using ground-based kinematics that have been corrected to the same characteristic mass. We model the cluster density profile using HST star counts for a sample that is also centered on the same characteristic mass (see Section 3.2). We make the simplifying assumption that there is no mass segregation among the luminous stars only so that we can extend the HST density profile to larger radii, where it can be tied to ground-based star-count data and integrated photometry that are based largely on giant stars (see Section 3.1). Since this approximation does not affect the very central region of the cluster, it should have little effect on our results for the central mass distribution. In fact, the multiplicative correction of discussed above also has little effect on this. We verified explicitly that none of the main conclusions of our paper change at a level that exceeds the formal uncertainties when the multiplicative correction is omitted.

### 5.2. HST Proper Motion Data

#### Sample selection

In Paper I we derived proper motions for stars observed in two fields observed by HST, a “central field” on the cluster center, and a “major axis field” positioned adjacent to it in non-overlapping fashion roughly along the cluster major axis. Here we use the “high-quality” sample presented in Paper I, which contains the non-saturated stars that are isolated enough and bright enough (-band instrumental magnitude ) for a particularly reliable proper motion measurement. This sample has 53382 stars in the central field and 19593 stars in the major axis field. In this sample, 95.1% of the stars have proper motion errors below . Our cut in proper motion accuracy is therefore not very different from that applied by vdV06. However, the median proper motion error per coordinate in the final sample of vdV06 is , whereas it is a factor two smaller, , for our sample. We note that our use of only the high-quality sample is rather conservative. Even the remainder of the sample of Paper I is quite accurate, and could well have been used to constrain the cluster dynamics.

The central HST field covers radii , while the major axis field covers radii . The central field has complete position angle coverage out to (the radius of the largest enclosed circle). For the coverage over the range of position angles becomes increasingly sparse. The major axis field is restricted at all radii to position angles within of the major axis. For the analysis of the present paper it is important to have coverage of all position angles at a given radius, so that average kinematical quantities over circular annuli can be calculated. For this reason we have used only the 25194 stars at to constrain the dynamical models.

In our analysis we ignore the 47781 stars with high-quality HST proper motions at radii . We do discuss relevant quantities derived from these data in Appendix C. More sophisticated axisymmetric modeling techniques (such as those in vdV06; van den Bosch et al. 2006; and Chaname, Kleyna, & van der Marel 2008) that might be applied in the future probably would not want to exclude these data in their analysis. By contrast, to include these data in our spherical models would require estimates of the ratios of the average of either or over a restricted range of position angles, to its average over a circle. We have experimented with such estimates, but concluded that the approximations that need to be made introduce sufficient uncertainties that these data do not really help to constrain our models. Conversely, we have found no evidence that inclusion of the HST data at , combined with reasonable estimates for the ratios , would alter in any way the conclusions draw below (see Appendix C).

To reject cluster non-members we created scatter plots of total two-dimensional proper motion versus radius , where is the proper motion vector of an individual star. Such a plot reveals a sparse population of field stars at large values of , in addition to the numerous cluster stars at low . We rejected those stars from the sample that reside at and for which . At the radii for which we have HST proper motions, this equation is a good approximation to the escape velocity curve for a spherical model with a canonical distance and mass-to-light ratio (, ; vdV06). With this criterion, only 27 non-member stars were identified at , the closest of which resides at