# Upper limits on the dark matter content in globular clusters

Javier Reynoso-Cordova  Marco Regis, and  Marco Taoso

###### Abstract

We present a systematic analysis on the possible presence of dark mass
components inside globular clusters (GCs). A spherical Jeans analysis is
applied to the stellar kinematics of 10 nearby GCs. On top of the mass
distribution provided by the luminous stellar component, we add either dark
matter (DM), described by an NFW or Burkert mass profile, or an intermediate
mass black-hole (IMBH), described by a point-like mass. Their existence would
have important implications in the context of indirect DM searches. After
profiling over the stellar parameters, we find no evidence neither for DM nor
for IMBH. Upper limits on the two components are reported.

##  1 Introduction

Globular clusters (GC) are old, very compact, gravitationally bound systems of
stars. They have a stellar mass of a few $\times\,10^{5}\,M_{\odot}$, similar
to the one of faint dwarf spheroidal (dSph) galaxies. On the other hand, GCs
have a much higher density, with the stellar mass confined within a region of
tidal radius of a few tens of parsec. DSph galaxies are instead significantly
more extended and they require the presence of a dominant dark matter (DM)
component to explain their dynamics. This is not the case for GCs, where there
is no compelling evidence, at the time of writing, that they host a dark
component different from dark remnants from stellar evolution. However, the
possible presence of a subdominant DM component is a very important question
to address. It is indeed crucial to understand GC formation, and, in
particular, whether GCs formed from a host halo which has been then
(partially) stripped away, or if GCs are simply big star clusters. Recent
simulations of the formation and evolution of GCs in DM halos can be found in
[1, 2, 3, 4, 5, 6] and references therein. The presence of a DM halo in GC
could have dramatic consequences for indirect searches of particle DM. Indeed,
several GCs in the Milky Way are very close to our location, and being very
compact, they could have a very high $J$-factor, which would make them a prime
target in the hunt for annihilating DM [7, 8, 9, 10, 11, 12, 13, 14].

In this work, we analyze the line-of-sight (LOS) velocity of stars belonging
to GCs, in order to infer the GC mass components. In the past years, several
works, involving different methods, have been devoted to the description of
the dynamics of GCs and their mass profiles, see e.g., the seminal work in
[15] and the discussion of the developments in [16, 17]. Recent measurements
are now offering the possibility to analyse several thousands of stars for
tens of GCs, greatly enhancing the statistical capability. Here we perform a
systematic study involving ten GCs and using most updated data, taken from
[16] and [18].

We investigate the possible presence of a DM mass component, with the spatial
density distributed following a NFW [19] or a Burkert [20] profile.
Additionally, we test the possibility that GCs host an intermediate mass black
hole (IMBH) at their center. This would constitute a dark “point-like” source
of mass, different from the extended distribution considered for DM. As for
the case of DM, the literature includes hints for a presence of IMBH in GCs
[21, 22, 23], as well as bounds [24, 25, 26, 27, 28, 29, 30], and this work
offers a systematic up-to-date address to the topic.

The paper is organized as follows. Section 2 describes the data on the LOS
velocity of stars employed in our work, and presents the Jeans analysis
adopted for the theoretical prediction. In Section 3, we describe the
statistical analysis and report the results in terms of bounds on the presence
of DM and IMBHs in GCs. We summarize consequences and conclusions in Section
4.

##  2 Data and methods

The analysis performed in this paper is based on two databases of LOS
velocities: the sample collected by [16] and recent measurements from the MUSE
collaboration [18]. The sample from [16] is a compilation of homogeneous
measures from spectra obtained at the Very Large Telescope and Keck telescope,
complemented by published measurements collected from the literature. The
reported velocities typically cover three decades in radius and have an
uncertainty around 0.5 km/s. Details on the sample selection and data analysis
can be found in [16] and [31], see also [32] for a recent update.

The dataset collected by the MUSE integral-field spectrograph [18] greatly
enlarges the number of stars acquired in each GC. Moreover, it provides LOS
velocities at radii typically closer to the GC centres, allowing a more
accurate derivation of the mass density in the innermost region.

GC | RA | DEC | D | $R_{\rm{half}}$ | Mass | # stars [B] | # stars [M]
---|---|---|---|---|---|---|---
47 Tuc | 6.02 | -72.08 | 4.52 | 4.10 | 8.95 | 3709 | 12666
NGC 1851 | 78.53 | -40.05 | 11.95 | 1.95 | 3.18 | 656 | 5559
NGC 2808 | 138.01 | -64.86 | 10.06 | 2.17 | 8.64 | 1523 | 3431
NGC 3201 | 154.4 | -46.41 | 4.74 | 3.33 | 1.60 | 2898 | 3283
$\omega$-cen | 201.70 | -47.48 | 5.43 | 7.93 | 36.40 | 2747 | 16282
M 80 | 244.26 | -22.97 | 10.34 | 2.01 | 3.38 | 434 | 1151
M 22 | 279.10 | -23.90 | 3.30 | 3.52 | 4.76 | 940 | 2465
NGC 6752 | 287.71 | -59.98 | 4.13 | 2.51 | 2.76 | 1578 | 4580
M 2 | 323.36 | -0.82 | 11.69 | 2.88 | 6.20 | 511 | 5238
M 30 | 325.09 | -23.18 | 8.46 | 2.9 | 1.40 | 800 | 4864
Table 1: GC parameters used in the analysis. I. Globular Cluster identifier.
II. Right ascension $[\deg]$. III. Declination $[\deg]$. IV. Distance [kpc].
V. Half-light radius as derived from the fit through pyGravSphere [pc]. VI.
Total Stellar Mass taken from [33] $[\times 10^{5}M_{\odot}]$. VII. Number of
analyzed stars from the dataset in [16]. VIII. Number of analyzed stars from
the MUSE dataset [18].

Among all the GCs surveyed in [16] and [18], we select the ones having more
than $10^{3}$ stars in both catalogues, in order to deal with large
statistics, and a reported mass-to-light ratio larger than 1.5, since we are
interested in testing the presence of a “dark” component on top of the
luminous one. This results into ten GCs. After this choice, we then select, in
each GC, only stars having a large membership probability (in details,
$P>0.01$ in [16] and $P>0.5$ in [18], following their definition of membership
probability). Binary stars can influence the velocity dispersion and therefore
bias the determination of the mass of the GC, see e.g. [34, 35]. In order to
mitigate this effect, we exclude binaries from the sample of stars for which
multi-epoch observations are available (we impose $P_{\rm SNGL}>0.05$ in [16]
and $P_{\rm variable}<0.8$ in [18], where $P_{\rm SNGL}$ and $P_{\rm
variable}$ are respectively the binary and single probabilities estimated in
those references). We find that this cut leads to a non-negligible reduction
of the velocity dispersions, and in general a better agreement between the two
samples we are considering. Finally, after applying these selection
requirements, for each GC we combine the two samples from [16] and [18] in an
single dataset, taking care of the possible double counting of stars included
in both surveys. The properties of all GCs composing our sample are reported
in Table 1.

We constrain the stellar surface density of the GCs using the results of [36],
where the radial number density profiles of a large sample of GCs have been
determined using data from Gaia DR2, from the Hubble Space Telescope, and
ground-based surface brightness measurements.

The data of the ten GCs of Table 1 are compared to a model built from a
spherical Jeans analysis. The assumptions we employ are that GCs can be
described through i) a collisionless Boltzmann equation, ii) in a steady-
state, iii) with a spherical phase-space distribution function, iv) without
rotation.
Clearly, none of these assumptions is exactly realized in nature. On the other
hand, there are indications that they can be good approximations for GCs. GCs
are very old stellar systems, with age much larger than relaxation time, i.e.,
with high number of gravitational encounters between stars. Deviations from
thermodynamic equilibrium in GCs are limited and, in this regime, the Jeans
equation provides a good approximate description [37], despite it is devised
for collisionless systems, whilst GCs are collisional. The high densities of
GC imply that they are more resistant to external tidal forces than extended
objects, like, e.g., dwarf spheroidal galaxies. However, there are indications
of tidal effects in the GC outskirt [38]. As described below, we will consider
relatively small radii where tidal forces are likely to be negligible.
Rotation has been measured in some GCs, e.g., $\omega$-Centauri [18, 39]. Its
contribution to the total kinetic energy of a GC is of the order of a few %
[39], and it is particularly negligible at the center, which is the most
relevant region to search for an IMBH or a highly concentrated DM profile. On
the other hand, since we are searching for subdominant components (i.e. DM and
IMBH), the impact on their final bounds might be sizable, and neglecting
rotation is a simplifying assumption of our method.

The above theoretical considerations are encoded in the spherical Jeans
equation [40, 41]:

$\frac{1}{\nu(r)}\frac{\partial}{\partial
r}\left(\nu(r)\sigma_{r}^{2}\right)+\frac{2\beta(r)\sigma_{r}^{2}}{r}=-\frac{GM(<r)}{r^{2}},\\\
$ (2.1)

where $\sigma_{r}$ is the radial velocity dispersion of stars, $M(<r)$ is the
cumulative GC mass as a function of radius $r$, $\nu$ is the stellar radial
density profile, and $\beta$ is the velocity anisotropy. We model the various
ingredients and solve Eq. (2.1) by using the non-parametric Jeans code
GravSphere [42] through the python implementation pyGravSphere [43]. For the
stellar tracer density, we sum three Plummer spheres [41]:

$\nu(r)=\sum_{j=1}^{3}\frac{3M_{j}}{4\pi
a_{j}^{3}}\left(1+\frac{r^{2}}{a_{j}^{2}}\right)^{-5/2}$ (2.2)

where $M_{j}$ and $a_{j}$ are free parameters, providing the mass and scale
length of each individual component. Let us mention that we do not identify
the Plummer spheres with specific stellar populations. Our approach is to
adopt a phenomenological effective function that provides a good description
of the measured total stellar density profile. In Appendix C we repeat our
analysis adopting a different stellar distribution, namely the King profile
[44].

The velocity anisotropy $\beta(r)$ describes the orbital structure of the
stellar system and it is modeled as:

$\beta(r)=\beta_{0}+\left(\beta_{\infty}-\beta_{0}\right)\frac{1}{1+\left(r_{a}/r\right)^{\eta}}$
(2.3)

where the free parameters are the “inner” anisotropy $\beta_{0}$, the “outer”
anisotropy $\beta_{\infty}$, the transition radius $r_{a}$, and the index
$\eta$ giving the sharpness of the transition.

The total mass is $M(r)=M_{*}(r)+M_{\rm{dark}}(r)$. We compute the stellar
mass $M_{*}$ assuming the same shape as the light profile given by Eq. (2.2)
but with a free parameter $M_{\star}$ fixing the overall normalization (which
corresponds to allowing a free mass-to-light ratio). In the case of DM, the
mass is computed by integrating the density distribution. We consider two
options: an NFW density profile [19]:

$\rho_{\rm{NFW}}(r)=\frac{\rho_{s}}{\left(\frac{r}{r_{s}}\right)\left(1+\frac{r}{r_{s}}\right)^{2}}\;,$
(2.4)

and a Burkert profile [20], describing a cored DM distribution:

$\rho_{\rm{Bur}}(r)=\frac{\rho_{s}}{\left(1+\frac{r}{r_{s}}\right)\left(1+\left(\frac{r}{r_{s}}\right)^{2}\right)}\;.$
(2.5)

In Eq. (2.4) and Eq. (2.5) the two free parameters $\rho_{s}$ and $r_{s}$ give
the normalization and the scale radius of the profile. In the case of IMBH, we
simply assume a point source of mass $M_{BH}$. From the radial velocity
dispersion derived by solving Eq. (2.1), we compute the LOS velocity
dispersion by means of

$\sigma_{\text{LOS}}^{2}(R)=\frac{2}{\Sigma_{*}(R)}\int_{R}^{\infty}\left(1-\beta\frac{R^{2}}{r^{2}}\right)\frac{\nu(r)\sigma_{r}^{2}(r)r}{\sqrt{r^{2}-R^{2}}}\mathrm{d}r,$
(2.6)

where $\Sigma_{*}(R)$ is the tracer surface mass density at the projected
radius $R$, derived from Eq. (2.2):

$\Sigma_{*}(R)=\sum_{j=1}^{3}\frac{M_{j}}{\pi
a_{j}^{2}}\left(1+\frac{R^{2}}{a_{j}^{2}}\right)^{-2}.$ (2.7)

In addition to $\sigma_{\text{LOS}}$, we also fit data concerning the surface
density profile of stars $\Sigma_{*}$, as well as the fourth-order moments of
the velocity distribution, called virial shape parameters (VSP) [42]

$\displaystyle v_{s1}$ $\displaystyle=$
$\displaystyle\frac{2}{5}\int_{0}^{\infty}G\,M\,\nu(r)\,[5-2\beta(r)]\,\sigma_{r}^{2}\,r\,dr=\int_{0}^{\infty}\Sigma_{*}(R)\,\langle
v_{\rm LOS}^{4}\rangle\,R\,dR$ (2.8) $\displaystyle v_{s2}$ $\displaystyle=$
$\displaystyle\frac{4}{35}\int_{0}^{\infty}G\,M\,\nu(r)\,[7-6\beta(r)]\,\sigma_{r}^{2}\,r^{3}\,dr=\int_{0}^{\infty}\Sigma_{*}(R)\,\langle
v_{\rm LOS}^{4}\rangle\,R^{3}\,dR\;,$ (2.9)

where the first term gives the model prediction, while the second one is the
estimator obtained from data.

The VSP can alleviate the degeneracy $\beta-\rho$ between the anisotropy of
the stellar component and the dark density, and they can be measured in this
context thanks to the large number of stars considered for each GC. One could
attempt to address the $\beta-\rho$ degeneracy also by considering tangential
velocities from proper motion data. However, the analysis of currently
available data can introduce more caveats than $\sigma_{\text{LOS}}$, in
particular a larger impact related to the description of GC rotation (that we
are not including in our model), see, e.g. [45]. For the sake of simplicity,
we restrict our analysis to LOS velocities.

##  3 Results

Figure 1:  Analysis of the $\omega$-Centauri observables in the case of a mass
model including an NFW DM profile. Left upper panel: line-of-sight velocity
dispersion. Blue points correspond to the kinematic data obtained from the
MUSE dataset [18] and from [16]. Right upper panel: the surface density
profile derived from [36] (orange points). In both panels the gray region
denotes the 95% interval of the posterior distribution. Lower panels:
measurements of the VSP parameters (red solid line) and the corresponding
2-$\sigma$ intervals (dashed red lines) compared to the posterior
distributions.

We compare the models with the observational data discussed in the previous
section considering the likelihood function $\mathcal{L}$ given by:

$\displaystyle-2\ln\mathcal{L}$ $\displaystyle=$ $\displaystyle\chi^{2}_{{\rm
LOS}}+\chi^{2}_{\Sigma_{*}}+\chi^{2}_{{\rm VSP1}}+\chi^{2}_{{\rm
VSP2}}=\chi^{2}.$ (3.1)

Our analysis includes measurements of the LOS velocity dispersion, the surface
density profile and virial shape parameters. The chi-squared $\chi^{2}$ for
the LOS velocity dispersion is obtained using data binned over the projected
radius $R$:

$\displaystyle\chi^{2}_{{\rm LOS}}=\sum_{j}\frac{\left(\sigma_{{\rm
LOS},j}-\bar{\sigma}_{{\rm LOS},j}\right)^{2}}{\sigma_{\sigma_{{\rm
LOS},j}}^{2}},$ (3.2)

where the quantity $\bar{\sigma}_{{\rm LOS},j}$ is the observed velocity
dispersion in the radial bin $j$. Following [42, 43], we use a number of bins
equal to $\sqrt{N_{*}},$ where $N_{*}$ is the number of stellar tracers. The
error in each bin, $\sigma_{\sigma_{{\rm LOS},j}}$ includes Poissonian and
sampling uncertainties, see [42, 43] for details. The dispersion velocity
predicted by the model, $\sigma_{{\rm LOS},j}$, is computed as explained in
the Section 2, see Eq. (2.6). Concerning the virial shape parameters, the chi-
squared $\chi^{2}_{{\rm VSP1}}$ and $\chi^{2}_{{\rm VSP2}}$, are obtained
analogously to Eq. (3.2), using the predicted values from the models computed
with Eq. (2.8) and Eq. (2.9), and the corresponding observational quantities
and estimated errors, obtained as detailed in [42, 43]:

$\displaystyle\chi^{2}_{{\rm
VSP1}}=\frac{\left(v_{s1}-\bar{v}_{s1}\right)^{2}}{\sigma_{v_{s1}}^{2}}\;\;\;,\;\;\;\chi^{2}_{{\rm
VSP2}}=\frac{\left(v_{s2}-\bar{v}_{s2}\right)^{2}}{\sigma_{v_{s2}}^{2}}.$
(3.3)

Finally, we fit the surface density profiles of [36]:

$\displaystyle\chi^{2}_{\Sigma_{*}}=\sum_{j}\frac{\left(\Sigma_{*,j}-\bar{\Sigma}_{*,j}\right)^{2}}{\sigma_{\Sigma_{*},j}^{2}},$
(3.4)

where $\bar{\Sigma}_{*,j}$ and $\sigma_{\Sigma_{*},j}$ are respectively the
surface number density in the radial bin $j$ and the corresponding error from
[36], while the prediction of the model $\Sigma_{*,j}$ is computed with Eq.
(2.7). In Eq. (3.4), for each GC, we include only measurements corresponding
to the radial distances probed by the LOS velocity dispersion in Eq. (3.2).
The reason is as follows. At the largest radii considered in [36], the
background contamination from field stars becomes dominant. This introduces a
systematic effect from a region which is not particularly relevant for our
purposes. On the other hand, in the innermost radii and for some GCs, there is
some level of mismatch between our model and the data. The same happens for
the models considered in [36]. This might be due to the incompleteness of the
sample of stars in these regions, or to the fact that the relative simple
models adopted are not tailored to fit simultaneously both the outskirts and
centres of the GCs. Again, to avoid the fit to be significantly affected by
data from outside the regions were we perform the dynamical analysis, we
exclude these internal regions from the surface density fit.

Figure 2: Same as Figure 1, but for the M2 GC.

The models outlined in Section 2 contain 11 parameters describing the stellar
system plus two (one) parameters for the DM (IMBH) component. We explore this
parameter space with a Markov Chain Monte Carlo (MCMC) analysis employing the
ensemble sampler Emcee [46]. We test different algorithms in order to ensure a
good coverage of the parameter space and to verify the robustness of the
results. For the plots presented here, we use 500 walkers and $10^{4}$ steps
with 5000 burn in.

To sample the parameters, we use flat priors within the range reported in
Table 2. We have 4 free parameters for the velocity anisotropy
($\beta_{\infty}$, $\beta_{0}$, $\eta$, $r_{a}$), a total of 7 parameters for
the stellar profile ($m_{1}$, $a_{1}$, $m_{2}$, $a_{2}$, $m_{3}$, $a_{3}$,
$M_{\star}$) and two extra parameters ($\rho_{s}$, $r_{s}$) in the DM case,
while one extra parameter ($\rho_{s}$) in the IMBH case. It is worth noticing
that, in the case of the stellar parameters, they are allowed to vary within a
range around a previously fitted value. Namely, prior the MCMC, pyGravSphere
fits the surface density through a minimization process determining the
associated best-fit values of the parameters of the stellar profile, and then,
during the MCMC, they are allowed to vary, around the surface density best-
fit, within the reported range. For the overall stellar mass $M_{\star}$ we
consider the value reported in Table 1 as the reference value, again allowing
some variation during the MCMC analysis. Also, since the anisotropy parameter
$\beta$ can take values from $-\infty$ to 1, GravSphere uses the symmetrized
anisotropy parameter [42] $\hat{\beta}=\beta/(2-\beta)$ which takes a value of
1 for a full radial anisotropy and $-1$ for tangential anisotropy. Note that,
for the sake of generality and to check the consistency between the DM and BH
findings, we allowed $r_{s}$ to take very small values, even though such
extremely concentrated DM halos are not expected. Clearly, we cannot
efficiently constrain profiles with $r_{s}$ much smaller than the radius of
the first data-point (or said in a different way, DM profiles with extremely
high concentration are allowed by data). The bounds on the DM mass are in any
case independent from choice of the priors, as we will discuss later on.

Parameters | min. value | max. value
---|---|---
$\log_{10}(\rho_{s}/(M_{\odot}/\rm{kpc}^{3}))$ | 0 | 30
$\log_{10}(r_{s}/\rm{kpc})$ | -8 | -1
$\beta_{\infty}$ | -1 | 1
$\beta_{0}$ | -1 | 1
$\eta$ | 1. | 10.
$r_{a}/\rm{kpc}$ | $R_{\rm{half}}/10$ | $10\cdot R_{\rm{half}}$
$M_{j}/M_{\odot}$ | $0.1(M_{j,\rm{bf}})$ | $1.9(M_{j,\rm{bf}})$
$a_{j}/\rm{kpc}$ | $0.1(a_{j,\rm{bf}})$ | $1.9(a_{j,\rm{bf}})$
$M_{\star}/M_{\odot}$ | $M_{\star}/2.$ | $2.5M_{\star}$
Table 2: Priors used in the MCMC. The parameters of the stellar profile can
vary around their best-fit values obtained by fitting the stellar surface
density profile only, before running the MCMC. The total stellar mass
$M_{\star}$ is allowed to vary around the value provided in Table 1.

We choose two illustrative examples to show the typical outcomes of the
analysis, focusing on the case with the DM component included in the fit. In
Fig. 1 we show the case of $\omega$-Centauri, which is among the brightest and
well studied GCs in the literature. We report the LOS velocity dispersion
profiles, the surface density profiles, and the estimates of the virial shape
parameters from the data, as well as the corresponding model predictions from
the fit. As an additional example we show the same plots but for the GC M2 in
Fig. 2. The LOS velocity dispersion profiles for all the others GCs are in the
Appendix A. In the Appendix B, we show triangle plots reporting the posterior
distributions of all the model parameters, including the stellar ones, for one
example ($\omega$-Centauri).

The main goal of this work is to derive upper limits on the presence of a DM
halo or an IMBH lying the center of GCs. In a Bayesian context, credible
intervals are obtained from the posterior distribution. However, given that we
find the parameters of the dark components, namely $\rho_{s}$ and
$M_{\rm{BH}}$, to be unconstrained from below, care must be taken to deal with
volume effects in the marginalization and to ensure that the results are
independent on the choice of the priors. For the sake of robustness, we decide
to adopt a frequentest approach, and derive our bounds performing a profile
likelihood analysis.

More specifically, we compute the likelihood ratio
$\lambda_{b}=\mathcal{L}(\vec{\theta}_{\rm b.f.},\vec{\Pi}_{\rm
b.f.})/\mathcal{L}(\vec{\theta},\vec{\Pi}^{\prime}_{\rm b.f.})$, where
$\vec{\theta}$ are the parameters of interest ($\rho_{s}$ and $r_{s}$ in the
DM case, and $M_{\rm{BH}}$ in the IMBH case), $\vec{\Pi}$ are the stellar
(nuisance) parameters, the suffix b.f. indicates the best-fit value, and the ′
symbol highlights that in the denominator we compute the best-fit values of
$\vec{\Pi}$ at any given value of $\vec{\theta}$, so they are in general
different from the global best-fit appearing in the numerator. It is easy to
see that
$-2\,\ln\lambda_{b}=\Delta\chi^{2}=\chi^{2}(\vec{\theta},\vec{\Pi}^{\prime}_{\rm
b.f.})-\chi^{2}(\vec{\theta}_{\rm b.f.},\vec{\Pi}_{\rm b.f.})$. The one-sided
95 % C.L. exclusion limits on the two parameters $\rho_{s}$ and $r_{s}$ (or on
the single parameter $M_{\rm{BH}}$) are obtained for $\Delta\chi^{2}>4.61$
($\Delta\chi^{2}>2.71).$

Results are shown in Figs. 3, 4, in Table 3 (DM), and in Table 4 (IMBH), and
discussed below. We note also that, for all the GCs analyzed, we do not find
any statistical preference for DM or IMBH, namely we find
$\mathcal{L}(\vec{\theta}_{\rm b.f.},\vec{\Pi}_{\rm
b.f.})\simeq\mathcal{L}(\vec{\theta}=0,\vec{\Pi}^{\prime}_{\rm b.f.})$. We
verified that the same conclusion can be drawn by looking at the marginalized
posteriors of $\rho_{s}$ and $M_{\rm{BH}}$, which are compatible with zero. In
Table 3, we display the number of data-points entering the statistical
analysis and report the best-fit values of the $\chi^{2}$ obtained from Eq.
3.1, for models including an NFW DM halo. Similar values are obtained for the
cases of a Burkert profile or an IMBH.

Figure 3: Upper limits at 95% C.L. on $\rho_{s}$ as a function of $r_{s}$ for
an NFW (top panel) and Burkert (bottom panel) DM profile. The bounds are
obtained following the profile likelihood procedure described in Section 3,
and the different panels show all the analyzed GCs. Figure 4: Upper limits at
95% C.L. on the DM mass enclosed within a radius $r$ for the different GCs
analyzed. Solid lines are the C.L. obtained using and NFW profile, dashed for
Burkert DM profile. Dots correspond to the minimum and maximum radii covered
by the measurements of the velocity dispersion.

###  3.1 Dark Matter

GC | $N_{\rm KIN}$ | $N_{\rm S.D.}$ | $\chi^{2}_{b.f.}$
---|---|---|---
47 Tuc | 127 | 115 | 262
NGC 1851 | 79 | 72 | 131
NGC 2808 | 70 | 60 | 141
$\omega$-cen | 138 | 33 | 215
M 80 | 40 | 134 | 120
M 22 | 56 | 59 | 85
NGC 6752 | 78 | 83 | 168
M 2 | 76 | 28 | 72
M 30 | 73 | 76 | 124
NGC 3201 | 71 | 51 | 126
Table 3: II. Number of kinematic data-points included in our statistical
analysis. III. Number of data points for the surface density. IV. $\chi^{2}$
corresponding to the best-fit for models including an NFW DM halo (similar
values are obtained for the Burkert profile and the IMBH case).

In the top panel of Fig. 3, we show the 95% C.L. upper limit on $\rho_{s}$ as
a function of $r_{s}$ obtained from the profile likelihood procedure described
above and assuming an NFW DM profile. Clearly, as the radius increases, the
bound on the normalization of the profile decreases, since it is the mass
distribution that is constrained by the kinematic analysis.

Fig. 4 reports the 95% C.L. upper bound on the DM mass enclosed within a given
GC radius $r$. It contains the same information as in Fig. 3, but shown in a
different way. We note that the increase in the mass bound is slower at low
radii and steeper at large radii. This can be understood as follows: for
radius $\lesssim$ first kinematic point (which is the dispersion velocity at
smallest observed radius), the profiles that can maximize the enclosed mass
are the ones with $r_{s}\lesssim$ first kinematic point. For an NFW profile,
the mass increases logarithmically for $r>r_{s}$ and so does the mass bound.
It is easy to see that for sufficiently small $r_{s}$ (again $\lesssim$ first
kinematic point) the bound is independent on $r_{s}$, so extending the prior
range to very low values has no impact on the constrained mass. When we
consider larger radii, DM profiles that can maximize the allowed mass have
larger $r_{s}$, and for $r_{s}\gtrsim r$, the mass of an NFW profile grows
more steeply than $\ln r$, up to $\propto r^{2}$ for $r_{s}\gg r$. Clearly,
the exact scaling also depends on the shape of the velocity dispersion profile
and on the stellar mass profile, that can be somewhat different in different
GCs. In Fig. 4, the minimum and maximum radius for the $\sigma_{\rm{LOS}}$
data points of each GC are reported.

Qualitatively, similar considerations hold for a Burkert DM profile. The upper
limits on $\rho_{s}$ as a function of $r_{s}$ are shown in the bottom panel of
Fig. 3. The bounds on the DM mass enclosed within a given GC radius, presented
with dashed lines in Fig. 4, are similar to the ones obtained for an NFW
distribution.

In Fig. 7 in the Appendix A, we show the posterior distributions of the Mass-
to-Light ratios derived from the analysis with an NFW profile and for the
different GCs. Very similar results are obtained considering a Burkert profile
or an IMBH, since these dark components provide only a small contribution to
the total mass. In the same figure, we also report the values of the Mass-to-
Light ratios presented in [47] (red lines). In general, our estimates are in
good agreement with those of [47], although a certain level of tension is
present in few cases. Similar mismatches (at the level of a few tens of %) are
found in the literature, see e.g. [48, 49, 50] and references therein, and are
likely due to systematic uncertainties related to the modeling and
measurements.

Recently, the possible presence of a DM halo in $\omega$-Centauri has been
investigated in [45, 13] on the basis of a Jeans analysis similar to the one
performed here. That work analyzed the LOS dispersion profiles from [16] as
well as measurements of the proper motion of stars. The upper bounds on the DM
content for an NFW distribution are comparable to the one in Fig. 4. However,
Ref. [45] finds significant evidence for the DM component. This different
conclusion is probably mainly arising from a different treatment of the
stellar component. In particular, our modeling include more free parameters to
describe the stellar profile and velocity anisotropy. Ref. [51] determined the
velocity dispersion profiles of several GCs using Gaia EDR3 data. For NGC 6752
(and NGC 6205) it is reported an hint of a rising velocity dispersion at radii
larger than $r\gtrsim 30-50$ pc, which correspond to distances beyond those
probed in our analysis. In principle, these results could be consistent with
the presence of a DM halo. However, further investigation is needed to
corroborate these findings, for instance to confirm the membership of the
stars to the cluster. Previous analyses in the literature have not find
evidence for DM in 47 Tucanae [52], in M22 and M30 [53], and in other GCs (not
analyzed in our work) as well [53, 54, 55, 56].

###  3.2 Intermediate Mass Black Holes

Table 4 reports the the 95% C.L. upper limit on $M_{\rm{BH}}$ obtained from
the profile likelihood analysis. The limits we find are of the order of (a few
times) $10^{4}\,M_{\odot}$.

Several methods are available to study IMBHs in GCs, namely to detect the
X-ray and radio emission from an IMBHs in the accretion phase, to measure the
acceleration of millisecond pulsars in the cluster, or to analyze the dynamics
of stars, as performed here. At present no solid evidence has been found from
these observations. A tentative detection of an IMBHs of $\sim 4\,\times
10^{4}\,M_{\odot}$ in $\omega$-Centauri has been suggested in [21, 22]. This
claim has been questioned by subsequent analyses [24, 25, 26], with ref. [24]
obtaining and upper limit on $M_{\rm{BH}}$ of $\sim 1.2\,\times
10^{4}\,M_{\odot.}$ Similarly, the presence for an IMBH of $\sim 2\times
10^{3}\,M_{\odot}$ in 47 Tucanae has been suggested in [23], while other works
concluded that current data do not favour such interpretation [27, 28, 29,
30]. An upper limits of $\sim 4\,\times 10^{3}\,M_{\odot}$ has been obtained
in [30]. A discussion on the searches for IMBHs in GCs can be found in the
recent review [57]. We notice that our bounds on $M_{\rm{BH}}$ in Table 3 are
weaker than the most stringent limits in the literature, which however are
derived employing different techniques and observations, more tailored to
probe the innermost region of GCs, rather than testing a full mass model as in
this work.

GC |  $M_{BH}^{\rm max}$ [$10^{4}\,M_{\odot}$]
---|---
47 Tuc | 2.9
NGC 1851 | 2.6
NGC 2808 | 4.0
$\omega$-cen | 6.0
M 80 | 3.3
M 22 | 2.2
NGC 6752 | 1.2
M 2 | 7.5
M 30 | 2.5
NGC 3201 | 1.4
Table 4: 95% C.L. upper limit on the mass of the IMBH.

##  4 Conclusions

We performed a spherical Jeans analysis of ten GCs of the Milky Way with the
goal of providing bounds on the mass components of these objects. The main
contribution to the mass budget is provided by the luminous stellar
distribution traced by the surface density profiles. Additional baryonic dark
components can consist of brown dwarfs, neutron stars and stellar-mass black
holes. We effectively include these contributions by making the simplifying
assumption that they follow the same spatial distribution of the luminous
stellar part, and by treating the total stellar mass as a free parameter,
i.e., allowing variations of the mass-to-light ratio. The main focus of the
work was in fact to test the possible presence of additional dark mass terms,
in the form of a dark halo, described by an NFW or Burkert mass profile, or an
IMBH, described by a point-like mass.

We described the mass distribution in each GC employing a model with 6
parameters pertinent to the stellar profile, 1 parameter providing the stellar
mass-to-light ratio, 4 parameters for the velocity anisotropy, and 2 (1)
parameters describing DM (IMBH). A different model for the stellar
distribution is adopted in Appendix C. We do not find any statistical
preference for DM or IMBH, and derived upper limits on their contribution to
the GC mass profile.

Results are shown in Figs. 3, 4, in Table 3 (DM), and in Table 4 (IMBH). The
DM mass is typically bounded to be $\lesssim 10^{5}\,M_{\odot}$ inside the GC
half-light radius. The upper limits on the IMBH mass are at the level of (a
few times) $10^{4}\,M_{\odot}$. The results presented in this work can be
employed to estimate the $J$\- and $D$-factors, i.e., to set the expected
emissions associated to indirect searches of DM.

## Acknowledgments

We warmly acknowledge Holger Baumgardt, Sebastian Kamann and Thomas de Boer
for providing the datasets used in this work and for useful suggestions. We
also thank Mariana Julio, Jarle Brinchmann, Mauro Valli, Eduardo Vitral and
Bas Zoutendijk for useful discussions. We acknowledge support by the PRIN
research grant “From Darklight to Dark Matter: understanding the galaxy/matter
connection to measure the Universe” No. 20179P3PKJ funded by MIUR, by the
research grant “The Dark Universe: A Synergic Multimessenger Approach No.
2017X7X85” funded by MIUR, by “Deciphering the high-energy sky via cross
correlation” and “Unveiling Dark Matter and missing baryons in the high-energy
sky” funded by the agreement ASI-INAF n. 2017-14-H.0, by the INFN grant
“LINDARK,” by the “Department of Excellence” grant awarded by the Italian
Ministry of Education, University and Research (MIUR), by the research grant
TAsP (Theoretical Astroparticle Physics) funded by Istituto Nazionale di
Fisica Nucleare (INFN).

##  Appendix A LOS velocity dispersion profiles

In this Appendix we present the LOS velocity dispersion profiles of all the
analyzed GCs, and the corresponding fits of the model including an NFW DM
profile (Fig. 5) or an IMBH (Fig. 6). Points are obtained from the datasets in
[18] and [16] as explained in Section 2. The gray shaded regions represent the
95%C.L. intervals obtained from the posterior distributions.

In Fig. 7 and for all the analyzed GCs, we show instead the posterior
distributions of the Mass-to-Light ratios derived from our analysis, compared
to the values reported in [47].

Figure 5: Line-of-sight velocity dispersion for all the GCs considered in this
work. Blue points correspond to the kinematic data obtained from the dataset
in [18] and [16] as explained in Section 2. The gray shaded regions show the
95% intervals of the posterior distributions for the analyses including an NFW
DM profile in the mass model.

Figure 6: Same as Fig. 5 but including an IMBH instead of an NFW DM profile in
the analysis. Figure 7: Posterior distributions of the Mass-to-Light ratios
[$M_{\odot}/L_{\odot}$] derived from our analysis compared to the values
reported in [47] (vertical solid red lines) with their associated 1-$\sigma$
uncertainties (vertical red dashed lines).

##  Appendix B Posterior distributions for one example: $\omega$-Centauri

The 1D and 2D posterior distributions of all the parameters considered in the
analyses are shown in Fig. 8 and Fig. 9, for the cases where the dark
component is in the form of an NFW DM profile and an IMBH, respectively. These
plots refer to one example, the analyses of the $\omega$-Centauri GC. Other
targets can be made available upon request to the authors.

Figure 8: Posterior distributions of the different parameters of the model. I.
$\log_{10}(\rho_{s}/[M_{\odot}\rm{kpc}^{-3}])$. II.
$\log_{10}(r_{s}/\rm{kpc}^{-3})$. III. $\beta_{0}$. IV. $\beta_{\infty}$. V.
$\log_{10}(r_{a}/\rm{kpc})$. VI. $\eta$. VII. $\log_{10}(m1/M_{\odot})$. VIII.
$a_{1}/\rm{kpc}$. IX. $\log_{10}(m2/M_{\odot})$. X. $a_{2}/\rm{kpc}$. XI.
$\log_{10}(m3/M_{\odot})$. XII. $a_{3}/\rm{kpc}$. XIII.
$\log_{10}(M_{\rm{star}}/M_{\odot})$ XIV. $-2\ln(\rm{prob})$ The analysis
refers to the $\omega$-Centauri GC and includes an NFW DM profile.  Figure 9:
Same as in Fig. 8 but including an IMBH instead of an NFW DM profile in the
analysis. Note that we report the value of
$\log_{10}(\rho_{s}/[M_{\odot}\rm{kpc}^{-3}])$ instead of the IMBH mass. In
the sampling we computed the IMBH mass using an NFW profile (Eq. 2.4) with
$r_{s}=10^{-9}$ kpc and truncated at $R_{\rm half}/100$ (with $R_{\rm half}$
taken from Table 1), i.e., $\rho_{{\rm NFW}}(r>R_{\rm half}/100)=0$. Then
$M_{BH}=4\pi\int_{0}^{R_{\rm half}/100}\,dr\,r^{2}\,\rho_{{\rm
NFW}}(r,r_{s}=10^{-9}\,{\rm kpc})$.

##  Appendix C Analysis with the King stellar profile

Figure 10: Upper limits at 95% C.L. on the DM mass enclosed within a radius
$r$ for the different GCs analyzed and using an NFW DM profile. Solid lines
show the bounds obtained using a three Plummer sphere density profile to model
the stellar distribution, while the dashed lines are the limits obtained by
using a King profile. Dots correspond to the minimum and maximum radii covered
by the measurements of the velocity dispersion.

In this section we repeat the analysis described in Section 2 adopting a
different model for the the stellar distribution. The goal is to study the
impact of this ingredient on our results. We substitute the three Plummer
spheres adopted in Section 2 with a King profile [44]. It constitutes a
simple, with only two parameters, and physically motivated distribution, often
used in the literature. The drawback of its simplicity is that it leads to
poor fits to the surface density profiles of the GCs that we have analyzed. In
fact, for all the targets, we find a best-fit $\chi^{2}_{\Sigma_{*}}$ in Eq.
(3.4) much larger than the one obtained employing three Plummer spheres, and
with a reduced $\chi^{2}_{\Sigma_{*}}$ value much larger than one. Similar
conclusions have been obtained in [36]. For this reason, namely, in order to
avoid a situation where with a global fit completely dominated by the surface
density part and with kinematic data having little impact, we slightly modify
the procedure discussed in Section 2. For each GC, the parameters of the King
model are fixed to the values which minimize $\chi^{2}_{\Sigma_{*}}$ and the
Mass-to-Light ratio is fixed to a physically motivated value (taken from the
best-fit of the analysis in the main text). Then, we exclude the surface
density data from the likelihood in Eq. (3.1), including only the measurements
of the LOS velocity dispersion and the virial shape parameters
($\chi^{2}=\chi^{2}_{\rm LOS}+\chi^{2}_{\rm VSP1}+\chi^{2}_{\rm VSP2})$. This
procedure is less general than the one of the main text, but offers an
alternative way to constrain the dark component using a physically motivated
model.

The results of the analysis including a NFW DM halo are presented in Fig. 10,
where we show the 95%C.L. upper limits on the DM mass enclosed within a given
GC radius $r$. We do not find statistical preference for a DM component in any
of the GCs analyzed. With respect to the bounds obtained in Section 3.1, by
modeling the stellar profile with three Plummer spheres, we see that, in the
case of the King profile, limits are typically similar or slightly stronger.
The latter can be ascribed to the fact that here we are adopting a more
constrained statistical procedure.

## References

  * [1] E. Ardi and H. Baumgardt.  Depletion of dark matter within globular clusters.  J. Phys. Conf. Ser., 1503(1):012023, 2020.
  * [2] Henriette Wirth, Kenji Bekki, and Kohei Hayashi.  Formation of massive globular clusters with dark matter and its implication on dark matter annihilation.  Mon. Not. Roy. Astron. Soc., 496(1):L70–L74, 2020, 2005.07128.
  * [3] Raymond G. Carlberg and Laura C. Keating.  Simulating Globular Clusters in Dark Matter Subhalos.  Astrophys. J., 924(2):77, 2022, 2105.13900.
  * [4] Pierre Boldrini, Roya Mohayaee, and Joseph Silk.  Embedding globular clusters in dark matter minihaloes solves the cusp–core and timing problems in the Fornax dwarf galaxy.  Mon. Not. Roy. Astron. Soc., 492(3):3169–3178, 2020, 1909.07404.
  * [5] Xiangcheng Ma, Michael Y. Grudić, Eliot Quataert, Philip F. Hopkins, Claude-André Faucher-Giguère, Michael Boylan-Kolchin, Andrew Wetzel, Ji-hoon Kim, Norman Murray, and Dušan Kereš.  Self-consistent proto-globular cluster formation in cosmological simulations of high-redshift galaxies.  Mon. Not. Roy. Astron. Soc., 493(3):4315–4332, 2020, 1906.11261.
  * [6] Piero Madau, Alessandro Lupi, Juerg Diemand, Andreas Burkert, and Douglas N. C. Lin.  Globular Cluster Formation from Colliding Substructure.  5 2019, 1905.08951.
  * [7] Gabrijela Zaharijas.  Implications of Intermediate Mass Black Hole in globular cluster G1 on Dark Matter detection.  Phys. Rev. D, 78:027301, 2008, 0712.4181.
  * [8] M. Wood, G. Blaylock, S. M. Bradbury, J. H. Buckley, K. L. Byrum, Y. C. K. Chow, W. Cui, I. de la Calle Perez, A. D. Falcone, S. J. Fegan, J. P. Finley, J. Grube, J. Hall, D. Hanna, J. Holder, D. Horan, T. B. Humensky, D. B. Kieda, J. Kildea, A. Konopelko, H. Krawczynski, F. Krennrich, M. J. Lang, S. LeBohec, T. Nagai, R. A. Ong, J. S. Perkins, M. Pohl, J. Quinn, H. J. Rose, G. H. Sembroski, V. V. Vassiliev, R. G. Wagner, S. P. Wakely, T. C. Weekes, and A. Weinstein.  A Search for Dark Matter Annihilation with the Whipple 10 m Telescope.  ApJ, 678(2):594–605, May 2008, 0801.1708.
  * [9] A. Abramowski, F. Acero, F. Aharonian, A. G. Akhperjanian, G. Anton, A. Balzer, A. Barnacka, U. Barres de Almeida, A. R. Bazer-Bachi, Y. Becherini, J. Becker, B. Behera, K. Bernlöhr, A. Bochow, C. Boisson, J. Bolmont, P. Bordas, V. Borrel, J. Brucker, F. Brun, P. Brun, T. Bulik, I. Büsching, S. Carrigan, S. Casanova, M. Cerruti, P. M. Chadwick, A. Charbonnier, R. C. G. Chaves, A. Cheesebrough, L. M. Chounet, A. C. Clapson, G. Coignet, J. Conrad, M. Dalton, M. K. Daniel, I. D. Davids, B. Degrange, C. Deil, H. J. Dickinson, A. Djannati-Ataï, W. Domainko, L. O’C. Drury, F. Dubois, G. Dubus, J. Dyks, M. Dyrda, K. Egberts, P. Eger, P. Espigat, L. Fallon, C. Farnier, S. Fegan, F. Feinstein, M. V. Fernandes, A. Fiasson, G. Fontaine, A. Förster, M. Füßling, Y. A. Gallant, H. Gast, L. Gérard, D. Gerbig, B. Giebels, J. F. Glicenstein, B. Glück, P. Goret, D. Göring, S. Häffner, J. D. Hague, D. Hampf, M. Hauser, S. Heinz, G. Heinzelmann, G. Henri, G. Hermann, J. A. Hinton, A. Hoffmann, W. Hofmann, P. Hofverberg, M. Holler, D. Horns, A. Jacholkowska, O. C. de Jager, C. Jahn, M. Jamrozy, I. Jung, M. A. Kastendieck, K. Katarzyński, U. Katz, S. Kaufmann, D. Keogh, D. Khangulyan, B. Khélifi, D. Klochkov, W. Kluźniak, T. Kneiske, Nu. Komin, K. Kosack, R. Kossakowski, H. Laffon, G. Lamanna, D. Lennarz, T. Lohse, A. Lopatin, C. C. Lu, V. Marandon, A. Marcowith, J. Masbou, D. Maurin, N. Maxted, T. J. L. McComb, M. C. Medina, J. Méhault, R. Moderski, E. Moulin, C. L. Naumann, M. Naumann-Godo, M. de Naurois, D. Nedbal, D. Nekrassov, N. Nguyen, B. Nicholas, J. Niemiec, S. J. Nolan, S. Ohm, J. P. Olive, E. de Oña Wilhelmi, B. Opitz, M. Ostrowski, M. Panter, M. Paz Arribas, G. Pedaletti, G. Pelletier, P. O. Petrucci, S. Pita, G. Pühlhofer, M. Punch, A. Quirrenbach, M. Raue, S. M. Rayner, A. Reimer, O. Reimer, M. Renaud, R. de los Reyes, F. Rieger, J. Ripken, L. Rob, S. Rosier-Lees, G. Rowell, B. Rudak, C. B. Rulten, J. Ruppel, F. Ryde, V. Sahakian, A. Santangelo, R. Schlickeiser, F. M. Schöck, A. Schulz, U. Schwanke, S. Schwarzburg, S. Schwemmer, M. Sikora, J. L. Skilton, H. Sol, G. Spengler, Ł. Stawarz, R. Steenkamp, C. Stegmann, F. Stinzing, K. Stycz, I. Sushch, A. Szostek, J. P. Tavernet, R. Terrier, O. Tibolla, M. Tluczykont, K. Valerius, C. van Eldik, G. Vasileiadis, C. Venter, J. P. Vialle, A. Viana, P. Vincent, M. Vivier, H. J. Völk, F. Volpe, S. Vorobiov, M. Vorster, S. J. Wagner, M. Ward, A. Wierzcholska, M. Zacharias, A. Zajczyk, A. A. Zdziarski, A. Zech, H. S. Zechlin, and H. E. S. S. Collaboration.  H.E.S.S. Observations of the Globular Clusters NGC 6388 and M15 and Search for a Dark Matter Signal.  ApJ, 735(1):12, July 2011, 1104.2548.
  * [10] Lei Feng, Qiang Yuan, Peng-Fei Yin, Xiao-Jun Bi, and Mingzhe Li.  Search for dark matter signals with Fermi-LAT observation of globular clusters NGC 6388 and M 15.  JCAP, 04:030, 2012, 1112.2438.
  * [11] Richard Bartels and Thomas Edwards.  Comment on “Understanding the $\gamma$-ray emission from the globular cluster 47 Tuc: Evidence for dark matter?”.  Phys. Rev. D, 100(6):068301, 2019, 1807.08800.
  * [12] Anthony M. Brown, Thomas Lacroix, Sheridan Lloyd, Céline Bœhm, and Paula Chadwick.  Understanding the $\gamma$-ray emission from the globular cluster 47 Tuc: evidence for dark matter?  Phys. Rev. D, 98(4):041301, 2018, 1806.01866.
  * [13] Anthony M. Brown, Richard Massey, Thomas Lacroix, Louis E. Strigari, Azadeh Fattahi, and Céline Bœhm.  The glow of annihilating dark matter in Omega Centauri.  arXiv e-prints, page arXiv:1907.08564, July 2019, 1907.08564.
  * [14] Javier Reynoso-Cordova, Oleg Burgueño, Alex Geringer-Sameth, Alma X. González-Morales, Stefano Profumo, and O. Valenzuela.  On the origin of the gamma-ray emission from Omega Centauri: millisecond pulsars and dark matter annihilation.  ”JCAP”, 2021(2):010, February 2021, 1907.06682.
  * [15] David Merritt, Georges Meylan, and Michel Mayor.  The stellar dynamics of omega centauri.  ”Astron. J.”, 114:1074–1086, September 1997, astro-ph/9612184.
  * [16] H. Baumgardt and M. Hilker.  A catalogue of masses, structural parameters, and velocity dispersion profiles of 112 Milky Way globular clusters.  MNRAS, 478(2):1520–1557, August 2018, 1804.08359.
  * [17] V Hé nault-Brunet, M Gieles, A Sollima, L L Watkins, A Zocchi, I Claydon, E Pancino, and H Baumgardt.  Mass modelling globular clusters in theGaiaera: a method comparison using mock data from anN-body simulation of m 4.  Monthly Notices of the Royal Astronomical Society, 483(1):1400–1425, nov 2018.
  * [18] S. Kamann, T. O. Husser, S. Dreizler, E. Emsellem, P. M. Weilbacher, S. Martens, R. Bacon, M. den Brok, B. Giesers, D. Krajnović, M. M. Roth, M. Wendt, and L. Wisotzki.  A stellar census in globular clusters with MUSE: The contribution of rotation to cluster dynamics studied with 200 000 stars.  MNRAS, 473(4):5591–5616, February 2018, 1710.07257.
  * [19] Julio F. Navarro, Carlos S. Frenk, and Simon D. M. White.  The Structure of Cold Dark Matter Halos.  ApJ, 462:563, May 1996, astro-ph/9508025.
  * [20] A. Burkert.  The Structure of dark matter halos in dwarf galaxies.  Astrophys. J. Lett., 447:L25, 1995, astro-ph/9504041.
  * [21] Eva Noyola, Karl Gebhardt, and Marcel Bergmann.  Gemini and Hubble Space Telescope Evidence for an Intermediate-Mass Black Hole in $\omega$ Centauri.  ApJ, 676(2):1008–1015, April 2008, 0801.2782.
  * [22] H. Baumgardt.  N -body modelling of globular clusters: masses, mass-to-light ratios and intermediate-mass black holes.  MNRAS, 464(2):2174–2202, January 2017, 1609.08794.
  * [23] Bülent Kızıltan, Holger Baumgardt, and Abraham Loeb.  An intermediate-mass black hole in the centre of the globular cluster 47 Tucanae.  Nature, 542(7640):203–205, February 2017, 1702.02149.
  * [24] Roeland P. van der Marel and Jay Anderson.  New Limits on an Intermediate Mass Black Hole in Omega Centauri: II. Dynamical Models.  Astrophys. J., 710:1063–1088, 2010, 0905.0638.
  * [25] Alice Zocchi, Mark Gieles, and Vincent Hénault-Brunet.  The effect of stellar-mass black holes on the central kinematics of $\omega$ Cen: a cautionary tale for IMBH interpretations.  MNRAS, 482(4):4713–4725, February 2019, 1806.02157.
  * [26] H. Baumgardt, C. He, S. M. Sweet, M. Drinkwater, A. Sollima, J. Hurley, C. Usher, S. Kamann, H. Dalgleish, S. Dreizler, and T. O. Husser.  No evidence for intermediate-mass black holes in the globular clusters $\omega$ Cen and NGC 6624.  MNRAS, 488(4):5340–5351, October 2019, 1907.10845.
  * [27] P. C. C. Freire et al.  Long-term observations of the pulsars in 47 Tucanae – II. Proper motions, accelerations and jerks.  Mon. Not. Roy. Astron. Soc., 471(1):857–876, 2017, 1706.04908.
  * [28] Christopher R. Mann, Harvey Richer, Jeremy Heyl, Jay Anderson, Jason Kalirai, Ilaria Caiazzo, Swantje Möhle, Alan Knee, and Holger Baumgardt.  A Multimass Velocity Dispersion Model of 47 Tucanae Indicates No Evidence for an Intermediate-mass Black Hole.  Astrophys. J., 875(1):1, 2019, 1807.03307.
  * [29] Federico Abbate, Mario Spera, and Monica Colpi.  Intermediate mass black holes in globular clusters: effects on jerks and jounces of millisecond pulsars.  Mon. Not. Roy. Astron. Soc., 487(1):769–781, 2019, 1905.05209.
  * [30] F. Abbate, A. Possenti, A. Ridolfi, P. C. C. Freire, F. Camilo, R. N. Manchester, and N. D’Amico.  Internal gas models and central black hole in 47 Tucanae using millisecond pulsars.  Monthly Notices of the Royal Astronomical Society, 481(1):627–638, November 2018, 1808.06621.
  * [31] H. Baumgardt, M. Hilker, A. Sollima, and A. Bellini.  Mean proper motions, space orbits, and velocity dispersion profiles of Galactic globular clusters derived from Gaia DR2 data.  MNRAS, 482(4):5138–5155, February 2019, 1811.01507.
  * [32] H. Baumgardt and E. Vasiliev.  Accurate distances to Galactic globular clusters through a combination of Gaia EDR3, HST, and literature data.  MNRAS, 505(4):5957–5977, August 2021, 2105.09526.
  * [33] Michael Hilker, Holger Baumgardt, Antonio Sollima, and Andrea Bellini.  Galactic Globular Clusters: A new catalog of masses, structural parameters, velocity dispersion profiles, proper motions and space orbits.  In Angela Bragaglia, Melvyn Davies, Alison Sills, and Enrico Vesperini, editors, Star Clusters: From the Milky Way to the Early Universe, volume 351, pages 451–454, January 2020, 1908.02778.
  * [34] Sara Rastello, Giovanni Carraro, and Roberto Capuzzo-Dolcetta.  Effect of Binarity in Star Cluster Dynamical Mass Determination.  ApJ, 896(2):152, June 2020, 2005.03037.
  * [35] Francisco I. Aros, Anna C. Sippel, Alessandra Mastrobuono-Battisti, Paolo Bianchini, Abbas Askar, and Glenn van de Ven.  Using binaries in globular clusters to catch sight of intermediate-mass black holes.  MNRAS, 508(3):4385–4398, December 2021, 2110.00590.
  * [36] T. J. L. de Boer, M. Gieles, E. Balbinot, V. Hénault-Brunet, A. Sollima, L. L. Watkins, and I. Claydon.  Globular cluster number density profiles using Gaia DR2.  MNRAS, 485(4):4906–4935, June 2019, 1901.08072.
  * [37] James Binney and Scott Tremaine.  Galactic Dynamics: Second Edition.  2008\.
  * [38] Andrés E. Piatti and Julio A. Carballo-Bello.  The tidal tails of Milky Way globular clusters.  A&A, 637:L2, May 2020, 2004.11747.
  * [39] A. Sollima, H. Baumgardt, and M. Hilker.  The eye of Gaia on globular clusters kinematics: internal rotation.  MNRAS, 485(1):1460–1476, May 2019, 1902.05895.
  * [40] J. Binney.  The radius-dependence of velocity dispersion in elliptical galaxies.  190:873–880, March 1980.
  * [41] James Binney and Gary A. Mamon.  M/L and velocity anisotropy from observations of spherical galaxies, or must M87 have a massive black hole?  Monthly Notices of the Royal Astronomical Society, 200(2):361–375, 09 1982, https://academic.oup.com/mnras/article-pdf/200/2/361/9403099/mnras200-0361.pdf.
  * [42] J. I. Read and P. Steger.  How to break the density-anisotropy degeneracy in spherical stellar systems.  MNRAS, 471(4):4541–4558, November 2017, 1701.04833.
  * [43] A. Genina, J. I. Read, C. S. Frenk, S. Cole, A. Benítez-Llambay, A. D. Ludlow, J. F. Navarro, K. A. Oman, and A. Robertson.  To $\beta$ or not to $\beta$: can higher order Jeans analysis break the mass-anisotropy degeneracy in simulated dwarfs?  MNRAS, 498(1):144–163, October 2020, 1911.09124.
  * [44] Ivan King.  The structure of star clusters. I. an empirical density law.  Astronomical Journal, 67:471, October 1962.
  * [45] Addy J. Evans, Louis E. Strigari, and Paul Zivick.  Stellar kinematics of the Omega Centauri globular cluster.  arXiv e-prints, page arXiv:2109.10998, September 2021, 2109.10998.
  * [46] Daniel Foreman-Mackey, David W. Hogg, Dustin Lang, and Jonathan Goodman.  emcee: The MCMC Hammer.  Publ. Astron. Soc. Pac., 125:306–312, 2013, 1202.3665.
  * [47] Michael Hilker, Holger Baumgardt, Antonio Sollima, and Andrea Bellini.  Galactic globular clusters: A new catalog of masses, structural parameters, velocity dispersion profiles, proper motions and space orbits.  Proceedings of the International Astronomical Union, 14(S351):451–454, 2019.
  * [48] Laura L. Watkins, Roeland P. van der Marel, Andrea Bellini, and Jay Anderson.  Hubble Space Telescope Proper Motion (HSTPROMO) Catalogs of Galactic Globular Clusters. III. Dynamical Distances and Mass-to-Light Ratios.  ApJ, 812(2):149, October 2015, 1509.00513.
  * [49] Dean E. McLaughlin and Roeland P. van der Marel.  Resolved Massive Star Clusters in the Milky Way and Its Satellites: Brightness Profiles and a Catalog of Fundamental Parameters.  ApJS, 161(2):304–360, December 2005, astro-ph/0605132.
  * [50] Richard R. Lane, László L. Kiss, Geraint F. Lewis, Rodrigo A. Ibata, Arnaud Siebert, Timothy R. Bedding, and Péter Székely.  Testing Newtonian gravity with AAOmega: mass-to-light profiles of four globular clusters.  MNRAS, 400(2):917–923, December 2009, 0908.0770.
  * [51] Raymond G. Carlberg and Carl J. Grillmair.  Testing for Dark Matter in the Outskirts of Globular Clusters.  Astrophys. J., 922(2):104, 2021, 2106.00751.
  * [52] Richard R. Lane, Laszlo L. Kiss, Geraint F. Lewis, Rodrigo A. Ibata, Arnaud Siebert, Timothy R. Bedding, and Peter Szekely.  Testing Newtonian Gravity with AAOmega: Mass-To-Light Profiles and Metallicity Calibrations From 47 Tuc and M55.  Mon. Not. Roy. Astron. Soc., 401:2521–2530, 2010, 0910.0882.
  * [53] Richard R. Lane, Laszlo L. Kiss, Geraint F. Lewis, Rodrigo A. Ibata, Arnaud Siebert, Timothy R. Bedding, and Peter Szekely.  Testing Newtonian Gravity with AAOmega: Mass-to-Light Profiles of Four Globular Clusters.  Mon. Not. Roy. Astron. Soc., 400:917–923, 2009, 0908.0770.
  * [54] Warren J. Hankey and Andrew A. Cole.  Radial Velocity and Metallicity of the Globular Cluster IC4499 Obtained with AAOmega.  Mon. Not. Roy. Astron. Soc., 411:1536, 2011, 1010.0469.
  * [55] H. Baumgardt, P. Cote, M. Hilker, M. Rejkuba, S. Mieske, S. G. Djorgovski, and P. Stetson.  The velocity dispersion and mass-to-light ratio of the remote halo globular cluster NGC 2419.  Mon. Not. Roy. Astron. Soc., 396:2051, 2009, 0904.3329.
  * [56] Rodrigo Ibata, Carlo Nipoti, Antonio Sollima, Michele Bellazzini, Scott Chapman, and Emanuele Dalessandro.  Do globular clusters possess Dark Matter halos? A case study in NGC 2419\.  Mon. Not. Roy. Astron. Soc., 428:3648, 2013, 1210.7787.
  * [57] Jenny E. Greene, Jay Strader, and Luis C. Ho.  Intermediate-Mass Black Holes.  Ann. Rev. Astron. Astrophys., 58:257–312, 2020, 1911.09678.