#  Coalescing black hole binaries from globular clusters: mass distributions
and comparison to gravitational wave data from GWTC–3

Fabio Antonini1, Mark Gieles2,3, Fani Dosopoulou4,5, and Debatri
Chattopadhyay1
1Gravity Exploration Institute, School of Physics and Astronomy, Cardiff
University, Cardiff, CF24 3AA, UK;
2 ICREA, Pg. Lluís Companys 23, E08010 Barcelona, Spain;
3Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB),
Martí Franquès 1, E08028 Barcelona, Spain
4 Princeton Center for Theoretical Science, Princeton University, Princeton,
NJ 08544, USA
5 Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA E-mail<EMAIL_ADDRESS>
###### Abstract

We use our cluster population model, cBHBd, to explore the mass distribution
of merging black hole binaries formed dynamically in globular clusters. We
include in our models the effect of mass growth through hierarchical mergers
and compare the resulting distributions to those inferred from the third
gravitational wave transient catalogue. We find that none of our models can
reproduce the peak at $m_{1}\simeq 10M_{\odot}$ in the primary black hole mass
distribution that is inferred from the data. This disfavours a scenario where
most of the sources are formed in globular clusters. On the other hand, a
globular cluster origin can account for the inferred secondary peak at
$m_{1}\simeq 35M_{\odot}$, which requires that the most massive clusters form
with half-mass densities $\rho_{\rm h,0}\gtrsim 10^{4}~{}{\rm M}_{\odot}\,{\rm
pc}^{-3}\ $. Finally, we find that the lack of a high mass cut–off in the
inferred mass distribution can be explained by the repopulation of an initial
mass gap through hierarchical mergers. Matching the inferred merger rate above
$\simeq 50M_{\odot}$ requires both initial cluster densities $\rho_{\rm
h,0}\gtrsim 10^{4}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}\ $, and that black holes
form with nearly zero spin. A hierarchical merger scenario makes specific
predictions for the appearance and position of multiple peaks in the black
hole mass distribution, which can be tested against future data.

###### keywords:

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

††pagerange:  Coalescing black hole binaries from globular clusters: mass
distributions and comparison to gravitational wave data from GWTC–3–A.2

##  1 Introduction

The analysis of gravitational wave (GW) observations has identified structures
in the mass distribution of the observed population (Tiwari & Fairhurst,
2021). Some of these structures already emerged from the analysis of the
second gravitational wave transient catalog (GWTC-2; Abbott et al., 2021b;
Abbott, 2020a). However, thanks to the increased number of events in the new
GWTC-3 (Abbott et al., 2021c, a), we are now more confident of their
statistical significance. In particular, three important features in the
underlining BH mass distribution have been uncovered: (i) the distribution of
primary BH masses has a strong peak at about $\simeq 10M_{\odot}$; (ii) there
is clear evidence for a secondary peak at $\simeq 35M_{\odot}$; and (iii)
there is no evidence for any mass gap above $\approx 40$–$60M_{\odot}$, which
is predicted by stellar evolution models due to pulsational pair–instability
and pair–instability in massive stars (e.g., Woosley, 2017; Spera & Mapelli,
2017; Olejak et al., 2022). In this article, we perform a large number of
cluster simulations to understand whether (i), (ii) and (iii) can be explained
by a globular cluster (GC) origin for the sources.

The formation of BH binary mergers, including those with components above the
upper mass gap, might be explained by several formation pathways. These
include binary stellar evolution (e.g., Dominik et al., 2012; Mink &
Belczynski, 2015; Mandel & De Mink, 2016; Farmer et al., 2020; Costa et al.,
2021), multiple star interactions (e.g., Silsbee & Tremaine, 2017; Antonini et
al., 2017; Fragione et al., 2020; Liu & Lai, 2021; Hamers et al., 2021;
Stegmann et al., 2022), stellar collisions in open clusters (e.g., Kremer et
al., 2020; Di Carlo et al., 2020; Dall’Amico et al., 2021; Banerjee, 2021a;
Chattopadhyay et al., 2022), primordial BHs (e.g., Ballesteros et al., 2018;
Gow et al., 2020; De Luca et al., 2021), and formation in active galactic
nuclei (e.g., Bartos et al., 2016; Stone et al., 2016; Yang et al., 2019;
Tagawa et al., 2021b). One widely discussed scenario is formation through
three body dynamical interactions in dense stellar environments such as
nuclear star clusters (e.g., O’Leary et al., 2009; Miller & Lauburg, 2009;
Antonini & Rasio, 2016; Antonini et al., 2019; Fragione et al., 2022) and GCs
(e.g., Sigurdsson & Hernquist, 1993; Kulkarni et al., 1993; Banerjee et al.,
2010; Rodriguez et al., 2015, 2016; Askar et al., 2017; Fragione & Kocsis,
2018; Chattopadhyay et al., 2022).

The mass distribution of coalescing BH binaries produced in GCs has been
investigated in several studies (e.g., Rodriguez et al., 2016; Askar et al.,
2017; Park et al., 2017; Antonini & Gieles, 2020a; Mapelli et al., 2022; Zevin
& Holz, 2022). Previous work suggests that GCs are an environment where BH
binaries can efficiently assemble and merge, providing one of the main
formation channels of BH binary coalescences in the Universe (Portegies Zwart
& McMillan, 2000). In particular, it has been argued that due to the high
escape velocities of GCs, BH mass growth can occur through consecutive
mergers, populating any mass gap created by stellar processes (Rodriguez et
al., 2018; Kimball et al., 2020; Rodriguez et al., 2020; Doctor et al., 2020;
Kimball et al., 2021; Tagawa et al., 2021a). In this scenario, a BH that is
formed from a previous merger and is retained inside the cluster, sinks back
to the cluster core where it dynamically couples with another BH and merges
with it after a series of binary-single encounters. If this process repeats
multiple times, significant mass growth can occur (Antonini & Rasio, 2016;
Fishbach et al., 2017; Gerosa & Berti, 2017). A direct comparison of model
predictions to data, however, are rare (Mould, Gerosa, & Taylor, 2022). It
remains therefore an open question whether a GC origin provides a plausible
explanation for the inconclusive evidence for an upper mass gap in the GW
data, and whether the other features of the inferred BH mass distribution can
also be reproduced. A putative successful GC model will then provide useful
constraints on the properties of GCs and their BHs at birth.

In this work, we adopt our new fast method for the evolution of star clusters
and their BBHs, cBHBd (Antonini & Gieles, 2020b), to study the mass
distribution of BHs produced dynamically in GCs, including the effect of
hierarchical mergers and a novel recipe for sampling masses of the BBH
components and the interlopers. Our efficient approach allows us to address
how model assumptions affect the final results, place error bars on merger
rate estimates, compare to the distributions inferred from the new GW data
catalog GWTC-3, and, finally, asses a hierarchical merger origin for the
formation of the most massive BHs detected by LIGO and Virgo.

In Section 2 we describe our methodology and approximations. In Section 3 we
describe our main results, and the importance of model assumptions. We
conclude and summarise our results in Section 4.

##  2 cluster models with hierarchical mergers

We simulate the evolution of BH binaries in star clusters using our code
cBHBd, which we modify in order to include hierarchical mergers. We define
here hierarchical mergers as binary mergers in which at least one of the two
BH components is a BH remnant that was formed from a previous merger.

Our method is based on Hénon’s principle (Hénon, 1975) which states that the
rate of heat generation in the core is a constant fraction of the total
cluster energy per half-mass relaxation time. Thus, the energy production rate
in the core, which we assume is produced by BH binaries, is regulated by the
energy demand of the entire system (Breen & Heggie, 2013). The lifetime and
the merger rate of BHs in the cluster can be linked to the evolution of the
cluster itself as described in details in Antonini & Gieles (2020b). Then,
three ingredients are needed in order to determine the formation of BH
binaries, their merger rate and their properties: (i) a model for the
evolution of the cluster global properties; (ii) a model for binary BH
dynamics; and (iii) a realistic set of initial BH masses.

We start by sampling the masses of the stellar progenitors of BHs from a
standard mass function, $\phi(m_{\star})\propto m_{\star}^{-2.3}$ (Salpeter,
1955; Kroupa, 2001), with masses in the range $20M_{\odot}$ to $130M_{\odot}$.
For a given cluster metallicity, $Z$, we evolve the stars to BHs using the
Single Stellar Evolution (SSE) package (Hurley et al., 2002), which we
modified to include updated prescriptions for stellar winds and mass loss
(following Vink et al., 2001), and for pair–instability in massive stars
(following Spera & Mapelli, 2017). We therefore evolve the BH progenitors as
single stars, assuming a zero binary fraction initially. At the end of this
phase, the total number of BHs in a cluster model is calculate by assuming a
Kroupa (2001) initial mass function in the mass range $0.1M_{\odot}$ to
$130M_{\odot}$. The value of the largest BH mass formed in the model depends
on metallicity and varies between $\simeq 25M_{\odot}$ for $Z=2\times 10^{-2}$
and $\simeq 55M_{\odot}$ for $Z=1\times 10^{-4}$. For each BH we compute a
natal kick velocity from a Maxwellian distribution with dispersion $265{\rm
km\ s^{-1}}$ (Hobbs et al., 2005), lowered by the amount of mass that falls
back into the forming compact object (Fryer et al., 2012). In most of our
models, we start the BHs all with the same value of the spin angular momentum,
$\chi$, where $\vec{\chi}=\vec{S}/m^{2}$ is the dimensionless spin of the BH
and $\vec{S}$ is the spin angular momentum in units of $m^{2}$. In one model
the initial value of $\chi$ is sampled from a distribution that is consistent
with that inferred from the GW data and is given by the median distribution
shown in figure 15 of Abbott et al. (2021c).

Then, we initialise and evolve the cluster model. The initial conditions are
determined by three parameters: the cluster density, $\rho_{\rm h,0}$; the
cluster mass, $M_{\rm 0}$; and the total mass in BHs, $M_{\rm BH,0}$. The
latter is set equal to the total mass in BHs obtained with SSE, assuming a
Kroupa initial mass function in the range $0.08-130\ M_{\odot}$ and taking
into account that a fraction of the BHs are ejected from the cluster by a
natal kick. The time evolution of the cluster properties is then obtained as
in Antonini & Gieles (2020a). Briefly, we integrate a set of first order
differential equations which determine the time evolution of $M$, $M_{\rm
BH}$, and the cluster half mass radius, $r_{\rm h}$. These models include
simple prescriptions for mass loss and expansion due to stellar evolution and
cluster ‘evaporation’, while BHs are assumed to be lost through dynamical
ejections.

Finally, we dynamically evolve the BH binaries that form via three-body
processes in the cluster core. Our treatment of binary BH formation and
evolution follows closely Antonini & Gieles (2020b). The first binary BH forms
after the cluster core-collapse time

$\tau_{\rm cc}=3.21t_{\rm rh,0}$ (1)

where $t_{\rm rh,0}$ is the initial cluster half-mass relaxation time (for the
definition, see equation 10 in Antonini & Gieles, 2020b). We assume that the
binary is formed with a semi-major axis at the soft-hard boundary, $a_{\rm
h}\simeq{G\mu/\sigma^{2}}$, with $\mu=m_{1}m_{2}/(m_{1}+m_{2})$, where $m_{1}$
and $m_{2}$ are the masses of the binary components, and $m_{1}>m_{2}$. The
expression of $a_{\rm h}$ above is only approximate, and valid under the
assumption of equal mass components. Later in Appendix A we introduce the
quantity $\beta$ and equipartition among BHs of different masses, then the
definition is $a_{\rm h}=0.5Gm_{1}m_{2}\beta$.

The pairing of BHs is done by sampling their masses from the set of BHs still
left inside the cluster. We first draw two mass values from the power law
probability distributions $p(m_{1})\propto m_{1}^{\alpha_{1}}$ and
$p(q)=q^{\alpha_{2}}$, with $\alpha_{1}=8+2\alpha$, $\alpha_{2}=3.5+\alpha$
and $q=m_{2}/m_{1}$. Here $\alpha$ is the power law index of the BH mass
function, which also evolves with time as the BH population is depleted. The
two BH components are then selected by choosing the two BHs that have the mass
closest to the values drawn from $p(m_{1})$ and $p(q)$ (or $p(m_{2})$).

Once selected, the binary is evolved through a sequence of binary-single
encounters. Similarly, we find the mass of the third BH interloper from the
power law distribution $p(m_{3})\propto m_{3}^{\alpha_{3}}$ with
$\alpha_{3}=\alpha+1/2$. The adopted expressions for $p(m_{1})$, $p(q)$ and
$p(m_{3})$ are motivated below in Appendix A. The power law exponent,
$\alpha$, is obtained at the start of the integration for each cluster from a
fit to the initial BH mass function after removing BHs that are ejected by
natal kicks. The value of $\alpha$ as well as the lower and upper bound of the
BH mass function are then recalculated after each time-step. Specifically, the
lower bound of the BH mass function is set equal to the mass of the lightest
BH in the cluster, and the upper bound is the mass of the most massive BH.
This procedure allows to take into account the evolution of the BH mass
function with time due to ejections and the growth of BHs through hierarchical
mergers.

Following Samsing (2018), we divide each binary-single encounter in a set of
$N_{\rm rs}=20$ resonant intermediate states and assume that the eccentricity
of the binary after each state is sampled from a thermal distribution
$N(<e)\propto e^{2}$. If

$\ \sqrt{1-e^{2}}<h\left(R_{\rm S}\over a\right)^{5/14}$ (2)

a merger occurs through a GW capture before the next intermediate binary-
single state is formed, where $R_{\rm S}=2G(m_{1}+m_{2})/c^{2}$ and $h$ is a
constant of order unity.

If the binary survives the 20 intermediate resonant states, we compute: (i)
the new binary semi-major axis, assuming that its binding energy decreases by
the fixed fraction $\Delta E/E=0.2$ (Samsing, 2018). (ii) the recoil kick due
to energy and angular momentum conservation experienced by the binary centre
of mass (Antonini & Rasio, 2016)

$v_{\rm bin}^{2}=0.2G{m_{1}m_{2}\over m_{1}+m_{2}+m_{3}}q_{3}/a$ (3)

with $q_{3}=m_{3}/(m_{1}+m_{2})$, and (iii) the recoil kick experienced by the
interloper:

$v_{\rm 3}=v_{\rm bin}/q_{3}\ .$ (4)

If $v_{\rm bin}>v_{\rm esc}$, the binary is ejected from the cluster; if
$v_{\rm 3}>v_{\rm esc}$, the interloper is also ejected from the cluster. If
the binary is ejected from the cluster, we compute its merger timescale due to
GW energy loss using the standard Peter’s formula (Peters, 1964).

If $v_{\rm bin}<v_{\rm esc}$, and the binary angular momentum at the end of
the triple interaction is such that (Antonini & Gieles, 2020b)

$\sqrt{1-e^{2}}<1.3\left[G^{4}(m_{1}m_{2})^{2}(m_{1}+m_{2})\over
c^{5}\dot{E}_{\rm bin}\right]^{1/7}a^{-5/7}\ ,$ (5)

then the BH binary merges before the next binary-single encounter takes place.
Binaries that undergo this type of evolution are often named ‘in-cluster
inspirals’ (Samsing, 2018; Rodriguez et al., 2018). We then assign the new
remnant BH a GW recoil kick, $v_{\rm GW}$, and compute its new spin and mass
following Rezzolla et al. (2008). If $v_{\rm GW}>v_{\rm esc}$ the remnant is
ejected from the cluster, otherwise we compute the dynamical friction
timescale to sink back to the cluster core

$\tau_{\rm df}\simeq 1.65\,r_{\rm in}^{2}{\sigma\over{\ln\Lambda Gm}}$ (6)

where (Antonini et al., 2019)

$r_{\rm in}=r_{\rm h}\sqrt{{v_{\rm esc}^{4}\over\left(v_{\rm esc}^{2}-v_{\rm
GW}^{2}\right)^{2}}-1}$ (7)

and only allow the BH to form a new binary after this time. If $v_{\rm
bin}<v_{\rm esc}$, but condition equation (5) is not satisfied, then a new
interloper is sampled from the given distribution and the binary is evolved
through a new binary-single interaction.

Each binary is evolved through a sequence of binary-single encounters until
either a merger occurs or it is ejected from the cluster. Then a new binary is
formed. The cluster is assumed to live in a state of balanced evolution in
which the binary disruption rate is equal to the binary formation rate. Under
this assumption, the lifetime of a binary, or the timescale until a new binary
is formed, is simply

$\tau_{\rm bin}={m_{\rm ej}\over\dot{M}_{\rm BH}},$ (8)

where $m_{\rm ej}$ is the total mass ejected by the binary, and $\dot{M}_{\rm
BH}$ is the BH mass loss rate given by the cluster model.

We continue selecting new binaries and evolve them through binary-single
encounters until either all BHs have been ejected from the cluster, or until a
maximum integration time of $t=13\rm Gyr$ has passed.

Figure 1: Local distributions of primary BH mass and mass ratio for merging BH
binaries in our best-case scenario where clusters start with high densities,
$\rho_{\rm h,0}=10^{5}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}\ $, and the BHs are
all initialised with zero dimensionless spin parameter $\chi=0$. Top, middle
and bottom panels correspond to the delayed supernova mechanism, the rapid
supernova mechanism, and the BH mass distribution of Belczynski et al. (2008),
respectively. Black lines show the corresponding distributions when
hierarchical mergers are not included in the calculation, i.e., it is assumed
that a BH formed from a previous merger is always ejected from the cluster.
Thick lines show median values of the merger rate value, and thin lines the
corresponding $90\%$ confidence intervals. The green lines and hatched regions
show the median and corresponding $90\%$ confidence regions inferred from the
GWTC-3(Abbott et al., 2021a).

###  2.1 Cluster initial mass function, formation time, and metallicity

In order to generate predictions for BH binary mergers, we need a GC initial
mass function, and a model for how the formation rate and metallicity of
clusters evolve with redshift. These ingredients of our models are described
below. For this we follow the approach of Antonini & Gieles (2020a).

The cluster initial mass function is obtained by fitting an evolved Schechter
mass function to the observed GC mass function in the Milky Way today (Jordán
et al., 2007)

$\phi_{\rm cl}=A(M+\Delta)^{-2}\exp\left(-\frac{M+\Delta}{M_{\rm c}}\right)\
.$ (9)

This gives the posterior distribution for the parameters $M_{\rm c}$ and
$\Delta$. Adopting a simple model for cluster evaporation and mass loss due to
stellar evolution, the corresponding initial GC mass function is given by:

$\phi_{\rm cl,0}=2AM_{0}^{-2}\exp\left(-\frac{M_{0}}{2M_{\rm c}}\right).$ (10)

The corresponding fractional mass loss due to evaporation and stellar
evolution is

$K=\frac{\rho_{\rm GC0}}{\rho_{\rm GC}}=\frac{\int_{M_{\rm
lo}}^{\infty}\phi_{\rm cl,0}M_{0}{\rm d}M_{0}}{\int_{M_{\rm
lo}}^{\infty}\phi_{\rm cl}M{\rm d}M}\ =32.5^{+86.9}_{-17.7}\ \rm(90\%\ cred.\
intervals).$ (11)

The spread in $K$ provides an estimate of the uncertainty in the fractional
mass loss from cluster until today, given the 156 Milky Way GC masses. The
cluster mass formed per unit volume integrated over all times is (Antonini &
Gieles, 2020a)

$\displaystyle\rho_{\rm GC0}=2.4^{+2.3}_{-1.2}\times 10^{16}~{}{\rm
M}_{\odot}\,{\rm Gpc}^{-3}\ \rm(90\%\ cred.\ intervals).$ (12)

The large error bars here imply that $\rho_{\rm GC0}$ is uncertain by a factor
of $\simeq 2$. In the next sections we include this uncertainty as well as the
uncertainty on $K$ in the predictions for the merger rate.

We obtain the distribution of cluster formation times from the semi-analytical
galaxy formation model of El-Badry et al. (2019). The resulting cluster
formation history peaks at a redshift of $\sim 4$, which is earlier than the
peak in the cosmic star formation history (redshift $\sim 2$, Madau &
Dickinson 2014). We sample the formation redshift of our cluster models from
the total cosmic cluster formation rate given by the fiducial model of El-
Badry et al. (2019) and integrated over all halo masses. This corresponds to
the formation rate per comoving volume of their Fig. 8 with their parameters
$\beta_{\Gamma}=1$ and $\beta_{\eta}=1/3$, where $\beta_{\Gamma}$ sets the
dependence of the cluster formation efficiency on surface density, and
$\beta_{\eta}$ the dependence of the star formation rate on the halo virial
mass. Here, we renormalise the cluster formation rate, $\phi_{z}(z)$, such
that $\int_{\infty}^{0}\phi_{z}{\rm d}z=1$. Thus, we only sample the cluster
formation time from El-Badry et al. (2019), and then rescale the cluster
formation rate such to reproduce the total mass density given by our equation
(12). We note that the cluster formation model has a negligible impact on the
local merger rate and properties of the merging binaries. In Antonini & Gieles
(2020a) we showed that unrealistic models where all clusters are assumed to
form at the same time (e.g., $z=3$) produce similar results than models in
which the cluster formation rate is varied with redshift.

For the cluster metallicity, we fit a quadratic polynomial to the observed
age-metallicity relation for the Milky Way GCs (VandenBerg et al., 2013), to
obtain the mean metallicity

$\log(Z_{\rm mean}/Z_{\odot})\simeq 0.42+0.046\left({t\over{\rm
Gyr}}\right)-0.017\left({t\over{\rm Gyr}}\right)^{2}\ .$ (13)

Given the cluster formation redshift, we then assume a log-normal distribution
of metallicity around the mean

$\phi_{Z}=\frac{\log(e)}{\sqrt{2\pi\sigma^{2}}}\exp\Big{\\{}-\frac{[\log(Z/{\rm
Z_{\odot}})-Z_{\rm mean}]^{2}}{2\sigma^{2}}\Big{\\}},$ (14)

with standard deviation $\sigma=0.25$ dex. This takes into account the large
spread found in the observed age-metallicity relation for the Milky Way GCs.

###  2.2 Merger rates and their error bars

Finally, we construct a library of cluster models over a 3-dimensional grid of
formation time, metallicity and cluster mass.

We sample the cluster formation redshift over the range $z\in[10;0]$ with
step-size $\Delta z=0.5$; at a given redshift, the metallicity of the cluster
is sampled in the range $Z\in[10^{-4};0.02]$ with logarithmic step size
$\Delta\log Z=0.1$; finally, for a given formation time and metallicity, the
initial mass of the cluster is varied in the range $M_{\rm
0}\in[10^{2};2\times 10^{7}]M_{\odot}$, with step size $\Delta\log
M_{0}/M_{\odot}=0.1$. The merger rate is then calculated over the grid of
cluster models as:

$\mathcal{R}(\tau)=K\rho_{\rm
GC}{\sum\limits_{z}\sum\limits_{Z}\sum\limits_{M_{0}}\dot{\mathcal{N}}(\tau;z,Z,M_{0})\phi_{z}\phi_{Z}\phi_{\rm
cl,0}M_{0}\over\sum\limits_{z}\phi_{z}\sum\limits_{Z}\phi_{Z}\sum\limits_{M_{0}}\phi_{\rm
cl,0}M_{\rm 0}^{2}},$ (15)

where $\dot{\mathcal{N}}(\tau;M_{0},r_{\rm h,0},Z)$ is the BH binary merger
rate at a look-back time $\tau$ corresponding to a cluster with an initial
mass $M_{0}$, metallicity $Z$ and that formed at a redshift $z$.

In order to take into account the uncertainties in the initial cluster mass
function, we sample 100 values over the posterior distributions of the
parameters $M_{\rm c}$ and $\Delta$ obtained from the MCMC fit to the Milky
Way GC mass function. We also take into account the uncertainty on the mass
density of GCs in the Universe, $\rho_{\rm GC}$. We assume that the parameter
$\rho_{\rm GC}$ follows a Gaussian distribution with mean $7.3\times
10^{14}M_{\odot}\ \rm Gpc^{-3}$ and dispersion $\sigma=2.6\times
10^{14}M_{\odot}\ \rm Gpc^{-3}$. We sample 100 values from this Gaussian
distribution and for each of them we use equation (15) to determine a merger
rate estimate for each of the [$M_{\rm c}$, $\Delta$] values, and thus obtain
a distribution of merger rate density values. Since in this work we are
interested in the mass distribution of local BH binary mergers, we consider
the differential merger rate in the local universe ${\rm
d}\mathcal{R}(z=0)/{\rm d}m_{1}$ and ${\rm d}\mathcal{R}(z=0)/{\rm d}q$, which
we compare to the distributions inferred from GWTC-3.

Figure 2: Differential local merger rate as a function of the initial cluster
mass. We also show the initial cluster mass function (in arbitrary units) for
our best fit value of Schechter mass, $\log M_{\rm c}/M_{\odot}=6.26$. The
delayed supernova mechanism has been adopted here.

Figure 3: Dependence of mass and mass ratio distributions on initial cluster
half-mass density, and initial BH spins. The delayed supernova model is
assumed here. Top panels use $\chi=0$ and the half-mass density of the cluster
is varied as indicated. In the bottom panel we take $\rho_{\rm
h,0}=10^{5}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}$ and change the initial value of
$\chi$; the black histograms show the results for a model where the initial
value of $\chi$ is sampled using the inferred distribution of BH spins shown
in Fig. 15 of Abbott et al. (2021c). In all the other models, the BHs all form
with the same value of $\chi$ as indicated.

##  3 Results

###  3.1 Primary BH mass and mass ratio distributions

In Antonini & Gieles (2020a) we found good agreement between model predictions
and the inferred distribution of primary BH masses within the range of values
$m_{1}=15M_{\odot}$ to $40M_{\odot}$. Outside this range, the binary BH merger
rate was found to be several orders of magnitude smaller than inferred. The
first question we address here is whether the inclusion of hierarchical
mergers can reduce the discrepancy at $m_{1}\gtrsim 40M_{\odot}$ between
models and the inferred astrophysical distributions.

In Fig. 1 we plot the distributions of $m_{1}$ and $q$ for three different
assumptions about the initial BH mass function. In the upper panels we use the
delayed supernova mechanism, in the middle panels the rapid supernova
mechanism (Fryer et al., 2012), and in the lower panels we use the BH mass
distribution from Belczynski et al. (2008). These prescriptions produce
somewhat different initial BH mass functions, and lead also to different natal
kick values. In these models, all clusters are initialised with the same half-
mass radius density of $\rho_{\rm h,0}=10^{5}~{}{\rm M}_{\odot}\,{\rm
pc}^{-3}\ $ and the BHs are all started with zero dimensionless spin
parameter, $\chi=0.$

In the left panels of Fig. 1, we see that the new models produce mergers above
the $\sim 50M_{\odot}$ threshold. These mergers are produced by BHs that grow
hierarchically through mergers – the vast majority of them are mergers between
a first generation BH and a second generation BH. When we include these
mergers in our calculation, we find good agreement between the models and the
inferred distributions at $m_{1}\gtrsim 50M_{\odot}$, in the sense that the
models give a merger rate that is comparable to the inferred value. However,
we also note that a simple power-law profile above this mass threshold is not
a good representation of the model distributions. Above $m_{1}\simeq
50M_{\odot}$, the model distributions are characterised by several peaks. Such
higher mass peaks are related to peaks in the BH mass distribution at lower
masses. From Fig. 1 we see that the merger rate between first generation BHs
peaks at $\simeq 30M_{\odot}$ and $40M_{\odot}$. Thus, mergers between first
and second generation BHs lead to additional peaks at
$\simeq(30+30)M_{\odot}=60M_{\odot}$, $(30+40)M_{\odot}=70M_{\odot}$ and
$(40+40)M_{\odot}=80M_{\odot}$. The presence of peaks within the
pair–instability mass gap and their relation to lower mass peaks in the BH
mass distribution is a clear prediction of a hierarchical merger model for the
origin of the binaries.

The black histograms in the left panels of Fig. 1 show the results from models
in which any remnant BH formed from a previous merger is ejected from the
cluster ‘by hand’. In these models, the distributions are sharply truncated at
$\sim 50M_{\odot}$ since BHs cannot grow hierarchically above this mass value.
The merger rate at $m_{1}\simeq 10M_{\odot}$ derived from all models is about
two orders of magnitude smaller than the inferred rate. This lower-mass peak
can be explained, however, through other scenarios, including formation in the
galactic field (e.g., Broekgaarden et al., 2022; van Son et al., 2022) and
formation in young and open star clusters because of their higher metallicity
(e.g., Banerjee, 2021b; Chattopadhyay et al., 2022). On the other hand, our
models reproduce the inferred merger rate near $m_{1}\simeq 30M_{\odot}$,
which can therefore be explained by a GC origin. This peak in the mass
distribution is due to mergers involving first-generation BHs, and it is not
related to hierarchical mergers.

By comparing the results in Fig. 1 with the models in Antonini & Gieles
(2020a), we find that the latter generated a merger rate at $m_{1}\lesssim
20~{}M_{\odot}$ higher by a factor $\lesssim 2$. The reason for this
difference is due to the adopted new recipe for sampling the black hole binary
components and the interloper masses. In Antonini & Gieles (2020a) we had
assumed that $m_{1}=m_{2}=m_{3}$=$m_{\rm max}$, where $m_{\rm max}$ is the
mass of the most massive BH in the cluster. The distributions in Appendix A
mean instead that in the current models $\langle m_{3}\rangle\ll\langle
m_{1}\rangle\simeq\langle m_{2}\rangle$. Thus, each binary ejects more low-
mass BH interloopers lowering the overall BH merger rate at low masses.

In the right panels of Fig. 1 we consider the distribution of the mass-ratio
$q$. The new models result in a significantly higher rate of merging binaries
with small mass ratio, $q\lesssim 0.5$, providing a better match to the
inferred distribution than models without hierarchical mergers. Most of these
additional low-$q$ systems are mergers between a first generation BH and a BH
that formed through a previous merger above the pair-instability mass limit.
At high values of $q$, instead, both models with and without hierarchical
mergers produce a similar merger rate, which, at $q\gtrsim 0.8$, is about one-
order of magnitude smaller than inferred from the data.

In Fig. 2 we show the differential contribution to the local merger rate with
respect to the initial cluster mass. This allows us to identify in which type
of clusters most of the mergers are formed. The contribution to the total
merger rate is nearly constant in log bins between $M_{0}=5\times
10^{4}M_{\odot}$ and $5\times 10^{6}M_{\odot}$, while it decrease
exponentially above this range because of the truncation of the initial GC
mass function at $\log M_{\rm c}/M_{\odot}\simeq 6.26$. We then show the same
cluster mass distribution, but only considering mergers with $m_{1}\geq
50M_{\odot}$ (red histogram). These mergers, involving a primary BH above the
mass gap limit, are mostly produced in clusters with relatively large masses,
between $\sim 2\times 10^{5}M_{\odot}$ and $5\times 10^{6}M_{\odot}$. The
fraction of these higher mass mergers to the total number of mergers in each
cluster mass bin increases with cluster mass. By comparing the blue and the
red histograms in the figure, we see that at $M_{0}\sim 10^{7}M_{\odot}$,
between $10$ to $30\%$ of mergers have a primary with $m_{1}>50M_{\odot}$. The
percentage goes down to $\sim 1\%$ in clusters with an initial mass lower than
$\simeq 10^{6}M_{\odot}$. Finally, we show the $M_{0}$ distributions for the
most massive mergers produced in our models, $m_{1}\geq 100M_{\odot}$. These
BHs originate from at least two previous mergers since their mass is larger
than twice the initial mass cut-off at $\simeq 50M_{\odot}$. These systems are
formed in the most massive GCs, with initial mass well above the Schechter
mass.

Based on the models of Fig. 2, we compute a local binary BH merger rate of
$4.1^{+2.2}_{-2.5}\rm Gpc^{-3}~{}yr^{-1}$ (delayed), $4.5^{+2.7}_{-2.9}\rm
Gpc^{-3}~{}yr^{-1}$ (rapid), and $6.0^{+3.6}_{-4.0}\rm Gpc^{-3}~{}yr^{-1}$ (BH
mass distribution from Belczynski et al., 2008), at $90\%$ confidence. The
binary BH merger rate inferred from the gravitational wave data is estimated
to be between $17.9\rm Gpc^{-3}~{}yr^{-1}$ and $44\rm Gpc^{-3}~{}yr^{-1}$
(Abbott et al., 2021a), and it is therefore a factor of $\simeq 2$ to $20$
larger than the rate computed from our models.

Figure 4:  Results for six alternative models. Left panels: blue histograms
are for a model in which the initial cluster half-mass radius is assumed to
scale with cluster mass as in equation (16); the blue histograms correspond to
a model in which the BH birth kicks are zero; the red histograms correspond to
a model where the recipes for pulsation pair instability (PPI) were switched
off. Right panels: blue and black histograms are the results obtained assuming
that the initial stellar mass function for massive stars scales as
$\phi(m_{\star})\propto m_{\star}^{-2}$ and $\phi(m_{\star})\propto
m_{\star}^{-2.6}$, respectively. In the red and orange histograms we have
multiplied our standard wind mass loss rate on the Wolf-Rayet stage by a
factor $f_{\rm WR}=0.1$ and $f_{\rm WR}=5$. We have used the delayed supernova
mechanism, and, unless otherwise specified, all the other model parameters are
the same as in Fig. 1.

###  3.2 Effect of initial cluster density, initial spins, and other model
variations

The number of heavier BHs produced by a cluster through hierarchical mergers
is affected by both the cluster density and the initial spin of the BHs. A
larger cluster density means a larger merger rate and escape velocity, and
therefore a larger probability that a remnant BH is retained inside the
cluster following a recoil kick. Similarly, if the BHs have negligible spins,
this translates into a smaller recoil velocity and higher retention
probability. While in the previous section we have looked at a somewhat
optimistic scenario in which clusters all form with high densities and the BHs
have initially zero spins, in this section we vary these assumptions and
investigate their effect on the BH binary merger rate and properties. We adopt
here the delayed supernova mechanism, but similar results are obtained with
the rapid supernova prescription and the Belczynski et al. (2008) mass
distribution.

In the upper panels of Fig. 3, we vary the initial cluster half-mass density
within the range $\rho_{\rm h,0}=10^{2}$ to $10^{5}~{}{\rm M}_{\odot}\,{\rm
pc}^{-3}\ $, and assume that the BHs have zero spins initially. The results
illustrate that although our models can in principle account for most mergers
above $m_{1}\gtrsim 20M_{\odot}$, this is only true under some specific
conditions. As we lower the initial cluster half-mass density the merger rate
goes down significantly at all values of mass and mass-ratio. The depletion is
more significant at masses above the cut-off mass of $50M_{\odot}$ and for
$q\lesssim 0.5$. Thus, a scenario where most merging BH binaries with
$m_{1}\gtrsim 20M_{\odot}$ form in GCs would imply a typical initial cluster
density $\rho_{0}\gtrsim 10^{4}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}$. It is
important to note that this condition would however only apply to clusters
with initial mass $M_{0}\gtrsim 5\times 10^{4}M_{\odot}$, where most of the
merging binaries are formed (see Fig. 2).

In the lower panels of Fig. 3 we show how the results change with changing the
initial BH spins. In these models we keep the initial density to the fixed
value $\rho_{0}=10^{5}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}$. We see that the
merger rate density distributions are not affected significantly for
$m_{1}\lesssim 50M_{\odot}$ and $q\gtrsim 0.5$. This is because the majority
of these binaries are made of first generation BHs. Hence, their merger rate
is not affected by the recoil kick velocity and by the initial choice of BH
spin. On the other hand, the number of BHs formed via hierarchical mergers
decreases significantly when higher initial spins are used due to the larger
recoil kicks. This leads to a lower merger rate at $m_{1}\gtrsim 50M_{\odot}$
when $\chi$ is increased. Even relatively modest initial spins, $\chi=0.1$,
lead to a distribution that does no longer match the inferred distribution.
The constrains on $\chi$ seems therefore quite strong as a hierarchical origin
for all mergers with $m_{1}\gtrsim 50M_{\odot}$ would require that BHs are
formed with nearly zero spin.

Finally, we consider six additional model realisations. In one model, we
assume that the BHs receive no kick at formation and that the initial density
is the same for all clusters, $\rho_{\rm h,0}=10^{5}~{}{\rm M}_{\odot}\,{\rm
pc}^{-3}\ $. In another model, we assume that the cluster half-mass radius
scales as

$\log\left({{r}_{\rm h,0}\over{\rm
pc}}\right)=-3.56+0.615\log\left({M_{0}\over M_{\odot}}\right)\ .$ (16)

This latter relation was derived by Gieles et al. (2010) from the results of
Hasegan et al. (2005) who fit this Faber-Jackson-like relation to ultra-
compact dwarf galaxies (UCDs) and elliptical galaxies. Gieles et al. (2010)
derived the initial mass-radius relation correcting for mass loss and
expansion by stellar evolution and correcting radii for projection. We
consider an additional model realisation where we did not include any
prescription for pair instability so that the initial BH mass function has no
upper gap and BHs can form above $50M_{\odot}$. Moreover, we consider two
models where the initial mass function above $0.5\ M_{\odot}$ is assumed to
scale as $\phi(m_{\star})\propto m_{\star}^{-2}$ (top-heavy) and
$\phi(m_{\star})\propto m_{\star}^{-2.6}$ (bottom heavy), respectively.
Finally, we evolve two additional models where our standard Wolf-Rayet winds
based on Hamann & Koesterke (1998) and Vink & de Koter (2005) are multiplied
by a factor $f_{\rm WR}=0.1$ and $f_{\rm WR}=5$ (e.g., Broekgaarden et al.,
2022). Unless otherwise specified, all the other model parameters are the same
as before, i.e., delayed supernova mechanism, $\chi=0$, fallback kicks, etc.

Fig. 4 shows that the mass properties of the BH binaries produced in the new
models without birth kicks and with the new $r_{\rm h}$-$M$ relation are
similar to those found previously in Section 3.1. The fact that adopting the
mass-radius relation equation (16) does not change significantly our results
is not surprising. Clusters with an initial mass $M_{0}\sim 10^{6}M_{\odot}$
contribute the most to the merger rate (see Fig.2). The initial half-mass
density of these clusters as derived from equation equation (16) is $\rho_{\rm
h,0}\simeq 5\times 10^{4}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}\ $. This is
comparable to the constant density value of $10^{5}~{}{\rm M}_{\odot}\,{\rm
pc}^{-3}$ adopted previously. Interestingly, the results of models with no
birth kicks show that assuming zero velocity kicks at birth increases slightly
the merger rate at the lower mass end of the $m_{1}$ distribution and the
number of merging binaries with asymmetric masses. On the other hand, the
shape and normalisation of the distributions at masses higher than
$m_{1}\gtrsim 20M_{\odot}$ remain virtually the same as in the fallback kick
model.

The model without pair instability physics leads to a mass distribution which
is significantly different from the other model realisations, showing how our
results can depend on the assumptions about stellar evolution and the adopted
prescriptions. In this case, the mass distribution still peaks at $m_{1}\sim
30~{}M_{\odot}$, while the other peak near $40~{}M_{\odot}$ is not longer
present. A secondary peak is found near $70M_{\odot}$. For masses larger than
this value, the merger rate drops and becomes much smaller than the rate
inferred from the GW data.

In the model with a top-heavy stellar mass function, the overall merger rate
is higher than for our standard models due to the larger number of BHs formed.
On the other hand, for a bottom-heavy mass function the total merger rate is
significantly reduced due to the fewer massive stars formed. Our model with
modified Wolf-Rayet wind mass-loss rate lead to results that are qualitatively
similar to those obtained under our more standard assumptions.

That our resulting mass distributions are sensitive to the initial BH mass
function, and therefore to the uncertain stellar evolution prescriptions is
expected. It is interesting, however, that most of our models share similar
properties. Specifically: (i) the inferred peak in the merger rate at
$10M_{\odot}$ is much lower than the one inferred from the data, and (ii) the
distribution of $m_{1}$ presents a main peak at near $35M_{\odot}$. The main
reason why there are so few mergers with small masses is because of the
relatively low number of light BHs in the initial mass function. This is due
to the low metallicity of GCs, which results in low wind mass loss and large
BH masses. The other reason why the mass distribution of merging binaries
peaks at relatively high values is dynamics. The masses of the binary
components tend to be sampled near the top end of the BH mass function, due to
the high value of the power law exponents that appear in the density
probability functions $p_{1}$ and $p_{2}$ (see Section A). On the other hand,
the flatter $p_{3}$ distribution means that the ejected BH interlopers will be
on average lighter than the binary BH components. These lighter BHs are
therefore no longer available for merging.

##  4 conclusions

In this work we have used our fast cluster evolution code, cBHBd, to
investigate the mass distribution of BH binaries produced dynamically in dense
GCs. We compared our results to the astrophysical distribution of BH binary
masses inferred from GWTC-3 to make inference about the astrophysical origin
of the sources. For the first time, we have included hierarchical mergers in
our models. This allowed us to address the question of whether a dynamical
formation scenario is a feasible explanation for the detected BHs within the
so called ‘upper mass gap’. Such a mass gap in the initial BH mass function is
predicted by stellar evolution theories, and in our models is located at
$\gtrsim 50M_{\odot}$. Because cBHBd is highly efficient compared to other
techniques (e.g., Monte Carlo, $N$-body), we were able to systematically
investigate the impact of model assumptions on or results. Our main
conclusions are summarised below:

  * i)

A purely GC formation scenario for the BH binaries detected by LIGO and Virgo
is inconsistent with the $\simeq 10M_{\odot}$ peak in the primary BH mass
distribution that is inferred from the data. This likely excludes a scenario
where the majority of the sources were formed in GCs.

  * ii)

A GC origin can easily account for the secondary mass peak at $m_{1}\simeq
35M_{\odot}$ inferred from the data. This requires that clusters form with
initial half-mass density $\gtrsim 10^{4}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}\ $.
Assumptions about the initial BH spins and the supernova mechanism have no
effect on this conclusion.

  * iii)

Dynamical formation in GCs can explain the inferred merger rate of all BH
binaries with $m_{1}\gtrsim 20M_{\odot}$ and $q\lesssim 0.8$, including
binaries with component masses lying above the assumed mass limit due to
pair–instability. For this to be true we require that both the most massive
GCs, $M_{0}\gtrsim 10^{5}M_{\odot}$, form with half-mass density $\gtrsim
10^{4}~{}{\rm M}_{\odot}\,{\rm pc}^{-3}\ $, and that the birth spins of BHs
are nearly zero. Even small deviations from this latter condition lead to a
merger rate above $50M_{\odot}$ that is orders of magnitude smaller than the
inferred rate.

  * iv)

A hierarchical merger scenario predicts the appearance of multiple peaks in
the primary BH mass distribution and within the upper mass gap due to a
pile–up of mergers between first and second generation BHs. Inter–generation
mergers lead to a simple relation between the mass value of any of such peaks
and that of peaks found at masses lower than the pair-instability limit. These
features can be tested against future GW data to place constrains on a GC
origin for the sources.

Additional constraints on the formation of BH mergers can be placed by
exploring correlations between binary parameters, which we have not considered
here, but plan to study in a future work. For example, a BH formed from a
previous merger will have a spin $\chi\simeq 0.7$. We expect therefore a
change in the value of the typical effective and precession spin parameters of
binaries with components within the upper mass gap (e.g., Tagawa et al.,
2021a; Baibhav et al., 2020) and an increase in spin magnitude for systems
with more unequal mass ratio. Binaries formed dynamically will also have
larger eccentricities, which can lead to a positive correlation between
eccentricity, spin and binary mass in the overall population. The analysis of
the data from GWTC-3 has shown marginal evidence that the spin distribution
broadens above $30M_{\odot}$, and that the mass ratio and spin are correlated
in the sense that spins are larger for more asymmetric binaries (Abbott et
al., 2021a). The evidence for these correlations remain weak, but it suggests
that future analysis based on larger data sets will soon be able to provide
more stringent constrains. The residual eccentricity of a binary is by itself
another potentially powerful tool for identifying sources formed in clusters.
Romero-Shaw et al. (2022) suggest that a significant fraction of the detected
GW sources in GWTC-3 show support for eccentricity $\gtrsim 0.1$ at $10$Hz.
Their results indicate that densely–populated star clusters may produce the
majority of the observed mergers.

Finally, it is worth noting that in our work we used the pair–instability
prescriptions from Spera & Mapelli (2017). This gives an upper limit in the
initial BH mass function of about $50M_{\odot}$. However, there are several
uncertainties in the modelling, and different assumptions can lead to
significantly different values for the high mass cut-off, generally in the
range $40M_{\odot}$ to $70M_{\odot}$ (Giacobbo et al., 2018; Farmer et al.,
2020; Costa et al., 2021; Fryer et al., 2022). Exploring the effect of these
uncertainties is beyond the scope of this paper, but should be considered in
future work.

## Acknowledgements

Some of the processing of the results has been done using the python
programming language and the following open source modules:
numpy111http://www.numpy.org, scipy222http://www.scipy.org,
matplotlib333http://matplotlib.sourceforge.net. We acknowledge the support of
the Supercomputing Wales project, which is part-funded by the European
Regional Development Fund (ERDF) via Welsh Government. FA acknowledges support
from a Rutherford fellowship (ST/P00492X/1) from the Science and Technology
Facilities Council. MG acknowledges support from the Ministry of Science and
Innovation (EUR2020-112157, PID2021-125485NB-C22, CEX2019-000918-M funded by
MCIN/AEI/10.13039/501100011033) and from AGAUR (SGR-2021-01069).

## Data availability

The data underlying this article will be shared on a reasonable request to the
corresponding author.

## References

  * Abbott (2020a) Abbott R. e. a., 2020a
  * Abbott et al. (2021a) Abbott et al., 2021a, arXiv e-prints, arXiv:2111.03606
  * Abbott et al. (2021b) Abbott et al., 2021b, ApJ, 913, L7
  * Abbott et al. (2021c) Abbott et al., 2021c, arXiv e-prints, arXiv:2111.03634
  * Antonini & Gieles (2020a) Antonini F., Gieles M., 2020a, Phys. Rev. D, 102, 123016
  * Antonini & Gieles (2020b) Antonini F., Gieles M., 2020b, MNRAS, 492, 2936
  * Antonini et al. (2019) Antonini F., Gieles M., Gualandris A., 2019, MNRAS, 486, 5008
  * Antonini & Rasio (2016) Antonini F., Rasio F. A., 2016, ApJ, 831, 187
  * Antonini et al. (2017) Antonini F., Toonen S., Hamers A. S., 2017, ApJ, 841, 77
  * Askar et al. (2017) Askar A., Szkudlarek M., Gondek-Rosińska D., Giersz M., Bulik T., 2017, MNRAS, 464, L36
  * Baibhav et al. (2020) Baibhav V., Gerosa D., Berti E., Wong K. W. K., Helfer T., Mould M., 2020, Phys. Rev. D, 102, 043002
  * Ballesteros et al. (2018) Ballesteros G., Serpico P. D., Taoso M., 2018, J. Cosmology Astropart. Phys., 2018, 043
  * Banerjee (2021a) Banerjee S., 2021a, MNRAS, 500, 3002
  * Banerjee (2021b) Banerjee S., 2021b, MNRAS, 503, 3371
  * Banerjee et al. (2010) Banerjee S., Baumgardt H., Kroupa P., 2010, Monthly Notices of the Royal Astronomical Society, 402, 371
  * Bartos et al. (2016) Bartos I., Kocsis B., Haiman Z., Márka S., 2016, 1
  * Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, The Astrophysical Journal Supplement Series, 174, 223
  * Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, Monthly Notices of the Royal Astronomical Society, 432, 2779
  * Broekgaarden et al. (2022) Broekgaarden F. S. et al., 2022, MNRAS
  * Chattopadhyay et al. (2022) Chattopadhyay D., Hurley J., Stevenson S., Raidani A., 2022, MNRAS, 513, 4527
  * Costa et al. (2021) Costa G., Bressan A., Mapelli M., Marigo P., Iorio G., Spera M., 2021, MNRAS, 501, 4514
  * Dall’Amico et al. (2021) Dall’Amico M., Mapelli M., Di Carlo U. N., Bouffanais Y., Rastello S., Santoliquido F., Ballone A., Arca Sedda M., 2021, MNRAS, 508, 3045
  * De Luca et al. (2021) De Luca V., Desjacques V., Franciolini G., Pani P., Riotto A., 2021, Phys. Rev. Lett., 126, 051101
  * Di Carlo et al. (2020) Di Carlo U. N., Mapelli M., Bouffanais Y., Giacobbo N., Santoliquido F., Bressan A., Spera M., Haardt F., 2020, MNRAS, 497, 1043
  * Doctor et al. (2020) Doctor Z., Wysocki D., O’Shaughnessy R., Holz D. E., Farr B., 2020, ApJ, 893, 35
  * Dominik et al. (2012) Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2012, DOUBLE COMPACT OBJECTS. I. THE SIGNIFICANCE OF THE COMMON ENVELOPE ON MERGER RATES
  * El-Badry et al. (2019) El-Badry K., Quataert E., Weisz D. R., Choksi N., Boylan-Kolchin M., 2019, MNRAS, 482, 4528
  * Farmer et al. (2020) Farmer R., Renzo M., de Mink S. E., Fishbach M., Justham S., 2020, ApJ, 902, L36
  * Fishbach et al. (2017) Fishbach M., Holz D. E., Farr B., 2017, ApJ, 840, L24
  * Fragione & Kocsis (2018) Fragione G., Kocsis B., 2018, Phys. Rev. Lett., 121, 161103
  * Fragione et al. (2022) Fragione G., Loeb A., Kocsis B., Rasio F. A., 2022, arXiv e-prints, arXiv:2204.03745
  * Fragione et al. (2020) Fragione G., Loeb A., Rasio F. A., 2020, ApJ, 895, L15
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, The Astrophysical Journal, 749, 91
  * Fryer et al. (2022) Fryer C. L., Olejak A., Belczynski K., 2022, ApJ, 931, 94
  * Gerosa & Berti (2017) Gerosa D., Berti E., 2017, Phys. Rev. D, 95, 124046
  * Giacobbo et al. (2018) Giacobbo N., Mapelli M., Spera M., 2018, MNRAS, 474, 2959
  * Gieles et al. (2010) Gieles M., Baumgardt H., Heggie D. C., Lamers H. J. G. L. M., 2010, MNRAS, 408, L16
  * Gow et al. (2020) Gow A. D., Byrnes C. T., Hall A., Peacock J. A., 2020, J. Cosmology Astropart. Phys., 2020, 031
  * Hamann & Koesterke (1998) Hamann W. R., Koesterke L., 1998, A&A, 335, 1003
  * Hamers et al. (2021) Hamers A. S., Fragione G., Neunteufel P., Kocsis B., 2021, MNRAS, 506, 5345
  * Hasegan et al. (2005) Hasegan M. et al., 2005, ApJ, 627, 203
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Hénon (1975) Hénon M., 1975, Dynamics of the Solar Systems
  * Hills & Day (1976) Hills J. G., Day C. A., 1976, Astrophys. Lett., 17, 87
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, Monthly Notices of the Royal Astronomical Society, 360, 974
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Jordán et al. (2007) Jordán A. et al., 2007, ApJS, 171, 101
  * Kimball et al. (2020) Kimball C., Talbot C., Berry C. P. L., Carney M., Zevin M., Thrane E., Kalogera V., 2020, ApJ, 900, 177
  * Kimball et al. (2021) Kimball C. et al., 2021, ApJ, 915, L35
  * Kremer et al. (2020) Kremer K. et al., 2020, ApJ, 903, 45
  * Kroupa (2001) Kroupa P., 2001, Monthly Notices of the Royal Astronomical Society, 322, 231
  * Kulkarni et al. (1993) Kulkarni S. R., Hut P., McMillan S. J., 1993, Nature, 364, 421
  * Liu & Lai (2021) Liu B., Lai D., 2021, MNRAS, 502, 2049
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Mandel & De Mink (2016) Mandel I., De Mink S. E., 2016, MNRAS, 458, 2634
  * Mapelli et al. (2022) Mapelli M., Bouffanais Y., Santoliquido F., Arca Sedda M., Artale M. C., 2022, MNRAS, 511, 5797
  * Miller & Lauburg (2009) Miller M. C., Lauburg V. M., 2009, ApJ, 692, 917
  * Mink & Belczynski (2015) Mink S. E. d., Belczynski K., 2015, The Astrophysical Journal, 814, 58
  * Mould et al. (2022) Mould M., Gerosa D., Taylor S. R., 2022, arXiv:2203.03651, arXiv:2203.03651
  * O’Leary et al. (2009) O’Leary R. M., Kocsis B., Loeb A., 2009, Monthly Notices of the Royal Astronomical Society, 395, 2127
  * Olejak et al. (2022) Olejak A., Fryer C. L., Belczynski K., Baibhav V., 2022, arXiv e-prints, arXiv:2204.09061
  * Park et al. (2017) Park D., Kim C., Lee H. M., Bae Y.-B., Belczynski K., 2017, MNRAS, 469, 4665
  * Peters (1964) Peters P., 1964, Physical Review, 136, B1224
  * Portegies Zwart et al. (2007) Portegies Zwart S., Gaburov E., Chen H.-C., Gürkan M. A., 2007, MNRAS, 378, L29
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Rezzolla et al. (2008) Rezzolla L., Barausse E., Dorband E. N., Pollney D., Reisswig C., Seiler J., Husa S., 2008, Phys. Rev. D, 78, 044002
  * Rodriguez et al. (2018) Rodriguez C. L., Amaro-Seoane P., Chatterjee S., Rasio F. A., 2018, Phys. Rev. Lett., 120, 151101
  * Rodriguez et al. (2016) Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Physical Review D, 93, 084029
  * Rodriguez et al. (2020) Rodriguez C. L. et al., 2020, ApJ, 896, L10
  * Rodriguez et al. (2015) Rodriguez C. L., Morscher M., Pattabiraman B., Chatterjee S., Haster C.-J., Rasio F. A., 2015, Physical Review Letters, 115, 051101
  * Romero-Shaw et al. (2022) Romero-Shaw I. M., Lasky P. D., Thrane E., 2022, arXiv e-prints, arXiv:2206.14695
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Samsing (2018) Samsing J., 2018, Phys. Rev. D, 97, 103014
  * Sigurdsson & Hernquist (1993) Sigurdsson S., Hernquist L., 1993, Nature, 364, 423
  * Silsbee & Tremaine (2017) Silsbee K., Tremaine S., 2017, The Astrophysical Journal, 836, 39
  * Spera & Mapelli (2017) Spera M., Mapelli M., 2017, MNRAS, 470, 4739
  * Stegmann et al. (2022) Stegmann J., Antonini F., Schneider F. R. N., Tiwari V., Chattopadhyay D., 2022, arXiv e-prints, arXiv:2203.16544
  * Stone et al. (2016) Stone N. C., Metzger B. D., Haiman Z., 2016
  * Tagawa et al. (2021a) Tagawa H., Haiman Z., Bartos I., Kocsis B., Omukai K., 2021a, MNRAS, 507, 3362
  * Tagawa et al. (2021b) Tagawa H., Kocsis B., Haiman Z., Bartos I., Omukai K., Samsing J., 2021b, ApJ, 908, 194
  * Tiwari & Fairhurst (2021) Tiwari V., Fairhurst S., 2021, ApJ, 913, L19
  * van Son et al. (2022) van Son L. A. C. et al., 2022, arXiv:2209.13609, arXiv:2209.13609
  * VandenBerg et al. (2013) VandenBerg D. A., Brogaard K., Leaman R., Casagrande L., 2013, The Astrophysical Journal, 775, 134
  * Vink & de Koter (2005) Vink J. S., de Koter A., 2005, A&A, 442, 587
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, Astronomy and Astrophysics, 369, 574
  * Woosley (2017) Woosley S. E., 2017, The Astrophysical Journal, 836, 244
  * Yang et al. (2019) Yang Y. et al., 2019, Phys. Rev. Lett., 123, 181101
  * Zevin & Holz (2022) Zevin M., Holz D. E., 2022, arXiv e-prints, arXiv:2205.08549

Figure 5: PDFs for $m_{1}$ (left) and $q$ (right) for different BH mass
function upper masses ($m_{\rm up}$) and logarithmic slopes ($\alpha$). In all
cases $m_{\rm lo}=5$. The results discussed in the text (equations 20 and 22)
are shown as full lines, while the dots with errors bars show Monte Carlo
realisations (as a check) obtained by drawing $10^{6}$ pairs from
$p_{I}(m_{I})$ (equation 19) and the black dashed lines show simple power-law
approximations.

##  Appendix A Mass sampling routines

###  A.1 Masses of 3-body binaries

We are interested in the probability density functions (PDFs) of the masses of
the two components of (BH) binaries that form in three-body interactions.
Following Heggie (1975), the formation rate of hard binaries per unit of
volume and energy is expressed as

$\Gamma_{\rm
3b}(m_{I},m_{II},m_{III},x)=n_{I}n_{II}n_{III}Q(m_{I},m_{II},m_{III},x),$ (17)

where $n_{i}=n(m_{i})$ is the number density of BHs with mass $m_{i}$, $m_{I}$
and $m_{II}$ are the masses of the stars ending up in the binary, $m_{III}$ is
the mass of the catalyst star and $x$ is the (positive) binding energy of the
binary. The rate function $Q$ is given by equation 4.14 in Heggie (1975) and
it is a function of the three masses, as well as their
$\beta_{i}=(m_{i}\sigma_{i}^{2})^{-1}$, where $\sigma_{i}$ is the one-
dimensional velocity dispersion of mass component $i$. Because binaries tend
to form from the most massive objects, for which energy equipartition is
established quickly, we assume $\beta_{I}=\beta_{II}=\beta_{III}=\beta$.
Integrating equation 4.14 from Heggie over all $x$, from the hard-soft
boundary (i.e. $x=\beta^{-1}$) to $\infty$ (i.e. all hard binaries), we find
that the formation rate of hard binaries per unit of volume is

$\Gamma_{\rm 3b}(m_{I},m_{II},m_{III})\propto
n_{I}n_{II}n_{III}\frac{(m_{I}m_{II})^{4}m_{III}^{5/2}}{\sqrt{(m_{I}+m_{II}+m_{III})(m_{I}+m_{II})}}\beta^{9/2}.$
(18)

For equal masses, this result reduces to the frequently used scaling
$\Gamma_{\rm 3b}(m)\propto n^{3}m^{5}\sigma^{-9}$. Equation (18) is symmetric
in $m_{I}$ and $m_{II}$, so the PDF for the mass of one of them is found from

$p_{I}(m_{I})=\int_{m_{\rm lo}}^{m_{\rm up}}\int_{m_{\rm lo}}^{m_{\rm up}}{\rm
d}m_{III}{\rm d}m_{II}\,\Gamma_{\rm 3b},$ (19)

and $p_{II}(m_{II})=p_{I}(m_{I})$. Here $m_{\rm lo}$ and $m_{\rm up}$ are the
lower and upper bound of the mass distribution, respectively.

We note here that Equation 18 includes the assumption of equipartition and
therefore takes into account the dependence of the velocity dispersion on
mass. On the other hand, we do not consider the change in the BH mass function
that is expected in the core due to mass segregation. Numerical simulations
have shown that the mass function in the core has a logarithmic slope that is
approximately only $+1$ steeper than the global mass function (e.g., Portegies
Zwart et al., 2007). Thus we expect the effect of mass segregation on our
results to be relatively small.

We now adopt the convention that $m_{1}$ and $m_{2}$ are the most massive and
least massive component, respectively, with corresponding PDFs $p_{1}(m_{1})$
and $p_{2}(m_{2})$. These correspond to the PDFs of the maximum and minimum
value, respectively, when a sample of two values are drawn from
$p_{I}(m_{I})$, and are given by

$\displaystyle p_{1}(m_{1})$ $\displaystyle=2P_{I}(m_{1})p_{I}(m_{1}),$ (20)
$\displaystyle p_{2}(m_{2})$
$\displaystyle=2\left[1-P_{I}(m_{2})\right]p_{I}(m_{2}),$ (21)

where $P_{I}(m_{I})=\int_{m_{\rm lo}}^{m_{I}}{\rm
d}m_{I}^{\prime}\,p_{I}(m_{I}^{\prime})$ is the cumulative density function
(CDF) of $p_{I}(m_{I})$.

The PDF of $q$ is a ratio distribution and can be found from the joint
distribution of the minimum and maximum values, which is given by
$p_{12}(m_{1},m_{2})=2p_{1}(m_{1})p_{2}(m_{2})$ and

$\displaystyle p_{q}(q)=\int_{m_{\rm lo}}^{m_{\rm up}}{\rm
d}m_{2}\,p_{12}(qm_{2},m_{2}).$ (22)

We then assume that the mass function is a power law such that $n_{i}\propto
m_{i}^{\alpha}$ between $m_{\rm lo}$ and $m_{\rm up}$. A value of $\alpha=0.5$
provides a good approximation of the mass function of BHs at low metallicities
($Z\lesssim 0.05Z_{\odot}$, see Fig. 4 of Antonini & Gieles, 2020b). In Fig. 5
we show the resulting $p_{1}(m_{1})$ and $p_{q}(q)$ for $m_{\rm
up}=[10,20,50]$ and $m_{\rm lo}=5$ and $\alpha=+0.5$ (approximate for metal-
poor GCs) and $\alpha=-1.5$ (approximate for metal-rich GCs). We find that
these PDFs can be reasonably well approximated by power-laws of the form:
$p_{1}(m_{1})\propto m_{1}^{8+2\alpha}$ and $p_{q}(q)\propto q^{3.5+\alpha}$,
for all values of $m_{\rm up}$ and $\alpha$. This approximation can be used to
sample $m_{1}$ and $m_{2}$ (via $q$).

###  A.2 Masses of interlopers

Assume a binary BH with mass $M_{12}=m_{1}+m_{2}$, moving in a field of BHs
with number density $n_{3}$. The rate of encounters between the BH binaries
and field BHs is (Hills & Day, 1976)

$\dot{N}_{3}=n_{3}\langle\Sigma v\rangle,$ (23)

where $v$ is the relative velocity between the binary BH and another BH and
$\Sigma$ is the cross section for an encounter, which in the gravitational
focusing regime is (Hills & Day, 1976)

$\displaystyle\Sigma\simeq\frac{2\pi GaM_{123}}{v^{2}},$ (24)

where $G$ is the gravitational constant, $M_{123}=M_{12}+m_{3}$ and $a$ is the
semi-major axis of the binary. We can find $\langle\Sigma v\rangle$ from
integrating over all velocities

$\displaystyle\langle\Sigma v\rangle$
$\displaystyle=\frac{4l^{3}}{\pi^{1/2}}\int_{0}^{\infty}\Sigma(v)v^{3}\exp(-l^{2}v^{2}){\rm
d}v,$ (25) $\displaystyle=4\pi^{1/2}GlM_{123}a.$ (26)

Here $l^{2}=\beta M_{12}m_{3}/(2M_{123})$ for our assumption of equipartition.
The semi-major axis is $a\propto Gm_{1}m_{2}\beta$ such that the interaction
rate scales with the masses as

$\dot{N}_{3}\propto m_{3}^{\alpha}\sqrt{m_{3}}\sqrt{M_{12}M_{123}}m_{1}m_{2}$
(27)

So, interactions with more massive BHs are slightly favoured wrt random draws
from the BH mass function. Because $p_{1}(m_{1})$ and $p_{2}(m_{2})$ are much
steeper than this distribution, we find that to good approximation
$p_{3}(m_{3})\propto m_{3}^{\alpha+1/2}$. It also means that exchanges are not
very important, because these happen when the intruder is more massive than
any of the binary members. Here we find for $m_{\rm up}/m_{\rm lo}=10$ and
$\alpha=+0.5$ that $\langle m_{1}\rangle\simeq 0.91m_{\rm up};\langle
m_{2}\rangle\simeq 0.76m_{\rm up}$ and $\langle m_{3}\rangle\simeq 0.55m_{\rm
up}$ and exchange interactions should therefore not be very common. Once the
width of the BH mass function has shrunk to $m_{\rm up}/m_{\rm lo}\simeq 2$,
$\langle m_{2}\rangle\simeq 0.78$ which is comparable to $\langle
m_{3}\rangle\simeq 0.75$ and exchange interactions (which we neglect) are more
relevant.