#  Multimass modelling of Milky Way globular clusters \- II. present-day black
hole populations

N. Dickson1, P.J. Smith1, V. Hénault-Brunet1, M. Gieles2,3, H. Baumgardt4
1Department of Astronomy and Physics, Saint Mary’s University, 923 Robie
Street, Halifax, NS B3H 3C3, Canada
2ICREA, Pg. Lluís Companys 23, E08010 Barcelona, Spain
3Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB),
Martí Franquès 1, E08028 Barcelona, Spain
4School of Mathematics and Physics, The University of Queensland, St Lucia,
QLD 4072, Australia E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Populations of stellar-mass black holes (BHs) in globular clusters (GCs)
influence their dynamical evolution and have important implications on one of
the main formation channels for gravitational wave sources. Inferring the size
of these populations remains difficult, however. In this work, multimass
models of 34 Milky Way GCs, first presented in Dickson et al., are used to
explore the present-day BH populations. Direct constraints on both the total
and visible mass components provided by several observables allow these models
to accurately determine the distribution of the dark mass (including BHs)
within clusters, as we demonstrate in a proof-of-concept fitting of the models
to mock observations extracted from Monte Carlo cluster models. New
constraints on the BH population retained to the present-day in each cluster
are inferred from our models. We find that BH mass fractions ranging from 0 to
1 per cent of the total mass are typically required to explain the
observations, except for $\omega$ Cen, for which we infer a mass fraction
above 5 per cent, in agreement with previous works. Relationships between the
dark remnant populations and other cluster parameters are examined,
demonstrating a clear anti-correlation between the amount of BHs and mass
segregation between visible stars, as well as a correlation between remnant
mass fractions and the dynamical age of clusters. Our inferred BH populations
are in good agreement overall with other recent studies using different
methodologies, but with notable discrepancies for individual clusters.

###### keywords:

galaxies: star clusters – globular clusters: general – stars: kinematics and
dynamics – stars: black holes

††pubyear: 2024††pagerange:  Multimass modelling of Milky Way globular
clusters \- II. present-day black hole populations –A

##  1 Introduction

It has been suggested that dynamical formation of close stellar-mass black
hole (BH) binaries in the dense cores of globular clusters (GCs) could be one
of the main formation channels for gravitational wave sources from BH-BH
mergers (e.g. Portegies Zwart & McMillan, 2000; Rodriguez et al., 2016; Abbott
et al., 2016; Antonini & Gieles, 2020a; Antonini et al., 2023). The amount of
dynamically formed BH-BH binaries in GCs and subsequent mergers, however,
depends on many uncertain physical ingredients, such as the initial mass and
number distribution of BHs (which itself is dependent on the stellar initial
mass function and stellar evolution; Spera et al., 2015, e.g.) and the
magnitude of natal kicks that BHs receive at the time of their formation in a
supernova explosion, which can eject them from the cluster (Chatterjee et al.,
2017), as well as the fraction of the initially retained BHs ejected due to
dynamical encounters throughout a cluster’s life, which depends on its initial
properties and dynamical evolution (e.g. Breen & Heggie, 2013a; Hurley et al.,
2016; Gieles et al., 2021).

The presence of massive BHs in GCs will have a large impact on their dynamical
evolution and present-day structure. As the most massive objects in the
system, any BHs which have not been ejected from the cluster by natal kicks
will rapidly segregate to the cluster centre, forming a concentrated
population of BHs in the core. While it was long argued that this dense
subsystem would dynamically decouple from the rest of the GC (e.g. Spitzer,
1969; Sigurdsson & Hernquist, 1993), and in short order eject all of the BHs
through strong dynamical interactions, recent work has shown that these BHs
can be retained on much longer timescales (e.g. Breen & Heggie, 2013a;
Morscher et al., 2013, 2015), and populations of BHs can be expected to
survive in most clusters to the present day, depending on the initial cluster
central densities. This theoretical work is complemented by recent detections
of BH candidates in binaries within GCs, based on radio/X-ray emission (from
accretion) or the radial velocity variations of a bright companion (Strader et
al., 2012; Miller-Jones et al., 2015; Giesers et al., 2018, 2019). Several
studies examining detailed evolutionary dynamical models (such as $N$-body or
Monte Carlo models) of clusters with an initial population of BHs have used
these expected impacts on observable quantities, such as the absence of mass
segregation among visible stars (e.g. Peuten et al., 2016; Alessandrini et
al., 2016; Weatherford et al., 2018, 2020), an elevated central mass-to-light
ratio (e.g. Baumgardt et al., 2019b) or a large effective radius of the
cluster (Torniamenti et al., 2023) and the presence of tidal tails (Gieles et
al., 2021), to argue for the presence of retained populations of BHs in
certain clusters at the present day.

Alternatively, equilibrium modelling approaches, such as Jeans modelling (e.g.
Kamann et al., 2014; Kamann et al., 2016; Abbate et al., 2019; Vitral & Mamon,
2021; Vitral et al., 2022) and multimass distribution-function (DF) based
models (e.g. Sollima et al., 2012, 2016; Zocchi et al., 2019; Hénault-Brunet
et al., 2019; Hénault-Brunet et al., 2020), have also been used to probe the
dark remnant content of GCs (including BHs) through fitting models to
observations of specific clusters while accounting for both visible and dark
mass components. Rapid and flexible equilibrium models, such as DF models,
allow for a more complete exploration of parameter space compared to more
computationally expensive evolutionary models, and, through the application of
statistical fitting techniques, the ability to very precisely reproduce the
kinematics and structure of real, observed GCs.

In Dickson et al. (2023, hereafter Paper I), we presented multimass DF models
fit to several observables (proper motions, line-of-sight velocities, stellar
densities and mass functions), for 34 Milky Way GCs. By inferring the global
stellar mass functions of these clusters and simultaneously constraining their
distributions of stellar remnants, the stellar initial mass function and its
possible dependence on metallicity was examined, in particular in the high-
mass regime ($\geq 1\,\mathrm{M}_{\odot}$) where stars in old GCs have evolved
into stellar remnants by the present day.

In this work, these models are used to now examine in detail the remnant
populations in this large sample of Milky Way GCs. In particular, the amount
and distribution of BHs in each cluster are inferred from our models. The
multimass limepy models (Gieles & Zocchi, 2015), mass function evolution
algorithm, observational datasets and model fitting procedures used are all
restated briefly, from Paper I, in Section 2. In Section 3, we provide a proof
of concept and show that our method is able to reliably recover the mass
fraction in BHs in GCs by fitting our models to mock data extracted from
snapshots of evolutionary models containing different amounts of BHs. The
overall distributions of BHs in our sample of GCs are presented in Section 4,
alongside an analysis of their relationships with other cluster parameters. In
Section 5 we discuss the implications of these results on the co-evolution of
GCs and their BHs, and provide comparisons between our results and the
inferred BH populations of other recent studies. Finally, we summarize our
results and conclude in Section 6.

##  2 Methods

In this section, the methodology used in Paper I to fit multimass dynamical
models to a number of observables is restated briefly, with emphasis on the
elements most crucial to inferring BH populations. For more detail on all
procedures summarized here, we refer the reader to Sections 2 through 4 of
Paper I.

###  2.1 Models

To model the mass distribution of the globular clusters we use the limepy
multimass distribution-function (DF) based models (Gieles & Zocchi, 2015). DF
based models are equilibrium models built around a distribution function which
describes the phase-space particle density of stars and satisfies the
collisionless Boltzmann equation. This DF is used to self-consistently solve
for the system’s potential ($\phi(r)$) using Poisson’s equation, and to derive
a variety of useful quantities for describing a GC, such as the projected
velocity dispersion, the projected surface density and the total mass.

Multimass DF models, defined by a sum of component DFs for individual mass
bins, allow for a more accurate description of real GCs, which are made up of
a spectrum of stellar and remnant masses, and are necessary in order to
account for the effects of mass segregation.

Our multimass models are defined by 10 free parameters dictating the mass
function and physical solution of the limepy DF. The overall structure of
these models is controlled by the (dimensionless) central potential parameter
$\hat{\phi}_{0}$, which defines how centrally concentrated the model is. To
mimic the effects of the host galaxy’s tides, the energy near the truncation
radius is reduced, lowering the escape velocity of stars and making it easier
to escape, with a sharpness of truncation defined by the parameter $g$ (lower
values resulting in a more abrupt truncation). The models can be expressed in
physical units by adopting relevant size and mass scales. In order to match
observations, we opt to scale the models using the parameters for total
cluster mass $M$ and 3D half-mass radius $r_{\mathrm{h}}$. The limepy models
allow for velocity anisotropy through an extra angular momentum term in the
DF, which produces an isotropic core followed by a degree of radial velocity
anisotropy at a distance from the centre defined by the anisotropy radius
parameter $r_{\mathrm{a}}$ and then returning to isotropy again near the
truncation radius. The multimass version of the limepy DF is defined in part
by the mass-dependent velocity scale $s_{j}$. This scaling captures the trend
towards kinetic energy equipartition among stars of different masses and
models the effects of mass segregation (Gieles & Zocchi, 2015; Peuten et al.,
2017; Hénault-Brunet et al., 2019), and is defined based on the parameter
$\delta$, such that $s_{j}\propto m_{j}^{-\delta}$.

The constituent discrete mass components which approximate the mass spectrum
of a GC are represented in the multimass models by the total ($M_{j}$) and
mean ($m_{j}$) masses of each mass bin. As DF-based models, such as limepy,
are equilibrium, instantaneous “snapshot” models, and do not directly simulate
any temporal astrophysical processes during their computation, we must instead
incorporate a separate prescription for stellar evolution from an initial mass
function, over the age of the cluster, to the present-day stellar and remnant
mass functions. In keeping with the formulation of canonical IMFs (e.g.
Kroupa, 2001), we use a 3-component broken power law, with power-law “slopes”
of each component given by the free parameters
$\alpha_{1},\,\alpha_{2},\,\alpha_{3}$ (with break masses at 0.5 and 1
$\mathrm{M}_{\odot}$ and bounded between 0.1 and 100 $\mathrm{M}_{\odot}$). To
evolve the population of stars to the present day we follow the algorithm
first described by Balbinot & Gieles (2018) and expanded upon in the ssptools
library111Available at https://github.com/SMU-clusters/ssptools and Paper I.
The amount of stars which evolve off the main sequence over the lifetime of
the clusters is dictated by a set of equations based on interpolated Dartmouth
Stellar Evolution Program models (Dotter et al., 2007; Dotter et al., 2008).
The types and masses of the stellar remnants formed by these evolved stars are
then determined based on their initial mass, metallicity and initial-final
mass relation (IFMR). The white dwarf (WD) IFMR, including the maximum initial
mass which will form a WD, is interpolated from the MIST 2018 isochrones
(Dotter, 2016; Choi et al., 2016). The BH IFMR, as well as the minimum initial
mass required to form a BH, is interpolated directly from a grid of stellar
evolution library (SSE) models (Banerjee et al., 2020), using the rapid
supernova scheme (Fryer et al., 2012). Stars with initial masses between the
WD and BH precursor masses are assumed to form neutron stars (NS) with a mass
of $1.4\,\mathrm{M}_{\odot}$.

The algorithm then must account for the loss of BHs through two different
channels. Firstly, the ejection of, primarily low-mass, BHs through natal
kicks is simulated. Beginning with the assumption that the kick velocity is
drawn from a Maxwellian distribution with a dispersion of
$265\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ (Hobbs et al., 2005) and
scaled down by $1-f_{\mathrm{fb}}$, where $f_{\mathrm{fb}}$ is the fallback
fraction, which we interpolated from the same grid of SSE models used for the
BH IFMR. The fraction of BHs retained, in each mass bin, is then found by
integrating the Maxwellian kick velocity distribution from 0 to the system
initial escape velocity, which we compute as
$v_{\mathrm{esc}}=2\sqrt{2\phi_{0}}$. BHs are also ejected over time from the
core of GCs due to dynamical interactions with one another (e.g. Breen &
Heggie, 2013a, b). This is simulated through the direct removal of BHs,
beginning with the heaviest mass bins (with larger gravitational interaction
cross-sections) through to the lightest (e.g. Morscher et al., 2015; Antonini
& Gieles, 2020b), until the combination of mass in BHs lost through both
ejection channels leads to a final retained mass in BHs equal to the
percentage of the initial mass in BHs for the given IMF specified by the BH
mass retention fraction parameter ($\mathrm{BH}_{\mathrm{ret}}$).

Finally, the heliocentric distance to the GCs $d$ is introduced as a free
parameter, to allow for the conversion between projected, linear model units
and the angular units of observations.

###  2.2 Fitting models to observations

In Paper I, best-fitting model parameters for 34 Milky Way GCs were determined
through the comparison of the phase-space distribution of stars in the models
to observations of GC structure and kinematics. These clusters were chosen
primarily to address the main hypothesis of Paper I (i.e. the possible
metallicity dependence of the IMF) and were selected based on the quantity and
quality of kinematic and mass function data available), however they also
allow us to conduct a census of BH populations in a significant sample of
Milky Way GCs (although our sample may be biased against low-mass, low-density
and metal-rich clusters, as discussed in section 5).

A variety of observational datasets were used to fit all chosen GCs, providing
direct constraints on the phase-space distribution of visible stars and the
overall mass of the cluster, and in turn providing indirect constraints on the
amount and distribution of dark mass (in both faint low-mass stars and dark
remnants). This is possible due to the fact that the distribution of the
different dark and visible mass components are arranged with limited
flexibility and linked, thanks to partial energy equipartition. Though BHs may
represent only up to a few per cent of the total mass of a GC, they can
actually dominate the mass density in the central regions of a cluster, with
significant dynamical effects, and thus can be probed. Details of the datasets
used for each cluster can be found in Appendix A of Paper I, but we summarize
the main sources below.

Radial profiles of proper motion (PM) and line-of-sight (LOS) velocity
dispersions of cluster stars are used to constrain the internal kinematics of
each cluster and thus its total (visible and dark) mass. PMs in both the
radial and tangential directions also provide constraints on the degree of
velocity anisotropy in the cluster, which is important given the degeneracy
between anisotropy and central dark mass (e.g. Zocchi et al., 2017). PM
dispersion (radial and tangential) profiles were computed in Paper I using
Gaia (DR3; Gaia Collaboration et al., 2022) proper motions for all clusters.
This data was supplemented with PM dispersion profiles in the cores of most
clusters from Hubble Space Telescope data (HST, Watkins et al., 2015;
Libralato et al., 2022). LOS velocity dispersion profiles are taken from
compilations of various ground-based (Baumgardt & Hilker, 2018; Kamann et al.,
2018; Dalgleish et al., 2020) and Gaia (Baumgardt et al., 2019a) datasets.

Radial profiles of the projected stellar number density in all GCs provide
vital constraints on the spatial distribution and concentration of the
clusters stars. The density profiles for all clusters are taken from de Boer
et al. (2019), consisting of combined Gaia star counts in the outskirts and
HST counts (Miocchi et al., 2013) or ground-based surface-brightness profiles
(SBPs, Trager et al., 1995) in the central regions.

Finally, constraints on the global present-day mass function of the clusters,
the degree of mass segregation and the total mass in visible stars are
provided by HST datasets, from which local stellar mass functions were
extracted. The mass function data for each cluster is taken from Baumgardt et
al. (2023), consisting of a stellar counts, based on a large amount of
archival HST data, in radial annuli and mass bins. The compilation of
photometry in each cluster is made up of several HST fields at varying
distances from the cluster centres, typically covering stars within a mass
range of $\sim$0.16\text{\,}\mathrm{-}$$0.8\text{\,}\mathrm{M}_{\odot}$$.

The models are constrained by these datasets in order to provide best-fitting
values and posterior probability distributions of the model parameters that
describe each cluster, determined through Bayesian parameter estimation
techniques222Latest model fits, cumulative mass profiles and black hole
distributions for all clusters are available online at
https://github.com/nmdickson/GCfit-results. The posterior probability
distributions of all parameters are sampled using dynamic nested sampling,
through the dynesty software package (Speagle, 2020). All fitting is carried
out using the software library and fitting pipeline GCfit333Available at
https://github.com/nmdickson/GCfit. For a discussion of the model fits and
parameter posterior distributions, see Paper I. It should be noted that the
latest version of GCfit, which includes various small improvements, such as to
the accuracy of the mass function likelihood, was used to refit all of the
models presented in this paper, however the conclusions of Paper I remain
unaffected.

The posterior probability distributions of the model-derived quantities used
throughout this work, such as the black hole mass fractions
($f_{\mathrm{BH}}=M_{\mathrm{BH}}/M_{\mathrm{cluster}}$), are constructed
based on the models representing the set of weighted posterior samples
retrieved from the nested sampler.

##  3 Validation of BH population inference

In order to test the reliability of our method in inferring BH populations, we
first apply it to simulated observations from Monte Carlo models with known BH
populations.

In order to explore a number of models with similar properties as real Milky
Way clusters, we select snapshots from the existing grid of Cluster Monte
Carlo (CMC; Rodriguez et al., 2022) models presented in Kremer et al.
(2020a)444Available at https://cmc.ciera.northwestern.edu/home. We select the
snapshots using the same methodology as Rui et al. (2021c) which we briefly
summarize here.

The selections are based on the SBPs of Trager et al. (1995) and the velocity
dispersion profiles (VDP) compiled by Baumgardt (2017), Baumgardt & Hilker
(2018) and Baumgardt et al. (2023)555Available at
https://people.smp.uq.edu.au/HolgerBaumgardt/globular . We search for
snapshots that are a good match to any clusters from the Harris (1996, 2010
edition) catalogue present in both the VDP and SBP compilations, leaving us
with about 100 clusters to match to snapshots. We first use the metallicities
from Harris (1996, 2010 edition) and the present-day galactocentric radii from
Baumgardt et al. (2019a) to select the subset of models which are closest to
the true values for each cluster.

From this subset, we then search every model for snapshots that match suitably
well to a given cluster’s observed SBP and VDP simultaneously 666See Rui et
al. (2021c) for details on how the model profiles are extracted from the CMC
snapshots for comparison with the observations.. In order to select a
snapshot, we adopt a threshold of
$s\equiv\max\left(\tilde{\chi}_{\mathrm{SBP}}^{2},\tilde{\chi}_{\mathrm{VDP}}^{2}\right)<10$
for the "fitting heuristic" $s$ of Rui et al. (2021c), which describes the
goodness-of-fit of a snapshot based on the $\tilde{\chi}^{2}$ statistic
between the observations and the interpolated model profiles. We have found
that a threshold of $s<10$ provides an acceptable balance between the number
and fit quality of the retained snapshots. While a number of snapshots passing
this filter have an apparently poor overall fit to one or both of the
observational profiles, we opt to still include these snapshots in the sample,
as our goal is not to select only snapshots which are perfect matches to
specific real clusters but instead to build a sample of snapshots that are
qualitatively similar to the Milky Way clusters examined in this work. For
each cluster covered by the observational datasets we select the single best-
fitting snapshot, where one exists.

###  3.1 Mock observations

The search described above results in a sample of 53 CMC model snapshots,
representative of Milky Way GCs. From these we next extract synthetic
observations designed to emulate the real observational data used to constrain
the models examined in this work.

We place each cluster at its respective heliocentric distance as reported by
Baumgardt & Vasiliev (2021) and then use the cmctoolkit library (Rui et al.,
2021c, a) to calculate projected positions and velocities as well as simulated
photometry for objects in each snapshot.

####  3.1.1 Number density profiles

We extract projected number density profiles from the models, designed to
emulate those of de Boer et al. (2019). We select all stars brighter than Gaia
$G=20$, sort them into 50 radial bins, with equal numbers of stars, and
calculate the number density in each radial bin. All densities are assigned a
Poisson counting error. The de Boer et al. (2019) profiles are combinations of
Gaia star counts in the outer regions and HST and archival SBPs in the inner
regions, where crowding becomes an issue, however we find that using the same
$G<20$ cut over the entire radial extent of the cluster results in well-
sampled profiles which cover a similar radial extent and have similar
uncertainties to the de Boer et al. (2019) profiles.

####  3.1.2 Proper motion dispersion profiles

We extract two sets of PM dispersion profiles for each snapshot, in order to
represent the performance of the two different sources of PM observations
used.

In the inner regions, we seek to emulate the performance of the HST based PM
dispersion profiles of Libralato et al. (2022). We select stars within the
central $100\text{\,}\mathrm{\SIUnitSymbolArcsecond}$ of the cluster to mimic
the footprint of an HST ACS footprint, and limit our selection to stars within
$15<V<18$. We split the stars into radial bins containing at least 120 stars
each, up to a maximum of five bins. This provides sufficient radial coverage
of the cluster while still allowing us to construct profiles for distant
clusters, where limited numbers of stars pass the magnitude cut. We assume a
typical uncertainty of
$0.1\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{/}\mathrm{yr}$ on all
stars. Within each bin we compute the mean velocity and velocity dispersion
along with their associated uncertainties, assuming the velocities are drawn
from a Gaussian error distribution, using MCMC. This is repeated for both the
radial and tangential components of PM. The median and $1\sigma$ values of the
dispersion in each bin are used going forward.

In the outer regions, we seek to emulate the Gaia DR3 based profiles of
Vasiliev & Baumgardt (2021). We base our magnitude cuts on their profiles,
selecting all stars in the $13<G<19$ range outside of the innermost
$100\text{\,}\mathrm{\SIUnitSymbolArcsecond}$, to avoid overlapping with the
HST profiles. We assign each star an uncertainty in proper motion based on its
$G$ band magnitude using the calibrations provided in Table 4 of Lindegren et
al. (2021), allowing us to replicate the performance of the Gaia DR3
catalogue. We again bin the stars using the same conditions as in the inner
HST profiles, and calculate the velocity dispersion in each bin using the same
method, again for both radial and tangential components.

####  3.1.3 Line-of-sight velocity dispersion profiles

In addition to the PM dispersion profiles, we also extract LOS velocity
dispersion profiles, designed to emulate those presented by Baumgardt (2017),
Baumgardt & Hilker (2018) and Baumgardt et al.
(2023)${}^{\ref{foot:baum_cat}}$. As the compilation of velocity dispersions
that make up these profiles consist of several different inhomogeneous
datasets, with varying precisions, we adopt the simplifying assumption of a
typical uncertainty of $1\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ on
all observed stars. We limit this dataset to only giants brighter than $V=17$,
which is typical of the datasets used in the observed compilations. We again
sort the stars into several radial bins, requiring at least 70 stars per bin,
up to a maximum of 10 bins, and compute the velocity dispersion for each in
the same way as for the PM profiles.

####  3.1.4 Stellar mass functions

We extract stellar mass function data for each of our snapshots, designed to
emulate the datasets presented in Baumgardt et al. (2023). These datasets
consists of star counts, binned by stellar mass, extracted from archival HST
observations. While the real HST fields are distributed somewhat randomly
around each GC according to the various goals of each proposal for which they
were originally observed, in general there is typically at least one exposure
centred on the cluster centre and a number of fields placed outside of the
central region. For simplicity, we opt to place one field over the centre of
the cluster, covering a range of $0\text{\,}\mathrm{\SIUnitSymbolArcminute}$
\- $1.6\text{\,}\mathrm{\SIUnitSymbolArcminute}$ in projected angular
separation from the centre, and two outer annuli at radial distances of
$2.5\text{\,}\mathrm{\SIUnitSymbolArcminute}$ and
$5\text{\,}\mathrm{\SIUnitSymbolArcminute}$, each sized such that they cover
the same area as the central field. The central field is split into 4 annuli,
giving us a total of 6 annular fields covering the central, intermediate and
outer regions of the cluster. In many of the more well-studied clusters in our
sample, such as NGC 104 and $\omega$ Cen, the large number of observed HST
fields actually provides much better coverage of the clusters than our fields
simulated here, and thus the results of this section may actually be
conservative. In each of these fields we extract stellar counts separately, in
bins of stellar mass with widths of .

In the real datasets, the faintest stars (lowest stellar masses) for which
stellar counts can be extracted with reasonable completeness ($>90$ per cent)
is a function of crowding. To replicate this effect in our synthetic mass
functions we construct an empirical relation between the surface number
density and the lowest observable mass within a field. We use NGC 104 as the
basis for this calibration because it covers a wide range of number densities
from its core to the outskirts and has a large number of HST fields for which
the mass function was extracted. We extract star counts in each field, down to
a lower mass limit calculated by the above relation and up to the main
sequence turn-off, replicating the performance of the observed HST stellar
counts.

We assign Poisson counting errors to our stellar counts, though to reflect the
scatter we typically see in the real data we also inflate these errors by a
factor of $F=3$ and re-sample each point within the errors, resulting in mass
functions that are very similar to those of Baumgardt et al. (2023) (see
Section 4.1.1.3 of Paper I for a description of this $F$ parameter and its
motivation).

###  3.2 Validation results

After extracting the synthetic datasets, we then directly apply our model
fitting method (as described in Section 2 and Paper I) and compare the
resulting inferred BH mass fractions to the known BH population in the CMC
models from which the mock data was extracted.

As in Paper I, we discard any obviously poor fits. We also discard any
snapshots which have much smaller datasets (mostly consisting of snapshots
matched to very distant clusters). This leaves us with 44 final snapshots with
datasets of similar quality to the Milky Way clusters we study in this work,
and with model fits which satisfyingly reproduce the mock observations and
recover the various cluster parameters well, such as the total mass, which we
generally recover within $\sim 10$ per cent and the half-mass radius, which we
generally recover within uncertainties. Our inferred values of
$f_{\mathrm{BH}}$, compared to the true values for our collection of
snapshots, are shown in Figure 1 (right panel), Figure 2, and Table 2.

Figure 1: Left panel: Validation snapshots plotted in the $\delta-
N_{\mathrm{relax}}$ plane. A selection of dynamically young clusters with high
inferred values of $\delta$ is shown by the grey shaded region
($N_{\mathrm{relax}}<3$, $\delta>0.4$). The points are coloured based on the
number of $\sigma_{f_{\mathrm{BH}}}$ the median inferred $f_{\mathrm{BH}}$ is
away from the true value (where $\sigma_{f_{\mathrm{BH}}}$ is the width of the
posterior for $f_{\mathrm{BH}}$). Under this colour scheme, snapshots for
which we underpredict (overpredict) $f_{\mathrm{BH}}$ are increasingly orange
(blue). Right panel: The $f_{\mathrm{BH}}$ values inferred based on the
multimass model fits to mock observations extracted from CMC models, against
the true values in those validation snapshots ($f_{\mathrm{BH,true}}$).
Snapshots which fall within the shaded region in the left panel are plotted
with unfilled markers in the right panel, demonstrating that they also largely
represent the snapshots for which we significantly underpredict
$f_{\mathrm{BH}}$.

While in general our fits satisfyingly recover the mass fraction in BHs, we
find a number of snapshots for which we underpredict the black hole mass
fraction. A common feature of these problematic snapshots is their dynamical
age. All of the mock clusters for which we significantly underpredict
$f_{\mathrm{BH}}$ have an age that is less than 3 times their present-day
half-mass relaxation time. In other words, they are dynamically very young.
There are also, however, many dynamically young clusters for which correctly
recover $f_{\mathrm{BH}}$.

To understand this behaviour, we plot in the left panel of Figure 1 all of the
mock clusters used in the validation of our models on the $\delta-
N_{\mathrm{relax}}$ plane (where $N_{\mathrm{relax}}$ is the ratio of cluster
age to present-day half-mass relaxation time), coloured by the distance of
their median inferred $f_{\mathrm{BH}}$ value from the true value. As $\delta$
is a measure of the degree of energy equipartition in our multimass models
(with typical values of $\delta\approx 0.5$ for mass-segregated, dynamically
evolved clusters), we should normally expect dynamically young clusters
($\lesssim 3\,N_{\mathrm{relax}}$) to have lower values of $\delta$. The mock
clusters for which we significantly underpredict $f_{\mathrm{BH}}$ are clearly
concentrated in the upper-left corner of the $\delta-N_{\mathrm{relax}}$
plane, highlighted with the grey-shaded rectangle in the left panel of Figure
1. Their value of $\delta$ is inconsistent with our expectations for these
systems.

Previous comparisons of the multimass limepy models used here with $N$-body
models (Peuten et al., 2017) have confirmed this intuition that dynamically
young clusters should have lower values of $\delta$ (with $\delta\sim 0.35$)
and have also shown that models with any significant black hole population
should have similarly low values of $\delta$ until they eject most of their
black holes (after several relaxation times). Therefore we should not expect
to find any clusters in the upper-left or lower-right sections of the $\delta-
N_{\mathrm{relax}}$ plane and any model fits which fall into these regions are
worthy of some suspicion, even though it is not immediately clear why the data
for some of these dynamically young clusters would prefer models with values
of $\delta$ close to 0.5. Simple test fits with $\delta$ forced to a lower
value ($\delta\leq 0.35$) have shown that a lower value of $\delta$ can result
in an increase in the inferred $f_{\mathrm{BH}}$ for these snapshots, as
generally expected.

Within our validation sample there are a handful of problematic clusters in
the upper-left, suspect region of the $\delta-N_{\mathrm{relax}}$ plane. In
the right panel of Figure 1, we show our inferred values of $f_{\mathrm{BH}}$
compared to the true values for our collection of snapshots, with the
snapshots that fall into this problematic region plotted with unfilled
markers. It can be clearly seen that, if these clusters are ignored, the
correlation between the inferred and true black hole mass fraction is even
stronger. In Figure 2, we plot the inferred values of $f_{\mathrm{BH}}$
compared to the true values but now ignoring these problematic snapshots and
showing the $1\sigma$ statistical error bars (in blue) for the retained mock
clusters.

With these caveats in mind we can turn to the sample of real clusters we
investigate in this work to look for clusters that fall into this suspect
region. In the real sample only three clusters out of the 34 fall into this
region of relatively high $\delta$ and low $N_{\mathrm{relax}}$ (compared to 9
of the 44 validation snapshots). The relative lack of problematic fits in the
real sample compared to the mock sample highlights the fact that our mock
sample is not a perfect imitation of the population of MW clusters we
investigate in this work and is limited by the constraints of the CMC grid.
One of these problematic clusters is NGC 3201, which is discussed further
below (Section 5.3.4). The other two are NGC 104 (47 Tuc) and NGC 5139
($\omega$ Cen), two of the most massive clusters in our sample, which still
possess notable BH populations despite their high $\delta$. Our results for
these two well-studied clusters concurs well with other literature values (see
the discussion of each in Section 5.3). One real cluster in our sample falls
into the lower-right portion of $\delta-N_{\mathrm{relax}}$. This is NGC 6624,
a core-collapsed cluster which is discussed further in Section 4.1, and for
which we argue that we incorrectly infer a non-zero number of black holes.

As may be expected, the statistical uncertainties derived solely from the
fitting procedure seem to underestimate the real uncertainties slightly. Our
fitting procedure operates under the assumption that our models are a perfect
representation of the data, and as such, in reality, may underestimate the
true errors. It has been shown that multimass DF models, such as those used
here, may underestimate the uncertainties when compared to more flexible
models, such as Jeans models (Hénault-Brunet et al., 2019), which could be
indicative of systematic errors not captured in the statistical uncertainties
and limitations in the ability of these models to perfectly reproduce the
data.

In an attempt to quantify this underestimation, we search for the factor by
which the statistical uncertainties on our inferred values of
$f_{\mathrm{BH}}$ need to be inflated to make them fully consistent with a
one-to-one relation with the true model values. We define a nuisance parameter
$\epsilon$ by which we multiply the statistical errors
$\sigma_{f_{\mathrm{BH}}}$ to determine the total, inflated error $\Delta$ on
each inferred $f_{\mathrm{BH}}$. We then fit a fixed, one-to-one (i.e. slope
of 1, intercept of 0) line through the points in Figure 2, assuming a Gaussian
likelihood and allowing $\epsilon$ to vary freely and inflate the
uncertainties on the inferred $f_{\mathrm{BH}}$. This fit results in a value
of $\epsilon=2.5\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$. We also show
the inflated $1\sigma$ error bars on the inferred $f_{\mathrm{BH}}$ in Figure
2 (orange), demonstrating the additional uncertainty needed to make our
inferred values consistent with the true values and illustrating the typical
uncertainties expected when applying our method to real data.

Figure 2: The $f_{\mathrm{BH}}$ values inferred based on the mock observations
extracted from CMC models, against the true values in those models
($f_{\mathrm{BH,true}}$). The one-to-one line is shown in grey, representing
perfect agreement. The median and $1\sigma$ values, based solely on the
statistical uncertainties from the fit, are shown in blue. The inflated
errors, based on the nuisance parameter $\epsilon$ ($\Delta=\epsilon\
\sigma_{f_{\mathrm{BH}}}$), are shown in orange.

Overall, this comparison with mock observations extracted from dynamical
simulations lends confidence in the ability of our methods to correctly
recover the mass fraction in BHs in GCs, with the important note that the
statistical uncertainties on these inferred values may underestimate the true
uncertainties by up to a factor of $\sim 2.5$, and with a reminder that
dynamically young clusters with inferred values of $\delta\sim 0.5$ should be
treated with caution.

##  4 Black Hole Populations

We will now explore the populations of black holes (and other dark remnants)
inferred from our best-fitting models. We will examine the distribution of the
total mass, mass fraction and amount of black holes found from these fits, and
explore possible correlations present between the remnant populations and
other cluster properties. We will also discuss in more detail the case of
‘core-collapsed’ clusters and how we model them.

Cluster | $f_{\mathrm{BH}}\ \left[\%\right]$ | ${M}_{\mathrm{BH}}\ \left[M_{\odot}\right]$ | ${N}_{\mathrm{BH}}$ | $\bar{m}_{\mathrm{BH}}\ \left[\mathrm{M}_{\odot}\right]$ | $f_{\mathrm{remn}}\ \left[\%\right]$
---|---|---|---|---|---
NGC104 | $0.046\begin{subarray}{c}+0.017\\\ -0.009\end{subarray}$ | $420\begin{subarray}{c}+150\\\ -80\end{subarray}$ | $60\begin{subarray}{c}+20\\\ -10\end{subarray}$ | $6.7$ | $46.9\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$
NGC288 | $0.09\begin{subarray}{c}+0.09\\\ -0.06\end{subarray}$ | $80\begin{subarray}{c}+80\\\ -50\end{subarray}$ | $6\begin{subarray}{c}+5\\\ -4\end{subarray}$ | $14.1$ | $56\begin{subarray}{c}+2\\\ -2\end{subarray}$
NGC362∗ | $0.070\begin{subarray}{c}+0.010\\\ -0.017\end{subarray}$ | $190\begin{subarray}{c}+30\\\ -50\end{subarray}$ | $14\begin{subarray}{c}+2\\\ -3\end{subarray}$ | $14.0$ | $45.1\begin{subarray}{c}+0.6\\\ -0.6\end{subarray}$
NGC1261 | $0.14\begin{subarray}{c}+0.07\\\ -0.07\end{subarray}$ | $300\begin{subarray}{c}+100\\\ -100\end{subarray}$ | $19\begin{subarray}{c}+10\\\ -9\end{subarray}$ | $14.0$ | $39\begin{subarray}{c}+2\\\ -2\end{subarray}$
NGC1851 | $0.060\begin{subarray}{c}+0.011\\\ -0.006\end{subarray}$ | $200\begin{subarray}{c}+40\\\ -20\end{subarray}$ | $15\begin{subarray}{c}+3\\\ -2\end{subarray}$ | $13.2$ | $41.0\begin{subarray}{c}+0.5\\\ -0.6\end{subarray}$
NGC2808 | $0.08\begin{subarray}{c}+0.01\\\ -0.02\end{subarray}$ | $700\begin{subarray}{c}+100\\\ -200\end{subarray}$ | $47\begin{subarray}{c}+10\\\ -10\end{subarray}$ | $15.7$ | $36.2\begin{subarray}{c}+0.5\\\ -0.6\end{subarray}$
NGC3201 | $0.008\begin{subarray}{c}+0.026\\\ -0.008\end{subarray}$ | $10\begin{subarray}{c}+50\\\ -10\end{subarray}$ | $0.9\begin{subarray}{c}+3.2\\\ -0.9\end{subarray}$ | $14.3$ | $46\begin{subarray}{c}+1\\\ -1\end{subarray}$
NGC5024 | $0.11\begin{subarray}{c}+0.15\\\ -0.08\end{subarray}$ | $600\begin{subarray}{c}+800\\\ -400\end{subarray}$ | $40\begin{subarray}{c}+60\\\ -30\end{subarray}$ | $14.1$ | $43\begin{subarray}{c}+1\\\ -2\end{subarray}$
NGC5139 | $5.9\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $182000\begin{subarray}{c}+5000\\\ -6000\end{subarray}$ | $12400\begin{subarray}{c}+700\\\ -300\end{subarray}$ | $14.8$ | $50.3\begin{subarray}{c}+0.9\\\ -0.3\end{subarray}$
NGC5272 | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $1500\begin{subarray}{c}+500\\\ -600\end{subarray}$ | $110\begin{subarray}{c}+40\\\ -40\end{subarray}$ | $13.4$ | $46.9\begin{subarray}{c}+0.8\\\ -0.7\end{subarray}$
NGC5904 | $0.05\begin{subarray}{c}+0.06\\\ -0.04\end{subarray}$ | $200\begin{subarray}{c}+200\\\ -200\end{subarray}$ | $20\begin{subarray}{c}+20\\\ -10\end{subarray}$ | $10.9$ | $55.7\begin{subarray}{c}+1.0\\\ -0.9\end{subarray}$
NGC5986 | $0.03\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$ | $80\begin{subarray}{c}+100\\\ -80\end{subarray}$ | $6\begin{subarray}{c}+7\\\ -6\end{subarray}$ | $13.8$ | $53\begin{subarray}{c}+1\\\ -1\end{subarray}$
NGC6093 | $0.62\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$ | $1800\begin{subarray}{c}+100\\\ -100\end{subarray}$ | $129\begin{subarray}{c}+9\\\ -9\end{subarray}$ | $14.1$ | $53.3\begin{subarray}{c}+0.8\\\ -0.8\end{subarray}$
NGC6121 | $0.09\begin{subarray}{c}+0.05\\\ -0.03\end{subarray}$ | $80\begin{subarray}{c}+50\\\ -30\end{subarray}$ | $6\begin{subarray}{c}+3\\\ -2\end{subarray}$ | $13.4$ | $64.6\begin{subarray}{c}+0.9\\\ -0.8\end{subarray}$
NGC6171 | $0.07\begin{subarray}{c}+0.04\\\ -0.04\end{subarray}$ | $40\begin{subarray}{c}+30\\\ -30\end{subarray}$ | $3\begin{subarray}{c}+2\\\ -2\end{subarray}$ | $13.9$ | $66.6\begin{subarray}{c}+1.0\\\ -0.9\end{subarray}$
NGC6205 | $0.4\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $1800\begin{subarray}{c}+800\\\ -700\end{subarray}$ | $130\begin{subarray}{c}+50\\\ -50\end{subarray}$ | $13.9$ | $55\begin{subarray}{c}+1\\\ -1\end{subarray}$
NGC6218 | $0.14\begin{subarray}{c}+0.01\\\ -0.02\end{subarray}$ | $140\begin{subarray}{c}+10\\\ -20\end{subarray}$ | $10.0\begin{subarray}{c}+0.8\\\ -1.4\end{subarray}$ | $14.1$ | $61.4\begin{subarray}{c}+0.7\\\ -0.7\end{subarray}$
NGC6254 | $0.09\begin{subarray}{c}+0.08\\\ -0.07\end{subarray}$ | $200\begin{subarray}{c}+200\\\ -100\end{subarray}$ | $14\begin{subarray}{c}+11\\\ -9\end{subarray}$ | $13.2$ | $57.0\begin{subarray}{c}+0.8\\\ -0.9\end{subarray}$
NGC6266∗ | $0.14\begin{subarray}{c}+0.01\\\ -0.01\end{subarray}$ | $1060\begin{subarray}{c}+80\\\ -70\end{subarray}$ | $130\begin{subarray}{c}+10\\\ -10\end{subarray}$ | $8.2$ | $52\begin{subarray}{c}+2\\\ -2\end{subarray}$
NGC6341 | $0.2\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $600\begin{subarray}{c}+400\\\ -300\end{subarray}$ | $50\begin{subarray}{c}+30\\\ -20\end{subarray}$ | $12.4$ | $52\begin{subarray}{c}+1\\\ -1\end{subarray}$
NGC6352 | $0.062\begin{subarray}{c}+0.008\\\ -0.010\end{subarray}$ | $57\begin{subarray}{c}+6\\\ -8\end{subarray}$ | $5.2\begin{subarray}{c}+0.6\\\ -0.7\end{subarray}$ | $10.9$ | $59\begin{subarray}{c}+2\\\ -2\end{subarray}$
NGC6362 | $0.03\begin{subarray}{c}+0.04\\\ -0.03\end{subarray}$ | $30\begin{subarray}{c}+50\\\ -30\end{subarray}$ | $2\begin{subarray}{c}+3\\\ -2\end{subarray}$ | $12.4$ | $65.1\begin{subarray}{c}+1.0\\\ -0.9\end{subarray}$
NGC6366 | $0.07\begin{subarray}{c}+0.08\\\ -0.05\end{subarray}$ | $20\begin{subarray}{c}+30\\\ -20\end{subarray}$ | $2\begin{subarray}{c}+2\\\ -1\end{subarray}$ | $13.0$ | $55\begin{subarray}{c}+2\\\ -2\end{subarray}$
NGC6397∗ | $0.019\begin{subarray}{c}+0.001\\\ -0.001\end{subarray}$ | $21\begin{subarray}{c}+2\\\ -2\end{subarray}$ | $3.1\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ | $6.7$ | $57.2\begin{subarray}{c}+0.8\\\ -0.8\end{subarray}$
NGC6541∗ | $0.054\begin{subarray}{c}+0.008\\\ -0.011\end{subarray}$ | $130\begin{subarray}{c}+20\\\ -30\end{subarray}$ | $12\begin{subarray}{c}+2\\\ -2\end{subarray}$ | $10.1$ | $53.3\begin{subarray}{c}+0.8\\\ -0.8\end{subarray}$
NGC6624∗ | $0.48\begin{subarray}{c}+0.02\\\ -0.03\end{subarray}$ | $490\begin{subarray}{c}+20\\\ -20\end{subarray}$ | $46\begin{subarray}{c}+2\\\ -2\end{subarray}$ | $10.8$ | $70.6\begin{subarray}{c}+0.7\\\ -0.7\end{subarray}$
NGC6681∗ | $0.056\begin{subarray}{c}+0.006\\\ -0.007\end{subarray}$ | $55\begin{subarray}{c}+6\\\ -6\end{subarray}$ | $6.5\begin{subarray}{c}+0.7\\\ -0.7\end{subarray}$ | $8.5$ | $72.8\begin{subarray}{c}+0.5\\\ -0.6\end{subarray}$
NGC6723 | $0.05\begin{subarray}{c}+0.07\\\ -0.04\end{subarray}$ | $90\begin{subarray}{c}+120\\\ -70\end{subarray}$ | $8\begin{subarray}{c}+9\\\ -6\end{subarray}$ | $12.2$ | $66\begin{subarray}{c}+1\\\ -1\end{subarray}$
NGC6752∗ | $0.45\begin{subarray}{c}+0.01\\\ -0.02\end{subarray}$ | $920\begin{subarray}{c}+30\\\ -40\end{subarray}$ | $83\begin{subarray}{c}+3\\\ -5\end{subarray}$ | $11.1$ | $67.5\begin{subarray}{c}+0.4\\\ -0.4\end{subarray}$
NGC6779 | $0.48\begin{subarray}{c}+0.06\\\ -0.04\end{subarray}$ | $770\begin{subarray}{c}+90\\\ -70\end{subarray}$ | $54\begin{subarray}{c}+6\\\ -5\end{subarray}$ | $14.3$ | $53.9\begin{subarray}{c}+0.9\\\ -1.0\end{subarray}$
NGC6809 | $0.05\begin{subarray}{c}+0.02\\\ -0.03\end{subarray}$ | $90\begin{subarray}{c}+40\\\ -50\end{subarray}$ | $6\begin{subarray}{c}+3\\\ -4\end{subarray}$ | $14.3$ | $57\begin{subarray}{c}+1\\\ -1\end{subarray}$
NGC7078∗ | $0.36\begin{subarray}{c}+0.05\\\ -0.04\end{subarray}$ | $2200\begin{subarray}{c}+400\\\ -300\end{subarray}$ | $290\begin{subarray}{c}+40\\\ -30\end{subarray}$ | $7.9$ | $50.4\begin{subarray}{c}+0.5\\\ -0.5\end{subarray}$
NGC7089 | $0.4\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$ | $2300\begin{subarray}{c}+700\\\ -700\end{subarray}$ | $130\begin{subarray}{c}+40\\\ -30\end{subarray}$ | $17.7$ | $40\begin{subarray}{c}+2\\\ -2\end{subarray}$
NGC7099∗ | $0.042\begin{subarray}{c}+0.003\\\ -0.002\end{subarray}$ | $56\begin{subarray}{c}+5\\\ -3\end{subarray}$ | $6.9\begin{subarray}{c}+0.6\\\ -0.4\end{subarray}$ | $8.2$ | $53.2\begin{subarray}{c}+0.7\\\ -0.6\end{subarray}$
Table 1: Median and $1\sigma$ credibility intervals of the BH mass fraction
($f_{\mathrm{BH}}$), total mass in BHs (${M}_{\mathrm{BH}}$), total number of
BHs (${N}_{\mathrm{BH}}$) and remnant mass fraction ($f_{\mathrm{remn}}$) in
all clusters in our sample, as well as the median mean BH mass
($\bar{m}_{\mathrm{BH}}$). All clusters classified as core-collapsed in Trager
et al. (1995) are denoted by an asterisk. Note that all uncertainties
presented here represent only the statistical uncertainties on the fits, and
likely underestimate the true uncertainties (see Section 3).  Figure 3: Violin
plots (in blue) of the posterior probability distribution of the mass fraction
in BHs (upper panel), the total mass in BHs (middle panel) and total number of
BHs (lower panel) in all clusters in our sample, except for NGC 5139, which
has a median total mass in BHs of approximately ($f_{\mathrm{BH}}\sim 5.9$ per
cent), and is excluded in order to better highlight the distributions of the
other clusters. The median and $1\sigma$ values are denoted by the horizontal
blue ticks within each distribution. These errors include only the statistical
uncertainties on our fits, and thus are likely underestimated (see Section 3).
Clusters are sorted based on the mass fraction in black holes. All clusters
classified as core-collapsed in Trager et al. (1995) are denoted by an
asterisk. The median and $1\sigma$ results are also shown for the
corresponding quantities in Weatherford et al. (2020) (red), Askar et al.
(2018) (orange) and Rui et al. (2021b) (purple), for all clusters in common
with our sample. Values from Weatherford et al. (2020) are computed using the
median clustercentric mass segregation parameter $\Delta_{r50}$ (Table 1 of
Weatherford et al. 2020), and any necessary conversions between total mass and
mass fraction are computed using our total cluster mass estimates.

Figure 3 shows the posterior probability distributions of the mass fraction in
BHs ($f_{\mathrm{BH}}$), the total mass in BHs ($M_{\mathrm{BH}}$), and total
number of BHs ($N_{\mathrm{BH}}$) inferred from our best-fitting models of
most clusters in our sample. The median and $1\sigma$ credibility intervals of
the distributions of these quantities, as well as the mean individual BH mass
($\bar{m}_{\mathrm{BH}}$) and the mass fraction in all dark remnants
($f_{\mathrm{remn}}$), are also presented in Table 1. NGC 5139 ($\omega$ Cen)
is not included in Figure 3, due to the very large inferred amount of black
holes ($5.9\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ per cent), but it is
discussed in more detail in Section 5.3.

A large number of the clusters are consistent, within $2\sigma$, with
harbouring little to no black holes, while the remainder possess, on average,
at most a few thousand $\mathrm{M}_{\odot}$ of stellar-mass black holes, with
constituent individual BH masses ranging between $\sim 5$ to $15\
\mathrm{M}_{\odot}$. Other than $\omega$ Cen, it is clear that none of our
models favour a very large population of black holes, with all clusters having
a mass fraction in black holes between 0 and 1 per cent at the present day.

It must be noted again that, as examined in Section 3, the errors on all of
these quantities represent only the statistical uncertainties, and in reality
the uncertainties on a quantity like $f_{\mathrm{BH}}$ may be underestimated
by a factor of $\sim 2.5$.

###  4.1 Core-collapsed clusters

After all BHs are dynamically ejected, the visible core ‘collapses’ (Breen &
Heggie, 2013a). Globular clusters having undergone core collapse are typically
defined based on the shape of their central density profiles, with core-
collapsed clusters showing a power-law density profile increasing all the way
to their centres, while non core-collapsed clusters possess larger, isothermal
cores with a flat central density profile (e.g. Djorgovski & King, 1986;
Trager et al., 1995). Core-collapsed clusters are expected to contain very
few, if any, black holes at the present day (Breen & Heggie, 2013a, b; Kremer
et al., 2020b). In GCs with a population of stellar-mass BHs, core collapse
occurs within the BH subsystem but, due to the efficient heat transfer from
BHs to stars, the visible core will actually remain large (relative to
$r_{\mathrm{h}}$). The presence of BHs in a cluster may thus play a large role
in explaining the observed population of core-collapsed Milky Way GCs, which,
given the ages of most clusters, is smaller than would be expected when
considering only stellar binaries as the sole mechanism delaying core
collapse. It is not until almost all BHs (and the last BH binary; Hurley,
2007) are ejected that a cluster core will collapse and exhibit the defining
power-law central density profile (Chatterjee et al., 2013; Kremer et al.,
2020b). Almost all GCs have likely reached a state of balanced evolution
(Hénon, 1961; Gieles et al., 2011) but, due to BHs, only a minority of GCs
($\sim 20\%$) appear to be post-collapse.

The nine clusters in our sample defined as core-collapsed in Trager et al.
(1995) (NGC 362, NGC 6266, NGC 6397, NGC 6541, NGC 6624, NGC 6681, NGC 6752,
NGC 7078, NGC 7099) are denoted in Figures 3 and 1 by an asterisk. Our best-
fitting models of some of these clusters clearly favour a non-negligible
amount of mass in black holes, however, as discussed below, we suspect the
results are unreliable for these few clusters.

Part of this discrepancy between the theoretical expectations and the inferred
BH populations for some core-collapsed clusters may arise simply due to the
limitations of the limepy models themselves. limepy models, by definition,
possess an isothermal core, characterized by a flat inner density profile,
which is incompatible with the central cusp of core-collapsed clusters (see
also Section 3.1.4 of Gieles & Zocchi 2015). As such, our models may struggle
to accurately capture the inner density profiles of these clusters right up to
the centre. Indeed this divergence can be seen in the profiles of the core-
collapsed clusters with for which we infer substantial BH populations from our
best-fitting models, which tend to underestimate the amount of stars within a
very small distance from the centre (typically $\sim$0.1\text{\,}\mathrm{pc}$$
in these clusters). While, in most core-collapsed clusters, the models may be
able to provide a satisfactory fit to the available density data simply by
having a sufficiently small core (below the radial reach of the data), in some
clusters the shape of other datasets (such as the mass functions) may require
a larger core, and cause the fitting procedure to sacrifice the quality of the
fit to the central density profiles.

In order to investigate these systems further, models were also fit to these
clusters, in the same fashion as before, but with the amount of retained black
holes at the present day now fixed to 0 (by fixing the
$\mathrm{BH}_{\mathrm{ret}}$ parameter to 0 per cent). As might be expected,
the most immediately noticeable change to the model fits is in the number
density profiles. Shown in Figure 4 are examples of the changes between the
sets of models for NGC 6624 and NGC 7078, the two core collapsed clusters in
our sample with the highest inferred mass fraction in BHs
($f_{\mathrm{BH}}\sim 0.5$ per cent) and the only two where significant
differences can be seen between the models with and without BHs. In the other
core-collapsed clusters, which favour smaller or negligible BH populations,
the fits do not change noticeably. Likely, as suggested above, the cores of
these models are small enough that even though the models are not truly core
collapsed, they are still able to represent the observed density profiles. In
the case of NGC 6624 and NGC 7078, Figure 4 clearly shows that models without
BHs are able to also reproduce the central density profiles, and would provide
a much better fit in that regime, however, the models containing BHs have a
higher overall likelihood, due to the fits to other datasets, which cannot be
reproduced as well by the models without BHs. Overall, it is clear that
caution must be applied when attempting to fit core-collapsed clusters using
limepy models.

Figure 4: The number density profiles of the best-fitting models of NGC 6624
and NGC 7078, with and without allowing for a population of BHs. The number
density data used to constrain the models is shown by the orange circles.
Background levels have been subtracted from all data. Inset frames show a
zoomed view on the model profiles near the cluster cores, to showcase the
differences between the sets of models.

The original models (with BHs) are used throughout this paper and in all
discussion of BHs, however the results for these core-collapsed clusters
should be regarded with great caution, especially when a large population of
BHs is inferred. All core-collapsed clusters are noted in the figures and
tables by an asterisk or a red outline. As the larger inferred BH populations
of some of these clusters are most likely artificial, the models introduced
here, fixed to 0 BHs, were used in Paper I to avoid any impacts on the
inferred mass function slopes (though these were found to be minimal).

###  4.2 Relationships between BH population and other parameters

We next examine how the population of black holes and other remnants in our
cluster models correlate with various related parameters.

Figure 5: Relations between the high-mass initial mass function exponent
($\alpha_{3}$), the black hole retention fraction parameter
($\mathrm{BH}_{\mathrm{ret}}$) and the mass fraction in black holes
($f_{\mathrm{BH}}$) for all clusters, except for NGC 5139, which has a high
$\mathrm{BH}_{\mathrm{ret}}$ of $\sim 35\begin{subarray}{c}+9\\\
-3\end{subarray}$ per cent and a more typical $\alpha_{3}$ value of $\sim
2.28\begin{subarray}{c}+0.08\\\ -0.02\end{subarray}$. All core-collapsed
clusters, whose inferred black hole populations may not be physical, are
highlighted by a red outline.

Figure 5 shows the relationships between the high-mass initial mass function
exponent $\alpha_{3}$ ($\geq$1\text{\,}\mathrm{M}_{\odot}$$) from Paper I, the
black hole retention fraction $\mathrm{BH}_{\mathrm{ret}}$ and the BH mass
fraction $f_{\mathrm{BH}}$. This serves to demonstrate the role of the
$\mathrm{BH}_{\mathrm{ret}}$ parameter, which is not directly proportional to
the number of BHs. At high values of $\alpha_{3}$ (i.e. steeper slopes), only
a small number of black holes can be formed initially from the IMF, and a
higher retention fraction is required to maintain any amount of black hole
mass at the present day. The top panel of Figure 5 shows that no clear
correlation is present between $\alpha_{3}$ and the mass fraction in BHs at
the present day. The visible pattern instead relates to the relationship with
$\mathrm{BH}_{\mathrm{ret}}$; to end up with similar mass fractions in BHs,
clusters with higher $\alpha_{3}$ values produce fewer BHs initially while
clusters with low $\alpha_{3}$ values produce more BHs initially but retain
few. This degeneracy can be clearly seen in the rightmost panel. This of
course does not imply a causal relationship between BH mass fraction and
either $\alpha_{3}$ or $\mathrm{BH}_{\mathrm{ret}}$, but merely helps to
explain the distribution of $\mathrm{BH}_{\mathrm{ret}}$. For a given cluster
there is generally not a direct degeneracy between $\alpha_{3}$ and
$\mathrm{BH}_{\mathrm{ret}}$ apparent in the posterior probability
distributions of the parameters, as the shape of the IMF and retention factor
also strongly impact the amount of other massive stars and remnants formed,
which may be strongly constrained by the observations.

Figure 6: Relation between the mass-dependant velocity scale ($\delta$) and
the mass fraction in black holes for all clusters, except NGC 5139. All core-
collapsed clusters, whose inferred black hole populations may not be physical,
are highlighted by a red outline.

We do find an interesting relationship between the mass fraction of BHs and
the parameter $\delta$, which sets the mass-dependence of the velocity scale
and acts as a proxy of mass segregation, as shown in Figure 6. Clusters with
little to no mass in BHs tend to converge near values of $\delta\sim
0.4\rm{-}0.5$, which is typical of evolved and mass-segregated clusters,
whereas the clusters with more substantial populations of black holes
congregate closer to the lower bound of 0.3. This is in agreement with the
study of Peuten et al. (2017), who find, by comparing limepy models against
$N$-body models with and without black holes, that the majority of mass-
segregated clusters should converge to a value of $\sim 0.5$, but also show
that, in models with a significant population of black holes, the degree of
mass segregation as traced by the parameter $\delta$ may be suppressed.

##  5 Discussion

###  5.1 The evolution of clusters and their BH populations

The co-evolution of a cluster and its BH population depends directly on the BH
mass fraction $f_{\mathrm{BH}}$ and the initial density of the cluster. Very
early on in the lifetime of a cluster, $f_{\mathrm{BH}}$ will rapidly
increase, as stellar mass is lost due to stellar evolution and BHs are formed.
This initial population of BHs, and $f_{\mathrm{BH}}$, depends primarily on
the cluster metallicity, the shape of the IMF and the velocity distribution of
natal kicks imparted on BHs when they form. The BHs initially retained (after
kicks) will rapidly segregate to the cluster core. During the proceeding
expansionary phase, BHs will be ejected from the core due to dynamical BH-BH
interactions, while the lowest mass stars, on the cluster periphery, will be
the most affected by tidal losses, leading to a decrease in $f_{\mathrm{BH}}$
(Breen & Heggie, 2013b; Gieles et al., 2021). The amount that
$f_{\mathrm{BH}}$ decreases during this stage is dependent on the initial
cluster density, relative to the tidal density. Higher density clusters will
eject a higher proportion of the initially formed BHs (Banerjee & Kroupa,
2011).

Once the cluster becomes tidally limited, the remaining $f_{\mathrm{BH}}$
determines the subsequent evolution of the cluster mass. Theory (Breen &
Heggie, 2013b) and $N$-body modelling (Gieles & Gnedin, 2023) suggests that
there exists a critical BH mass fraction $f_{\mathrm{crit}}$, at which the
dynamical losses of BHs (relative to the total cluster mass loss) will exactly
equal $f_{\mathrm{BH}}$, causing $f_{\mathrm{BH}}$ to remain constant until
the dissolution of the cluster. Given an $f_{\mathrm{BH}}$ greater(less) than
this critical fraction, tidal stellar mass losses will be greater(less) than
the dynamical ejections of BHs from the core. Clusters above the critical
fraction will progress towards 100 per cent BHs, while clusters below will
eject most of the BHs and begin to evolve similarly to a 0 BH system. This
critical fraction was determined to be $\sim$ 10 per cent (Breen & Heggie,
2013b) for idealized single-mass clusters filling their tidal radius, however
recent $N$-body modelling placed it closer to 2.5 per cent for clusters with a
full mass spectrum (Gieles et al., 2021).

In our sample, while there is a range of values between individual clusters,
we find that nearly all clusters have a $f_{\mathrm{BH}}$ below 1 per cent
(except for $\omega$ Cen). Gieles & Gnedin (2023) modelled the Milky Way GC
population using a parametrization for the mass-loss rate from GCs with BHs,
assuming that the initial GC half-mass density is 30 times higher than Roche
filling. Under this assumption, they predicted that the majority of GCs should
now have $f_{\mathrm{BH}}\sim 2$ per cent. Our empirical results are lower
than this, which could indicate that their initial densities are not high
enough. However, the authors do note that the constant density (relative to
the tidal density) they have used for all GCs is only meant to approximate the
average density of GCs. In reality, a distribution of ratios of initial
densities over the tidal densities must exist, which would lead to both a
population of GCs that can eject all their BHs (the highest density GCs) and a
population of extended GCs (like the Palomar GCs) which fill the Roche volume
(the lowest density GCs) (Baumgardt et al., 2010). Lower density GCs would
also have a higher mass-loss rate, and therefore there may exist a survival
bias in the clusters visible at the present-day, biasing our sample. It must
also be noted that the initial $f_{\mathrm{BH}}$ is highly dependent on the
IMF of the clusters. The IMF of our clusters is likely bottom-light compared
to a Kroupa IMF (see Paper I and Baumgardt et al., 2023). This would lead to
an increased initial $f_{\mathrm{BH}}$, and thus require an even higher
density in order to reach the critical fraction.

As the $f_{\mathrm{BH}}$ we see in our clusters is well below the critical
mass fraction noted by Gieles et al. (2021), it is clear that the
$f_{\mathrm{BH}}$ in these clusters must have been below this critical amount
if/when they became tidally limited. The clusters can thus be expected to
continue to evolve towards BH free clusters.

Early theoretical expectations of the behaviour of GCs and their compact BH
subsystem suggested that the subsystem would become dynamically decoupled from
the rest of the stars in the cluster, succumbing to the Spitzer Instability
(Spitzer, 1969, 1987), and rapidly ejecting nearly all BHs. However, more
recent modelling work has shown that clusters do not necessarily become
Spitzer unstable (Morscher et al., 2013, 2015). In our models, a large
fraction ($24/34$) of the clusters remain Spitzer stable (according to the
classification given in Spitzer 1987), consisting mostly of clusters with
$f_{\mathrm{BH}}\lesssim 0.1$ per cent. Breen & Heggie (2013a), examining an
idealized, two-component model of clusters, suggested a relationship between
the mass fraction in BHs and the relative size of the BH subsystem,
proportional to $f_{\mathrm{BH}}^{3/2}$ for Spitzer unstable systems and
$f_{\mathrm{BH}}$ for stable systems. While it is difficult to compare
directly with the fits on this relation by their simple models, they do
predict a BH subsystem size of $r_{h,\mathrm{BH}}/r_{h}\sim 0.04$ near
$f_{\mathrm{BH}}\sim 0.1$ per cent, which is similar to the bulk of our
clusters.

Baumgardt et al. (2023), examining a suite of $N$-body models compared to
observations, demonstrated that the trend seen between the global low-mass
stellar mass function slopes, as derived from observations, and the dynamical
age of clusters could be reproduced only by models with a maximum BH retention
fraction of 30 per cent (immediately after natal kicks), ruling out a high
initial retention rate. Their models with higher initial BH retention
fractions (which suppresses mass segregation and the preferential loss of low-
mass stars) cannot reproduce the trend since they are not able to produce
clusters strongly depleted in low-mass stars. After the subsequent BH
hardening and dynamical ejections, they estimate an average surviving BH mass
fraction of $f_{\mathrm{BH}}\sim 0.03-0.1$ per cent at the present day, in
generally good agreement with many of our clusters.

It should be noted that any analysis of our results with respect to the
complete population of Milky Way GCs could be biased by the choice of clusters
examined. The sample of clusters chosen was limited primarily by the
availability of good quality data, especially mass function depth and radial
coverage, requiring adequate deep HST photometry. These criteria bias our
sample somewhat towards more massive and nearby clusters. In comparison to the
overall population of GCs (as given in Harris, 1996, 2010 edition), our chosen
clusters tend to be slightly more massive, have a smaller core radius and a
lower galactocentric radius.This could indeed bias our sample slightly towards
clusters with lower $f_{\mathrm{BH}}$, as they do not contain many low-density
or ‘fluffy’ outer halo GCs, which may have much higher $f_{\mathrm{BH}}$
(Gieles et al., 2021). Correspondingly, as was noted in Paper I, our sample is
deficient in the most metal-rich clusters, which are expected to produce
notably different remnant populations.

Figure 7: Relations between mass fraction in BHs (top panel) and mass fraction
in all dark remnants (bottom panel) with the remaining mass fraction for all
clusters, except for NGC 5139 ($M_{\rm{today}}/M_{\rm{initial}}=0.50$,
$f_{\mathrm{BH}}=5.9$ per cent, $f_{\mathrm{remn}}=50.3$ per cent). All core-
collapsed clusters, whose inferred BH populations may be unreliable, are
highlighted in red.

We can also examine how the populations of remnants in our clusters relate to
their dynamical age. Figure 7 shows the relationship between the fraction of
mass in BHs and the fraction of the cluster mass in all remnants (WD, NS and
BHs), against the dynamical age of the clusters, estimated based on the
remaining mass fraction, as was described in Paper I:

$\frac{M_{\mathrm{today}}}{M_{\mathrm{initial}}}=0.55\times\left(1-\frac{\mathrm{Age}}{\tau_{\mathrm{life}}}\right),$
(1)

where the factor 0.55 reflects the typically assumed mass loss from stellar
evolution of $\sim 45$ per cent of the initial cluster mass in the first Gyr
of a cluster’s evolution and the dissolution time $\tau_{\mathrm{life}}$
represents the estimated total lifetime of the cluster and is taken from
Baumgardt et al. (2023).

While we can see here that clusters with substantial populations of BHs tend
to be less evolved, there is no strong correlation between the BH mass
fraction and the dynamical age of the clusters, likely indicative of their
different initial conditions and the effects they would have on the evolution
of the BH population over time, as discussed above.

The evolution of the remnant mass fraction, which includes all types of
stellar remnants, shows a stronger relationship with the dynamical age of the
clusters, as might be expected, and has been previously reported by Sollima &
Baumgardt (2017). As a cluster evolves and loses mass, the mass lost is
preferentially in the form of lower-mass stars from the outer parts of the
cluster, rather than the heavier remnants, and as such the fraction of mass in
remnants should increase as the cluster’s low-mass stars are depleted.
Interestingly, some of the most dynamically evolved clusters have around 70
per cent of their mass in dark remnants at the present day, something worth
bearing in mind when interpreting the mass-to-light ratios and inferring
masses of unresolved GCs in distant galaxies.

These relatively high remnant mass fractions are in good agreement with the
results of Sollima & Baumgardt (2017), for the clusters in common, although it
should be noted that these authors adopt a simpler prescription for the mass
function of remnants and a fixed, canonical high-mass IMF, unlike the
flexibility allowed in our models for both of these quantities. The $N$-body
models of Baumgardt & Makino (2003) also showed that evolved clusters can
consist of up to nearly 70 per cent WDs, by mass, in line with our results. As
detailed in Paper I and Baumgardt et al. (2023), our dynamical masses are also
in excellent agreement with other recent works, and therefore the cluster
mass-to-light ratios implied by our results are in keeping with the recent
literature (e.g. Baumgardt et al., 2020). These mass-to-light ratios are
confined to a narrow range near $M/L_{V}\sim 2$ and have been shown to be
consistent with the the mass-to-light ratios predicted by stellar population
models once the depletion of low-mass stars is taken into account (Baumgardt
et al., 2020). Interestingly, the higher remnant mass fractions compensate the
lack of low-mass stars to maintain the mass-to-light ratios within this narrow
range, a finding consistent with the results of Sollima et al. (2012).

###  5.2 Comparison with literature results

Also shown in Figure 3 are comparisons with the distribution of BH mass
fraction, total mass in BHs, and/or number of BHs in our models with those of
Askar et al. (2018), Weatherford et al. (2020) and Rui et al. (2021c). In
order to estimate the BH mass fraction in a number of Milky Way GCs,
Weatherford et al. (2020) compared the amount of visible mass segregation in
these clusters to the anti-correlation found in Weatherford et al. (2018)
between the degree of mass segregation in a cluster and its BH population in
the Cluster Monte Carlo (CMC) catalogue of models. In similar fashion to their
analysis, we also scale their computed estimates of $f_{\mathrm{BH}}$ (based
on the median clustercentric mass segregation parameter $\Delta_{r50}$) by the
total cluster mass determined by our models in order to compare the total mass
in BHs. Askar et al. (2018) predicted the amount of BHs in a number of Milky
Way GCs based on the correlations found in Arca Sedda et al. (2018) between
the density of inner BH-subsystems and the central surface brightness of the
clusters in the Monte Carlo Cluster Simulator (MOCCA) survey database. A
somewhat analogous analysis may be found in Rui et al. (2021c), who matched
the surface brightness and velocity dispersion profiles of 26 Milky Way GCs to
a grid of CMC models and explored seven in more detail (see Section 3 for more
information). The number of BHs reported for the three clusters in common with
Rui et al. (2021c) is shown in Figure 3. Further analysis and comparison with
other studies of interesting individual clusters is presented below in Section
5.3

Figure 8: The (error weighted) residuals of the mass fraction in BHs between
literature sources (Weatherford et al. (2020) in red and Askar et al. (2018)
in orange) and our models, with respect to the mass fraction in BHs of our
models, for all clusters overlapping with our sample, excepting any core-
collapsed clusters. Weatherford et al. (2020) values are computed using the
median clustercentric mass segregation parameter $\Delta_{r50}$ (Table 1 of
Weatherford et al. 2020). Conversions between the total BH masses presented by
Askar et al. (2018) and BH mass fractions are computed using our total cluster
mass estimates.

The majority of our clusters agree well, within $2\sigma$, with the amount of
mass in BHs estimated in these studies, however notable discrepancies can be
seen between individual clusters.

In Figure 8, the residuals (normalized by the combined $1\sigma$
uncertainties) between our BH mass fraction results and those of the
literature sources discussed above are shown. We can see a clear trend showing
that for clusters where we infer small BH populations, we tend to predict
fewer BHs than other studies, while for clusters where we infer larger BH
populations, we tend to predict more BHs than previously reported in the
literature. This is especially pronounced in the comparisons with Askar et al.
(2018), where the clusters with smaller BH populations do not even agree
within $2\sigma$. In other words, the distributions of total BH masses between
clusters predicted by Askar et al. (2018) and Weatherford et al. (2020) are
somewhat flatter (i.e. less variation in BH mass fraction across the sample)
than our results. Despite the differences between specific clusters in our
samples, we do agree with the overall conclusion that, in general, the mass
fraction of BHs retained in clusters at the present day is small, between 0
and 1 per cent.

Our analysis of the BH populations of individual clusters may be more robust
than many of these literature results, which rely on general correlations
between models with only a few varied initial parameters, and are fit on only
a single observed property (mass segregation between three stellar populations
for Weatherford et al. (2020) and the central surface luminosity for Askar et
al. (2018)), whereas we self-consistently include the effect of BHs in our
fits of numerous cluster observables with many free parameters. In particular,
as noted in Weatherford et al. (2020), the correlations of Askar et al. (2018)
rely on a number of chained parametric fits, which may bias the final values,
and the models used to construct the correlations in Arca Sedda et al. (2018)
exclude those with $N_{\mathrm{BH}}<15$, which may lead to overpredictions in
their inferred numbers of BHs. Rui et al. (2021c) matched available CMC model
snapshots to profiles of surface brightness and velocity dispersion
observations, but this study is based on a limited grid of models with only a
few varied parameters. We are thus able to, in most clusters, place tighter
constraints on the mass in BHs through our fits. It should, however, be noted
again that our uncertainties account solely for the statistical uncertainties
on the parameter fits and could thus be underestimated by a factor of up to
$\sim 2.5$ (see Section 3). In addition, there are a number of astrophysical
processes that could possibly have a small effect on the inferred amount of
which we do not model, such as binaries and cluster rotation.

Another possible major source of differences between these results is the
(initial) mass function formulation, which was identified by Weatherford et
al. (2020) and Rui et al. (2021c) as a potential source of unexplored
uncertainty in their analysis. The freedom in the shape of our (initial) mass
function, on a per-cluster basis, allows us to best explore the population of
BHs and other heavy remnants, as well as their relative abundance compared to
lower-mass stars. The generally more bottom-light, low-mass mass function
slopes found in our fits (see Paper I) and in similar analyses (Baumgardt et
al., 2023), compared to the canonical Kroupa IMF assumed by Weatherford et al.
(2020) and others, may for example affect the mass segregation correlation
found in Weatherford et al. (2018) and thus the amount of mass in BHs found by
Weatherford et al. (2020). However, the exact effects of such a bottom-light
IMF on cluster evolution, and the evolution of BHs within the clusters, remain
to be further explored.

We can also place our BH results in the context of a number of studies that
have searched for the presence of an intermediate-mass black hole (IMBH) in
various Milky Way GCs, based on a number of indirect dynamical inferences
(e.g. Noyola et al., 2008; McNamara et al., 2012; Kızıltan et al., 2017;
Perera et al., 2017). All of these reported IMBH candidates have been
controversial, and in many cases have been rebutted, either by the
introduction of improved data, or by other, more plausible, physical
interpretations of the data (e.g. van der Marel & Anderson, 2010; Zocchi et
al., 2017, 2019; Gieles et al., 2018; Baumgardt et al., 2019b). Our models are
able to reproduce well the observables of all clusters in our sample without
the need for any central IMBH. It should, however, be noted that it is
currently not possible to self-consistently include an IMBH in our models to
compare directly and to explore any partial degeneracy possible between a
central IMBH and a central concentration of stellar-mass BHs (e.g. Lützgendorf
et al., 2013). It has also been shown that an IMBH with a mass fraction below
$\sim 2$ per cent of the cluster mass may produce similar dynamical effects as
a population of stellar-mass BHs (Aros & Vesperini, 2023). Therefore we cannot
preclude their existence with certainty with our models, given these
degeneracies. Some of the clusters previously claimed to host an IMBH are
discussed in more detail in Section 5.3 below.

###  5.3 Comments on individual clusters

We now compare our results with other dynamical studies of BHs in GCs, based
on a number of different methods, in certain particularly interesting
clusters.

####  5.3.1 NGC 5139

NGC 5139, or $\omega$ Cen, is the most massive Milky Way GC, and stands apart
from the population of Milky Way GCs due to its mass and stellar populations
(e.g. Harris, 1996, 2010 edition; Baumgardt et al., 2019a). It has been
suggested that $\omega$ Cen may not be a classical globular cluster, but
rather the possible remnant nuclear star cluster of an accreted and disrupted
dwarf galaxy (e.g. Meza et al., 2005). It has also been hypothesized to
harbour an elusive IMBH (Noyola et al., 2008; van der Marel & Anderson, 2010).

While our models are able to fit the large amount of available data very well,
it does not appear in many of the figures above given its significantly larger
inferred mass fraction in BHs. Our fits favour a mass fraction in BHs of
$5.9\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$ per cent
($182\,000\begin{subarray}{c}+5000\\\ -6000\end{subarray}\
\mathrm{M}_{\odot}$), largely concentrated within the central regions of the
cluster and driven, not by a top-heavy IMF producing more BHs initially
($\alpha_{3}=2.28\begin{subarray}{c}+0.08\\\ -0.02\end{subarray}$, close to
Salpeter; Paper I), but by a very large BH retention fraction
($35\begin{subarray}{c}+9\\\ -3\end{subarray}$ per cent). This amount of BHs
is substantially higher than any other cluster in our sample, but is in good
agreement with the results of other studies (e.g. Zocchi et al., 2019;
Baumgardt et al., 2019b).

Zocchi et al. (2019) modelled $\omega$ Cen using two-component (one
representing stellar-mass BHs and one capturing all other lower-mass remnants
and visible stars) limepy models. Our agreement with their results is
interesting, given our inclusion of the full mass spectrum, and our fitting of
the visible mass function, and reinforces the assertion of Zocchi et al.
(2019) that a two-component model is a reasonable approximation when modelling
$\omega$ Cen, given its large amount of BHs, long two-body relaxation time,
and young dynamical age.

Many claims have been made for the presence of an IMBH in the centre of
$\omega$ Cen. As in the studies of Zocchi et al. (2019) and Baumgardt et al.
(2019b), our models suggest that an IMBH is not needed to match the data,
however we are also limited by the extent of the kinematical data available in
the very centre of the cluster. As was also noted in Zocchi et al. (2019,
Figure 5), our models would be discrepant with the velocity dispersion
profiles of the different IMBH-containing models presented by Noyola et al.
(2008), van der Marel & Anderson (2010) and Baumgardt (2017) mostly within the
innermost few arcseconds of the cluster, where data is currently lacking. As
such we cannot say for certain whether some of the dark mass we find may
actually be in the form of a central IMBH, given the degeneracy between the
effects produced by such an IMBH and a central concentration of smaller BHs.

There is one caveat to our results; the Gaia proper motion anisotropy profile
shows that $\omega$ Cen transitions, at about
$20\text{\,}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{m}\mathrm{i}\mathrm{n}$ from
the centre, from being radially anisotropic to being slightly tangentially
anisotropic. Our limepy models are unable to reproduce any amount of
tangential anisotropy (Gieles & Zocchi, 2015), and thus cannot match this
feature. Instead, when tangentially biased anisotropy is present in our data.
the models will favour a mostly isotropic fit as a compromise between the
radial and tangential regimes (Peuten et al., 2017). There is a degeneracy
present between the degree of radial anisotropy in a cluster and its mass in
black holes (Zocchi et al., 2017), however the difference in the BH mass
fraction between the isotropic and radially anisotropic models of Zocchi et
al. (2019) is only on the order of $\sim 0.7$ per cent. While further
exploration of the effects of tangential anisotropy on mass models of $\omega$
Cen would be interesting, given our excellent fits of all other datasets, this
should have a negligible impact on the results presented here.

####  5.3.2 NGC 104

NGC 104, or 47 Tuc, is one of the nearest and most massive Milky Way GCs, and
as such has been extensively studied in the past. Recent modelling efforts
using Monte Carlo cluster models (Weatherford et al., 2020; Ye et al., 2022),
action-integral based DF models (Della Croce et al., 2023) and multimass
limepy models (Hénault-Brunet et al., 2020) have provided predictions on the
amount of BHs in the cluster, while candidate stellar and intermediate mass
BHs (Miller-Jones et al., 2015; Paduano et al., 2024) have been detected
within the cluster through combined radio and X-ray observations. As shown in
Figure 9, our models tend to favour a similar amount of mass in BHs
($f_{\mathrm{BH}}=0.046\begin{subarray}{c}+0.017\\\ -0.009\end{subarray}$ per
cent) to other modelling studies, and are consistent within $2\sigma$ with
all. It is notable that we are able to not only constrain an upper limit on
the mass in BHs (such as was presented in Hénault-Brunet et al. 2020), but
also now place clear and tight bounds on it, thanks to the updated treatment
of BHs in our multimass models and the updated mass function data we are
using.

Figure 9: Posterior probability distributions of the BH mass fraction, total
mass and number of BHs in 47 Tuc. The results (median and $1\sigma$) of
various recently inferred values from the literature are shown in the bottom
panels.

It was postulated by Kızıltan et al. (2017) that 47 Tuc may host an IMBH of
around , based on the analysis of the accelerations of millisecond pulsars in
the cluster and comparisons with $N$-body simulations. However follow-up
studies using equilibrium models fit to various cluster observables (Hénault-
Brunet et al., 2020; Mann et al., 2019, although see Mann et al., 2020)
determined that there was no need for an IMBH to explain the observations, and
that a central concentration of less-massive dark remnants could explain the
data. Della Croce et al. (2023) placed an upper limit on the dark mass in the
centre of 47 Tuc at , based on the 3D kinematics of inner stars. These
conclusions are again reinforced by our results, which favour a smaller
central concentration of stellar-mass black holes, alongside a population of
white dwarfs and neutron stars.

####  5.3.3 NGC 6397

NGC 6397 is a metal-poor, core-collapsed Milky Way GC at a very short
heliocentric distance ($\sim$2.4\text{\,}\rm\mathrm{k}\mathrm{p}\mathrm{c}$$;
Harris, 1996, 2010 edition), which has been well studied in the past. Kamann
et al. (2016) first showed that models including an IMBH or very centrally
concentrated cluster of stellar-mass BHs of
$\sim$600\text{\,}\mathrm{M}_{\odot}$$ could best reproduce the central
kinematics of this cluster. Vitral & Mamon (2021) showed, in turn, that Jeans
models with more robust proper motion fitting disfavoured an IMBH, and instead
proposed an inner subcluster of unresolved dark remnants measuring $\sim
1000\rm{-}2000\ \mathrm{M}_{\odot}$, which they suggested is dominated by
stellar-mass BHs. However, Rui et al. (2021b, c) demonstrated, through fits of
CMC models, that no BHs were required to explain the kinematics of NGC 6397
and that the inner density profile of the cluster argues against the presence
of BHs. This is in line with the core-collapsed nature of NGC 6397 and
reinforced by the mass segregation based estimates of Weatherford et al.
(2020). These results suggest instead that the central dark subcluster could
be made up largely of white dwarfs (Kremer et al., 2021). A subsequent re-
examination of the Jeans modelling of NGC 6397 by Vitral et al. (2022), with
updated proper motion datasets, lowered the claimed mass of the central “dark”
cluster to $\sim$800\text{\,}\mathrm{M}_{\odot}$$, and concurred with a
subcluster dominated by white dwarfs, instead of stellar-mass BHs.

Our best-fitting models of NGC 6397, despite the caveats of modelling core-
collapsed clusters discussed in Section 4.1, favour a negligible population of
black holes ($f_{\mathrm{BH}}=0.019\begin{subarray}{c}+0.001\\\
-0.001\end{subarray}$ per cent), consistent with the results of Weatherford et
al. (2020) and Rui et al. (2021b). Our models also favour a large population
of white dwarfs and neutron stars dominating the core of the cluster, and it
is clear that they concur with the general consensus that NGC 6397 hosts a
massive central concentration of WDs, and little to no BHs, nor an IMBH.

####  5.3.4 NGC 3201

NGC 3201 is a nearby Milky Way GC which has a notably low and flat core
density profile (i.e. far from core-collapsed), and is the host of three
confirmed stellar-mass black hole candidates in detached binaries (Giesers et
al., 2018, 2019).

CMC models of NGC 3201 (Kremer et al., 2018, 2019) suggested that models with
$\sim 120$ BHs were best able to recreate the velocity dispersion and surface
brightness profiles, in general agreement with the results of Askar et al.
(2018) and the inner subcluster of dark remnants found by Vitral et al.
(2022). Weatherford et al. (2020) in turn favoured a slightly lower, but still
consistent, $\sim 44$ BHs ($f_{\mathrm{BH}}\sim 0.6$ per cent), with large
uncertainties. In contrast, our best-fitting models of NGC 3201 favour a
remarkably small amount of BHs, with the distribution peaking at 0 BH (95%
probability of containing less than 10 BHs,
$f_{\mathrm{BH}}=0.007\begin{subarray}{c}+0.02\\\ -0.007\end{subarray}$ per
cent). This is somewhat surprising, given the literature results and the shape
of the cluster density profile, but is technically in agreement (within
errors) with the results from the literature, and follows the trend in our
results of predicting fewer BHs than other studies in this regime (see Figure
8).

As was noted in Section 3, NGC 3201 stands out in our sample as a cluster with
a long relaxation time and young dynamical age which was best fit by a high
$\delta$ value (and thus low number of BHs). As discussed above, this is
unexpected for dynamically young clusters, and is likely indicative that this
cluster falls within a regime our models may struggle to capture. Therefore
the results for this cluster based on our multimass model fit should be taken
with extra caution. NGC 3201 is the cluster in our sample for which this issue
is potentially the most severe.

####  5.3.5 NGC 6121

NGC 6121, or M 4, is the nearest GC ($d\sim$1.85\text{\,}\mathrm{kpc}$$;
Baumgardt & Vasiliev, 2021) and as such has been extensively observed.
Recently, Vitral et al. (2023) utilized Jeans modelling of M 4, followed up by
a comparison with Monte Carlo models, to fit on HST and Gaia kinematic data
and suggest an excess of dark mass of around
$800\pm$300\text{\,}\mathrm{M}_{\odot}$$, concentrated within the inner
0.016-$0.034\text{\,}\mathrm{pc}$ of the cluster. They also explore the
possibility that this mass is concentrated in a single IMBH. Our models have a
comparable amount of dark mass within the central regions (dominated almost
entirely by white dwarfs) but in a less concentrated mass profile, reaching a
similar cumulative mass of around near $0.1\text{\,}\mathrm{pc}$ instead. Our
models also suggest a small population of BHs, totalling around
$80\begin{subarray}{c}+50\\\ -30\end{subarray}\ \mathrm{M}_{\odot}$
($f_{\mathrm{BH}}=0.09\begin{subarray}{c}+0.05\\\ -0.03\end{subarray}$ per
cent). This mass in BHs is significantly smaller than the dark mass suggested
by the Jeans models of Vitral et al. (2023), however, notably, this is
actually in good agreement with the best matching CMC model they also
presented.

One key difference between our analyses lies in the data used. Vitral et al.
(2023) introduce new HST proper motion data which reaches deeper into the core
of the cluster than the Libralato et al. (2022) data we use and their velocity
dispersion profile seems to increase toward the centre in the inner few
arcseconds, albeit with very large uncertainties. However, given these large
uncertainties and our very good fits to all other datasets, it is unlikely
that our results would change significantly with the inclusion of these new
datapoints.

##  6 Conclusions

In this work, we have utilized the best-fitting multimass models of 34 GCs,
first presented in Paper I, to explore the BH and remnant populations of a
large sample of Milky Way clusters, yielding a number of important
conclusions:

  1. 1.

The models allow us to infer best-fitting, posterior probability distributions
for the total mass, number and mass fraction of BHs in our sample of clusters.
These results indicate that a large number of the GCs are consistent with
hosting little to no BHs, with the largest BH populations reaching masses in
BHs up to a few thousand $\mathrm{M}_{\odot}$ and mass fractions of around 1
per cent (save for $\omega$ Cen, for which $f_{\mathrm{BH}}\sim 5$ per cent).

  2. 2.

We find an anti-correlation between the BH mass fraction and the $\delta$
parameter, a proxy of mass segregation, with clusters having little mass in
BHs congregating around $\sim 0.5$, while $\delta$ is lower (closer to 0.3)
and mass segregation is increasingly suppressed in clusters with more
substantial BH populations, in agreement with the findings of Peuten et al.
(2017).

  3. 3.

As the $f_{\mathrm{BH}}$ we see in the clusters of our sample are well below
the critical mass fraction of $\sim 2.5$ per cent noted by Gieles et al.
(2021), these clusters are expected to continue to evolve towards a point when
they will have ejected all their BHs. The inferred present-day
$f_{\mathrm{BH}}$ encode information about the dynamical evolution and initial
density of GCs that can be used in future work to infer the initial conditions
of the population of Milky Way GCs.

  4. 4.

A clear correlation is also found between the dynamical age of the clusters
and the overall remnant mass fraction, which increases as clusters evolve and
lose low-mass stars. Our results show that the most evolved GCs in our sample
are made up of around 70 per cent dark remnants, by mass, at the present day.

  5. 5.

We find typically good agreement overall, within uncertainties, between our
results and those of other studies inferring BH populations in GCs in the
literature (e.g. Askar et al., 2018; Weatherford et al., 2020), but with
notable discrepancies between individual clusters. Our inferred masses in BHs
are, generally, slightly smaller (larger) than these studies in clusters with
small (large) BH populations.

  6. 6.

Closer inspection of a number of interesting clusters with previous claims of
hosting an elusive IMBH reveal no need for such an object to explain the large
amount of data used in our model fitting.

## Acknowledgements

ND is grateful for the support of the Durland Scholarship in Graduate
Research. VHB acknowledges the support of the Natural Sciences and Engineering
Research Council of Canada (NSERC) through grant RGPIN-2020-05990. MG
acknowledges support from the Ministry of Science and Innovation
(EUR2020-112157, PID2021-125485NB-C22, CEX2019-000918-M funded by
MCIN/AEI/10.13039/501100011033) and from AGAUR (SGR-2021-01069).

This research was enabled in part by support provided by ACENET (www.ace-
net.ca) and the Digital Research Alliance of Canada (https://alliancecan.ca).

This work has also benefited from a variety of Python packages including
astropy (Astropy Collaboration et al., 2013, 2018), dynesty (Speagle, 2020),
emcee (Foreman-Mackey, 2016), h5py (Collette et al., 2022), cmctoolkit (Rui et
al., 2021a), JAX (Bradbury et al., 2018), blackjax (Lao & Louf, 2020),
matplotlib (Hunter, 2007), numpy (Harris et al., 2020), scipy (Virtanen et
al., 2020) and shapely (Gillies et al., 2022).

## Data Availability

The data underlying this article are available at
https://github.com/nmdickson/GCfit-results.

## References

  * Abbate et al. (2019) Abbate F., Possenti A., Colpi M., Spera M., 2019, ApJ, 884, L9
  * Abbott et al. (2016) Abbott B. P., et al., 2016, ApJ, 818, L22
  * Alessandrini et al. (2016) Alessandrini E., Lanzoni B., Ferraro F. R., Miocchi P., Vesperini E., 2016, ApJ, 833, 252
  * Antonini & Gieles (2020a) Antonini F., Gieles M., 2020a, Phys. Rev. D, 102, 123016
  * Antonini & Gieles (2020b) Antonini F., Gieles M., 2020b, MNRAS, 492, 2936
  * Antonini et al. (2023) Antonini F., Gieles M., Dosopoulou F., Chattopadhyay D., 2023, MNRAS, 522, 466
  * Arca Sedda et al. (2018) Arca Sedda M., Askar A., Giersz M., 2018, MNRAS, 479, 4652
  * Aros & Vesperini (2023) Aros F. I., Vesperini E., 2023, arXiv e-prints, p. arXiv:2308.03845
  * Askar et al. (2018) Askar A., Arca Sedda M., Giersz M., 2018, MNRAS, 478, 1844
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Balbinot & Gieles (2018) Balbinot E., Gieles M., 2018, MNRAS, 474, 2479
  * Banerjee & Kroupa (2011) Banerjee S., Kroupa P., 2011, ApJ, 741, L12
  * Banerjee et al. (2020) Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R., Spurzem R., Wang L., 2020, A&A, 639, A41
  * Baumgardt (2017) Baumgardt H., 2017, MNRAS, 464, 2174
  * Baumgardt & Hilker (2018) Baumgardt H., Hilker M., 2018, MNRAS, 478, 1520
  * Baumgardt & Makino (2003) Baumgardt H., Makino J., 2003, MNRAS, 340, 227
  * Baumgardt & Vasiliev (2021) Baumgardt H., Vasiliev E., 2021, MNRAS, 505, 5957
  * Baumgardt et al. (2010) Baumgardt H., Parmentier G., Gieles M., Vesperini E., 2010, MNRAS, 401, 1832
  * Baumgardt et al. (2019a) Baumgardt H., Hilker M., Sollima A., Bellini A., 2019a, MNRAS, 482, 5138
  * Baumgardt et al. (2019b) Baumgardt H., et al., 2019b, MNRAS, 488, 5340
  * Baumgardt et al. (2020) Baumgardt H., Sollima A., Hilker M., 2020, Publ. Astron. Soc. Australia, 37, e046
  * Baumgardt et al. (2023) Baumgardt H., Hénault-Brunet V., Dickson N., Sollima A., 2023, MNRAS, 521, 3991
  * Bradbury et al. (2018) Bradbury J., et al., 2018, JAX: composable transformations of Python+NumPy programs, http://github.com/google/jax
  * Breen & Heggie (2013a) Breen P. G., Heggie D. C., 2013a, MNRAS, 432, 2779
  * Breen & Heggie (2013b) Breen P. G., Heggie D. C., 2013b, MNRAS, 436, 584
  * Chatterjee et al. (2013) Chatterjee S., Umbreit S., Fregeau J. M., Rasio F. A., 2013, MNRAS, 429, 2881
  * Chatterjee et al. (2017) Chatterjee S., Rodriguez C. L., Rasio F. A., 2017, ApJ, 834, 68
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Collette et al. (2022) Collette A., et al., 2022, h5py: 3.7.0, doi:10.5281/zenodo.6575970
  * Dalgleish et al. (2020) Dalgleish H., et al., 2020, MNRAS, 492, 3859
  * Della Croce et al. (2023) Della Croce A., Pascale R., Giunchi E., Nipoti C., Cignoni M., Dalessandro E., 2023, arXiv e-prints, p. arXiv:2310.15221
  * Dickson et al. (2023) Dickson N., Hénault-Brunet V., Baumgardt H., Gieles M., Smith P. J., 2023, MNRAS, 522, 5320
  * Djorgovski & King (1986) Djorgovski S., King I. R., 1986, ApJ, 305, L61
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Dotter et al. (2007) Dotter A., Chaboyer B., Jevremović D., Baron E., Ferguson J. W., Sarajedini A., Anderson J., 2007, AJ, 134, 376
  * Dotter et al. (2008) Dotter A., Chaboyer B., Jevremović D., Kostov V., Baron E., Ferguson J. W., 2008, ApJS, 178, 89
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, The Journal of Open Source Software, 1, 24
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2208.00211
  * Gieles & Gnedin (2023) Gieles M., Gnedin O. Y., 2023, MNRAS, 522, 5340
  * Gieles & Zocchi (2015) Gieles M., Zocchi A., 2015, mnras, 454, 576
  * Gieles et al. (2011) Gieles M., Heggie D. C., Zhao H., 2011, MNRAS, 413, 2509
  * Gieles et al. (2018) Gieles M., Balbinot E., Yaaqib R. I. S. M., Hénault-Brunet V., Zocchi A., Peuten M., Jonker P. G., 2018, MNRAS, 473, 4832
  * Gieles et al. (2021) Gieles M., Erkal D., Antonini F., Balbinot E., Peñarrubia J., 2021, Nature Astronomy, 5, 957
  * Giesers et al. (2018) Giesers B., et al., 2018, MNRAS, 475, L15
  * Giesers et al. (2019) Giesers B., et al., 2019, A&A, 632, A3
  * Gillies et al. (2022) Gillies S., van der Wel C., Van den Bossche J., Taves M. W., Arnott J., Ward B. C., et al., 2022, Shapely, doi:10.5281/zenodo.7263102
  * Harris (1996) Harris W. E., 1996, AJ, 112, 1487
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hénault-Brunet et al. (2019) Hénault-Brunet V., Gieles M., Sollima A., Watkins L. L., Zocchi A., Claydon I., Pancino E., Baumgardt H., 2019, mnras, 483, 1400
  * Hénault-Brunet et al. (2020) Hénault-Brunet V., Gieles M., Strader J., Peuten M., Balbinot E., Douglas K. E. K., 2020, MNRAS, 491, 113
  * Hénon (1961) Hénon M., 1961, Ann. Astrophys., 24, 369; English trans.: arXiv:1103.3499
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Hurley (2007) Hurley J. R., 2007, MNRAS, 379, 93
  * Hurley et al. (2016) Hurley J. R., Sippel A. C., Tout C. A., Aarseth S. J., 2016, Publ. Astron. Soc. Australia, 33, e036
  * Kamann et al. (2014) Kamann S., Wisotzki L., Roth M. M., Gerssen J., Husser T. O., Sandin C., Weilbacher P., 2014, A&A, 566, A58
  * Kamann et al. (2016) Kamann S., et al., 2016, A&A, 588, A149
  * Kamann et al. (2018) Kamann S., et al., 2018, MNRAS, 473, 5591
  * Kızıltan et al. (2017) Kızıltan B., Baumgardt H., Loeb A., 2017, Nature, 542, 203
  * Kremer et al. (2018) Kremer K., Ye C. S., Chatterjee S., Rodriguez C. L., Rasio F. A., 2018, ApJ, 855, L15
  * Kremer et al. (2019) Kremer K., Chatterjee S., Ye C. S., Rodriguez C. L., Rasio F. A., 2019, ApJ, 871, 38
  * Kremer et al. (2020a) Kremer K., et al., 2020a, ApJS, 247, 48
  * Kremer et al. (2020b) Kremer K., Ye C. S., Chatterjee S., Rodriguez C. L., Rasio F. A., 2020b, in Bragaglia A., Davies M., Sills A., Vesperini E., eds, Vol. 351, Star Clusters: From the Milky Way to the Early Universe. pp 357–366, doi:10.1017/S1743921319007269
  * Kremer et al. (2021) Kremer K., Rui N. Z., Weatherford N. C., Chatterjee S., Fragione G., Rasio F. A., Rodriguez C. L., Ye C. S., 2021, ApJ, 917, 28
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Lao & Louf (2020) Lao J., Louf R., 2020, Blackjax: A sampling library for JAX, http://github.com/blackjax-devs/blackjax
  * Libralato et al. (2022) Libralato M., et al., 2022, ApJ, 934, 150
  * Lindegren et al. (2021) Lindegren L., et al., 2021, A&A, 649, A2
  * Lützgendorf et al. (2013) Lützgendorf N., et al., 2013, A&A, 552, A49
  * Mann et al. (2019) Mann C. R., et al., 2019, ApJ, 875, 1
  * Mann et al. (2020) Mann C. R., et al., 2020, ApJ, 893, 86
  * McNamara et al. (2012) McNamara B. J., Harrison T. E., Baumgardt H., Khalaj P., 2012, ApJ, 745, 175
  * Meza et al. (2005) Meza A., Navarro J. F., Abadi M. G., Steinmetz M., 2005, MNRAS, 359, 93
  * Miller-Jones et al. (2015) Miller-Jones J. C. A., et al., 2015, MNRAS, 453, 3918
  * Miocchi et al. (2013) Miocchi P., et al., 2013, ApJ, 774, 151
  * Morscher et al. (2013) Morscher M., Umbreit S., Farr W. M., Rasio F. A., 2013, ApJ, 763, L15
  * Morscher et al. (2015) Morscher M., Pattabiraman B., Rodriguez C., Rasio F. A., Umbreit S., 2015, ApJ, 800, 9
  * Noyola et al. (2008) Noyola E., Gebhardt K., Bergmann M., 2008, ApJ, 676, 1008
  * Paduano et al. (2024) Paduano A., et al., 2024, ApJ, 961, 54
  * Perera et al. (2017) Perera B. B. P., et al., 2017, MNRAS, 468, 2114
  * Peuten et al. (2016) Peuten M., Zocchi A., Gieles M., Gualandris A., Hénault-Brunet V., 2016, MNRAS, 462, 2333
  * Peuten et al. (2017) Peuten M., Zocchi A., Gieles M., Hénault-Brunet V., 2017, MNRAS, 470, 2736
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Rodriguez et al. (2016) Rodriguez C. L., Morscher M., Wang L., Chatterjee S., Rasio F. A., Spurzem R., 2016, MNRAS, 463, 2109
  * Rodriguez et al. (2022) Rodriguez C. L., et al., 2022, ApJS, 258, 22
  * Rui et al. (2021a) Rui N. Z., Kremer K., Weatherford N. C., Chatterjee S., Rasio F. A., Rodriguez C. L., Ye C. S., 2021a, NicholasRui/cmctoolkit: First release, doi:10.5281/zenodo.4579951, https://doi.org/10.5281/zenodo.4579951
  * Rui et al. (2021b) Rui N. Z., Weatherford N. C., Kremer K., Chatterjee S., Fragione G., Rasio F. A., Rodriguez C. L., Ye C. S., 2021b, Research Notes of the American Astronomical Society, 5, 47
  * Rui et al. (2021c) Rui N. Z., Kremer K., Weatherford N. C., Chatterjee S., Rasio F. A., Rodriguez C. L., Ye C. S., 2021c, ApJ, 912, 102
  * Sigurdsson & Hernquist (1993) Sigurdsson S., Hernquist L., 1993, Nature, 364, 423
  * Sollima & Baumgardt (2017) Sollima A., Baumgardt H., 2017, MNRAS, 471, 3668
  * Sollima et al. (2012) Sollima A., Bellazzini M., Lee J. W., 2012, ApJ, 755, 156
  * Sollima et al. (2016) Sollima A., et al., 2016, MNRAS, 462, 1937
  * Speagle (2020) Speagle J. S., 2020, MNRAS, 493, 3132
  * Spera et al. (2015) Spera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086
  * Spitzer (1969) Spitzer Lyman J., 1969, ApJ, 158, L139
  * Spitzer (1987) Spitzer L. S., 1987, Dynamical Evolution of Globular Clusters. Princeton University Press
  * Strader et al. (2012) Strader J., Chomiuk L., Maccarone T. J., Miller-Jones J. C. A., Seth A. C., 2012, Nature, 490, 71
  * Torniamenti et al. (2023) Torniamenti S., Gieles M., Penoyre Z., Jerabkova T., Wang L., Anders F., 2023, MNRAS, 524, 1965
  * Trager et al. (1995) Trager S. C., King I. R., Djorgovski S., 1995, AJ, 109, 218
  * Vasiliev & Baumgardt (2021) Vasiliev E., Baumgardt H., 2021, MNRAS, 505, 5978
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vitral & Mamon (2021) Vitral E., Mamon G. A., 2021, A&A, 646, A63
  * Vitral et al. (2022) Vitral E., Kremer K., Libralato M., Mamon G. A., Bellini A., 2022, MNRAS, 514, 806
  * Vitral et al. (2023) Vitral E., Libralato M., Kremer K., Mamon G. A., Bellini A., Bedin L. R., Anderson J., 2023, MNRAS, 522, 5740
  * Watkins et al. (2015) Watkins L. L., van der Marel R. P., Bellini A., Anderson J., 2015, ApJ, 803, 29
  * Weatherford et al. (2018) Weatherford N. C., Chatterjee S., Rodriguez C. L., Rasio F. A., 2018, ApJ, 864, 13
  * Weatherford et al. (2020) Weatherford N. C., Chatterjee S., Kremer K., Rasio F. A., 2020, ApJ, 898, 162
  * Ye et al. (2022) Ye C. S., Kremer K., Rodriguez C. L., Rui N. Z., Weatherford N. C., Chatterjee S., Rasio F. A., Fragione G., 2022, ApJ, 931, 84
  * Zocchi et al. (2017) Zocchi A., Gieles M., Hénault-Brunet V., 2017, MNRAS, 468, 4429
  * Zocchi et al. (2019) Zocchi A., Gieles M., Hénault-Brunet V., 2019, MNRAS, 482, 4713
  * de Boer et al. (2019) de Boer T. J. L., Gieles M., Balbinot E., Hénault-Brunet V., Sollima A., Watkins L. L., Claydon I., 2019, MNRAS, 485, 4906
  * van der Marel & Anderson (2010) van der Marel R. P., Anderson J., 2010, ApJ, 710, 1063

##  Appendix A Results of the validation of BH population inference

CMC Model | Snapshot |  Closest Milky Way Analogue | $f_{\mathrm{BH,true}}\ \left[\%\right]$ | $f_{\mathrm{BH}}\ \left[\%\right]$
---|---|---|---|---
N4e5-rv1-rg2-Z0.0002 | 248 | NGC 6558 | 0.0 | $0.000\begin{subarray}{c}+0.003\\\ -0.000\end{subarray}$
N4e5-rv1-rg8-Z0.0002 | 217 | NGC 5946 | 0.008 | $0\begin{subarray}{c}+0\\\ -0\end{subarray}$
N4e5-rv1-rg8-Z0.0002 | 217 | NGC 5986 | 0.008 | $0.03\begin{subarray}{c}+0.02\\\ -0.02\end{subarray}$
N4e5-rv1-rg8-Z0.0002 | 252 | NGC 6544 | 0.015 | $0.004\begin{subarray}{c}+0.007\\\ -0.004\end{subarray}$
N4e5-rv2-rg2-Z0.002 | 165 | NGC 6352 | 0.152 | $0.15\begin{subarray}{c}+0.08\\\ -0.07\end{subarray}$
N4e5-rv2-rg8-Z0.0002∗ | 160 | NGC 6584 | 0.222 | $0.13\begin{subarray}{c}+0.06\\\ -0.05\end{subarray}$
N4e5-rv2-rg8-Z0.0002 | 165 | NGC 6981 | 0.147 | $0.13\begin{subarray}{c}+0.08\\\ -0.07\end{subarray}$
N4e5-rv4-rg8-Z0.0002∗ | 128 | NGC 5897 | 0.786 | $0.3\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$
N4e5-rv4-rg8-Z0.0002 | 134 | NGC 288 | 0.624 | $0.6\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$
N8e5-rv0.5-rg2-Z0.0002 | 608 | NGC 6355 | 0.0 | $0.019\begin{subarray}{c}+0.011\\\ -0.008\end{subarray}$
N8e5-rv0.5-rg2-Z0.0002 | 608 | NGC 6681 | 0.0 | $0.02\begin{subarray}{c}+0.01\\\ -0.01\end{subarray}$
N8e5-rv0.5-rg2-Z0.002 | 516 | NGC 6624 | 0.0 | $0\begin{subarray}{c}+0\\\ -0\end{subarray}$
N8e5-rv1-rg2-Z0.0002 | 375 | NGC 6293 | 0.064 | $0.10\begin{subarray}{c}+0.05\\\ -0.03\end{subarray}$
N8e5-rv1-rg2-Z0.0002 | 375 | NGC 6453 | 0.064 | $0.07\begin{subarray}{c}+0.01\\\ -0.02\end{subarray}$
N8e5-rv1-rg2-Z0.0002 | 375 | NGC 6522 | 0.064 | $0.24\begin{subarray}{c}+0.06\\\ -0.06\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 340 | NGC 6316 | 0.229 | $0.09\begin{subarray}{c}+0.02\\\ -0.03\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 340 | NGC 6539 | 0.229 | $0.13\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 345 | NGC 6569 | 0.18 | $0.09\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 345 | NGC 6637 | 0.18 | $0.10\begin{subarray}{c}+0.03\\\ -0.03\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 345 | NGC 6642 | 0.18 | $0.25\begin{subarray}{c}+0.10\\\ -0.09\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 354 | NGC 6638 | 0.179 | $0.09\begin{subarray}{c}+0.02\\\ -0.03\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 358 | NGC 6401 | 0.175 | $0.09\begin{subarray}{c}+0.02\\\ -0.02\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 362 | NGC 6256 | 0.155 | $0.09\begin{subarray}{c}+0.03\\\ -0.02\end{subarray}$
N8e5-rv1-rg2-Z0.002 | 362 | NGC 6304 | 0.155 | $0.09\begin{subarray}{c}+0.02\\\ -0.02\end{subarray}$
N8e5-rv1-rg20-Z0.0002 | 394 | NGC 1904 | 0.154 | $0.08\begin{subarray}{c}+0.02\\\ -0.04\end{subarray}$
N8e5-rv1-rg8-Z0.0002 | 379 | NGC 6284 | 0.188 | $0.09\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$
N8e5-rv1-rg8-Z0.0002 | 386 | NGC 6779 | 0.153 | $0.08\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$
N8e5-rv1-rg8-Z0.0002 | 386 | NGC 6934 | 0.153 | $0.09\begin{subarray}{c}+0.03\\\ -0.04\end{subarray}$
N8e5-rv1-rg8-Z0.002∗ | 366 | NGC 6760 | 0.234 | $0.18\begin{subarray}{c}+0.05\\\ -0.04\end{subarray}$
N8e5-rv2-rg2-Z0.0002 | 315 | NGC 6218 | 0.456 | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
N8e5-rv2-rg2-Z0.002∗ | 273 | NGC 6723 | 0.977 | $0.3\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
N8e5-rv2-rg2-Z0.002∗ | 286 | NGC 6712 | 0.798 | $0.17\begin{subarray}{c}+0.08\\\ -0.10\end{subarray}$
N8e5-rv2-rg20-Z0.0002∗ | 296 | NGC 1261 | 0.735 | $0.09\begin{subarray}{c}+0.07\\\ -0.07\end{subarray}$
N8e5-rv2-rg8-Z0.0002 | 299 | NGC 3201 | 0.694 | $0.5\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
N8e5-rv4-rg8-Z0.0002∗ | 194 | NGC 4372 | 2.218 | $0.3\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$
N8e5-rv4-rg8-Z0.0002∗ | 194 | NGC 6101 | 2.218 | $0.4\begin{subarray}{c}+0.3\\\ -0.3\end{subarray}$
N1.6e6-rv0.5-rg2-Z0.002 | 875 | NGC 6440 | 0.133 | $0.32\begin{subarray}{c}+0.04\\\ -0.05\end{subarray}$
N1.6e6-rv1-rg2-Z0.0002 | 694 | NGC 6402 | 0.898 | $0.8\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$
N1.6e6-rv1-rg2-Z0.0002 | 727 | NGC 6656 | 0.734 | $1.3\begin{subarray}{c}+0.2\\\ -0.2\end{subarray}$
N1.6e6-rv1-rg2-Z0.0002 | 737 | NGC 6541 | 0.679 | $0.8\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
N1.6e6-rv1-rg2-Z0.0002 | 755 | NGC 6333 | 0.584 | $0.9\begin{subarray}{c}+0.1\\\ -0.1\end{subarray}$
N1.6e6-rv1-rg2-Z0.0002 | 755 | NGC 6342 | 0.584 | $1.7\begin{subarray}{c}+0.4\\\ -0.3\end{subarray}$
N1.6e6-rv1-rg8-Z0.002∗ | 698 | NGC 6356 | 0.715 | $0.10\begin{subarray}{c}+0.06\\\ -0.08\end{subarray}$
Table 2: Results of the validation fits described in Section 3. The “True”
model values ($f_{\mathrm{BH,true}}$) and median and $1\sigma$ credibility
intervals inferred by the fits to mock observations ($f_{\mathrm{BH,infer}}$)
are shown. The model names (as given in the CMC cluster catalogue) indicate
the initial number of particles, virial radii, galactocentric radii and
metallicity of each CMC model. The Milky Way cluster to which each snapshot
was matched is also given. All dynamically young snapshots which fall into the
suspect region identified in Section 3 are denoted by an asterisk.