Testing for Dark Matter in the Outskirts of Globular Clusters

The proper motions of stars in the outskirts of globular clusters are used to estimate cluster velocity dispersion profiles as far as possible within their tidal radii. We use individual color-magnitude diagrams to select high probability cluster stars for 25 metal-poor globular clusters within 20 kpc of the sun, 19 of which have substantial numbers of stars at large radii. Of the 19, 11 clusters have a falling velocity dispersion in the 3-6 half mass radii range, 6 are flat, and 2 plausibly have a rising velocity dispersion. The profiles are all in the range expected from simulated clusters started at high redshift in a zoom-in cosmological simulation. The 11 clusters with falling velocity dispersion profiles are consistent with no dark matter above the Galactic background. The 6 clusters with approximately flat velocity dispersion profiles could have local dark matter, but are ambiguous. The 2 clusters with rising velocity dispersion profiles are consistent with a remnant local dark matter halo, but need membership confirmation and detailed orbital modeling to further test these preliminary results.


INTRODUCTION
Whether or not globular clusters are surrounded by locally bound dark matter sub-halos as they orbit within the galaxy is an important observational question that will test theories of globular cluster origins and to some degree the nature of dark matter (Peebles 1984). Star clusters forming in a self-gravitating galactic gas disk rotating within a dark matter halo are not expected to have an excess of dark matter above the galactic background dark matter. The disk population of globular clusters (Searle & Zinn 1978) is relatively metal rich, and younger than the halo clusters (VandenBerg et al. 2013). However, even the metal poor halo clusters may have formed as the high mass, high gas pressure, regime of low metal gas disk star cluster formation at high redshift (Portegies Zwart et al. 2010;Sun et al. 2018). The halo clusters are very old (VandenBerg et al. 2013), likely forming in pre-galactic sub-halos which later merged into the Milky Way's halo (Searle & Zinn 1978).
One possibility is that some halo clusters formed at or near the center of small dark matter halos. Tidal fields in an assembling galaxy later remove the outer dark matter, but clusters can retain dark matter within their tidal radii to the present day (Ibata et al. 2013;Boldrini raymond.carlberg@utoronto.ca carl@ipac.caltech.edu et al. 2020;Carlberg & Keating 2021). An example is the cluster NGC 6715 (M54), located at the center of the Sagittarius dwarf galaxy, an accreted dwarf galaxy being dispersed into the Milky Way (Ibata et al. 1997;Vasiliev & Belokurov 2020). The velocity dispersion of M54 begins to rise at about 5 half mass radii (Bellazzini et al. 2008;Ibata et al. 2009) likely due to the dwarf galaxy's remnant dark matter halo. M54 is not typical, since it qualifies as a nuclear star cluster, having both an extended old stellar population and younger stars in the inner region (Alfaro-Cuello et al. 2020). More generally, the lower mass nuclear star clusters are considered to be a globular cluster (or the merger of several) in the centers of dark halos Neumayer et al. (2020). That is, some massive old globular clusters clearly were formed in the central regions of dark halos. On the other hand, the well studied Milky Way halo globular clusters are well fit with purely stellar models (Baumgardt & Hilker 2018;Vasiliev & Baumgardt 2021) which provide the dynamically expected mass segregation which explains the modest rise in mass-to-light in the outer regions. Nevertheless a significant amount of dark matter is still allowed in the outskirts, though limits on the mass remain poorly constrained (Ibata et al. 2013).
The proper motions of individual globular cluster stars is coming within the range of observational capabilities (Baumgardt 2017;Baumgardt & Hilker 2018;Baumgardt et al. 2019;Sollima 2020;Wan et al. 2021). Current observational results show no evidence for any dark matter above the Galactic background, though these results are generally limited to 5 half-mass radii or less. A primary issue is to be able to discriminate cluster members from background stars. A second issue is the development of dynamical models that predict stellar velocities to large cluster radii.
Numerical modeling shows that extended stellar density profiles arise as a result of internally driven expansion and the heating from external galactic tides (Spitzer 1987;Fukushige & Heggie 2000;Daniel et al. 2017;Tiongco et al. 2016;Peñarrubia et al. 2017). Beyond the tidal radius the outer region loses stars into tidal streams (Odenkirchen et al. 2003;Grillmair & Dionatos 2006;Kuzma et al. 2018). Developing a model which describes the entire range and its orbital phase dependence remains an active area of research. This paper analyzes the available Gaia EDR3 data (Gaia Collaboration et al. 2016 for metal poor halo globular clusters sufficiently nearby to have good internal velocities from proper motions. The main innovations are color-magnitude weighting of the data and comparison to n-body simulation results that form globular clusters at various locations of high redshift sub-halos, including the sub-halo centers. This paper is a preliminary investigation using currently available data, straightforward measurement methods and a simple comparison of models to the data.

GLOBULAR CLUSTER SELECTION
We select Milky Way clusters that have [Fe/H] below -1 and are within 20 kpc of the sun as promising targets for dark matter in the outskirts. The more distant clusters will have proportionally larger proper motion uncertainties, but we include them to determine where the falloff in internal kinematic measurement precision occurs for the currently available data. Setting a Galactic latitude limit of at least 20 degrees minimizes extinction and extinction variations with position. Colormagnitude filtering requires accurate colors and brightnesses, which can be significantly compromised by high and variable extinction. Higher latitudes also greatly reduce the surface density of contaminating field stars.
The cluster sample is listed in Table 1 along with some of their properties as listed in Baumgardt et al. (2019) and available at https://people.smp.uq.edu.au/-HolgerBaumgardt/globular/. The fifth column is a qualitative assessment of the shape of the outer velocity dispersion-radius relation in the range of 3-6 half-mass radii in Figures 3-12 below. The last two columns give the current orbital radius relative to the apocenter and the orbital eccentricity. The clusters in the table are in the same order as the plots below. Stars are selected using a matched filter technique, which depends on how close a star is to a colormagnitude locus at the distance and metallicity of the cluster (Rockosi et al. 2002;Grillmair 2009). Specifically, we created a color-magnitude locus for each cluster using the observed Gaia G, G BP − G RP measurements, extinction corrected using the reddening maps of Schlegel et al. (1998), themselves corrected using the prescription of Schlafly & Finkbeiner (2011). We employed theoretical isochrones with appropriate metallicities from http://stev.oapd.inaf.it/cgi-bin/cmd (Girardi et al. 2004), though we adjusted these isochrones to better match the observed red giant branches, where the photometric uncertainties are particularly small. The color offset from the color-magnitude relation normal- Figure 1. The color-magnitude diagram of the extinction corrected stars in the NGC 6752 cluster with less than 1σ color offset used here are shown as blue points. The kinematic mixture modeling (Vasiliev & Baumgardt 2021) 90% probability cluster members, with no extinction correction, are shown as red points. Our color-magnitude selection does not include horizontal branch stars.
ized to the quoted photometric uncertainties of every star less than 2σ from the locus is included as a measurement. Only stars brighter than G = 20 mag are included in our sample. Figure 1 compares our sample of stars for NGC 6752 to the 90% probability cluster stars as found in a kinematically-based Gaussian mixture model (Vasiliev & Baumgardt 2021) for stars between projected radii of 5 and 100 pc of the cluster center. Stellar crowding in the inner 5 pc increases the photometric errors and reduces the limiting magnitude (Pancino et al. 2017;Riello et al. 2021), but these innermost stars are largely irrelevant to the goals of this paper. The horizontal branch stars can not be sufficiently narrowly defined to be usefully included in the photometric model of the cluster, but their numbers are sufficiently small that little is lost. The effectiveness of our photometric selection technique in reducing the number of unrelated field stars is shown in Figure 2 where we plot the ratio of the total number of stars in a radial bin the number of stars photometrically classified as cluster members at 1σ.
Every star in the sample has an estimated error in each of the two components of proper motion, which we combine in quadrature for a total proper motion error. Only stars with proper motion errors less 0.15 mas/year are included in the analysis. The proper motion variance is also used with weighted measures of the velocity dispersion. The color offset is used to construct a Gaussian weight exp [−(c/σ c ) 2 /2], where σ c is set at 1.5 which gives a weight at the color edge of 0.41. The same analysis with σ c = 1 gives very similar results. For weighted velocity dispersion measurements, the color weight and the proper motion weight are multiplied together to give a total weight for each star. The radial component of the proper motion velocity is corrected with the viewing angle velocity, ∆v r = V r · r/D , where V r is the radial velocity of cluster measured from the sun at distance D and r is the projected radial distance of the star from the cluster center. The correction has no effect on the velocity dispersion.

VELOCITY DISPERSION PROFILE ANALYSIS
The Gaia proper motions of the selected stars are converted to cluster-centric velocities using the cluster distances and center of mass velocity components of Baumgardt et al. (2019). For simplicity, the two components of the velocity are added together to produce a total velocity in the plane of the cluster. While there is useful information in the velocity anisotropy at large distances, its measurement beyond about 3 half-mass radii is very uncertain with the current data. These clusters have been previously shown to have nearly isotropic planar velocities and relatively little rotation (Vasiliev & Baumgardt 2021). All velocities are measured with respect to the cluster center and the published velocity of the clus-  ter. Our analysis assumes that all the color-magnitude selected stars are at the same distance. The mean radial and tangential velocities are 1-2 km s −1 , generally consistent with zero and small enough to have no dynamical significance.
Potential issues in estimating the velocity dispersion are the presence of unbound background stars and velocity errors that depend on brightness. Although the clusters do exhibit mass segregation (Ebrahimi et al. 2020) any variation of velocity dispersion with stellar mass should be small, since the cluster stars bright enough to have high precision proper motions are generally near the main sequence turn-off and above, so have comparable masses even though their luminosities vary substantially. The maximum velocity for inclusion in the cluster has no radial dependence and is chosen in the range of 20 to 50 km s −1 based on the velocity dispersion of the cluster, with the velocity indicated as the maximum of the velocity range of the scatter plots shown in the left-hand panels of the velocity plots. The velocity cuts    were set to be at least 4 times the velocity dispersion at about 30 pc. The velocity cuts are at least 2.35 times the half-mass circular velocity, GM c /(2r h ), which is above the escape velocity in the outskirts of the cluster. The color-magnitude selection and the velocity cut reduce the number of contaminating field stars in the clusters to very low levels.
The two clusters with the largest remnant background counts, NGC 6254 and NGC 6752, have 6 background stars between 159 and 200 pc of the center. The background can be assumed to have a uniform surface density distribution of stars, which implies that in the radial range from 40 to 50 pc there should be 3/8 of a background star. We actually find 6 stars meeting our selection criteria, a factor of 16 above the estimated background. For the bin extending from 32 to 40 pc we would expect 0.24 background stars. The 16 stars meeting our criteria therefore exceed the estimated background level by a factor of 66. Nevertheless, the higher velocity stars within the tidal radius should be prime targets for spectroscopic membership tests. Inset in the velocity dispersion plot is the surface density profile of the selected stars within the velocity cut and proper motion accuracy cut. The surface density of cluster stars within about half the tidal radius is at least an order of magnitude above the background at large radius. The surface density plot also shows that crowding diminishes the numbers of selected stars at small radius, though the accuracy of the velocity dispersions at small radii is not important for this paper.
The simplest approach to measuring the velocity dispersion is simply to calculate the variance of the velocities, with and without weighting of the stars. The error in the velocity dispersion σ is estimated as σ/ √ 2N . We tried a double Gaussian fit to the cluster stars and the background, but there are so few background stars predicted within half the tidal radius that the results were essentially identical to the two simple velocity dispersion estimates. The velocity dispersion measurements are displayed in Figures 3 to 12. The plots show the data with the unweighted velocity dispersion in the left panel and both measures of the velocity dispersion profile along with model curves in the right panel. The last column of Table 1 gives our assessment of whether the velocity dispersion is rising, falling, flat, or no stars at a projected distance in the range extending to around 6 times the half-mass radius. Of the 25 clusters there are two with rising velocity dispersion profiles, six with no stars at large radius, eleven with falling dispersion profiles, and six that are approximately flat.
For comparison we plot projected velocity dispersion profiles for a W=7 King model (King 1966), with veloc-ities calculated from points generated with the McCluster model-maker (Küpper et al. 2011). The King model gives densities and velocities that go to zero at a finite radius, hence is not expected to be a good representation of the outer density or velocity profile of a cluster that is undergoing slow mass loss. We also plot two models from a simulation of globular clusters in a cosmological model (Carlberg & Keating 2021). The green line shows the relation for a cluster with no dark matter other than the galactic background. The blue line is for a typical cluster at the center of its small dark matter halo. The dark matter mass density is nearly constant with radius over this range so there is essentially no detectable dark matter within the stellar cluster, but it begins to dominate at about 10 half mass radii. Beyond the tidal radius the dark matter is subject to tidal removal. The cluster simulations were projected on the sky and measured in the same as the observational data. The simulations find essentially a continuous range of possibilities for the large radius velocity dispersion profile. The model profiles are normalized to the velocity dispersion at the half-mass radius or nearby. There is no fitting of the models to the data at this stage.
The velocity dispersion profiles of NGC 6752 and NGC 6205 exhibit a rise in the velocity dispersion in the range 30-50 pc. In both cases the velocity anisotropy becomes radial at the onset of the rise, agreeing with Baumgardt et al. (2019) for the radii in common. There is no clear indication in the current data that the rise in velocity dispersion in the 30-50 pc range is the result of background stars. For instance, for NGC 6205 the rise in velocity dispersion at 30-40 pc is dependent on three stars with velocities between 19 and 31 km s −1 . However the density of all background stars in this radial range is a factor of 80 lower, making these stars fairly high probability members. Furthermore, the background has no stars in this range of velocities. A concern is that the rising velocity dispersion at large radius is not present in stars brighter than G = 17 mag as shown in Figure 13, however this is entirely because there are very few bright stars at larger radii, likely a result of mass segregation (Ebrahimi et al. 2020). That is, the bright stars do not sample the large radius velocity dispersion in these clusters.

DISCUSSION
Testing for the presence of dark matter around a globular cluster requires pushing kinematic measurements to as close to the tidal radius as possible where the dark matter density of a cosmological halo may begin to dominate over the stellar density (Ibata et al. 2013;Carlberg & Keating 2021). Many globular clusters with high qual- ity Gaia globular cluster proper motion have been modeled with great success inside 3-5 half mass radii, with no need to include a significant distributed dark matter component (Baumgardt et al. 2019;Vasiliev & Baumgardt 2021). However, for a typical cluster with a tidal radius of 100 pc, conclusively ruling out a local dark matter halo requires kinematic measurements at larger radii, where the surface density of cluster stars drops below the that of Galactic field stars. The approach here is to use a color magnitude diagram carefully tuned to each cluster to identify likely cluster stars, which are further culled of background stars using a velocity cut that eliminates stars in the outskirts that are not bound to the cluster.
The outer velocity dispersion profile depends on the orbit of the cluster in the galaxy. Clusters on noncircular orbits are subject to regularly varying tides, which lead to a regular orbital phase dependence of the velocity dispersion at large radius, whether the cluster is embedded in dark matter halo or not. That is, as a cluster passes the orbital pericenter, stars are quickly accelerated, with the velocity boost increasing outward in proportion to their distance from the cluster center. The stars move outwards settling into a new equilibrium after apocenter (Binney & Tremaine 2008), with the unbound stars leaving in a tidal stream.
The n-body simulations found that the velocity dispersion anisotropy is an indicator of the presence or absence of dark matter (Carlberg & Keating 2021). That is, the velocity dispersion ellipsoid of cluster without dark matter tends to be quite radial in the outer region (Bianchini et al. 2017), but nearly isotropic if dark matter is present. Figure 14 displays the 2D version of the velocity anisotropy parameter, β 2D = 1 − σ 2 t /σ 2 r . The clusters identified as rising profile are conspicuously near zero velocity anisotropy, whereas the other clusters have a significantly larger spread in β 2D , with no clear indication of an offset to positive β 2D , (radial anisotropy). All of the clusters in our sample have been previously studied at smaller radii using a different measure of ve- Figure 14. The velocity anisotropy parameter for the different classes of velocity dispersion profile plotted as a function of the radial distance normalized to the half mass radius circular velocity. Only radial bins with at least 20 velocities are plotted.
Clusters that have velocity dispersions in their outskirts that rise above the expectations of a selfgravitating cluster of stars are not new. Baumgardt & Hilker (2018) show results for 59 clusters with more than 100 stars with kinematics. Of those, 8 have velocity dispersion profiles at large radius that are significantly above the self-gravitating model, and NGC 6715/M54, has a rising velocity dispersion profile from 5 to 8 halfmass radii, a consequence of being the at the center of the Sagittarius dwarf (Ibata et al. 1997(Ibata et al. , 2013Vasiliev & Belokurov 2020).
We selected 25 nearby halo clusters, of which 6 have effectively no stars at large radii, hence with no evidence for or against a local dark matter halo. Of the remaining 19 clusters, 11 have a falling velocity dispersion with cluster radius. These clusters are consistent with having their origin in an environment with no dark matter above the galactic background at a radius of a few kpc in a high redshift dark matter halo. The 6 velocity dispersion profiles that are approximately flat at large radius are ambiguous at present. The Carlberg & Keating (2021) n-body simulations find that both clusters centered in their dark halos as well as significantly offset clusters lead to flat velocity dispersion profiles at large radius.
We find that two clusters, NGC 6205 and 6752, plausibly have a rising velocity dispersion at large radius. The rise is similar to what is seen in n-body simula-tions of star clusters with local dark matter halos and is tentative evidence of local dark matter halos for these two clusters. The rise in velocity dispersion profiles is subject to the measurement problem of interloper background stars. For these two clusters the velocity dispersion profile falls for stars brighter than G = 17 mag, but rises for stars between G = 17 and 20, which are more radially extended as expected in these mass segregated clusters.
NGC 6205 and 6752 are at Galacto-centric distances, R GC , of 8.3 and 5.2 kpc, respectively, with estimated peri-Galacticons of 1.55 and 3.23 kpc, and apo-Galacticons of 8.32 and 5.17 kpc, respectively (Baumgardt et al. 2019). Estimating the tidal radius for a cluster of mass M c as r t = R GC [GM c /(2V 2 c )] 1/3 (Webb et al. 2013) where V c = 220 km s −1 is the assumed inner galaxy circular velocity (Bovy et al. 2012), gives tidal radii of 113 and 65 pc for NGC 6205 and 6752, respectively. These tidal radii are well beyond the radii at which the velocity dispersion profiles begin their rise. Cluster membership criteria will be strengthened with additional photometric and/or spectroscopic data to confirm that all the stars share the chemical abundance signatures of the cluster stars. Future Gaia data releases will also greatly reduce the kinematic errors. This increased accuracy will also enable more distant clusters to be examined for evidence of local dark halos.
The cluster sample has a mean radial position relative to their apocenter of 0.78 and a mean orbital eccentricity of 0.64. The orbital distribution of the sample is almost certainly missing a few members closer to pericenter where they are below our Galactic latitude cut. Both NGC 6205 and 6794 have a negative radial velocity with respect to the Galactic center so they are just past their orbital apocenters, in which case we expect there to be minimal or no radial flow within the tidal radius, which is observed., although stars beyond the tidal radius NGC 6752 have an outward flow of as much as 5 km s −1 . Near pericenter a rise in velocities at large radius is much more likely than just past apocenter, where both of these clusters are located. Nevertheless, the complications of orbital phase effects and the possibility that background stars contaminate the outer part of the velocity dispersion measurement means it is not yet possible to say whether the rising velocity dispersion profiles are the result of a local dark matter sub-halos or are the response of the outer stars in an eccentric galactic orbit.
A local dark matter sub-halo would deepen the gravitational potential of a cluster and could have some bearing on the multiple stellar populations of globular clusters (Gratton et al. 2012). Detailed abundance studies do show some differences from field stars in NGC 6752 (Norris et al. 1981) but no clear differences for NGC 6205 (Carretta & Bragaglia 2021). NGC 6752 stands out in a comparison of widths of the asymptotic giant branch and the red giant branch, but NGC 6205 does not (Lagioia et al. 2021). We conclude that there is no clear correlation of cluster metallicity trends with velocity dispersion profile within our small sample of clusters. In summary, the velocity dispersion measurements presented here provide consistent, but not compelling, evidence for the presence of remnant dark matter within the tidal radii of 2 of the 19 clusters, with a further 6 clusters with flat profiles being ambiguous.

CONCLUSIONS
Determining which stars are globular cluster members is the central issue of measuring their velocity dispersion at large radii. The approach here is to use a color-magnitude relation tuned to the high precision Gaia photometry, which has the advantage that the primary membership criterion is completely independent of cluster kinematics. The method has been proven to work well in uncovering faint stellar streams (Grillmair 2009), but is even more powerful for globular clusters with known distances and large numbers of stars in the central region which can be used to accurately calibrate the color-magnitude diagram. The color-magnitude cut is augmented with a velocity cut that is at least the 2.35 times the scale velocity of the cluster, the circular velocity at the half mass radius, GM c /(2r h ). The colormagnitude cut reduces the numbers of stars at least a factor of 2 and more typically a factor of 5 in the radial range of interest, as shown in Figure 2. Of the 25 clusters in our sample, 19 have sufficient numbers and accuracy of proper motion velocities at large radii to give a velocity dispersion profile to 5 half mass radii or more. Of the 19, 11 have velocity dispersion profiles that are falling between 3 and 6 half mass radii, with a further 6 with larger errors that are flat or falling. Therefore, one result of our analysis is that 11/19 (58%) of the sample clusters show no evidence for significant local dark matter based on their falling velocity dispersion profiles at large radius. The 2 clusters with rising velocity dispersion profiles are consistent with local dark matter halos, but are not conclusive evidence. Those two, along with the 6 flat dispersion profile clusters are interesting targets for further observations and modeling. This research was supported by NSERC of Canada.