Merger rate of black hole binaries from globular clusters: theoretical error bars and comparison to gravitational wave data from GWTC-2
MMerger rate of black hole binaries from globular clusters: theoretical error bars andcomparison to gravitational wave data from GWTC-2
Fabio Antonini ∗ and Mark Gieles
2, 3, † Gravity Exploration Institute, School of Physics and Astronomy,Cardiff University, Cardiff, CF24 3AA, United Kingdom ICREA, Pg. Llu´ıs Companys 23, E08010 Barcelona, Spain Institut de Ci`encies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB),Mart´ı Franqu`es 1, E08028 Barcelona, Spain (Dated: November 10, 2020)Black hole binaries formed dynamically in globular clusters are believed to be one of the mainsources of gravitational waves in the Universe. Here, we use our new population synthesis code, cBHBd , to determine the redshift evolution of the merger rate density and masses of black holebinaries formed in globular clusters. We simulate ∼ R ( z ) = R (1 + z ) κ at z <
2, and marginalizing over uncertainties, we find: R = 7 . +21 . − . Gpc − yr − and κ = 1 . +0 . − . (90% credibility). The rate parameters for binaries that merge inside the clusters are R , in =1 . +1 . − . Gpc − yr − and κ in = 2 . +1 . − . ; ∼
20% of these form as the result of a gravitational-wavecapture, implying that eccentric mergers from globular clusters contribute (cid:46) . − yr − to thelocal rate. A comparison to the merger rate reported by LIGO-Virgo shows that a scenario in whichmost of the detected black hole mergers are formed in globular clusters is consistent with currentconstraints, and requires initial cluster half-mass densities (cid:38) M (cid:12) pc − . Interestingly, thesemodels also reproduce the inferred primary black hole mass distribution in the range 13 − M (cid:12) .However, all models under-predict the data outside this range, suggesting that other mechanismsmight be responsible for the formation of these sources. I. INTRODUCTION
Several black hole (BH) binaries have been detectedby the Advanced LIGO and Virgo interferometers [1–10].The recently released Second Gravitational-Wave Tran-sient Catalog (GWTC-2), includes a total of 44 confidentBH binary (BHB) events [11, 12]. While the astrophys-ical origin of these sources is still unknown, one widelydiscussed possibility is that they formed in the dense coreof globular clusters (GCs) through dynamical three-bodyinteractions [13–17].The realistic modeling of the dynamical evolution ofBHs in the core of a GC represents a complex computa-tional challenge requiring an enormous dynamical rangein both space and time. For this reason, it is only very re-cently, thanks to major improvements in computationalmethods and hardware, that it became possible to makerobust predictions about the numbers and physical prop-erties of BH binary (BHB) mergers produced in GCs [18–22]. Thanks to these past efforts it is now clear thata large fraction of the sources detected by LIGO-Virgocould have been dynamically assembled in GCs. How-ever, as discussed below, a full characterisation of the ∗ antoninif@cardiff.ac.uk † [email protected] model uncertainties related to the BHB merger rate fromthe GC channel is still missing.Several recent studies only considered the contributionfrom clusters that have survived to the present day [e.g.,19, 23]. These studies found that the present-day popu-lation of GCs produces BHB mergers at a local rate of ≈ − yr − . This represents a lower limit to theactual merger rate as there likely existed a population ofclusters which did not survive to the present, but thatcontributed significantly to the local merger rate [24]. Infact, it is believed that the GC mass function (GCMF)today is the result of an initial GCMF that was shapedby dynamical processes [e.g., 25–30]. These processes,e.g., relaxation driven evaporation and tidal shocking,are particularly efficient at destroying low-mass clusters.A key uncertainty in estimating a merger rate from all GCs is that the amount of such disrupted clusters is notknown.Previous estimates for the BHB merger rate ignoredthe fact that the fractional mass that has been lost fromthe GC population by the present time, K (see equa-tion 6 below), is very uncertain as it cannot be tightlyconstrained from the present-day properties of the GCpopulation. We will show that once this uncertainty istaken properly into account, it becomes one of the dom-inant factors in setting the error bars on local mergerrate estimates from the GC channel. For example, both a r X i v : . [ a s t r o - ph . H E ] N ov Fragione and Kocsis [24] and Rodriguez and Loeb [31]used a single value for K which was derived under oneassumption for the initial GCMF. Although Rodriguezand Loeb [31] considered the effect of different GCMFson their results, they neglected that K should be relatedto the choice of initial GCMF and that its value can beconstrained by the present-day GCMF once evaporationmass loss is taken into account. A different initial GCMFnot only changes the mass of the GCs which make theBHBs, but it also sets the amount of mass that is lostfrom the GC system; and it turns out that the BHBmerger rate is quite sensitive to both effects.The properties and merger rate of BHBs depend onseveral other physical processes, many of which lackstrong observational constraints [32]. For example, thedistribution of the natal kicks controls the number of BHsthat are ejected from the GC upon formation, as well asthe fraction of BHs that retain their binary companionafter supernova. Different assumptions about the earlystages of BH formation will also reflect on the evolution ofthe host cluster, affecting its total lifetime and final prop-erties. Moreover, merger rates are expected to be sensi-tive to the assumed density and the related mass-radiusrelation of GCs at formation which is also unconstrainedobservationally [33]. The full implications of these un-certainties is still not fully explored. The main reasonfor this is that standard numerical techniques such as N -body and Monte Carlo simulations are still too slow toallow a full parameter space exploration. This is why inthis study we employ our new population synthesis code clusterBHBdynamics (hereafter cBHBd ) [34] to sys-tematically vary assumptions made for the model param-eters and over the full range of initial conditions relevantto real GCs. Thus, we examine the effect of these initialassumptions on the number and properties of mergingBHBs using a suite of about 20 million cluster models.In summary, the merger rate of BHBs produced dy-namically in GCs has been studied by multiple teams[e.g., 16, 19, 20, 35]. Here we build on former studiesin two ways which allow us to place error bars on the-oretical estimates for the BHB merger rate density andon its redshift evolution: (i) we constrain the fractionalmass that has been lost from GCs over cosmic time by fit-ting an evolved Schechter mass function to the observedGCMF in the Milky Way today, and using a simple modelfor cluster evaporation. (ii) we employ our new popu-lation synthesis code cBHBd to explore how the BHBmerger rate depends on uncertain parameters in the mod-els (e.g., initial cluster densities, BH formation recipes,natal kicks), and explore the parameter space that is rel-evant to real GCs and over all mass scales.The paper is organized as follows. In Section II wecompute the GC formation rate density as a function oftime using constraints from the present-day GCMF. InSection III we describe our population synthesis modeland detail the modifications we made to it with respectto the version used in [34]. Section IV describes our mainresults. We discuss the implications of our results and conclude in Section V. II. CLUSTER FORMATION RATE
In order to compute a BHB merger rate we need thecluster formation rate density (i.e. per unit of volume) asa function of time: ˙ ρ GC ( t ). We do this by imposing thatin our model: (i) the present-day GC mass density in theUniverse, ρ GC , is consistent with its empirically inferredvalue and (ii) the present-day GCMF is consistent withthe observed mass function of the Milky Way GCs. A. Globular clusters density in the Universe
To derive ρ GC we use the same approach as [18], whouse the empirically established relation between the totalmass of a GC population ( M GCs ) and the dark matterhalo mass of the host galaxy ( M h ). The ratio of thesetwo quantities is remarkable constant over a large rangeof halo masses (10 (cid:46) M h /M (cid:12) (cid:46) ) and for differentgalaxy types: η ≡ M GCs /M h (cid:39) (3 − × − [36–40].We can also estimate this ratio for the Milky Way: thetotal luminosity of all Milky Way GCs from the Harriscatalogue [41, 42] is 1 . × L V, (cid:12) . Adopting a mass-to-light ratio in the V -band of Υ V = 2 M (cid:12) /L V, (cid:12) and avirial mass of the Milky Way of M h = 1 . × M (cid:12) [43]we find η = 2 . × − for our Galaxy. Table 1 in [39]summarises 8 results from different studies. We use the7 studies that include at least 25 galaxies and combinethis with the result of η = 2 . × − by [40] that waspublished after this summary. We also add the MilkyWay estimate from above to find a mean value of (cid:104) η (cid:105) = (4 . ± . × − . (1)We determine the dark matter halo mass function fromsimulations of large scale structure formation by [44] us-ing the HMFcalc tool [45]. The total mass density in darkmatter halos with individual masses M h ≥ M (cid:12) /h is ρ DM = 3 . × h M (cid:12) Gpc − . Combined with ourvalue for η from equation (1) and h = 0 .
674 from thePlanck Collaboration [46] we find ρ GC = (cid:104) η (cid:105) ρ DM = (7 . ± . × M (cid:12) Gpc − . (2)The relation between M GCs and M h may hold down todwarf galaxy masses of M h (cid:39) M (cid:12) [47], and includingthese low-mass galaxies would increase ρ DM and ρ GC byabout 15%, but because of the uncertain GC occupationfraction below M h (cid:39) M (cid:12) , we continue with theresult of equation (2) . For an average GC mass of (cid:104) M (cid:105) = 3 × M (cid:12) this mass densityimplies a number density of n GC = 2 . ± . − and insection V we discuss how this compares to other studies. B. Globular cluster mass function
For our population synthesis model of the next sec-tion, we need the initial mass density of GCs in the Uni-verse ( ρ GC0 ). To find the relation between ρ GC0 and ρ GC from the previous section, we adopt a simple model forthe mass evolution of GCs. We assume that the initialGCMF is a Schechter-type function [48], i.e. a power-lawwith index − M c , as is found for young massiveclusters in the Local Universe [49]. We assume that thisinitial GCMF is universal throughout the Universe andacross cosmic time. There are arguments for a flatter ini-tial GCMF (i.e. fewer low-mass GCs) in dwarf galaxiesat high redshift and low metallicity [50], but we proceedwith the assumption of a Universal initial GCMF. Wewill discuss the effect of flatter initial GCMFs on theBHB merger rate in section V.To find an expression for the GCMF today, resultingfrom the initial GCMF we follow a similar approach as[27, 30] and assume that all GCs have lost an amountof mass ∆ = | ˙ M | t , where ˙ M is the mass loss rate fromescaping stars and t is the age of the GCs. We do notspecify the escape mechanism and let ∆ be constrainedby the Milky Way GCMF. Details of the various pro-cesses can be found in literature: relaxation driven evap-oration [51]; disc and bulge shocks [52, 53]; interactionswith molecular gas clouds (at young ages) [54–56] andcombinations of the various effects [25, 27, 57, 58]. Fromhere on we refer to the mechanism responsible for ˙ M , re-gardless of what the underlying physical process may be,as ‘evaporation’.We then assume that ∆ is a constant, i.e. independentof GC mass, host galaxy, orbit and formation epoch. Thisis clearly not realistic, because ˙ M depends on the (time-dependent) tidal field and the GC orbit within theirgalaxy [51]. However, this exercise is merely meant toarrive at an order of magnitude estimate of how muchmass GCs lose between formation and now, rather thandeveloping a realistic description of GC evolution. Thepresent-day GCMF, φ cl , defined as the number of GCsper unit volume ( n GC ) in the mass range [ M, M + d M ],is given by the ‘evolved Schechter function’ [30] φ cl = A ( M + ∆) − exp (cid:18) − M + ∆ M c (cid:19) . (3)At low masses, where the GCMF is affected by massloss ( M (cid:46) ∆ (cid:46) M c ), this function approaches a constant φ cl (cid:39) A/ ∆ . In fact, any initial GCMF evolves towardsa uniform φ cl at low masses if ˙ M is constant [59]. TheGCMF is often plotted as the number of GCs in logarith-mic mass bins ( ∝ d N/ d log M ), which increases linearlywith M at low masses and peaks at M peak (cid:39) ∆ (for∆ (cid:46) M c ). The simple functional form of equation (3)provides a good description for the Milky Way GCMFand the luminosity function of GCs in external galaxies[30]. The constant of proportionality A is found from the constraint that all GCs must add up to the present-dayGC mass density in the Universe: (cid:82) ∞ M lo φ cl M d M = ρ GC ,with ρ GC from equation (2) and M lo = 100 M (cid:12) .The GC evolution model we use in the next sectionalso considers mass loss by stellar evolution, which mostlyhappens in the first few 100 Myr. The fraction of massthat clusters lose as a result of stellar evolution dependson metallicity, the stellar initial mass function, stellarevolution details and on whether BHs are ejected, or not.For the cluster evolution model of the next section weneed the remaining mass in stars and white dwarfs ( M (cid:63) ).We use SSE to compute M (cid:63) at 11 Gyr for a Kroupa IMFin the range 0 . − M (cid:12) . We find that for metallici-ties of [0.01, 0.1, 1] Solar, the remaining mass fractionis M (cid:63) (11 Gyr) /M (cid:63) (0) (cid:39) [0 . , . , . M (cid:63) is in-dependent of metallicity and when excluding mass loss byBH ejections a cluster loses approximately half of its ini-tial mass by stellar evolution. We therefore assume thatclusters lose half their mass by stellar evolution alone.Next, we assume that stellar evolution and escape affectthe GCMF sequentially (i.e. first stellar mass loss andthen escape). We can then write M = 2( M + ∆) and wecan find the initial GCMF from the continuity equation[27] φ cl , ≡ d n GC d M = φ cl ( M ) (cid:12)(cid:12)(cid:12)(cid:12) ∂M∂M (cid:12)(cid:12)(cid:12)(cid:12) . (4)Because ∂M/∂M = 0 . φ cl ( M ) = A ( M / − exp [ − M / (2 M c )], the initial GCMF thatcorresponds to the present-day GCMF of equation (3) isgiven by φ cl , = 2 AM − exp (cid:18) − M M c (cid:19) . (5)We note that the Schechter mass of the initial GCMF is2 M c , where M c is derived from the present-day GCMF.We then introduce a factor K for the ratio ρ GC0 over ρ GC , i.e. K = ρ GC0 ρ GC = (cid:82) ∞ M lo φ cl , M d M (cid:82) ∞ M lo φ cl M d M . (6)In the next section we include the contribution fromclusters of all masses, and it is therefore important tounderstand the exact value of K , or better the distri-bution of K . To find a posterior distribution for K ,we fit the evolved Schechter functions from equation (3)to the Milky Way GCs and then derive K using equa-tions (5) and (6). We use the V -band luminosities ofthe 156 GCs in the 2010 edition of the Harris catalogue[41, 42] and then assume a constant mass-to-light ratioof Υ V = 2 M (cid:12) /L V, (cid:12) to convert luminosities to masses.A histogram of the resulting mass function is shown inFig. 1. The Milky Way values are binned in bins with 15GCs each, with the highest mass bin containing 6 GCs.The black dots are the average masses of the GCs in eachbin, while the horizontal error bars show the bin rangeand the vertical error bars show the Poisson errors.We then use the normalised evolved Schechter func-tion of equation (3) as a likelihood function to find∆ and M c . We use the Markov Chain Monte Carlo(MCMC) code emcee [60] to maximise the log-likelihoodand vary log ∆ and log M c , assuming flat priors in therange 3 ≤ log(∆ /M (cid:12) ) ≤ ≤ log( M c /M (cid:12) ) ≤ walker positionsof the converged chains we compute the GCMF (equa-tion 3) and the initial GCMF (equation 5) from ∆ and M c and at each mass we determine the [5% , , . ± . M c = 5 . ± .
1. Note that Jord´an et al. did not con-sider stellar evolution mass loss, so their initial GCMFwas truncated at M c , while ours is truncated at 2 M c .We also compute K with equation (6) for these 10 walker positions and find K = 32 . +86 . − . (90% credible).The spread in K provides an estimate of the uncertaintyin K , given the 156 Milky Way GC masses. The mergerrate will not increase by the same factor of K . This isfirstly because half of the value of K is due to stellarmass loss. The decrease in GC population mass by evap-oration is K ∆ = K/ . +43 . − . . Using the approxima-tion for the number of mergers in the observable redshiftrange from [34], N merge ∝ M . r − . , , we can estimatethe fractional increase in the merger rate ( K merge ) as afunction of K ∆ . The relation between K merge and K ∆ depends on the adopted mass-radius relation. If we pa-rameterise this as r h , ∝ M µ , then we find that for µ = 0(i.e. a constant initial radius) that then K merge (cid:39) . µ = 1 / K merge (cid:39) µ = 0 . K merge (cid:39)
8. The reason that K merge increases with µ , is because for large µ , the low-mass clusters are denser and produce more BHB mergers.Because we do not know the initial mass-radius relation,the value of K merge is thus in the range 2 . − .
6, corre-sponding the range in K ∆ mentioned above: 10 . − . K = 2 . K merge = K . This value for K merge is on the lower bound-ary of our estimated range of K ∆ for clusters with a con-stant radius, corresponding to a factor of ∼ K . For our upperboundary of K ∆ for µ = 0 . K merge is afactor of 6 . .
4) lower than our K ∆ ( K ). M / M )1234 l o g [ d N / d l o g ( M / M )] log( / M ) = 5.33 +0.240.18 log(2 M c / M ) = 6.26 +0.350.25 GCMFInitialGCMF
FIG. 1. GCMF of 156 Milky Way GCs from the Harris cat-alogue (black dots with error bars). The blue line shows anevolved Schechter function fit (equation 3). The resulting ini-tial GCMF, corrected for mass loss by stellar evolution (factorof 2) and evaporation (∆), is shown as a green dashed line.The shaded regions and uncertainties of quoted values indi-cate the 90% credible intervals. The inferred K value impliesthat the total mass of the Milky Way GC population was 32.5times higher initially. Half of this mass reduction is becauseof stellar mass loss and the remaining factor of 16 is due toevaporation. C. The GC formation rate
Next, we compute the cluster formation rate ˙ ρ GC ≡ ˙ ρ GC ( τ ), where τ is lookback time, for a set of modelassumptions. We do this in terms of a normalised GCformation rate R ≡ R ( τ ), such that (cid:82) ∞ R d τ = 1. Inthe next section we will derive R from a model for GCformation across cosmic time. For a given present-day ρ GC (equation 2), and no mass loss, the GC formation isthen found from ˙ ρ GC = ρ GC R . We now show how ˙ ρ GC can be easily derived for a population of GCs with a givenpresent-day mass function that have lost mass by stellarevolution and/or escape. Thus, after we determine φ cl , ,we can find ˙ ρ GC from imposing˙ ρ GC = Kρ GC R. (7)The cluster mass formed per unit volume integratedover all times is ρ GC0 = (cid:90) ∞ ˙ ρ GC d τ = Kρ GC , (8)= 2 . +2 . − . × M (cid:12) Gpc . (9)The large error bars are because of the uncertainty in K and (cid:104) η (cid:105) , and imply that ρ GC0 is uncertain by a factor of2. In the next section we include this uncertainty in thepredictions for the merger rate.
III. METHODOLOGY
The evolution of the BHBs in our cluster models iscomputed using the fast code cBHBd . While the detailsof this method are described in Antonini and Gieles [34],here we give a brief summary of the model philosophy, in-cluding the full set of differential equations that are usedto compute the secular evolution of the cluster modelsand the merging BHBs they produce. A. clusterBH We assume that the cluster consists of two types ofmembers: BHs and all the other members (i.e. otherstellar remnants and stars). Each contribute a total massof M BH and M (cid:63) , respectively, such that the total clustermass is M cl = M (cid:63) + M BH . We assume that after sev-eral relaxation time-scales the cluster reaches a state ofbalanced evolution [63, 64], so that the heat generatedin the core by the BHBs and the evolution of the clusterglobal properties are related as [59, 64, 65]:˙ E = ζ | E | t rh , (10)where E (cid:39) − . GM /r h is the total energy of the clus-ter, with M cl the total cluster mass and r h the half-massradius. The constant ζ (cid:39) . t rh is the averagerelaxation time-scale within r h which is given by [e.g., 66] t rh = 0 . (cid:114) M cl r G (cid:104) m all (cid:105) ψ ln κ . (11)Here (cid:104) m all (cid:105) is the mean mass of the stars and remnants(initially (cid:104) m all (cid:105) = 0 . M (cid:12) ), and ln κ is the Coulomblogarithm, which varies slowly with N , but we fix it toln κ = 10. The quantity ψ depends on the mass spectrumwithin r h , for which we adopt the following form: δ = 1 + a f BH , (12)where f BH = M BH /M cl is the fraction of mass in BHsto the total cluster mass and a is a constant that wasdetermined from a comparison to N -body models (seebelow). We define the start of the balanced evolution as t cc = N rh t rh , . (13)Under the above assumptions, the set of coupled or-dinary differential equations given below are integratedforward in time to obtain solutions for M BH ( t ), M cl ( t )and r h ( t ). The mass loss rate of BHs is coupled to the energygeneration rate, which itself is coupled to the total E and t rh of the cluster (equation 10), such that [64]˙ M BH = , t < t cc or M BH = 0 , − β M cl t rh , t ≥ t cc and M BH > . (14)The cluster mass loss due to stellar evolution is˙ M (cid:63), sev = , t < t sev , − ν M (cid:63) t , t ≥ t sev , (15)with t sev (cid:39) M (cid:63), ev = − ∆ (cid:104) t (cid:105) , (16)with (cid:104) t (cid:105) (cid:39)
10 Gyr the averaged cluster formation time.The total mass loss rate of the cluster is then˙ M = ˙ M (cid:63), sev + ˙ M (cid:63), ev + ˙ M BH . (17)In balanced evolution, the expansion rate of the clusterradius as the result of relaxation is˙ r h , rlx = ζ r h t rh + 2 ˙ MM cl r h . (18)Before balanced evolution, the cluster radius expandsadiabatically as the result of stellar mass loss at a rate˙ r h , sev = − ˙ M (cid:63), sev M cl r h . (19)The final expression for the half-mass radius evolution is˙ r h = (cid:40) ˙ r h , sev , t < t cc , ˙ r h , sev + ˙ r h , rlx , t ≥ t cc . (20)The remaining parameters were obtained in Antoniniand Gieles [34] by fitting the results of N -body simula-tions: N rh = 3 . β = 2 . × − , ν = 8 . × − and a = 1 . × . B. BHBdynamics
The initial contraction of the cluster core due to two-body relaxation leads to high central densities of BHswhich favor the formation of binaries through three-bodyprocesses. The energy produced by the BHBs reverts thecontraction process of the core and powers the subsequentexpansion of the cluster as described by equation (18).We can then relate the BHB hardening rate to the rateof energy generation ˙ E bin = − ˙ E , (21)where ˙ E bin is the hardening rate of all core binaries.Equation (21) allows us to couple in a simple way theevolution of the BHBs to the evolution of the clustermodel. It is important to stress that the hardening rateequation (21) depends neither on the number of binariespresent in the cluster core nor on the exact mechanismleading to their formation.In order to compute a merger rate and the binaryproperties from equation (21), we need to further spec-ify the dynamical processes that lead to the hardeningand merger of the binaries. We are interested in merg-ers that occur through (strong) binary-single dynamicalencounters in the cluster core [e.g., 67, 68]. Thus weconsider: (i) mergers that occur in between binary-singleencounters while the binary is still bound to its parentcluster (in-cluster inspirals); (ii) mergers that occur dur-ing a binary-single (resonant) encounter as two BHs aredriven to a short separation such that gravitational wave(GW) radiation will lead to their merger (GW captures);and (iii) mergers that occur after the binary is ejectedfrom its parent cluster. We use BHBdynamics to de-termine the rate and masses of the BH binary mergersproduced by these three dynamical channels.In balanced evolution, after a binary is ejected ormerges a new binary must quickly form to meet the en-ergy demand from the cluster. Under such conditions,the binary formation rate nearly equals the binary ejec-tion rate and it is given therefore by the BH mass ejectionrate equation (14) divided by the total mass ejected byeach binary Γ bin (cid:39) − ˙ M BH m ej , (22)where m ej was computed using equation (38) in [34], andit is approximately a fixed number m ej (cid:39) m . The num-ber of BHBs that merge before a time t from the forma-tion of the cluster is [34] N ( < t ) = (cid:90) t Γ bin [ P in + P ej ( t − t (cid:48) )] d t (cid:48) , (23)where P ej ( t − t ej ) is the probability that a binary ejecteddynamically from the cluster at a time t ej will mergedue to GW emission within a time t from the forma-tion of the cluster, and P in is the probability that abinary merges inside the cluster. Specifically, P in isthe sum of the probability that that a binary mergesthrough an in-cluster inspiral and the probability of aGW capture. As described in [34], the probability thata binary merges in between two binary-single interac-tions is given by integrating the differential merger prob-ability per binary-single encounter over the total numberof binary-single interactions experienced by the binary.Similarly, the probability that a binary merges througha GW capture is obtained by dividing each binary-singleencounter into 20 intermediate resonant states as in [68],and by integrating the differential merger probability per resonant encounter over all encounters experienced by abinary. The merger probabilities are computed by as-suming that the eccentricity of the BHBs follows thatof a so-called thermal distribution N ( e ) ∝ e [e.g., 69],and their full expressions can be found in Antonini andGieles [34]. Moreover, when evaluating P ej and P in weset E bin = − Gm m / a , with a the binary semi-majoraxis and m and m the mass of the BH components,i.e., we have assumed that only one BHB is responsiblefor all the heating at any given time. However, becausethe dependence is weak, e.g. P in ∝ ˙ E − / [34], and thenumber of hard binaries is expected to be of order unityin the type of clusters we consider [70], this simplificationis reasonable. C. Black hole mass function and natal kicks
In order to calculate the merger rate through equa-tion (23) we need a physically motivated model for theBH mass function and its time evolution.We sample 100 stellar progenitor masses from the ini-tial mass-function φ (cid:63) ( m (cid:63) ) ∝ m − . (cid:63) [71] between m (cid:63), lo and 100 M (cid:12) , with m (cid:63), lo (cid:39) M (cid:12) the stellar mass abovewhich a BH forms. The resulting masses of the BHs arethen obtained using the fast stellar evolution code SSE[72] which we modified to include up to date prescrip-tions for stellar wind driven mass loss [73], compact-object formation and supernova kicks [74, 75], and wealso include prescriptions to account for pulsational-pairinstabilities and pair-instability supernovae [76]. The ini-tial mass fraction in BHs is set equal to the total mass inBHs divided by the total mass in stars between 0 . m (cid:63), lo for a Kroupa [71] initial mass function, and rangesfrom f bh (cid:39) .
04 to 0 .
07 depending on the metallicityand the adopted prescription for compact-object forma-tion. These mass fractions first increase as the result ofstellar evolution mass loss and then they reduce due tothe ejection of BHs caused by natal kicks and dynam-ical ejections. The BH natal kicks are computed usinga standard fallback model in which the BHs receive akick drawn from a Maxwellian distribution with disper-sion σ = 265 km s − [77], lowered by the fraction of theejected supernova mass that falls back into the compact-object. The fallback fraction and remnant masses aredetermined according to the chosen remnant-mass pre-scription. We adopt here the rapid supernova prescrip-tion described in Section 4 of Fryer et al. [75], in whichthe explosion is assumed to occur within the first 250msafter bounce. But later in Section IV B we also exploreother choices for the compact-object formation recipe.The cluster dynamically processes its BH populationsuch that the mass of the merging BHBs progressivelydecreases with time because the most massive BHs arethe first to reach the cluster core, form hard binariesand merge [e.g., 78]. Assuming that the merger prod-ucts of BHB mergers are ejected, simulations of densestar clusters also show that the merging BHBs have adistribution of mass-ratio that is strongly peaked aroundone [e.g., 20]. Thus, we assume that the BHs taking partin the dynamical interactions have the same mass, andthat at any given time this mass is equal to that of thelargest BH in the cluster, i.e., m = m = m = m max .The value of m max at a given time, can be easily linkedto the time evolution of the total mass in BHs given by clusterBH . For a generic BH mass function φ • , we usethe fact that (cid:90) m max m lo φ • m d m = M BH , (24)where the integral in the left hand side is simply thecumulative distribution of BH masses computed withSSE. We then invert numerically this relation to find m max ( M BH ). D. Cluster formation and initial properties
We sample the initial cluster masses from the Schechtermass function equation (5), i.e., we assume that bothevaporation and stellar mass loss are important. The val-ues [∆, M c ] needed to compute φ cl , and K are sampledfrom their posterior distributions derived in Section II B.An important property of balanced evolution is thatthe value of a cluster half-mass radius today is largelyindependent of its initial value [33, 59]. It is there-fore not possible to infer the initial density of GCsfrom their properties today. We therefore consider threechoices for the initial stellar mass density within the half-mass radius, ρ = 3 M / (8 πr , ), which we set equalto ρ = 10 M (cid:12) pc − for our canonical Mod1, and in-crease(decrease) by a factor 10 in Mod2(Mod3) to explorethe effect of initial cluster density.For the cluster metallicity, we fit a quadratic polyno-mial to the observed age-metallicity relation for the MilkyWay GCs [79], to obtain the mean metallicitylog( Z mean /Z (cid:12) ) (cid:39) .
42 + 0 . (cid:18) t Gyr (cid:19) − . (cid:18) t Gyr (cid:19) . (25)Given the cluster age, t , we then assume a log-normal dis-tribution of metallicity around the mean, with standarddeviation σ = 0 . et al. [80]. The same model has also been used inrecent work [e.g., 31, 81], allowing a direct comparisonof our results to literature. El-Badry et al. [80] describethe process of GC formation as resulting from star forma-tion activity in the high-density disks of gas-rich galax-ies. Motivated by the results of simulations of molecular cloud collapse, they assume that massive bound clustersform preferentially when the gas surface density exceedsa certain threshold. Applying this Ansatz to a semi-analytic gas model built on dark matter merger trees,they make predictions for the cosmological formation rateof GCs. The resulting cluster formation history peaksat a redshift of ∼
4, which is earlier than the peak in thecosmic star formation history (redshift ∼
2, [82]). Wesample the formation redshift of our cluster models fromthe total cosmic cluster formation rate given by the fidu-cial model of El-Badry et al. [80] and integrated over allhalo masses. This corresponds to the formation rate percomoving volume of their Figure 8 with their parameters β Γ = 1 and β η = 1 /
3, where β Γ sets the dependence ofthe cluster formation efficiency on surface density, and β η the dependence of the star formation rate on the halovirial mass. We then normalize the GC formation tounit total number to obtain R (equation 7). Thus, weonly sample the cluster formation redshift from the El-Badry et al. [80] model, while the total cluster formationrate is given by our equation (7). For this model, ap-proximately 25% of the cosmic star formation [82] is instar clusters at redshifts (cid:38) K = 32 .
5, implying that K = 130 is an upper limit toensure that the cluster formation rate is below the starformation rate. This limit corresponds to 2 . σ in our K distribution and hence it is unlikely to occur. Later,in order to determine the importance of our assumptionabout the cluster formation history, we will consider an-other class of models with different values for β Γ and β η .Given the initial cluster mass, radius, metallicity, andformation time, the BHB merger rate, ˙ N , is obtainedfrom cBHBd . IV. BINARY BLACK HOLE MERGER RATE
The merger rate density of BHBs at a lookback time τ is R ( τ ) = (cid:90) (cid:90) (cid:90) φ cl , ( M ) R ( τ (cid:48) ) s ( Z )˙ N ( τ (cid:48) − τ ; M , r h , , Z )d τ (cid:48) d M d Z (26)where ˙ N ( t ; M , r h , , Z ) is the BHB merger rate corre-sponding to a cluster with an initial mass M , half-massradius r h , and metallicity Z at a time t after its forma-tion; R ( τ ) is the normalized cluster formation rate and s ( Z ) is the normalized formation rate of clusters with ametallicity Z at a time τ , (cid:82) s ( Z ; τ )d Z = 1, which can becalculated given a model for the time evolution of metal-licity, e.g., equation (25).In practice, for each model assumption in Table I wesample 100 values over the posterior distributions of theparameters M c and ∆ obtained from the MCMC fit tothe observed Milky Way GCMF. Then, for each pair [ M c ,∆], we evolve N cl = 600 models with masses sampled overa grid of constant logarithmic step size, δ log M/M (cid:12) = Model Density Z Cluster formation BH masses Natal ˙ M (cid:63), ev R κ R , in − cluster M (cid:12) pc − kicks [Gpc − yr − ] [Gpc − yr − ]Mod1 10 Eq. (25) β Γ = 1; β η = 1 / . +12 . − . . +0 . − . . +1 . − . Mod2 10 - - - - - 12 . +32 . − . . +0 . − . . +1 . − . Mod3 10 - - - - - 3 . +5 . − . . +0 . − . . +1 . − . Mod4 10 (cid:12) - - - - 6 . +13 . − . . +0 . − . . +1 . − . Mod5 - 0.1Z (cid:12) - - - - 7 . +16 . − . . +0 . − . . +2 . − . Mod6 - Z (cid:12) - - - - 6 . +15 . − . . +0 . − . . +2 . − . Mod7 - Eq. (25) All form at z=3 - - - 7 . +15 . − . . +0 . − . . +1 . − . Mod8 - - β Γ = 0; β η = 1 / . +16 . − . . +0 . − . . +2 . − . Mod9 - - β Γ = 1; β η = 1 / . +13 . − . . +0 . − . . +1 . − . Mod10 - - β Γ = 1; β η = 1 / . +13 . − . . +0 . − . . +1 . − . Mod11 - - - [86] - - 7 . +16 . − . . +0 . − . . +1 . − . Mod12 - - - Delayed - - 7 . +14 . − . . +0 . − . . +1 . − . Mod13 - - - Rapid No kicks - 10 . +21 . − . . +0 . − . . +2 . − . Mod14 - - - - momem. - 8 . +14 . − . . +0 . − . . +2 . − . Mod15 - - - - fallback 0 2 . +1 . − . . +0 . − . . +0 . − . Mod16 10 (cid:12) - - - Eq. (16) 12 . +29 . − . . +0 . − . . +1 . − . Mod17 - 0.1Z (cid:12) - - - - 12 . +32 . − . . +0 . − . . +1 . − . Mod18 - Z (cid:12) - - - - 10 . +24 . − . . +0 . − . . +2 . − . Mod19 - Eq. (25) All form at z=3 - - - 11 . +27 . − . . +0 . − . . +1 . − . Mod20 - - β Γ = 0; β η = 1 / . +29 . − . . +0 . − . . +2 . − . Mod21 - - β Γ = 1; β η = 1 / . +28 . − . . +0 . − . . +1 . − . Mod22 - - β Γ = 1; β η = 1 / . +34 . − . . +0 . − . . +1 . − . Mod23 - - - [86] - - 14 . +36 . − . . +0 . − . . +1 . − . Mod24 - - - Delayed - - 11 . +27 . − . . +0 . − . . +1 . − . Mod25 - - - Rapid No kicks - 18 . +36 . − . . +0 . − . . +2 . − . Mod26 - - - - momem. - 14 . +32 . − . . +0 . − . . +2 . − . Mod27 - - - - fallback 0 3 . +2 . − . . +0 . − . . +0 . − . Mod28 10 (cid:12) - - - Eq. (16) 2 . +4 . − . . +0 . − . . +0 . − . Mod29 - 0.1Z (cid:12) - - - - 3 . +5 . − . . +0 . − . . +1 . − . Mod30 - Z (cid:12) - - - - 3 . +7 . − . . +0 . − . . +2 . − . Mod31 - Eq. (25) All form at z=3 - - - 2 . +5 . − . . +0 . − . . +1 . − . Mod32 - - β Γ = 0; β η = 1 / . +5 . − . . +0 . − . . +1 . − . Mod33 - - β Γ = 1; β η = 1 / . +4 . − . . +0 . − . . +1 . − . Mod34 - - β Γ = 1; β η = 1 / . +4 . − . . +0 . − . . +0 . − . Mod35 - - - [86] - - 2 . +4 . − . . +0 . − . . +0 . − . Mod36 - - - Delayed - - 2 . +5 . − . . +0 . − . . +1 . − . Mod37 - - - - No kicks - 3 . +6 . − . . +0 . − . . +1 . − . Mod38 - - - - momem. - 2 . +3 . − . . +0 . − . . +1 . − . Mod39 - - - - fallback 0 0 . +0 . − . . +0 . − . . +0 . − . TABLE I. Model parameters used in the calculations. Here β Γ and β η refer to the parameters of the cosmological models inEl-Badry et al. [80] that are used to sample the cluster ages in our simulations. The rightmost three columns give the localmerger rate density of BHBs, the rate evolution parameter κ , and the local merger rate of in-cluster mergers (including GWcaptures). .
01, in the range 10 − M (cid:12) , and use that: R ( z ) (cid:39) Kρ GC N cl (cid:80) i =1 ˙ N ( z ; M ,i ) φ cl , ( M ,i ) M ,iN cl (cid:80) i =1 φ cl , ( M ,i ) M , i , (27) where the formation time of each cluster is randomlysampled from the corresponding R ( τ ) distribution; wethen use equation (25) to compute the mean metallic-ity, Z mean , that corresponds to that formation time, andthus find the cluster metallicity by drawing from a log-normal distribution around Z mean . The half-mass radiusof the cluster is obtained from the cluster mass given FIG. 2. Median of the merger rate distribution (solid lines)for different initial cluster half-mass densities correspondingto Mod1, Mod2 and Mod3 in Table I. The dashed lines con-tain 90% of all model realizations (between the 5 and 95 per-centiles). The upper panel gives the merger rate as a functionof redshift. The middle and lower panels show the differentiallocal merger rate as a function of the the primary BH massand initial cluster mass, respectively. We compare our resultsto the median (solid), and the 90% credible intervals (hatchedregions) inferred from the GWTC-2 catalogue in [12]. In themiddle panel, we have used their ‘Power Law & Peak’ modeland the thick-red line gives the BH initial mass function (in ar-bitrary units). In the lower panel we show the initial GMF (inarbitrary units) four our best fit value log(2 M c /M (cid:12) ) = 6 . M (cid:12) BHBs [83]. The observable horizonsof the Einstein telescope [84] and the Cosmic Explorer [83]extend to the very early Universe and are both to the rightof the x -axis range in the figure. the assumed half-mass density. Note that because eachcluster has its own metallicity, each time we generate anew BH population using SSE; the fraction of clusterswith Z < . Z (cid:12) is (cid:39) ρ GC . We assume that the parameter ρ GC followsa Gaussian distribution with mean 7 . × M (cid:12) Gpc − and σ = 2 . × M (cid:12) Gpc − . We sample 1000 val- FIG. 3. Distribution of the rate parameter κ , and the localmerger rate, R , for each of the three models of Fig. 2. Colorsare as in Fig.2. ues from this latter distribution and for each of them weuse equation (27) to determine a merger rate estimate foreach of the [ M c , ∆] values, and thus obtain a distribution of merger rate density values.Because our results turn out to be more sensitive tothe cluster initial density than to other parameters, wefirst focus on Mod1, Mod2 and Mod3 in Table I. In thesemodels we vary ρ in a range that is relevant to real glob-ular clusters, while keeping fixed all the other parametersas given in the table. This allows us to bracket a plausi-ble range of values for the local merger rate density andits redshift evolution. In Fig. 2 we plot the median of themerger rate distribution and credible intervals as a func-tion of redshift, and the primary BH mass distributionof binaries merging at redshifts z < ∼
10 for any density assumption. Thisis due to the fact that
R ∝ K at a very good approx-imation, and, as we discussed above, tight constraintson K cannot be placed from the present day Milky WayGCMF. Moreover, rates are not too sensitive to the ini-tial cluster density – two orders of magnitude differencein the initial density leads to a factor of ∼ R . From this we conclude that the ini-tial density uncertainty is as important as the unknowninitial GCMF.For each initial GCMF, corresponding to new valuesof [ M c , ∆] and (cid:104) η (cid:105) , we fit the redshift distribution of themerger rate density at z < R ( z ) = R (1 + z ) κ , (28)to derive a distribution of values for the parameters R and κ for each of our three density assumptions. In thisanalysis we neglect the uncertainties associated to eachfit because their standard deviations are much smallerthan the variation in the inferred parameters across thedifferent models. The initial parameters for each of the0three densities and the corresponding median values anduncertainties (5 and 95 percentiles) of R and κ aregiven in Table I (Mod1, Mod2 and Mod3). The distri-butions obtained from this analysis are shown in Fig. 3.The local BH merger rate density from GCs varies inthe range R (cid:39) − yr − . A comparison tothe local merger rate inferred from the GW detections,23 . +14 . − . Gpc − yr − (for their redshift independent re-sults and 19 . +16 . − . Gpc − yr − when the merger rate isallowed to evolve with redshift) [12], shows that BHBsformed dynamically in GCs are likely to explain a sig-nificant fraction of the BHB mergers detected by LIGO-Virgo. Note, however, that if GCs are formed with highdensities, ∼ M (cid:12) pc − , then our merger rate estimatesare consistent with the LIGO-Virgo merger rates withinuncertainties. We note that although this high densityis preferred to explain the overall rates, the rates in themass range 13 − M (cid:12) are somewhat better reproducedby Mod1 (10 M (cid:12) pc − , see Fig. 2). Combining thethree models together, i.e., assuming a universe in whichone third of the clusters form as in Mod1, one third as inMod2 and the remaining as in Mod3, we find R = 7 . +21 . − . Gpc − yr − ; κ = 1 . +0 . − . , (29)where uncertainties refer to the 90% credible intervals.The middle panel of Fig. 2 shows the primary BH massdistribution normalized to the BH merger rate densityin the local Universe ( z < − M (cid:12) . The BH mass distributionappears to be sensitive to the initial cluster density, inthe sense that higher densities lead to a higher fractionof lower mass BHs among the merging population. Itis interesting that the higher density models, Mod1 andMod2, provide a better match to the inferred BH massfunction and the rates. The additional low-mass BHs inhigh density GCs is due to the higher retention fractionof lower mass BHs in higher density models after a natalkick due to the high escape velocities from these clus-ters, and to the fact that denser clusters process theirBH populations faster, thereby ‘eating’ away their BHmass function more.Given the number of complex features that can be seenin the BH mass distributions we do not attempt here aparametrization over the full range of BH masses. More-over, as we will show later in Section IV B, these distri-butions are sensitive to the uncertain prescriptions forBH formation, natal kicks and metallicity. We insteadconsider the mass range 13 M (cid:12) to 30 M (cid:12) where a sim-ple power law model, dR/dm ∝ m α , does a reasonablejob. In this mass range we find α = 0 . +0 . − . (Mod3), α = − . +0 . − . (Mod1), and α = − . +0 . − . (Mod2), wherethe reported values are the median of the distributionsobtained by fitting each of the 100 BH mass distributionscorresponding to different [ M c , ∆], and the uncertaintiesrefer to the 5 and 95 percentiles. As before we neglect theuncertainties associated to each fit because their standard errors are much smaller than the variation of α across thedifferent models. By combining the three density modelstogether, we find α = − . +1 . − . . (30)Negative values of α are preferred, though positive val-ues are also acceptable. The value of the power lawindex found by us is broadly consistent with the valueof α = − . +0 . − . reported by [12] for their low-mass( < M (cid:12) ) slope of the “Broken Power Law” model. Forcomparison, the BH initial mass function integrated overall metallicities is shown in the middle panel of Fig. 2.Within the same BH mass range, m = 13 M (cid:12) to 30 M (cid:12) ,the best fit power law model to the initial mass func-tion had a spectral index α ≈ − .
8. The most strikingfeature, however, is that all mass distributions in Fig. 2are strongly depleted at m (cid:46) M (cid:12) . The fraction ofBHs below this mass decreased by more than one orderof magnitude with respect to the BH initial mass func-tion. Due to this, all our models underpredict the num-ber of BHBs at m ∼ M (cid:12) compared to the mass distri-bution inferred from the LIGO-Virgo detections. More-over, we note some features that are common to the threemodels considered here. All three models show peaks at m (cid:39) M (cid:12) , 20 M (cid:12) and 30 M (cid:12) . Above m = 30 M (cid:12) theBH mass distribution starts to decline rapidly with massuntil a break at m (cid:39) M (cid:12) above which the declinebecomes much steeper. All models show essentially noBHs with mass above 40 M (cid:12) or below 5 M (cid:12) . The lowmerger rate value at (cid:38) M (cid:12) is a consequence of thestellar mass loss prior to the formation of the BHs be-cause a down-turn above 30 M (cid:12) is also seen in the BHIMF and we find it even in models that do not include anyprescription for pair instability [87]. Our pulsational-pairinstabilities and pair-instability supernovae prescriptionsare taken from [76], and for the maximum initial stellarmass we considered, 100 M (cid:12) , they have little or no effecton the resulting BH masses. Note also that we do notconsider hierarchical mergers [88]. Their contribution tothe merger rate is sensitive to the distribution of BH na-tal spins, which is poorly constrained. Assuming thatBHs are formed with no spin, Rodriguez et al. [89] findsthat ∼
10% of BHB mergers come from previous mergers;when the BH dimensionless spin parameter is increasedabove 0 .
1, the contribution drops to only a few per centor less. In addition, in the discussion (Section V C) weshow that with our adopted mass-radius relation, less 2ndgeneration mergers are expected compared to Rodriguez et al. [90]. Thus, including hierarchical mergers is notexpected to significantly change our integrated mergerrate estimates. However, the high mass BHs resultingfrom multiple mergers can partly fill up the mass distri-butions above m ∼ M (cid:12) where the merger rates fromour models are nearly zero.The bottom panel of Fig. 2 shows the differential con-tribution of clusters with different masses to the localmerger rate ( z < M (cid:12) , above which the rate1starts to rapidly decrease because of the exponentialtruncation of the initial GCMF above M c . The con-tribution of clusters with masses larger than 10 M (cid:12) orsmaller than 10 M (cid:12) is negligible. In our models, clustersthat have a mass M < × M (cid:12) are fully disrupted bythe present time. These clusters have been neglected insome previous work [e.g., 18, 20, 23, 78], but we find thatthey contribute a significant fraction of the local mergerrate: ≈ .
33, 0 .
48 and 0 .
30 for Mod1, Mod2 and Mod3respectively.
A. In-cluster vs. ejected binaries
The merger of a BHB in our models can occur eitherafter the binary has been ejected dynamically from thecluster, or within the cluster itself. In-cluster mergersare relevant because they can lead to the formation ofBHs with mass above the pair-instability gap at ≈ M (cid:12) [88, 90, 91], and the observational implications of thishave been discussed in a number of previous papers, e.g.[92–94]. Among all in-cluster mergers, GW captures arealso particularly important because a fraction of them areexpected to have a finite eccentricity at the moment theyfirst chirp within the the LIGO frequency band above 10Hz [68, 90, 95–97]. Thus, they could be identified amongother binaries due to their unique eccentric signature.In Fig. 4 we show separately the rate evolution ofBH mergers occurring among the ejected binaries, thoseforming inside the cluster and GW captures, as well asthe mass distributions of mergers at z < z (cid:38) (cid:39) − yr − and that of GW cap-tures is (cid:39) . − yr − , with little dependence on theinitial density assumed. However, the fractional contri-bution of in-cluster mergers does depend quite stronglyon the cluster initial conditions, in the sense that higherdensities lead to lower fractions — for ρ = 10 M (cid:12) pc − ,nearly 40% of the local mergers are formed in-cluster,while for the highest initial densities, ρ = 10 M (cid:12) pc − ,in-cluster mergers only contribute (cid:39)
15% of the to-tal. The percentage of all in-cluster mergers that occurthrough a GW capture is ≈
20% in the local Universeand increases smoothly with redshift, reaching ≈ (cid:38) .
1) above 10 Hzfrequency, we conclude that eccentric mergers from glob-ular clusters contribute (cid:46) . − yr − to the mergerrate in the local universe. This low rate is consistent withthe non detection of eccentric binaries in current searches[98].Our models Mod1 and Mod2 show a decrement forthe local fraction of in-cluster mergers over previous es-timates which bracketed this between 30% and 50% of FIG. 4. As Fig.2 but we now show separately the median ofthe merger rate distribution of all in-cluster mergers (dashedlines), only GW captures (dot-dashed lines), and mergersamong the ejected binaries (solid lines). the total rate [e.g., 31]. This difference arises from thefact that this previous work only considered clusters withmass (cid:38) M (cid:12) for which about half of the overall mergerrate is due to in-cluster binaries. However, as shown inFig. 4, lower mass clusters contribute significantly to thelocal rate, although the mergers they produce only oc-cur among the ejected population. The reason for thisis that their BHs have all been ejected by z = 1. Thus,including these systems leads to an overall reduction ofthe contribution of in-cluster mergers and also affects theredshift evolution of the merger rate density.Fig. 5 shows the distributions of R and κ obtained byfitting equation (28) to the merger rate evolution of in-cluster and ejected binaries separately. The merger rateof in-cluster binaries evolves steeply with redshift, κ ≈ κ ≈
1. If, as before, we as-sign equal probability to each of our density assumptionsand fit the total merger rate density of in-cluster binariesusing equation (28) we find R , in = 1 . +1 . − . Gpc − yr − ; κ in = 2 . +1 . − . , (31)2while for ejected binaries R , ej = 5 . +21 . − . Gpc − yr − ; κ ej = 1 . +0 . − . . (32)We now use a simplified analytical model to gain somephysical insights on this result.From equation (23), the merger rate for in-cluster bi-naries is d N d t (cid:12)(cid:12)(cid:12) in = Γ bin P in ( t ) , (33)while for ejected binaries we haved N d t (cid:12)(cid:12)(cid:12) ej = dd t (cid:90) t Γ bin P ej ( t − t (cid:48) )d t (cid:48) . (34)The merger probabilities that enter in the integral equa-tions above can be linked to the evolution of the clus-ter properties in a simple way under some simplifyingassumptions. If we neglect cluster mass loss and thatthe BHs have a range of masses – both have little ef-fect on the merger rate evolution (see Secion IV B) –we can write t rh ( t ) = t rh , (cid:0) ζt/t rh , (cid:1) , and ρ ( t ) = ρ (cid:0) ζt/t rh , (cid:1) − [91, 99]. If r h ∝ t / and M =constant, then ρ ∝ t − , v esc ∝ t − / and there-fore v esc ∝ ρ / . Then from [34] we know P in ∝ v / ,hence P in ∝ ρ / (neglecting captures) and P ej ∝ ( t − t ej ) / ρ ( t ) / [34]. We can then determine theredshift dependence of the merger rate through equa-tions (33) and (34).At times t (cid:29) t rh , /ζ , for in-cluster binaries we haved N d t (cid:12)(cid:12)(cid:12) in ∝ (1 + z ) . , (35)where we used that t ( z ) ∝ (1 + z ) − / in order to converttime into redshift. For ejected binaries we have d N d t (cid:12)(cid:12)(cid:12) ej ∝ t − / , or d N d t (cid:12)(cid:12)(cid:12) ej ∝ (1 + z ) . . (36)Although we have neglected some important ingredients(e.g., mass loss, BH mass function), the expected valueof κ for the two populations is consistent with the onesfound above and, as expected, it is much steeper for in-cluster mergers. This fits in the view that the rate atwhich the merging BHBs are produced by a cluster iscontrolled by the relaxation process within the clusteritself, providing a physical interpretation to our results.Moreover it implies that most of the merging BHBs atproduced by clusters that are still in the expansion phase.Another result of our analysis is that the local rate ofin-cluster inspirals and GW captures are nearly indepen-dent of the initial density assumed. In general, we findthat other model variations also have little effect on themerger rate of in-cluster binaries. This result is becausecluster evolution during balanced evolution, i.e. at late FIG. 5. Distribution of the rate parameter κ , and the localmerger rate, R , for each of the three models of Fig. 2, and forthe in-cluster mergers and ejected binary mergers separately.Colors are as in Fig.2.FIG. 6. Example of the BHB merger rate evolution for twocluster masses, and separately for in-cluster and ejected bina-ries. The middle panels gives the cluster relaxation time, andthe lower panels the total BH mass in units of the initial value.Different lines correspond to different initial half-mass density, ρ = 0 . , . , . × M (cid:12) pc − . The initial densityincreases with line thickness. Here we set ∆ = 3 × M (cid:12) andthe simulations are terminated either after 13 Gyr of evolutionor after all BHs have been ejected. FIG. 7. Median of the merger rate distribution (solid lines) and 90% confidence intervals (dashed lines) for the models Mod4to Mod15 in Table I where the initial half-mass density is set to 10 M (cid:12) pc − . Middle panels give the distribution of primaryBH masses for mergers at z <
1. The lower panels show the mass distribution of clusters where these merging binaries wereformed. stages when in-cluster mergers that we can observe occur,is insensitive to the initial conditions.Due to the expansion powered by the BHBs, all clus-ters evolve asymptotically to (approximately) approachthe same value of half-mass relaxation time. Hence,after some time, the merger rate of in-cluster binariesmust also become approximately the same for all clus-ters. This concept is illustrated in Fig. 6 where we showthe evolution of a set of cluster models with the samemass but different ρ . For M = 3 × M (cid:12) (left pan-els), the BHB merger rate at (cid:38)
10 Gyr only varies by afactor of ∼ t/t rh increases until t rh ∼ t rh , /ζ , after which t rh ∝ t [33].In the right panel of Fig. 6 we consider the evolution ofclusters with an initial lower mass M = 3 × M (cid:12) . Inthese models all BHs are ejected at t (cid:46) M (cid:38) M (cid:12) (e.g., lower panel of Fig. 4). Moreover, the binary mergerrate near the end of the simulations shows a larger varia-tion among different models than in the high mass clus-ter case. This simply reflects the large difference in thecluster density at early times when these binaries wereformed and ejected. B. Dependence on model parameters
In the previous section we consider the merger ratedensity evolution for three choices of initial cluster den-sity. Here we discuss the results for a larger set of mod-els in which for each of the three density assumptionswe vary the prescription for the cluster metallicity, thecosmological cluster formation model, the BH formationmechanism, the BH natal kicks, and the cluster mass lossrate. All the model parameters we considered are listed4
FIG. 8. Same as Fig. 7 but for ρ = 10 M (cid:12) pc − , i.e., Mod16 to Mod27 in Table I. in Table I (Mod4 to Mod39), together with the corre-sponding values of merger rate evolution parameters anduncertainties. We stress that not all the models anal-ysed here are realistic representations of a globular clus-ter population. They are nevertheless useful in order tounderstand the impact of different model parameters onthe merger rate evolution and BH mass distribution. Ourmain message here is that variations in other model as-sumptions have little effect on the local value of the BHBmerger rate density and its redshift evolution. Thus, weconclude this section by presenting the results from anadditional set of models where ρ is varied over a widerrange of values than in Section IV to more systematicallyexplore its effect on the BHB merger rate. Metallicity . In order to explore the dependence of themerger rate on metallicity, we consider models where theclusters all have the same metallicity which we set to Z = 0 .
01, 0.1 or 1 × Z (cid:12) . Since mass loss due to stellarwinds is less effective in metal-poor stars, the formingmerger remnant mass increases with decreasing metal-licity. At Solar metallicity, the mass distribution of thefinal merger products spans from a few solar masses up toabout 30 M (cid:12) and peaks near 10 M (cid:12) . At lower metallic-ity, Z = 0 .
01 and 0.1, the distribution of remnant masses is much wider with its maximum at ∼ M (cid:12) . This hasan obvious effect on the mass distribution of the merg-ing BHBs as can be seen in Fig. 7, 8 and 9. The valueof metallicity affects also what type of clusters make theBH mergers in the local universe, with their mass dis-tribution being skewed towards higher values for Solarmetallicities. The important result here, however, is thatthe evolution of the merger rate density is largely un-affected by the choice of metallicity and its dependenceon cluster age. Even in the unrealistic case in which allclusters are formed at Solar metallicity, the merger ratedensity only starts to deviate significantly from the othermodels at z >
2. Such lower merger rate at early timesis expected and it is a consequence of the longer t rh dueto the lower initial BH mass fraction.We conclude that a detailed knowledge of the metal-licity distribution of GCs and its dependence on time isnot necessary in order to determine a BHB merger rate,although it has an important effect on their mass distri-bution. Cluster ages.
We implemented two additional choicesfor the parameters in the cosmological model of El-Badry et al. [80] which determine the distribution of clusterages: [ β Γ = 1 , β η = 1 / β Γ = 0 , β η = 1 / FIG. 9. Same as Fig. 7 but for ρ = 10 M (cid:12) pc − , i.e., Mod28 to Mod39 in Table I. two models are shown in Fig. 8 of El-Badry et al. [80].Moreover, we consider an additional case of a burst-likecluster formation history in which all clusters are formedat z = 3.From Fig. 7, 8, and 9 we can see that, for a given initialdensity, our results at z (cid:46) z = 3 leads to a merger rate density and BH massdistribution that are consistent with those obtained fromthe full cosmological models. At z >
2, however, the red-shift evolution of the merger rate is clearly affected withits peak coinciding with the peak of cluster formationactivity in each model.
BH formation.
We consider three more recipes to com-puting the BH mass distribution based on different core-collapse/supernova models. We use the delayed model inwhich the supernova explosion is allowed to occur over amuch longer timescale than in the previously employedrapid model [75]. We then use the compact-object massprescriptions from [85] and [86]. These two latter modelsuse slightly different recipes for the proto-compact objectmasses while adopting the same formulae to determine the amount of fallback material. We note that the effectof the BH formation recipe is two folds as it influencesboth the mass distribution of the BHs as well as theirnatal kicks. Apart from the effect on the BH mass func-tion, however, there is very little change of the mergerrate evolution among the various prescriptions, with thedelayed model leading to a slightly lower merger rate atall redshifts than the others.
Natal kicks.
Two additional assumptions about theBH natal kicks are explored. In one the BHs are formedwith no kick, and in the other the BHs receive the samemomentum kick as neutron stars, meaning that their kickvelocities are drawn from a Maxwellian distribution withdispersion σ = 265 km s − [77] and then reduced by theneutron star to BH mass ratio, 1 . M (cid:12) /m . Among themodel variations considered in this section, the BH natalkick prescription has the largest (but still mild) impacton our results.The zero kick and the momentum kick prescriptionslead, respectively, to a larger and smaller retention frac-tion of BHs compared to the fallback prescription [100].The difference becomes especially important in clusterswith initial mass M (cid:46) M (cid:12) because of their lower es-6cape velocities. In these clusters, virtually no BHs are leftafter the momentum kicks have been imparted, which isreflected in the mass distribution of useful clusters shownin the bottom-right panels of Fig. 7, 8, and 9. Cluster evaporation.
Our mass-independent and orbit-independent mass loss rate for cluster evaporation is cer-tainly a simplified one. To understand its effect on thecluster and BHB evolution, we computed three additionalmodels with exactly the same initial conditions as inMod1, Mod2 and Mod3 but with ˙ M (cid:63), ev = 0. Here westill compute the initial GCMF from equation (5) anduse the [ M c , ∆] values obtained from the MCMC anal-ysis above, but we do not include any prescription formass loss when evolving the clusters. Thus, this exerciseis only meant to determine the importance of the massloss effect on the secular evolution of the clusters and theBHBs they produce. We find that in these new mod-els, the local value of the merger rate density and of κ ,as well as the BH mass and progenitor cluster mass dis-tributions are consistent with those found in the modelswith cluster evaporation included. For the same ini-tial conditions as in Mod1, Mod2 and Mod3 the medianvalues of the local merger rate are R = 6 . − yr − ,14 . − yr − and 3 . − yr − , respectively. Thisshows that cluster evaporation has a small effect on thedynamics of the BHBs.In our models, however, tidal mass loss must becomeimportant at some point, e.g., for high enough ∆, GCswill evaporate before they can produce BHBs. We nowquantify how high ∆ needs to be in order to changethe BH dynamics significantly. To do this we comparethe tidal mass loss timescale, t ev ≡ M / | ˙ M (cid:63), ev | , to thetimescale after which the BHs have been nearly depletedby dynamical ejections, which we define to be t BH ≡ M BH , / | ˙ M BH | . We should expect that for t BH < t ev mostBHBs will have formed already before the cluster masshas changed significantly due to evaporation. This willhappen if ∆ is smaller than the critical value∆ c (cid:39) (cid:104) t (cid:105) t rh , βf BH M . (37)For ρ = 10 M (cid:12) pc − and f BH = 0 .
05, we find ∆ c ≈ M (cid:12) independent of the initial cluster mass; for ρ =10 M (cid:12) pc − , we have ∆ c ≈ M (cid:12) . These values arelarger than any value of ∆ used in our models (see Fig. 1),explaining the small impact of cluster evaporation on theresults.While the models discussed above show that the im-pact of cluster evaporation on the BHB dynamics issmall, they do not asses its effect on the merger rate.Thus, we consider three new models with ˙ M (cid:63), ev = 0 butnow use an initial GCMF that only accounts for massloss due to stellar evolution. If only stellar evolution isincluded, the initial GCMF that gives rise to the present-day GCMF shown in Fig. 1 becomes: φ (cid:48) cl , = 0 . A ( M / − exp [ − ( M / /M c ] , (38) FIG. 10. The upper panel shows the median of the mergerrate density distribution as a function of redshift; the middlepanel gives the distribution of primary BH mass for z < and K (cid:39)
2. These new models provide us with a safelower limit on the BHB merger rate for each density as-sumption; they are Mod15, Mod27 and Mod39 in Ta-ble I and Fig. 7, 8 and 9. From these results we see thatthe merger rate in models without evaporation are aboutthree times smaller than in models where the effect ofcluster evaporation is included.
Cluster density.
The model variations explored aboveshow that for a given initial GCMF, the initial clusterdensity is clearly the most important parameter for set-ting the BHB merger rate density and its redshift evo-lution. Thus, here we perform a more systematic explo-ration of such dependence by running an additional setof models where the initial cluster density is varied inthe range ρ = 10 M (cid:12) pc − to 10 M (cid:12) pc − . All other7 FIG. 11. Merger rate parameters as a function of the initial cluster half-mass density for the models of Fig. 10. The blackpoints represent median values, while the lower and upper error bars give the 5 and 95 percentiles of the distributions. model parameters are set as in Mod1 of Table I. The re-sults from these additional models are shown in Fig. 10and Fig. 11.From Fig. 10 we see that the peak of the merger rateand the redshift at which it occurs in each model increasewith ρ , and vary from (cid:39) − yr − at z = 3 for ρ = 10 M (cid:12) pc − , to (cid:39) Gpc − yr − at z = 4 . ρ = 10 M (cid:12) pc − . The situation is different, however,when we look at the merger rate in the local Universe.In Fig. 11 we see that the median value of R has amaximum value of (cid:39)
20 Gpc − yr − at ρ (cid:39) M (cid:12) pc − .This is an important result as it shows that the localBHB merger rate density from the GC channel has arobust upper limit of (cid:39)
50 Gpc − yr − – the upper errorbar estimate for the ρ = 10 M (cid:12) pc − model.The reason why R decreases with ρ above a cer-tain density can be understood from the lower panel inFig. 10. This plot shows the mass distribution of clustersfrom which the BHBs that merge in the local universeare formed. For initial densities above 10 M (cid:12) pc − thecontribution from clusters in the mass range 10 M (cid:12) (cid:46) M (cid:46) M (cid:12) gradually decreases as the distribution ofcluster masses contributing to the mergers becomes bi-modal. Such narrowing of the range of cluster massesthat can produce local mergers explains the relativelylow merger rate in the higher density models. It also af-fects the distribution of the primary BH masses as shownin the middle panel of Fig. 10.We looked into the density dependence of the clustermass distribution shown in the lower panel of Fig. 10in more details. We find that the lower mass peak seenin the cluster mass distribution for ρ = 10 M (cid:12) pc − and 10 M (cid:12) pc − is only due to ejected binaries while the higher mass peak is only due to in-cluster mergers. Thus,we can explain the depletion of BHBs that come from in-termediate mass clusters by considering the behaviour ofthe two merging populations when varying ρ . Abovea certain initial cluster mass, v esc becomes large enough( (cid:38) − ) that all BHBs merge inside the cluster.But, because v esc ∝ M / ρ / , the value of initial clustermass that still allows for dynamical ejections to occurgoes down as ρ increases. This explains why the up-per end of the mass distribution of clusters that producemergers from ejected BHBs decreases as ρ increases.Clusters with an initial mass larger than this value onlyproduce in-cluster mergers. Such clusters, however, canonly produce mergers in the local universe if they stillcontain BHs at the present time. Because the BHs areprocessed at a rate t − ∝ √ ρ/M , the value of the ini-tial cluster mass above which BHs are still present in thelocal universe increases with density. This explains whythe lower end of the mass distribution of systems thatproduce in-cluster mergers moves towards larger massesas ρ increases; clusters with a mass smaller than thisvalue get rid of all their BHs by z = 1.Fig. 10 shows that only models in which GCs startwith an initially high density ρ (cid:38) M (cid:12) pc − can ac-count for a large fraction of the LIGO-Virgo BHB merg-ers. Interestingly, these models also give a better fit tothe inferred BH mass function above m (cid:38) M (cid:12) as seenin Fig 2. Future GW observations will reduce the errorbars associated with the merger rate estimates and theBH mass distribution, providing important clues on thethe initial densities of GCs.Finally, we consider two additional model realizations.In one we evolve the same initial conditions as in [31] and8 FIG. 12. Merger rate evolution, primary BH mass of mergersat z < r h , =0.8 pc and the other half have r h , =1.6pc similar to [31](black lines), and for a model where r h , ∝ M . (blue lines). [62] where half of the clusters have r h , = 0 . r h , =1.6 pc; in the other model, thecluster half-mass radius scales aslog (cid:18) r h , pc (cid:19) = − .
56 + 0 .
615 log (cid:18) M M (cid:12) (cid:19) . (39)This latter relation was derived by [33] from the resultsof Has , egan et al. [61] who fit this Faber-Jackson-like re-lation to ultra-compact dwarf galaxies (UCDs) and el-liptical galaxies. Gieles et al. [33] derived the initialmass-radius relation correcting for mass loss and expan-sion by stellar evolution and correcting radii for projec-tion. All the other model parameters are the same as inMod1 of Table I. The results of these two new models areshown in Fig. 12. Interestingly, both give a local mergerrate, (cid:39)
10 Gpc − yr − , which is similar to the maximummerger rate value we obtained before. Moreover, theseresults show how the choice of initial half-mass radiusrelation has a significant effect on both the BH mass and initial cluster mass distributions. For r h , ∝ M . ,the cluster mass distribution becomes nearly flat so that,roughly speaking, all clusters with initial mass in therange 10 − M (cid:12) contribute equally to the local mergerrate. V. DISCUSSION AND CONCLUSIONS
Our models are based on some assumptions and ap-proximations. These are discussed in the following sec-tions, which also present a more detailed comparison tothe literature. We end the paper with a brief summaryof our main results.
A. Present-day GC mass density
Our value of ρ GC is larger than what was found by[18]: if we adopt their assumption of an average GCmass of 3 × M (cid:12) , we find that equation (2) corre-sponds to a GC number density of n GC = 2 . − ,which is 3 . n GC = 0 .
72 Mpc − ,but similar to the value used in [15] and [23]. Partof this difference is because we adopted a larger valueof η : if we use their mild M h -dependent η from [39],we find n GC = 1 .
50 Mpc − , which corresponds to ourlower error bar (cid:104) η (cid:105) . However, this is still a factor of 2.1higher than what was found by [18]. We are not surewhat causes this remaining difference, but we note that n GC = 1 .
50 Mpc − is about a factor h − larger than n GC = 0 .
72 Mpc − (Carl Rodriguez, private communi-cation). B. Initial GC density in the Universe
To derive ρ GC0 , a different approach was adopted by[31] and [62]. They use the total mass density of GCsforming in the semi-analytical galaxy formation modelof El-Badry et al. [80]. They approximate the numer-ical results with analytical functions and find a total ρ GC0 = 5 . × M (cid:12) Gpc − , about 15% higher thanEl-Badry et al. [80] and 20% lower than our adopted ρ GC (equation 2). They then assume that the initial massesof all GCs were a factor of 2.6 higher (from [101]) be-cause of stellar mass loss and evaporation and find thatinitial mass density of GCs more massive than 10 M (cid:12) is ρ GC0 ( M > M (cid:12) ) (cid:39) . × M (cid:12) Gpc − . Thisis a factor of ∼ M > M (cid:12) that lost (only) a factor of 2.6in mass implies a mass loss rate that is much lower in ourmodels. In our models the present-day ρ GC0 is made ofGCs with M (cid:38) × M (cid:12) .In addition, we do consider the contribution to themerger rate of lower mass GCs with M < M (cid:12) . We9Use cBHBd to evolve the same initial conditions as in[31] and [62] where half of the clusters have r h , = 0 . r h , =1.6pc but extendedthe initial GCMF down to M lo = 100 M (cid:12) (see Fig. 12).We find that (cid:39)
10% of the local mergers come from GCswith M < M (cid:12) , and therefore conclude that for theexact same initial GCMF, our models would lead to alocal merger rate that is still (cid:39) M c (cid:39) × M (cid:12) , so wecompare to the results from [31] for mass functions with M c = (2 . − × M (cid:12) . For these, the rates in [31]are 5 to 10Gpc − yr − . The median of the merger ratedistribution computed from our models is 10Gpc − yr − .Rodriguez and Loeb [31] use a fit to the results of a setof Monte Carlo simulations to determine the number ofmergers produced by a cluster as a function of time. Be-cause these fitting formulae are not public, it is currentlydifficult to establish the reason why our rates are only 1to 2 times, and not 4 times, those in [31]. We note inpassing that we compared our models to the number ofmergers from the two examples shown in Fig. 2B of [31]and found very good agreement. C. O-star ejections and IMBH formation
Our cBHBd model makes the simplifying assumptionthat all BHs are in place when the cluster forms. Becausethe typical timescale of GC evolution (i.e. 100 Myr - Gyr)is much longer than the timescale of BH formation (i.e.10 Myr), this is fine in most cases. However, for verydense, low-mass clusters, the relaxation time is so shortthat O-stars are ejected as ‘runaway stars’ before theyform BHs [102, 103]. As a result, the initial BH fractionis lower in these clusters than what we assume in ourmodels, possibly affecting the merger rate and propertiesof the mergers. To quantify this, we use the fact that thedynamical process that ejects O-stars is the same as theone that ejects BHs at a later stage. We therefore adopt clusterBH and replace the BH population by a massivestar population between 10 − M (cid:12) , with a logarithmicslope of − .
3, and a mass fraction of 18%, as appropri-ate for a Kroupa IMF. We then determine for a grid ofinitial cluster masses and half-mass radii the maximummass of massive stars that form BHs inside the cluster. InFig. 13 we show contours for 20 M (cid:12) (the minimum massof an O-star to produce a BH), 35 M (cid:12) (approximatelyhalf of the mass in BHs is produced by stars more mas-sive than this) and 100 M (cid:12) (the upper limit of our IMF).We also overplot the 3 initial cluster densities adopted inthe previous section. From this we see that clusters with M (cid:46) M (cid:12) are affected by O-star ejections, whichaffects about half of the mass in the initial GCMF. How-ever, these low-mass GCs are only responsible for ∼ cBHBd is re-peated mergers of BHs. After a merger, the BH mergerremnant receives a general relativistic momentum kickof several 100 km/s, and if this is smaller than the es-cape velocity from the center of the cluster, it can be in-volved in subsequent mergers [88–90] possibly forming anintermediate-mass BH (IMBH) [91]. This can only occurfor an initial escape velocity (cid:38)
300 km/s and in Fig. 13we show that only in our densest ( (cid:38) M (cid:12) / pc ), mostmassive ( (cid:38) M (cid:12) ) models this could happen. Ignor-ing the effect of IMBH formation via dry BH mergers istherefore not affecting our results. Although the forma-tion of IMBHs through repeated mergers is unlikely, wenote that hierarchical mergers can still contribute to theBH merger rate. As discussed above, hierarchical merg-ers represent only ten percent or less of the total numberof BHB mergers expected from GCs [89, 90]. Thus, theywill not affect significantly our integrated merger rate es-timates. On the other hand, second-generation mergerscan produce BHs with a mass higher than predicted bystellar evolution alone, broadening the BH mass distribu-tions we derived and populating them above ∼ M (cid:12) . D. Cluster mass loss and initial GCMF
Our models adopt a constant mass loss for all clus-ters of ∆ (cid:39) × M (cid:12) . This is what is required toevolve the initial GCMF with a power-law slope of − N -body simulations of tidally limited clusters show that˙ M ∝ M / [51], rather than ˙ M ∝ M . Including thismass dependence in ˙ M and maintaining the constraintthat all GCs formed with the same Universal initial massfunction, implies that clusters need to lose more massfor the turn-over in the GCMF to move to 2 × M (cid:12) [104], resulting in a twice as large value of K (cid:39)
64 [105]as we found for a mass independent ˙ M . Secondly, themodels of [51] show that ∆ for a typical Milky Way GCis smaller and depends on the apocenter and eccentric-ity of the Galactic orbit. For the median Galactocentricdistance of Miky Way GCs ( ∼ (cid:39) × M (cid:12) , i.e. a factor of 5smaller than what we assumed. These simulations con-sidered the secular evolution of clusters in a static tidalfield and therefore underestimate mass loss of clusters ifadditional disruption processes are important. For ex-ample, interactions with massive gas clouds in the earlyevolution (first Gyr) can be disruptive [54, 55], have asimilar mass dependence as relaxation driven evapora-tion in a static tidal field [58] and lead to a turn over in0 v e s c = k m / s = M / p c M / M l o g r h / p c r h M . ( G i e l e s + ) m a x ( m ) = M FIG. 13. GC initial mass-radius diagram showing the 3 ini-tial cluster densities adopted in this work with dashed lines.The magenta, full line shows the initial mass-radius relationthat describes the most massive GCs ( (cid:38) M (cid:12) ) and ultra-compact dwarf galaxies ( (cid:38) M (cid:12) ) [33, 61]. Full lines showthe maximum mass of O-stars that form BHs inside the clus-ter. Clusters with M (cid:46) M (cid:12) will eject some O-stars be-fore they become BHs, and these clusters will therefore haveslightly lower BHB mergers than in our model. The red,full line shows an initial escape velocity of 300 km/s, which isthe minimum escape velocity required for IMBH formation tooccur [91]. This process is not playing a role in our adoptedinitial conditions. the GCMF [56, 106]. If this is the cause for the value of∆, than | ˙ M | is much higher in the early evolution, whichwould affect the resulting merger rate. Because the re-laxation time decreases if the mass reduces, includingthis type of mass evolution will lead to a higher mergerrate than in our models with an ˙ M that is constant intime. The models of [51] also do not contain BHs and ithas been shown that retaining a BH population signifi-cantly increases the escape rate of stars [107, 108]. TheBH population can increase | ˙ M | by an order magnitude(Gieles et al., in prep), especially towards the end of theevolution. This implies a relatively low(high) | ˙ M | ( t rh ) inthe early evolution compared to our models, leading to areduction of the merger rate. In addition, for ˙ M ∝ M γ ,with γ <
0, the required K to get the turn-over at theright mass is lower than for γ = 0. If BHs are responsi-ble for the value of ∆, our merger rates could thereforealso be slightly overestimated for this reason. However,we do not expect this effect to be important for denseclusters ( (cid:38) M (cid:12) / pc ) because their BHB mergers areproduced when the clusters are still unaffected by theGalactic tides. We plan to include the effect of relax- ation driven escape in a tidal field in a future version of cBHBd to address this issue.Finally, we have assumed that all GC masses are drawnfrom an initial GCMF that is constrained by the shape ofthe Milky Way GCs. Although the present-day GCMFis remarkably universal across galaxies, variations in theinferred M c and ∆ values of a factor of ∼ M c and lower ∆ values are found in brightergalaxies. Although variations in M c and ∆ are partiallycaptured by the uncertainties in M c and ∆, this accountsonly for up to a factor of ∼
2. We may therefore under-populate the most massive clusters.
E. Primordial binaries
The effect of binaries that form in the star formationprocess and undergo stellar evolution in the first stagesof cluster evolution has not been discussed so far. We ar-gue here that primordial binaries have a negligible effecton the merger rate and the distribution of the BH masseswe derived. Because of H´enon principle, the energy gen-eration rate by binaries is determined by the relaxationprocess in the cluster as a whole. Whether dynamicallyactive binaries form in three-body encounters from sin-gle BHs, or in encounters involving BHBs that formedfrom primordial binaries will therefore result in a centralbinary with the same properties. However, primordialbinaries might affect the initial BH mass function due tobinary evolution processes. But, because BHB mergersfrom primordial binaries in GCs are a subdominant pop-ulation at low redshifts [see Fig. 2 in 18], the effect onthe local BHB mass distribution is also expected to besmall.
F. Conclusions
In this paper we have considered the dynamical for-mation of BHB mergers in GCs. Using our new pop-ulation synthesis code cBHBd we have evolved a largenumber of models covering a much wider set of initialconditions than explored in the literature. This allowedus to place robust error bars on the merger rate and massdistributions of the merging BHBs. We find that the GCchannel produces BHB mergers in the local universe ata rate of 7 . +21 . − . Gpc − yr − , where the error bars aremostly set by the unknown initial GC mass function andinitial cluster density. By comparing to the merger rateinferred by LIGO-Virgo, our results imply that a modelin which most of the detected mergers come from GCsis consistent with current constraints. This would re-quire, however, that GCs form with half-mass densitieslarger than (cid:38) M (cid:12) pc − , and suppression of other for-mation mechanisms. All our models show a drop in themerger rate of binary with primary BH mass outside therange (cid:39) − M (cid:12) , for which there is no evidence in1the gravitational wave data. This might suggest thatanother mechanism is responsible for the production ofthese sources.Our results have a number of implications for the for-mation of BHB mergers and GCs. The dependence ofthe merger rate and BHB properties (e.g., eccentricity,mass) on the model parameters suggests that a directcomparison to the gravitational wave data will allow usto place constraints on the initial conditions of GCs andtheir evolution. Our models will also help to understandother uncertain parameters that control the formation ofBHs and their natal kicks. While these latter param-eters have little effect on the merger rate, they have asignificantly impact on the masses of the merging BHBs.Thus, useful constraints could be placed once the numberof gravitational wave detections will be large enough toto allow for a statistically significant comparison to theinferred BH mass function.In the future, we plan to consider other type of clus-ters such as open and nuclear star clusters which are also believed to be efficient factories of gravitational wavesources [88, 91, 109–111]. The study of these systems willrequire us to add additional physics to cBHBd . VI. ACKNOWLEDGEMENTS
We thank the internal referee of the LVC, MichelaMapelli, for her suggestions, which helped us improvethis work. MG thanks Duncan Forbes for helpful dis-cussions on the relation between halo mass and GCpopulation mass. We thank the referee for their con-structive comments that helped to improved the paper.FA acknowledges support from a Rutherford fellowship(ST/P00492X/1) from the Science and Technology Fa-cilities Council. We acknowledge the support of the Su-percomputing Wales project, which is part-funded bythe European Regional Development Fund (ERDF) viaWelsh Government. [1] LIGO Scientific Collaboration, Classical and QuantumGravity , 074001 (2015), arXiv:1411.4547 [gr-qc].[2] Virgo Collaboration, Classical and Quantum Gravity , 024001 (2015), arXiv:1408.3978 [gr-qc].[3] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Aber-nathy, F. Acernese, K. Ackley, C. Adams, T. Adams,P. Addesso, and R. X. Adhikari, Phys. Rev. Lett. ,061102 (2016), arXiv:1602.03837 [gr-qc].[4] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Aber-nathy, F. Acernese, K. Ackley, C. Adams, T. Adams,P. Addesso, and R. X. Adhikari, Physical Review X ,041015 (2016), arXiv:1606.04856 [gr-qc].[5] LIGO Scientific Collaboration and Virgo Collaboration,The Astrophysical Journal Letters , L24 (2019),arXiv:1811.12940 [astro-ph.HE].[6] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham,F. Acernese, K. Ackley, C. Adams, R. X. Adhikari, V. B.Adya, C. Affeldt, and et al., Physical Review X ,031040 (2019), arXiv:1811.12907 [astro-ph.HE].[7] LIGO Scientific Collaboration and Virgo Collaboration,Phys. Rev. D , 064064 (2019), arXiv:1907.09384[astro-ph.HE].[8] LIGO Scientific Collaboration and Virgo Collaboration,The Astrophysical Journal Letters , L44 (2020),arXiv:2006.12611 [astro-ph.HE].[9] The LIGO Scientific Collaboration and the Virgo Col-laboration, arXiv e-prints , arXiv:2004.08342 (2020),arXiv:2004.08342 [astro-ph.HE].[10] LIGO Scientific Collaboration and Virgo Collaboration,Phys. Rev. Lett. , 101102 (2020), arXiv:2009.01075[gr-qc].[11] R. Abbott et al. , (2020a), arXiv:2010.14527 [gr-qc].[12] R. Abbott et al. (LIGO Scientific, Virgo), (2020),arXiv:2010.14533 [astro-ph.HE].[13] S. Sigurdsson and L. Hernquist, Nature , 423 (1993).[14] S. R. Kulkarni, P. Hut, and S. J. McMillan, Nature , 421 (1993). [15] S. F. Portegies Zwart and S. L. W. McMillan, TheAstrophysical Journal Letters , L17 (2000), astro-ph/9910061.[16] S. Banerjee, H. Baumgardt, and P. Kroupa, MonthlyNotices of the Royal Astronomical Society , 371(2010).[17] J. M. B. Downing, M. J. Benacquista, M. Giersz, andR. Spurzem, Monthly Notices of the Royal AstronomicalSociety , no (2011).[18] C. L. Rodriguez, M. Morscher, B. Pattabiraman,S. Chatterjee, C.-J. Haster, and F. A. Rasio, Physi-cal Review Letters , 051101 (2015).[19] A. Askar, M. Szkudlarek, D. Gondek-Rosi´nska,M. Giersz, and T. Bulik, MNRAS , L36 (2017),arXiv:1608.02520 [astro-ph.HE].[20] D. Park, C. Kim, H. M. Lee, Y.-B. Bae, and K. Bel-czynski, MNRAS , 4665 (2017), arXiv:1703.01568[astro-ph.HE].[21] J. Samsing, A. Askar, and M. Giersz, Astrophys. J. , 124 (2018), arXiv:1712.06186 [astro-ph.HE].[22] M. Zevin, J. Samsing, C. Rodriguez, C.-J. Haster,and E. Ramirez-Ruiz, Astrophys. J. , 91 (2019),arXiv:1810.00901 [astro-ph.HE].[23] K. Kremer, C. S. Ye, N. Z. Rui, N. C. Weatherford,S. Chatterjee, G. Fragione, C. L. Rodriguez, M. Spera,and F. A. Rasio, The Astrophysical Journal SupplementSeries , 48 (2020), arXiv:1911.00018 [astro-ph.HE].[24] G. Fragione and B. Kocsis, Phys. Rev. Lett. , 161103(2018), arXiv:1806.02351 [astro-ph.GA].[25] O. Y. Gnedin and J. P. Ostriker, Astrophys. J. , 223(1997).[26] E. Vesperini and D. C. Heggie, MNRAS , 898(1997), arXiv:astro-ph/9705073 [astro-ph].[27] S. M. Fall and Q. Zhang, Astrophys. J. , 751 (2001),arXiv:astro-ph/0107298 [astro-ph].[28] H. Baumgardt and J. Makino, Monthly Notices of theRoyal Astronomical Society , 227 (2003). [29] M. Gieles and H. Baumgardt, MNRAS , L28 (2008),arXiv:0806.2327.[30] A. Jord´an, D. E. McLaughlin, P. Cˆot´e, L. Ferrarese,E. W. Peng, S. Mei, D. Villegas, D. Merritt, J. L. Tonry,and M. J. West, The Astrophysical Journal SupplementSeries , 101 (2007), arXiv:astro-ph/0702496.[31] C. L. Rodriguez and A. Loeb, Astrophys. J. , L5(2018), arXiv:1809.01152 [astro-ph.HE].[32] S. Chatterjee, C. L. Rodriguez, and F. A. Rasio, As-trophys. J. , 68 (2017), arXiv:1603.00884.[33] M. Gieles, H. Baumgardt, D. C. Heggie, andH. J. G. L. M. Lamers, MNRAS , L16 (2010),arXiv:1007.2333 [astro-ph.GA].[34] F. Antonini and M. Gieles, MNRAS , 2936 (2020),arXiv:1906.11855 [astro-ph.HE].[35] C. L. Rodriguez, S. Chatterjee, and F. A. Rasio, Phys.Rev. D , 084029 (2016), arXiv:1602.02444 [astro-ph.HE].[36] L. R. Spitler and D. A. Forbes, MNRAS , L1 (2009),arXiv:0809.5057 [astro-ph].[37] I. Y. Georgiev, T. H. Puzia, P. Goudfrooij, andM. Hilker, MNRAS , 1967 (2010), arXiv:1004.2039[astro-ph.CO].[38] W. E. Harris, G. L. H. Harris, and M. Alessi, TheAstrophysical Journal , 82 (2013).[39] W. E. Harris, G. L. Harris, and M. J. Hudson, As-trophys. J. , 36 (2015), arXiv:1504.03199 [astro-ph.GA].[40] W. E. Harris, J. P. Blakeslee, and G. L. H. Harris,Astrophys. J. , 67 (2017), arXiv:1701.04845 [astro-ph.GA].[41] W. E. Harris, Astronomical Journal , 1487 (1996).[42] W. E. Harris, arXiv:1012.3224 (2010), arXiv:1012.3224[astro-ph.GA].[43] L. L. Watkins, R. P. van der Marel, S. T. Sohn,and N. W. Evans, Astrophys. J. , 118 (2019),arXiv:1804.11348 [astro-ph.GA].[44] J. L. Tinker, B. E. Robertson, A. V. Kravtsov,A. Klypin, M. S. Warren, G. Yepes, and S. Gottl¨ober,The Astrophysical Journal , 878 (2010).[45] S. G. Murray, C. Power, and A. S. G. Robotham, As-tronomy and Computing , 23 (2013), arXiv:1306.6721[astro-ph.CO].[46] Planck Collaboration, arXiv:1807.06209 ,arXiv:1807.06209 (2018), arXiv:1807.06209 [astro-ph.CO].[47] D. A. Forbes, J. I. Read, M. Gieles, and M. L. M.Collins, MNRAS , 5592 (2018), arXiv:1809.07831[astro-ph.GA].[48] P. Schechter, “An analytic expression for the luminosityfunction for galaxies.” (1976).[49] S. F. Portegies Zwart, S. L. McMillan, and M. Gieles,Annual Review of Astronomy and Astrophysics , 431(2010).[50] V. Bromm and C. J. Clarke, The Astrophysical JournalLetters , L1 (2002), arXiv:astro-ph/0201066.[51] H. Baumgardt and J. Makino, MNRAS , 227 (2003),arXiv:astro-ph/0211471 [astro-ph].[52] J. P. Ostriker, L. J. Spitzer, and R. A. Chevalier, TheAstrophysical Journal Letters , L51 (1972).[53] O. Y. Gnedin, L. Hernquist, and J. P. Ostriker, Astro-phys. J. , 109 (1999).[54] L. J. Spitzer, Astrophys. J. , 17 (1958). [55] M. Gieles, S. F. Portegies Zwart, H. Baumgardt,E. Athanassoula, H. J. G. L. M. Lamers, M. Sip-ior, and J. Leenaarts, MNRAS , 793 (2006), astro-ph/0606451.[56] B. G. Elmegreen, The Astrophysical Journal Letters , L184 (2010), arXiv:1003.0798.[57] J. M. D. Kruijssen, Monthly Notices of the Royal As-tronomical Society , 1658 (2015).[58] M. Gieles and F. Renaud, MNRAS , L103 (2016),arXiv:1605.05940.[59] M. H´enon, Annales d’Astrophysique , 369 (1961).[60] D. Foreman-Mackey, D. W. Hogg, D. Lang, andJ. Goodman, arXiv preprint arXiv: . . . , 1 (2012).[61] M. Has , egan, A. Jord´an, P. Cˆot´e, S. G. Djorgovski, D. E.McLaughlin, J. P. Blakeslee, S. Mei, M. J. West, E. W.Peng, and L. Ferrarese, Astrophys. J. , 203 (2005),arXiv:astro-ph/0503566 [astro-ph].[62] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee,K. Kremer, F. A. Rasio, J. Samsing, C. S. Ye,and M. Zevin, Phys. Rev. D , 123005 (2018),arXiv:1811.04926 [astro-ph.HE].[63] M. H´enon, in Dynamics of the Solar Systems , IAU Sym-posium, Vol. 69, edited by A. Hayli (1975) p. 133.[64] P. G. Breen and D. C. Heggie, Monthly Notices of theRoyal Astronomical Society , 2779 (2013).[65] M. Gieles, D. C. Heggie, and H. Zhao, MNRAS ,2509 (2011), arXiv:1101.1821 [astro-ph.GA].[66] L. Spitzer, Jr. and M. H. Hart, Astrophys. J. , 399(1971).[67] K. G¨ultekin, M. C. Miller, and D. P. Hamilton, Astro-phys. J. , 221 (2004), astro-ph/0402532.[68] J. Samsing, Phys. Rev. D , 103014 (2018),arXiv:1711.07452 [astro-ph.HE].[69] D. C. Heggie, Monthly Notices of the Royal Astronom-ical Society , 729 (1975).[70] J. Binney and S. Tremaine, Galactic Dynamics: SecondEdition, by James Binney and Scott Tremaine. ISBN978-0-691-13026-2 (HB). Published by Princeton Uni-versity Press, Princeton, NJ USA, 2008. (PrincetonUniversity Press, 2008).[71] P. Kroupa, MNRAS , 231 (2001).[72] J. R. Hurley, O. R. Pols, and C. A. Tout, MNRAS ,543 (2000), arXiv:astro-ph/0001295.[73] J. S. Vink, A. de Koter, and H. J. G. L. M. Lamers,Astronomy and Astrophysics , 574 (2001).[74] K. Belczynski, T. Bulik, C. L. Fryer, A. Ruiter,F. Valsecchi, J. S. Vink, and J. R. Hurley, Astrophys.J. , 1217 (2010), arXiv:0904.2784 [astro-ph.SR].[75] C. L. Fryer, K. Belczynski, G. Wiktorowicz, M. Do-minik, V. Kalogera, and D. E. Holz, The AstrophysicalJournal , 91 (2012).[76] K. Belczynski, A. Heger, W. Gladysz, A. J. Ruiter,S. Woosley, G. Wiktorowicz, H. Y. Chen, T. Bulik,R. O’Shaughnessy, D. E. Holz, C. L. Fryer, andE. Berti, Astronomy & Astrophysics , A97 (2016),arXiv:1607.03116 [astro-ph.HE].[77] G. Hobbs, D. R. Lorimer, A. G. Lyne, and M. Kramer,Monthly Notices of the Royal Astronomical Society ,974 (2005).[78] C. L. Rodriguez, S. Chatterjee, and F. A. Rasio, Phys-ical Review D , 084029 (2016).[79] D. A. VandenBerg, K. Brogaard, R. Leaman, andL. Casagrande, The Astrophysical Journal , 134(2013). [80] K. El-Badry, E. Quataert, D. R. Weisz, N. Choksi,and M. Boylan-Kolchin, MNRAS , 4528 (2019),arXiv:1805.03652 [astro-ph.GA].[81] J. Samsing, T. Venumadhav, L. Dai, I. Martinez,A. Batta, M. Lopez, E. Ramirez-Ruiz, and K. Kremer,Phys. Rev. D , 043009 (2019), arXiv:1901.02889[astro-ph.HE].[82] P. Madau and M. Dickinson, Annu Rev Astron Astr ,415 (2014), arXiv:1403.0007 [astro-ph.CO].[83] LIGO Scientific Collaboration and J. Harms, Clas-sical and Quantum Gravity , 044001 (2017),arXiv:1607.08697 [astro-ph.IM].[84] M. P. et al., Classical and Quantum Gravity , 194002(2010).[85] K. Belczynski, V. Kalogera, and T. Bulik, The Astro-physical Journal , 407 (2002).[86] K. Belczynski, V. Kalogera, F. A. Rasio, R. E. Taam,A. Zezas, T. Bulik, T. J. Maccarone, and N. Ivanova,The Astrophysical Journal Supplement Series , 223(2008).[87] M. Spera, M. Mapelli, and A. Bressan, Monthly Noticesof the Royal Astronomical Society , 4086 (2015).[88] F. Antonini and F. A. Rasio, Astrophys. J. , 187(2016), arXiv:1606.04889 [astro-ph.HE].[89] C. L. Rodriguez, M. Zevin, P. Amaro-Seoane, S. Chat-terjee, K. Kremer, F. A. Rasio, and C. S. Ye, Phys. Rev.D , 043027 (2019), arXiv:1906.10260 [astro-ph.HE].[90] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, andF. A. Rasio, Phys. Rev. Lett. , 151101 (2018),arXiv:1712.04937 [astro-ph.HE].[91] F. Antonini, M. Gieles, and A. Gualandris, MNRAS , 5008 (2019), arXiv:1811.03640 [astro-ph.HE].[92] M. Fishbach, D. E. Holz, and B. Farr, The Astrophys-ical Journal Letters , L24 (2017), arXiv:1703.06869[astro-ph.HE].[93] D. Gerosa and E. Berti, Phys. Rev. D , 124046 (2017),arXiv:1703.06223 [gr-qc].[94] C. Kimball, C. Talbot, C. P. L. Berry, M. Car-ney, M. Zevin, E. Thrane, and V. Kalogera, arXive-prints , arXiv:2005.00023 (2020), arXiv:2005.00023[astro-ph.HE].[95] F. Antonini, N. Murray, and S. Mikkola, The Astro-physical Journal , 45 (2014). [96] J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, TheAstrophysical Journal , 71 (2014).[97] F. Antonini, S. Chatterjee, C. L. Rodriguez,M. Morscher, B. Pattabiraman, V. Kalogera, andF. A. Rasio, The Astrophysical Journal , 65 (2016).[98] LIGO Scientific Collaboration and Virgo Collaboration,Astrophys. J. , 149 (2019).[99] M. H´enon, Annales d’Astrophysique , 62 (1965).[100] S. Banerjee, K. Belczynski, C. L. Fryer, P. Berczik,J. R. Hurley, R. Spurzem, and L. Wang, arXive-prints , arXiv:1902.07718 (2019), arXiv:1902.07718[astro-ph.SR].[101] N. Choksi, M. Volonteri, M. Colpi, O. Y. Gnedin, andH. Li, Astrophys. J. , 100 (2019), arXiv:1809.01164[astro-ph.GA].[102] A. Poveda, J. Ruiz, and C. Allen, Boletin de los Ob-servatorios Tonantzintla y Tacubaya , 86 (1967).[103] M. S. Fujii and S. Portegies Zwart, Science , 1380(2011), arXiv:1111.3644 [astro-ph.GA].[104] M. Gieles, MNRAS , 2113 (2009), 0901.0830.[105] J. M. D. Kruijssen and S. F. Portegies Zwart,The Astrophysical Journal Letters , L158 (2009),arXiv:0905.3744 [astro-ph.GA].[106] J. M. D. Kruijssen, MNRAS , 1658 (2015),arXiv:1509.02163.[107] M. Giersz, A. Askar, L. Wang, A. Hypki, A. Lev-eque, and R. Spurzem, MNRAS , 2412 (2019),arXiv:1904.01227 [astro-ph.GA].[108] L. Wang, MNRAS , 2413 (2020), arXiv:1911.05077[astro-ph.GA].[109] B. M. Ziosi, M. Mapelli, M. Branchesi, and G. Tormen,Monthly Notices of the Royal Astronomical Society ,3703 (2014).[110] U. N. Di Carlo, N. Giacobbo, M. Mapelli, M. Pasquato,M. Spera, L. Wang, and F. Haardt, MNRAS , 2947(2019), arXiv:1901.00863 [astro-ph.HE].[111] S. Rastello, P. Amaro-Seoane, M. Arca-Sedda,R. Capuzzo-Dolcetta, G. Fragione, and I. Tosta e Melo,MNRAS483