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

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 𝜔 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.


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; e.g.Spera et al. 2015) 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 ★ E-mail: nolan.dickson@smu.caand 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. 2013Morscher et al. , 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. 2018Giesers et al. , 2019)).Several stud-ies examining detailed evolutionary dynamical models (such as 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. 2018Weatherford et al. , 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. 2014Kamann et al. , 2016;;Abbate et al. 2019;Vitral & Mamon 2021;Vitral et al. 2022) and multimass distributionfunction (DF) based models (e.g.Sollima et al. 2012Sollima et al. , 2016;;Zocchi et al. 2019;Hénault-Brunet et al. 2019, 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-ofsight 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 (≥ 1 M ⊙ ) 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.

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.

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 (()) 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 φ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  (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  and 3D half-mass radius  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  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   .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 , such that   ∝  −   .The constituent discrete mass components which approximate the mass spectrum of a GC are represented in the multimass models by the total (  ) and mean (  ) 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  1 ,  2 ,  3 (with break masses at 0.5 and 1 M ⊙ and bounded between 0.1 and 100 M ⊙ ).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 library1 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 M ⊙ .
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 km s −1 (Hobbs et al. 2005) and scaled down by 1 −  fb , where  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  esc = 2 √︁ 2 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 (BH ret ).
Finally, the heliocentric distance to the GCs  is introduced as a free parameter, to allow for the conversion between projected, linear model units and the angular units of observations.

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 ∼ 0.16 −0.8 M ⊙ .
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 techniques2 .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 GCfit3 .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 (  BH =  BH / cluster ), are constructed based on the models representing the set of weighted posterior samples retrieved from the nested sampler.

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) 4 .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) 5 .We search for snapshots that are a good match to any clusters from the Harris (1996Harris ( , 2010 edition) 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 (1996Harris ( , 2010 edition) 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 simultaneously6 .In order to select a snapshot, we adopt a threshold of  ≡ max χ2 SBP , χ2 VDP < 10 for the "fitting heuristic"  of Rui et al. (2021c), which describes the goodness-of-fit of a snapshot based on the χ2 statistic between the observations and the interpolated model profiles.We have found that a threshold of  < 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.

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;Rui et al. 2021a) to calculate projected positions and velocities as well as simulated photometry for objects in each snapshot.

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  = 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  < 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.
from the CMC snapshots for comparison with the observations.

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 ′′ of the cluster to mimic the footprint of an HST ACS footprint, and limit our selection to stars within 15 <  < 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 mas/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 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 <  < 19 range outside of the innermost 100 ′′ , to avoid overlapping with the HST profiles.We assign each star an uncertainty in proper motion based on its  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.

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)  5 .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 km s −1 on all observed stars.We limit this dataset to only giants brighter than  = 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.

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 ′ -1.6 ′ in projected angular separation from the centre, and two outer annuli at radial distances of 2.5 ′ and 5 ′ , 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  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 0.1 M ⊙ .
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  = 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  parameter and its motivation).

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 ∼ 10 per cent and the half-mass radius, which we generally recover within uncertainties.Our inferred values of  BH , compared to the true values for our collection of snapshots, are shown in Figure 1 (right panel), Figure 2, and Table A1.
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  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  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  −  relax plane (where  relax is the ratio of cluster age to present-day half-mass relaxation time), coloured by the distance of their median inferred  BH value from the true value.As  is a measure of the degree of energy equipartition in our multimass models (with typical values of  ≈ 0.5 for mass-segregated, dynamically evolved clusters), we should normally expect dynamically young clusters (≲ 3  relax ) to have lower values of .The mock clusters for which we significantly underpredict  BH are clearly concentrated in the upper-left corner of the  −  relax plane, highlighted with the grey-shaded rectangle in the left panel of Figure 1.Their value of  is inconsistent with our expectations for these systems.
Previous comparisons of the multimass limepy models used here with -body models (Peuten et al. 2017) have confirmed this intuition that dynamically young clusters should have lower values of  (with  ∼ 0.35) and have also shown that models with any significant black hole population should have similarly low values of  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  −  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  close to 0.5.Simple test fits with  forced to a lower value ( ≤ 0.35) have shown that a lower value of  can result in an increase in the inferred  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  −  relax plane.In the right panel of Figure 1, we show our inferred values of  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  BH compared to the true values but now ignoring these problematic snapshots and showing the 1 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  and low  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 ( Cen), two of the most massive clusters in our sample, which still possess notable BH populations despite their high .Our results for these two wellstudied 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  −  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  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  by which we multiply the statistical errors   BH to determine the total, inflated error Δ on each inferred  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  to vary freely and inflate the uncertainties on the inferred  BH .This fit results in a value of  = 2.5 +0.3 −0.3 .We also show the inflated 1 error bars on the inferred  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.
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 ∼ 2.5, and with a reminder that dynamically young clusters with inferred values of  ∼ 0.5 should be treated with caution.

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.2020) are computed using the median clustercentric mass segregation parameter Δ 50 (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.MNRAS 000, 1-17 (2024) Figure 3 shows the posterior probability distributions of the mass fraction in BHs (  BH ), the total mass in BHs ( BH ), and total number of BHs ( BH ) inferred from our best-fitting models of most clusters in our sample.The median and 1 credibility intervals of the distributions of these quantities, as well as the mean individual BH mass ( mBH ) and the mass fraction in all dark remnants (  remn ), are also presented in Table 1.NGC 5139 ( Cen) is not included in Figure 3, due to the very large inferred amount of black holes (5.9 +0.2 −0.2 per cent), but it is discussed in more detail in Section 5.3.A large number of the clusters are consistent, within 2, with harbouring little to no black holes, while the remainder possess, on average, at most a few thousand M ⊙ of stellar-mass black holes, with constituent individual BH masses ranging between ∼ 5 to 15 M ⊙ .Other than  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  BH may be underestimated by a factor of ∼ 2.5.

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 corecollapsed 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  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 (∼ 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 Figure 3 and table 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.4of 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 ∼ 0.1 pc in these clusters).While, in most corecollapsed 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 BH 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 (  BH ∼ 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.
The original models (with BHs) are used throughout this paper and in all discussion of BHs, however the results for these corecollapsed 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).

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 shows the relationships between the high-mass initial mass function exponent  3 (≥ 1 M ⊙ ) from Paper I, the black hole retention fraction BH ret and the BH mass fraction  BH .This serves to demonstrate the role of the BH ret parameter, which is not directly proportional to the number of BHs.At high values of  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  3 and the mass fraction in BHs at the present day.The visible pattern instead relates to the relationship with BH ret ; to end up with similar mass fractions in BHs, clusters with higher  3 values produce fewer BHs initially while clusters with low  3 values produce more 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.
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  3 or BH ret , but merely helps to explain the distribution of BH ret .For a given cluster there is generally not a direct degeneracy between  3 and BH 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.We do find an interesting relationship between the mass fraction of BHs and the parameter , 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  ∼ 0.4−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 -body models with and without black holes, that the majority of mass-segregated clusters should converge to a value of ∼ 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  may be suppressed.

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  BH and the initial density of the cluster.Very early on in the lifetime of a cluster,  BH will rapidly increase, as stellar mass is lost due to stellar evolution and BHs are formed.This initial population of BHs, and  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  BH (Breen & Heggie 2013b;Gieles et al. 2021).The amount that  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  BH determines the subsequent evolution of the cluster mass.Theory (Breen & Heggie 2013b) and -body modelling (Gieles & Gnedin 2023) suggests that there exists a critical BH mass fraction  crit , at which the dynamical losses of BHs (relative to the total cluster mass loss) will exactly equal  BH , causing  BH to remain constant until the dissolution of the cluster.Given an  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 ∼ 10 per cent (Breen & Heggie 2013b) for idealized single-mass clusters filling their tidal radius, however recent -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  BH below 1 per cent (except for  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  BH ∼ 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  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  BH , and thus require an even higher density in order to reach the critical fraction.
As the  BH we see in our clusters is well below the critical mass fraction noted by Gieles et al. (2021), it is clear that the  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;Spitzer 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(Morscher et al. , 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  BH ≲ 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  3/2 BH for Spitzer unstable systems and  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  ℎ,BH / ℎ ∼ 0.04 near  BH ∼ 0.1 per cent, which is similar to the bulk of our clusters.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  BH ∼ 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 1996Harris , 2010 edition) 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  BH , as they do not contain many low-density or 'fluffy' outer halo GCs, which may have much higher  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.
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: where the factor 0.55 reflects the typically assumed mass loss from stellar evolution of ∼ 45 per cent of the initial cluster mass in the first Gyr of a cluster's evolution and the dissolution time  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 -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 /  ∼ 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).

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 anticorrelation 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  BH (based on the median clustercentric mass segregation parameter Δ 50 ) by the total cluster mass determined by our models in order to compare 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Weatherford et al., 2020Askar et al., 2018 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 Δ 50 (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 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 The majority of our clusters agree well, within 2, 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 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.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  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 ∼ 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. 2017Zocchi et al. , 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 ∼ 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.

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.

NGC 5139
NGC 5139, or  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 1996Harris , 2010 edition; edition;Baumgardt et al. 2019a).It has been suggested that  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 +0.2 −0.2 per cent (182 000 +5000 −6000 M ⊙ ), largely concentrated within the central regions of the cluster and driven, not by a top-heavy IMF producing more BHs initially ( 3 = 2.28 +0.08 −0.02 , close to Salpeter; Paper I), but by a very large BH retention fraction (35 +9 −3 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  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  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  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  Cen transitions, at about 20 arcmin 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 ∼ 0.7 per cent.While further exploration of the effects of tangential anisotropy on mass models of  Cen would be interesting, given our excellent fits of all other datasets, this should have a negligible impact on the results presented here.

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 (  BH = 0.046 +0.017 −0.009 per cent) to other modelling studies, and are consistent within 2 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.
It was postulated by Kızıltan et al. (2017) that 47 Tuc may host an IMBH of around 2300 M ⊙ , based on the analysis of the accelerations of millisecond pulsars in the cluster and comparisons with -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 578 M ⊙ , 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.

NGC 6397
NGC 6397 is a metal-poor, core-collapsed Milky Way GC at a very short heliocentric distance (∼ 2.4 kpc; Harris 1996Harris , 2010 edition) 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 ∼ 600 M ⊙ 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 ∼ 1000−2000 M ⊙ , 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 ∼ 800 M ⊙ , 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 (  BH = 0.019 +0.001 −0.001 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.

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(Giesers et al. , 2019)).
CMC models of NGC 3201 (Kremer et al. 2018(Kremer et al. , 2019) ) suggested that models with ∼ 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, ∼ 44 BHs (  BH ∼ 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,  BH = 0.007 +0.02 −0.007 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  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.

NGC 6121
NGC 6121, or M 4, is the nearest GC ( ∼ 1.85 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 ± 300 M ⊙ , concentrated within the inner 0.016-0.034pc 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 800 M ⊙ near 0.1 pc instead.Our models also suggest a small population of BHs, totalling around 80 +50 −30 M ⊙ (  BH = 0.09 +0.05 −0.03 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.

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: (i) 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 M ⊙ and mass fractions of around 1 per cent (save for  Cen, for which  BH ∼ 5 per cent).
(ii) We find an anti-correlation between the BH mass fraction and the  parameter, a proxy of mass segregation, with clusters having little mass in BHs congregating around ∼ 0.5, while  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).
(iii) As the  BH we see in the clusters of our sample are well below the critical mass fraction of ∼ 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  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.
(iv) 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.
(v) 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 indi-vidual clusters.Our inferred masses in BHs are, generally, slightly smaller (larger) than these studies in clusters with small (large) BH populations.
(vi) 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.

Figure 1 .
Figure 1.Left panel: Validation snapshots plotted in the  −  relax plane.A selection of dynamically young clusters with high inferred values of  is shown by the grey shaded region ( relax < 3,  > 0.4).The points are coloured based on the number of   BH the median inferred  BH is away from the true value (where   BH is the width of the posterior for  BH ).Under this colour scheme, snapshots for which we underpredict (overpredict)  BH are increasingly orange (blue).Right panel: The  BH values inferred based on the multimass model fits to mock observations extracted from CMC models, against the true values in those validation snapshots (  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  BH .

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

Figure 3 .
Figure3.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 1.82 × 10 5 M ⊙ (  BH ∼ 5.9 per cent), and is excluded in order to better highlight the distributions of the other clusters.The median and 1 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 inTrager et al. (1995) are denoted by an asterisk.The median and 1 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 fromWeatherford et al. (2020) are computed using the median clustercentric mass segregation parameter Δ 50 (Table1ofWeatherford et al. 2020), and any necessary conversions between total mass and mass fraction are computed using our total cluster mass estimates.MNRAS 000,1-17 (2024)

Figure 4 .
Figure4.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.

Figure 5 .Figure 6 .
Figure5.Relations between the high-mass initial mass function exponent ( 3 ), the black hole retention fraction parameter (BH ret ) and the mass fraction in black holes (  BH ) for all clusters, except for NGC 5139, which has a high BH ret of ∼ 35 +9 −3 per cent and a more typical  3 value of ∼ 2.28 +0.08 −0.02 .All core-collapsed clusters, whose inferred black hole populations may not be physical, are highlighted by a red outline.

Figure 7 .
Figure7.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 ( today / initial = 0.50,  BH = 5.9 per cent,  remn = 50.3per cent).All core-collapsed clusters, whose inferred BH populations may be unreliable, are highlighted in red.

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

Table 1 .
Trager et al. (1995)lity intervals of the BH mass fraction (  BH ), total mass in BHs ( BH ), total number of BHs ( BH ) and remnant mass fraction (  remn ) in all clusters in our sample, as well as the median mean BH mass ( mBH ).All clusters classified as core-collapsed inTrager 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).