# Gravitational wave background from primordial black holes in globular
clusters

Eleonora Vanzan  Sarah Libanore  Lorenzo Valbusa Dall’Armi  Nicola Bellomo
Alvise Raccanelli

###### Abstract

Primordial black holes still represent a viable candidate for a significant
fraction, if not for the totality, of dark matter. If these compact objects
have masses of order tens of solar masses, their coalescence can be observed
by current and future ground-based gravitational wave detectors. Therefore,
finding new gravitational wave signatures associated with this dark matter
candidate can either lead to their detection or help constraining their
abundance. In this work we consider the phenomenology of primordial black
holes in dense environments, in particular globular clusters. We model the
internal structure of globular clusters in a semi-analytical fashion, and we
derive the expected merger rate. We show that, if primordial black holes are
present in globular clusters, their contribution to the GW background can be
comparable to other well-known channels, such as early- and late-time
binaries, thus enhancing the detectability prospects of primordial black holes
and demonstrating that this contribution needs to be taken into account.

##  1 Introduction

Gravitational wave (GW) astronomy opened the door to the exploration of the
Universe through gravitational radiation [1, 2]. With already
$\mathcal{O}(100)$ GW events detected from the LIGO-Virgo-KAGRA Collaboration
[3, 4, 5], we started constraining both astrophysical properties of the
sources [6, 7] and cosmological properties of the Universe [8, 9, 10]. For
instance, a single binary neutron star event with an electromagnetic
counterpart has been remarkably effective in ruling out entire classes of
modified gravity theories [11, 12, 13]. Given the great potential of the next
observational runs, it is paramount to find novel ways to exploit this new
probe for cosmological purposes, especially because GWs have the potential to
answer some of the open questions in cosmology.

One very exciting possibility is related to shedding light on the nature of
dark matter [14]. So far, many viable candidates have been proposed, see,
e.g., refs. [15, 16], including Weakly Interacting Massive Particles (WIMPs)
[17, 18], axions [19, 20, 21], sterile neutrinos [22], and Primordial Black
Holes (PBHs). Despite being initially proposed in the early ’70s [23, 24, 25,
26], PBHs received a newly found interest after the first LIGO detection [27,
28, 29]. In contrast to ordinary astrophysical black holes (ABHs) created at
the end of the stellar cycle of very massive stars, PBHs form (in the most
standard scenario) in the early Universe, during the radiation-dominated era.
Moreover, since their formation mechanism is not related to any astrophysical
process, they can in principle have any mass.

One guaranteed signature of this dark matter candidate is GW emission from
dynamical processes involving PBHs, such as binaries’ inspirals and mergers,
see, e.g., refs. [30, 31]. In fact, current constraints on PBH abundance come
from a variety of different probes, see, e.g., ref. [32, 33] for recent
reviews, including GW experiments. Depending on the nature of the process and
the mass of PBHs, the characteristic GW emission can peak over a broad
frequency range: from the nHz band probed by pulsar timing arrays [34, 35, 36,
37, 38], to the mHz-dHz band probed by future space-borne interferometers
(e.g., LISA [39], Big Bang Observer [40], and DECIGO [41, 42]), up to the Hz-
kHz band targeted by current and future ground-based GW observatories (e.g.,
LIGO, Virgo, KAGRA, Einstein Telescope [43, 44], and Cosmic Explorer [45,
46]). Therefore, PBHs are expected to leave an imprint on top of the
gravitational wave background (GWB) sourced by other known astrophysical
sources, such as astrophysical black holes of any mass, neutron stars and
white dwarfs.111In this work, we do not consider the so called cosmological
GWB, sourced in the very early Universe, see e.g., [33] for a review. In that
case, the process of PBH formation can further contribute to GW emission. It
then becomes extremely timely and important to characterize PBH GW signatures,
and possibly find novel techniques to isolate PBH contributions from the
expected astrophysical ones.

Several authors already analyzed detectability prospects of PBH binaries [47,
48, 49, 50, 51, 52, 53] that formed either in the early Universe before
matter-radiation equality [29, 54, 55], or at late times in small dark matter
halos [27, 28]. However, limited attention has been payed to PBH phenomenology
in dense astrophysical environments, such as stellar clusters [56, 57]. Since
it is widely recognized that compact object binaries in different kinds of
cluster strongly contribute to GW emission due to enhanced interaction rates
[57], it is also reasonable to expect a strong PBHs contribution coming from
the very same type of environments. If so, current estimates of the GWB are
missing a potentially significant contribution coming not only from PBH-PBH
and PBH-ABH binaries, but also from GW bursts generated by compact objects
flybys. Several studies in recent years highlighted the role of ABHs in
globular [58, 59, 60], nuclear [61, 62, 63, 64, 65] and young [66, 67, 68, 69,
70] clusters. Despite all these types of star clusters being gravitationally
bound systems, they have substantially different formation histories and
properties. In particular, globular clusters (GCs) represent the most
interesting candidate to study the impact of PBHs on dense environments,
because of their early formation time and large central density.

This work aims to provide, for the first time, a comprehensive modeling of the
GWB produced by PBH dynamical encounters in dense environments. We build a
semi-analytical model that accounts for PBHs in GCs, and we derive from first
principles all the relevant probability distribution functions that regulate
gravitational interactions between different kinds of compact objects in the
cluster. Then we estimate the contribution to the GWB generated not only by
binaries formed via direct capture and three-body processes, but also by
hyperbolic encounters. The former contribution peaks in the frequency band
probed by the ground-based GW observatories, and we show that this binary
formation channel is potentially detectable for next generation instruments.
Moreover, we show that this novel PBH GW signature can be more prominent than
those of other more established channels, such as early- and late-time PBH
binaries. On the other hand, despite the large interaction rate, the
hyperbolic encounters contribution peaking in the mHz band will not be
detectable by future space-borne GW detectors.

The paper is organized as follows. In section 2 we describe PBH binary
formation channels, and we start introducing a coherent model of GCs that
incorporates PBHs. The GC model is fully characterized in section 3. In
section 4 we compute the PBH contribution to the GWB, including different
kinds of dynamical interaction in different frequency bands, and we compare it
to the more established contributions, both primordial and astrophysical.
Finally, we conclude in section 5. Appendices contain additional information
regarding typical timescales associated with the GC model (appendix A),
relative velocity distributions of compact objects inside GCs (appendix B),
hyperbolic encounters dynamics (appendix C) and probability distribution
functions (appendix D).

Fiducial Cosmology. In this work we assume a spatially flat $\Lambda$CDM
cosmology [71] with present-day Hubble expansion rate
$H_{0}=67\,\mathrm{km/s/Mpc}$, present-day matter and radiation relative
density $\Omega_{\mathrm{m}0}=0.307$ and $\Omega_{\mathrm{r}0}=9.061\times
10^{-5}$, respectively.

##  2 Primordial black holes

In the most common scenarios, PBHs form before Big Bang Nucleosynthesis from
the gravitational collapse of large density perturbations; see, e.g., [31, 30]
for extensive reviews of different formation mechanisms. However, regardless
of their origin, only PBHs with masses larger than $M_{\mathrm{PBH}}\gtrsim
10^{-18}\,M_{\odot}$ evaporate via Hawking radiation on a timescale larger
than a Hubble time, and are considered as a viable dark matter candidate.
Their abundance is parameterized by the relative abundance parameter
$f_{\mathrm{PBH}}=\bar{\rho}_{\mathrm{PBH}}/\bar{\rho}_{\mathrm{dm}}$, where
$\bar{\rho}_{\mathrm{PBH}}$ and $\bar{\rho}_{\mathrm{dm}}$ are the background
PBH and dark matter energy densities. The mass distribution of PBHs is
described as [72]

$\frac{df_{\mathrm{PBH}}}{dM_{\mathrm{PBH}}}=f_{\mathrm{PBH}}\frac{d\Phi_{\mathrm{PBH}}}{dM_{\mathrm{PBH}}}\,,$
(2.1)

where $d\Phi_{\mathrm{PBH}}/dM_{\mathrm{PBH}}$ encodes the shape of the PBH
mass distribution, and it is normalized to unity. Different formation
mechanisms typically produce different mass distributions, however in the
following we focus exclusively on the so-called monochromatic case, widely
used in the literature, where
$d\Phi_{\mathrm{PBH}}/dM_{\mathrm{PBH}}=\delta^{D}(M_{\mathrm{PBH}}-M^{\star}_{\mathrm{PBH}})$
is a Dirac delta centered at mass $M^{\star}_{\mathrm{PBH}}$. Our results can
be converted so to account for an extended mass function e.g., by following
the approach in ref. [72].

While current abundance constraints appear to rule out PBHs as a significant
constituent of dark matter [31, 30] on some mass ranges, these upper bounds on
$f_{\mathrm{PBH}}$ suffer from many kinds of uncertainty. In particular,
constraints in the $\mathcal{O}(10)\ M_{\odot}$ mass range have been showed to
not be robust under changes in the PBH mass distribution, see, e.g., ref.
[72], or under changes in the modeling of the observable, see, e.g., refs.
[73, 74, 75]. Since in the next sections we are mainly interested in GW
emission in the frequency band probed by current and future ground-based GW
observatories, we assume as benchmark values for PBH masses $M_{\rm
PBH}=1,10,100\,M_{\odot}$. Because of the aforementioned uncertainties, we
test different values for $f_{\mathrm{PBH}}$, including the extreme case
$f_{\mathrm{PBH}}=1$ in which PBHs constitute the totality of dark matter.
Finally, we assume PBHs to be spinless.

In this work we investigate three main PBH contributions to the GWB: early-
and late-time binaries, in sections 2.1 and 2.2, respectively, and the novel
PBH contribution sourced in dense environments, which we introduce in section
2.3. We discuss additional contributions that will be the subject of future
work in section 5.

###  2.1 Early-time primordial black hole binaries

Even during the radiation-dominated era, pairs of PBHs can form bound binary
systems if they decouple from the background expansion [76, 77, 29]. Tidal
forces generated by surrounding bodies and by the overall dark matter field
(in the case PBHs do not constitute the totality of dark matter) avoid heads-
on collisions and delay the merger, making the signal potentially observable
today [55, 54, 78]. Following the formalism outlined in ref. [52], we
parameterize the merger rate of early PBH binaries for a monochromatic PBH
mass function as

$R_{\rm EPBH}(z)=f_{\rm
PBH}^{53/37}\mathcal{A}_{m}\left(\frac{t_{0}}{t(z)}\right)^{34/37}\left(\frac{M_{\rm
PBH}}{30\ M_{\odot}}\right)^{-32/37}\,,$ (2.2)

where $t(z)$ is the cosmic time at the time of the merger, $t_{0}=t(0)$ is the
age of the Universe, and $\mathcal{A}_{m}$ is an amplitude parameter. There
are significant uncertainties on the amplitude, and thus in the early-time
binary merger rate, since this parameter accounts for the properties of the
binary environment at formation and during binary evolution across cosmic
times [54, 79, 80, 81, 82].

###  2.2 Late-time primordial black hole binaries

Another well-known PBH binary formation channel is given by direct capture
processes in low-mass dark matter halos [27, 28]. In this kind of
environments, where relative velocities between halo objects are relatively
low, unbounded PBHs can form binaries via the emission of GWs. Also in this
case, we follow ref. [52] and we parameterize the late-time binary merger rate
as

$R_{\rm LPBH}(z)=f_{\rm PBH}^{2}\int dt_{\rm d}p(t_{\rm d})\int dM_{\rm
h}\frac{dn_{\rm h}}{dM_{\rm h}}\mathcal{R}_{{\rm BF,h}}(M_{\rm h},z)\,,$ (2.3)

where $t_{\rm d}$ is the time delay between the binary formation and merger,
$p(t_{\rm d})$ is the time delay probability distribution function, $M_{\rm
h}$ is the mass of the host DM halo, $dn_{\rm h}/dM_{\rm h}$ the halo mass
function, and $\mathcal{R}_{{\rm BF,h}}(M_{\rm h},z)$ is the binary formation
rate per halo due to a direct capture process.

###  2.3 Primordial black holes in globular clusters

The dark matter content of GCs is largely debated. Even though numerical
simulations usually include dark matter at initial time [83, 84, 85, 86, 87],
this component is either retained or tidally stripped by the host galaxy
depending on the nature of dark matter and on which astrophysical
prescriptions are included in the simulations. The details of these processes
are currently highly uncertain and, depending on the parameters of the
simulation, the final outcome can be significantly different, see, e.g., refs.
[88, 89, 90, 91, 92]. Despite this underlying uncertainty, it is commonly
accepted that, if dark matter is retained inside GCs, it is more likely to be
found in their cores [93, 92, 94]. If PBHs account for at least part of the
dark matter and remain bounded to the clusters across their evolution, the
possibility of finding them in the core is even more likely, due to the
process of mass segregation, described in the next section. In this case, we
expect GC formation and evolution dynamics to be considerably influenced by
the presence of PBHs in their core, contributing not only to shape the
structure of such systems, but also to the total GW emission generated by the
cluster components.

In the following, we model the PBH abundance in GCs by introducing the cluster
abundance parameter

$f^{\mathrm{cl}}_{\mathrm{PBH}}=\frac{M^{\mathrm{tot}}_{\mathrm{PBH}}}{M_{\mathrm{b}}}\,,$
(2.4)

where $M_{\rm PBH}^{\rm tot}$ and $M_{\mathrm{b}}$ are the total PBH and
baryonic222As described in the next section, by baryonic mass we refer to
stars and ABHs, while we neglect the presence of subdominant contributions,
e.g., neutron stars, white dwarfs or interstellar gas. mass of the cluster,
respectively. Therefore, the total mass of a cluster is given by
$M_{\mathrm{cl}}=M_{\mathrm{b}}(1+f^{\mathrm{cl}}_{\mathrm{PBH}})$, and the
number of PBHs per cluster for a monochromatic mass distribution reads as
$N_{\rm PBH}=M_{\rm b}f_{\rm PBH}^{\rm cl}/M_{\rm PBH}$. In general,
$f^{\mathrm{cl}}_{\mathrm{PBH}}$ is different from the relative abundance
parameter $f_{\mathrm{PBH}}$ introduced in equation (2.1). In the following
sections we consider the benchmark cases for which
$f^{\mathrm{cl}}_{\mathrm{PBH}}=0.1,1$, and we argue about a possible link
between this parameter and $f_{\mathrm{PBH}}$ in section 3.3.

##  3 Globular clusters

Globular clusters are among the oldest structures in the Universe, they have
an almost-spherical shape, large mass and a dense compact core. In this work,
we assume GCs are made of three different components, namely stars, ABHs, and
PBHs, and we develop a semi-analytical model to account for their mass,
velocity, and density distributions. Hereafter, we refer to stars and ABHs as
baryonic mass components, and we account for them in the $M_{\rm b}$ parameter
introduced in equation (2.4). Our formalism can be generalized to other types
of cluster, for instance young or nuclear clusters, by changing the
characteristic mass, radius and redshift distributions. While an accurate
description of these environments can be obtained only via numerical
simulations, our modeling represents a useful tool to understand the relevance
of PBH GW signatures.

###  3.1 Stars and astrophysical black holes

We assume that GCs at formation time are made only of small-intermediate
$(M_{\star}\leq 8\,M_{\odot})$ and massive $(M_{\star}\geq 8\,M_{\odot})$
stars. We model the star distribution at initial time as a Kroupa mass
function $dN_{\star}/dM_{\star}=A_{\star}M^{-\alpha}_{\star}$ [95], where the
normalization constant reads as

$A_{\star}=\frac{M^{\mathrm{ini}}_{\mathrm{b}}}{\displaystyle\int^{M^{\mathrm{max}}_{\star}}_{M_{\star}^{\mathrm{min}}}dM_{\star}M^{1-\alpha}_{\star}}\,,$
(3.1)

$M^{\mathrm{ini}}_{\mathrm{b}}$ is the initial baryonic mass of the cluster,
the scaling exponent is $\alpha=2.3$, and the initial minimum and maximum star
masses are taken to be $M_{\star}^{\mathrm{min}}=0.5\,M_{\odot}$ and
$M_{\star}^{\mathrm{max}}=40\,M_{\odot}$, respectively. Since observations
seem to suggest that interstellar gas does not play a crucial role in GCs, we
neglect its contribution.

The small-intermediate star population is characterized by a long lifetime and
negligible mass loss during its evolution; hence, we assume that the low-mass
tail of the mass distribution does not evolve in time and we do not account
for small mass compact objects. Conversely, stars in the massive population
have shorter lifetimes and form BHs at the end of their stellar cycle.

We assume a negligible mass loss during the stellar lifetime and during the
gravitational collapse that leads to the ABH formation, i.e.,
$M_{\mathrm{ABH}}\simeq M_{\star}$, and we expect a population of ABHs to form
on a timescale of order $\mathcal{O}(10)\ \mathrm{Myr}$ (see also appendix A)
in the mass range $M_{\mathrm{ABH}}^{\mathrm{min}}=8\,M_{\odot}\leq
M_{\mathrm{ABH}}\leq 40\,M_{\odot}=M_{\mathrm{ABH}}^{\mathrm{max}}$.333 Stars
with masses of order $100\,M_{\odot}$ can also form in GCs. However, those
stars are expected to suffer from a strong stellar mass loss before entering
the BH formation stage [85]. Hence, even choosing a larger maximum star mass,
we still produce BHs in the same mass range we are considering here. Given
that the relative number of very massive stars is rather low because of the
steepness of the star mass function, we choose
$M_{\star}^{\mathrm{max}}=40\,M_{\odot}$, avoiding any modeling related to
stellar mass loss. After all massive stars have disappeared, the ABH mass
function can only evolve due to hierarchical mergers or cluster evaporation.
Regarding the former, simulations [96] suggest that the evolution of the mass
function is weak in the mass range of interest, thus we neglect it in our
parametrization. On the other hand, evaporation can be prompted by large kick
velocities acquired by compact objects during the supernova process. We
estimate this effect by introducing a retention fraction parameter $f_{\rm
ret}$ to rescale the initial baryonic mass in ABHs and account for such mass
loss. Typical values for the retention fraction are $f_{\mathrm{ret}}\lesssim
1/2$, but exact values are still largely debated in the literature [97, 98].
Following ref. [97] we define

$f_{\rm ret}=\int_{0}^{v_{\rm esc}}dv_{\rm kick}\,p(v_{\rm kick})\,,$ (3.2)

where $v_{\rm kick}$ is the natal kick velocity after the supernova event and
$v_{\rm esc}$ is the GC escape velocity. We model the natal kick probability
distribution function as a Maxwell-Boltzmann distribution

$p(v_{\rm kick})=\frac{4\pi v_{\rm
kick}^{2}}{\left(2\pi\sigma^{2}_{\mathrm{kick}}\right)^{3/2}}e^{-v_{\rm
kick}^{2}/2\sigma_{\rm kick}^{2}}\,,$ (3.3)

and we verify that, for a typical velocity dispersion of $\sigma_{\rm
kick}\simeq 10\,{\rm km/s}$ [97], GCs with baryonic mass $M_{\mathrm{b}}\simeq
10^{5}\,M_{\odot}$ have $f_{\rm ret}\simeq 0.3$. Thus we set
$f_{\mathrm{ret}}=1/3$ as conservative fiducial value for ABHs in all GCs,
independently on the cluster mass. Because of the mass loss due to natal
kicks, the (final) baryonic mass of the cluster is lower than the initial one,
i.e., $M_{\mathrm{b}}<M^{\mathrm{ini}}_{\mathrm{b}}$. We assume that after
ABHs with large natal kick velocities have left the cluster the mass functions
of small-intermediate stars and of the remaining ABHs do not evolve further,
therefore $M_{\rm b}$ can be effectively considered as a constant.

In summary, star and ABH populations are characterized by a number of objects
and a total mass that read as

$N_{\star}=A_{\star}\int^{M^{\mathrm{min}}_{\mathrm{ABH}}}_{M^{\mathrm{min}}_{\star}}dM_{\star}M^{-\alpha}_{\star}\,,\qquad
M^{\mathrm{tot}}_{\star}=A_{\star}\int^{M^{\mathrm{min}}_{\mathrm{ABH}}}_{M^{\mathrm{min}}_{\star}}dM_{\star}M^{1-\alpha}_{\star}\,,$
(3.4)

and

$N_{\mathrm{ABH}}=A_{\star}\int_{M^{\mathrm{min}}_{\mathrm{ABH}}}^{M^{\mathrm{max}}_{\mathrm{ABH}}}dM_{\star}M^{-\alpha}_{\star}f_{\mathrm{ret}}\,,\qquad
M^{\mathrm{tot}}_{\mathrm{ABH}}=A_{\star}\int_{M^{\mathrm{min}}_{\mathrm{ABH}}}^{M^{\mathrm{max}}_{\mathrm{ABH}}}dM_{\star}M^{1-\alpha}_{\star}f_{\mathrm{ret}}\,,$
(3.5)

respectively, explicitly defining the baryonic mass of the cluster
$M_{\mathrm{b}}=M_{\star}^{\mathrm{tot}}+M^{\mathrm{tot}}_{\mathrm{ABH}}$.

###  3.2 Globular cluster internal structure

The evolution of the internal structure of a GC depends on many physical
effects, for instance stellar evolution, $N$-body interactions, the effect of
the tidal field of the hosting galaxy [99]. Despite this intrinsic complexity,
observational properties can be recovered under a modest set of assumptions
that capture the main features of these systems [100, 101].

The spatial and velocity distribution of objects inside a cluster can be well
represented by a distribution function $f(v,r)$, where $v$ is the velocity and
$r$ the distance from the cluster center. The energy of an object in the
cluster is defined as $E(v,r)=\frac{1}{2}v^{2}-\Psi(r)$, where $\Psi(r)$ is
the effective gravitational potential, defined in such a way that $\Psi(r_{\rm
t})=0$ at the truncation radius $r_{\rm t}$, where objects are effectively
stripped away from the cluster. In this work, we consider an isotropic multi-
mass King model [102, 103]444 Generalizations of this model, for instance
including anisotropies, are also employed in the literature, see, e.g., ref.
[104]; however, the King profile describes the cluster dynamics well enough
for our purposes. of the form

$f_{j}(E)=\frac{\rho_{0j}}{(2\pi\sigma^{2}_{j})^{3/2}\mathcal{I}_{0j}}\left(e^{-E(v,r)/\sigma^{2}_{j}}-1\right)\Theta_{H}\left(-E(v,r)\right)\,,$
(3.6)

where $j=\\{\rm stars,\,ABH,\,PBH\\}$ are the GC components, $\Theta_{H}(x)$
is the Heaviside Theta function and $\sigma_{j}$ is a velocity scale, which
approaches the value of the 1-dimensional velocity dispersion in the limit of
very concentrated clusters, i.e., when $\sqrt{\Psi}/\sigma_{j}\gg 1$. The
normalization constant $\mathcal{I}_{0j}$ is defined to recover the central
density of each component, which is labelled as
$\rho_{j}(\Psi_{0})=\rho_{0j}$. Velocity and density profiles are indirectly
linked via the gravitational potential through

$\rho_{j}(\Psi)=4\pi\int_{0}^{\sqrt{2\Psi}}dv\
v^{2}f_{j}(E)=\frac{\rho_{0j}}{\mathcal{I}_{0j}}\left[e^{\Psi/\sigma^{2}_{j}}{\rm
erf}\left(\frac{\sqrt{\Psi}}{\sigma_{j}}\right)-\sqrt{\frac{4\Psi}{\pi\sigma^{2}_{j}}}\left(1+\frac{2\Psi}{3\sigma^{2}_{j}}\right)\right]\,,$
(3.7)

where the radial dependence on $r$ is left implicit in $\Psi$ and the
integration upper bound is the escape velocity
$v_{\mathrm{esc}}(r)=\sqrt{2\Psi(r)}$. The total density in the center of the
cluster, $\rho_{0}=\sum_{j}\rho_{0j}$, is determined by the Poisson equation

$\frac{d}{dr}\left(r^{2}\frac{d\Psi}{dr}\right)=-4\pi
Gr^{2}\sum_{j}\rho_{j}(\Psi)\,,$ (3.8)

where $G$ is the Newton constant. The velocity scale for each component
$\sigma_{j}$ is set by energy equipartition:555 More sophisticated models of
the cluster internal structure account for the evolution in the density
profile driven by interactions and energy exchange between different bodies.
It has been shown that more massive bodies increase their velocities through
these interactions, which in turn can lead to evaporation from the cluster and
to the dynamical evolution of the overall structure across cosmic time. For
this reason, clusters with a more complex dynamical model never reach the
condition of energy equipartition [105, 106, 107, 108]. A careful modeling of
these processes would require numerical simulations and goes beyond the scope
of this work. because of two-body gravitational interactions, energy is
redistributed between objects in the cluster, tending towards an equilibrium
state [109]. Under the (partial) equipartition condition, the velocities of
two classes of objects with masses $m_{i}$ and $m_{j}$ scale as
$\sigma_{i}m^{\delta}_{i}=\sigma_{j}m^{\delta}_{j}$, where $\delta\in[0,1/2]$,
and the case $\delta=1/2$ corresponds to a perfect energy equipartition
scenario.

We use the public code limepy666https://github.com/mgieles/limepy [110] to
compute density and velocity profiles of GCs. In limepy, we are required to
specify the concentration parameter $\Psi_{0}/\sigma^{2}$, where
$\sigma=(m_{j}/\bar{m})^{\delta}\sigma_{j}$ is the velocity scale
corresponding to the central density weighted mean mass
$\bar{m}=\sum_{j}m_{j}\rho_{0j}/\rho_{0}$, and the half-mass radius $r_{\rm
h}$, i.e., the radius that encloses half of the mass of the cluster. We
consider a benchmark value of $\Psi_{0}/\sigma^{2}=10$ for the concentration
parameter [100, 111], and of

$r_{\rm
h}=3\left(\frac{M_{\mathrm{cl}}}{{10^{5}\,M_{\odot}}}\right)^{1/3}\,\mathrm{pc}\,,$
(3.9)

for the half-mass radius. These values are compatible with existing fit to
observations, see, e.g., ref. [112].

Figure 1: Star (left panels), ABH (central panels) and PBH (right panels)
density (top panels) and average squared velocity (bottom panels) radial
profiles for different value of the PBH cluster abundance parameter
$f^{\mathrm{cl}}_{\mathrm{PBH}}$. We show profiles for clusters with
$M_{\mathrm{b}}=10^{5}\,M_{\odot}$ and $M_{\mathrm{PBH}}=100\,M_{\odot}$.

As a result of gravitational interactions, more massive objects experience a
dynamical friction effect, migrate towards the core of the GC, and become
effectively segregated. The typical timescale of the equipartition process is
of the order of the half-mass relaxation time [113, 114],

$t_{\mathrm{rel}}=0.17\frac{N}{\log(\lambda
N)}\sqrt{\frac{r^{3}_{h}}{GM_{\mathrm{cl}}}}\,,$ (3.10)

where $N$ is the total number of objects in the cluster, and $\lambda=0.11$ is
a coefficient calibrated on numerical simulations [115]. We show in figure 1
the stable configuration reached by clusters with different PBH abundance for
$M_{\mathrm{PBH}}=100\ M_{\odot}$ after a typical timescale of order
$t_{\mathrm{rel}}(M^{\mathrm{ini}}_{\mathrm{b}}=10^{5}\ M_{\odot})\simeq 300\
\mathrm{Myr}$. As we see from the bottom panels, increasing the PBH abundance
consistently boosts the typical average squared velocity

$\displaystyle\left\langle v^{2}\right\rangle_{j}$
$\displaystyle=\dfrac{\displaystyle 4\pi\int_{0}^{\sqrt{2\Psi}}dv\
v^{4}f_{j}(E)}{\rho_{j}(\Psi)}$ (3.11)
$\displaystyle=3\sigma^{2}_{j}\dfrac{e^{\Psi/\sigma_{j}^{2}}{\rm
erf}\left(\sqrt{\Psi}/\sigma_{j}\right)-\sqrt{4\Psi/\pi\sigma_{j}^{2}}\left(1+2\Psi/3\sigma_{j}^{2}+4\Psi^{2}/15\sigma_{j}^{4}\right)}{e^{\Psi/\sigma^{2}_{j}}{\rm
erf}\left(\sqrt{\Psi}/\sigma_{j}\right)-\sqrt{4\Psi/\pi\sigma_{j}^{2}}\left(1+2\Psi/3\sigma_{j}^{2}\right)}\,.$

of each component, which more massive components (PBHs in the case displayed
here) having a smaller average velocity because of energy equipartition.
Density profiles are affected by the presence of massive PBHs in a non-trivial
way: on the one hand, objects tend to be more concentrated because of the
dense core formed by PBHs, on the other the increased average velocity tends
to dilute objects lighter than PBH. While these two effects appear to average
out in the case of stars, we observe that the latter tend to dominate in the
case of ABHs, which present a lower core density than in the
$f^{\mathrm{cl}}_{\mathrm{PBH}}=0$ case. Historically, the size of the GC core
has been described by the King radius $r^{2}_{0}=9\sigma^{2}/(4\pi
G\rho_{0})$. In our multi-mass King model we use the King radii
$r_{0j}=(\bar{m}/m_{j})^{\delta}r_{0}$ as a measure of each component “core”,
where mass densities and velocities can be considered spatially uniform.
Outside this internal regions, density and velocity profiles rapidly decay, as
can be seen in figure 1. In the following section we will refer to these King
radii as “segregation radii”.

###  3.3 Globular cluster number density

Assuming that GC mass distribution scales as
$(M^{\mathrm{ini}}_{\mathrm{b}})^{-2}$ [116, 117, 118, 86], we have that

$\frac{dn_{\mathrm{cl}}}{dM^{\mathrm{ini}}_{\mathrm{b}}}=\frac{m_{\mathrm{b,tot}}(z)}{\left(M^{\mathrm{ini}}_{\mathrm{b}}\right)^{2}\log(M^{\mathrm{ini}}_{\mathrm{b,max}}/M^{\mathrm{ini}}_{\mathrm{b,min}})}\,,$
(3.12)

where we consider $M^{\mathrm{ini}}_{\rm b,min}=5\times 10^{4}\,M_{\odot}$ and
$M^{\mathrm{ini}}_{\rm b,max}=10^{7}\,M_{\odot}$ as the minimum and maximum
initial baryonic mass of a cluster, and $m_{\mathrm{b,tot}}(z)$ is the total
baryonic mass density in GCs at redshift $z$. The latter is defined as
$\displaystyle m_{\mathrm{b,tot}}(z)=\int^{t(z)}_{0}dt\ \dot{m}_{\mathrm{b}}$,
where the GC formation density rate is fitted to numerical simulations as [86,
119]

$\dot{m}_{\mathrm{b}}(z)=\frac{a}{\sigma\sqrt{2\pi}}\exp\left[-(z-\mu)^{2}/(2\sigma^{2})\right]\left[1+{\rm
erf}\left(\frac{\alpha(z-\mu)}{\sqrt{2}\sigma}\right)\right]\,,$ (3.13)

and the fitting parameters read as
$\\{\sigma,\mu,\alpha,a\\}=\\{2.7,2.7,1.5,9.0\times
10^{5}\,M_{\odot}\,\mathrm{Gpc^{-3}\,Myr^{-1}}\\}$. We find
$n_{\mathrm{cl}}(z=0)\simeq 3\ {\rm Mpc}^{-3}$, in agreement with current
literature [120, 121, 122]. The implementation of a time-evolving cluster mass
function is left for future work. Notice that the cluster mass function is
expressed in terms of the initial mass $M_{\rm b}^{\rm ini}$, which differs
from the baryonic mass $M_{\rm b}$ by the amount of ABHs that have been
ejected from the cluster, as described in section 3.1.

In section 2 we introduced two abundance parameters to describe PBHs in the
Universe, $f_{\mathrm{PBH}}$, and in clusters,
$f^{\mathrm{cl}}_{\mathrm{PBH}}$. While on the one hand we already showed in
sections 2.1 and 2.2 that GW emission from early- and late-time binaries
scales with powers of $f_{\mathrm{PBH}}$, on the other hand we expect PBHs
contribution from clusters to depend on $f^{\mathrm{cl}}_{\mathrm{PBH}}$. If
we compare the total energy density in PBH
($\bar{\rho}_{\mathrm{PBH}}=f_{\mathrm{PBH}}\bar{\rho}_{\mathrm{dm}}$) to the
total PBH energy density in GCs, which reads as

$\displaystyle\rho^{\mathrm{cl}}_{\mathrm{PBH}}(z)$ $\displaystyle=\int
dM^{\mathrm{ini}}_{\mathrm{b}}\frac{dn_{\mathrm{cl}}}{dM^{\mathrm{ini}}_{\mathrm{b}}}f^{\mathrm{cl}}_{\mathrm{PBH}}M_{\mathrm{b}}$
(3.14)
$\displaystyle=\left[1-(1-f_{\mathrm{ret}})\frac{\int_{M^{\mathrm{min}}_{\mathrm{ABH}}}^{M^{\mathrm{max}}_{\mathrm{ABH}}}dM_{\star}M_{\star}^{1-\alpha}}{\int_{M^{\mathrm{min}}_{\star}}^{M^{\mathrm{max}}_{\star}}dM_{\star}M_{\star}^{1-\alpha}}\right]f^{\mathrm{cl}}_{\mathrm{PBH}}m_{\mathrm{b,tot}}(z)\,,$

we observe that, for instance, at redshift $z=0$, we can have the extreme case
where GCs have equal parts of PBHs and baryons (i.e.,
$f^{\mathrm{cl}}_{\mathrm{PBH}}=1$), while having PBH relative abundance
parameter as low as $f_{\mathrm{PBH}}\simeq 2\times 10^{-5}$. Even if this
example is clearly extreme, it showcases that there are scenarios where PBH GW
emission in certain channels is heavily suppressed while GW emission from
dense environments is not.

##  4 Primordial black hole gravitational wave background

In full generality, the total GWB in any frequency band can be written as the
sum of independent GWBs generated by different sources and/or processes, i.e.,
$\Omega_{\mathrm{GWB}}=\sum_{\beta,i,j}\Omega^{(\beta)}_{ij}$, where $\beta$
labels the different the process, while $i,j$ indicate the sources (either ABH
or PBH). Not all combinations of dynamical processes $(\beta)$ and sources
$(i,j)$ are necessarily possible; in the following, we only discuss the ones
that are relevant for our work.

The individual backgrounds are given by [123, 124]

$\Omega^{(\beta)}_{ij}=\frac{f_{\rm
obs}}{\rho_{0c}c^{2}}\int\frac{dz\,d\boldsymbol{\theta}\,p(\boldsymbol{\theta},z)}{(1+z)H(z)}\mathcal{R}_{ij}^{(\beta)}(\boldsymbol{\theta},z)\frac{dE^{(\beta)}_{ij}}{df}(\boldsymbol{\theta},z,f_{\rm
obs})\,,$ (4.1)

where $f_{\rm obs}=f/(1+z)$ and $f$ are the observed and emitted GW frequency,
respectively, $z$ is the GW emission redshift, $\rho_{0c}$ is the present-day
critical density, $c$ is the speed of light, $H(z)$ the Hubble expansion rate.
We label with $\boldsymbol{\theta}$ the set of parameters that describe both
the process and the sources, and we introduce $p(\boldsymbol{\theta},z)$ to
characterize their probability distribution functions. These are used to
compute the dynamical process density rate $\mathcal{R}_{ij}^{(\beta)}$, and
the GW energy spectrum $dE^{(\beta)}_{ij}/df$. The latter is commonly written
as [125]

$\frac{dE^{(\beta)}_{ij}}{df}=\frac{2\pi^{2}D^{2}_{L}c^{3}f^{2}_{\rm
obs}}{G(1+z)^{2}}\left|h^{(\beta)}_{ij}\right|^{2}\,,$ (4.2)

where $D_{L}$ is the luminosity distance and
$h^{(\beta)}_{ij}(\boldsymbol{\theta},z,f_{\rm obs})$ is the emitted GW
strain.

The dimensionality of the parameter space in equation (4.1) rapidly increases
when accurate descriptions of dynamical processes are actually implemented, as
in this work. Apart from the emission redshift $z$, we are interested in
integrating over the GC initial baryonic masses
$M^{\mathrm{ini}}_{\mathrm{b}}$, the masses of the compact objects
$M_{i},M_{j}$, and any other astrophysical parameter that describes the
dynamical process at hand. In order to make the problem numerically tractable,
we choose to treat ABH and PBH mass distributions on equal grounds, i.e., as
monochromatic mass functions. While for PBHs this is already the case by
assumption, we compute an effective ABH mass $M_{\mathrm{ABH}}^{\mathrm{eff}}$
by requiring that the GWB produced by such monochromatic ABH mass distribution
matches the true GWB produced by the real ABH mass distribution in the GC
core. In practice, the ABH effective mass for each process and combination of
sources (in which at least one of them is an ABH) is computed by solving,
either analytically or numerically, the equation

$\Omega^{(\beta)}_{ij}\left(\frac{dn_{\mathrm{ABH}}}{dM_{\mathrm{ABH}}}=\frac{\rho_{0\rm
ABH}}{M_{\mathrm{ABH}}^{\mathrm{eff}}}\delta^{D}(M_{\mathrm{ABH}}-M_{\mathrm{ABH}}^{\mathrm{eff}})\right)=\Omega^{(\beta)}_{ij}\left(\frac{dn_{\mathrm{ABH}}}{dM_{\mathrm{ABH}}}=\frac{(2-\alpha)M^{-\alpha}_{\mathrm{ABH}}\rho_{0\rm
ABH}}{(M^{\mathrm{max}}_{\mathrm{ABH}})^{2-\alpha}-(M^{\mathrm{min}}_{\mathrm{ABH}})^{2-\alpha}}\right)\,,$
(4.3)

where on the right hand side we adopted the Kroupa mass function defined in
section 3.2.

Equation (4.1) accounts for all GW emissions produced by the $(i,j)$ objects
in the process $(\beta)$; depending on the signal-to-noise ratios, these will
be either resolved events or part of the background. Since the aim of this
work is not to characterize detection prospects of the different GWB presented
in this work, we do not distinguish between sources that are effectively
resolved or not by a specific GW detector network. Instead, all the GWBs
reported hereafter refer to the total contribution generated by all the
sources in any given processes. We leave accurate detectability forecasts for
future works.

###  4.1 Early- and late-time primordial black hole binaries

Following ref. [52], we compute the GWB generated by early- and late-time
binaries using the external modules of CLASS_GWB [126]. The external modules
create a catalog of GW binaries by sampling the properties of the GW sources
from the characteristic probability distribution functions of the parameters
that describe the merging event. For instance, in the case of early- and late-
time PBH binaries, we sampled the merging redshift starting from the merger
rates in equations (2.2) and (2.3), their masses, spins, sky localization,
inclination angle and polarization angle. The GWB is then computed by summing
over the contribution of all the sources in the catalog. More detail on this
procedure can be found in ref. [126].

###  4.2 Astrophysical black hole binaries

Another useful term of comparison for our results is given by the GWB
generated by a generic population of astrophysical BHs with properties
compatible with data from the latest LIGO-Virgo-KAGRA collaboration results
[6, 7]. In this case, we do not specify the process that originated each
binary, even if there are already some attempt to characterize different
astrophysical formation channels [63]. Also in this case, the GWB sourced by
this population is computed using the external modules of CLASS_GWB following
the model and the procedure described in ref. [126].

###  4.3 Dense environment binaries

The formation of binary BHs in dense stellar environments is driven by several
different physical mechanisms, see, e.g., ref. [127] and references therein.
The dominant formation channel is expected to be the two-body direct capture
($\beta={\rm DC}$ hereafter) process [128, 129], where pairs of unbounded
compact objects end up in a bound system due to GW emission during their
flyby. However, other dynamical interactions are also possible; in this work,
we account also for three-body interactions ($\beta=3{\rm B}$) [127], where
the formation of the bound system is catalyzed by gravitational interaction
with a third body. Other processes such as multiple encounters, hardening of
the binary, disruption and exchanges of binary components with objects in the
environment, are expected to provide a subdominant contribution [127], thus we
neglect them in our analysis.

For each process and source pair, we define the interaction rate per compact
object mass per cluster as

$\Gamma^{(\beta)}_{ij}=\left(1-\frac{1}{2}\delta^{K}_{ij}\right)\frac{dn_{i}}{dM_{i}}\frac{dn_{j}}{dM_{j}}V_{\rm
seg}\sigma^{(\beta)}v_{\rm rel}\,,$ (4.4)

where the Kronecker delta $\delta^{K}_{ij}$ avoids double counting for
identical sources, i.e., when $i=j$. $V_{\rm seg}={\rm min}[r_{0i},r_{0j}]$ is
the volume defined by the smallest segregation radius, inside which number
densities $dn_{i,j}/dM_{i,j}$ are spatially uniform, $v_{\rm rel}$ is the
value of the relative velocity between the objects $i$ and $j$ and
$\sigma^{(\beta)}$ is the cross section of the dynamical process. In the case
of three body encounters, we estimate $\sigma^{(}\beta)$ by marginalizing over
the properties of the third body that catalyze the interaction. Therefore, the
total interaction rate inside each cluster is given by

$\Gamma_{ij,\mathrm{tot}}^{(\beta)}=\left(1-\frac{1}{2}\delta^{K}_{ij}\right)\int
dM_{i}dM_{j}\frac{dn_{i}}{dM_{i}}\frac{dn_{j}}{dM_{j}}V_{\rm
seg}\left\langle\sigma^{(\beta)}v_{\rm rel}\right\rangle\,,$ (4.5)

where we sum over the masses of the compact objects and we average the
velocity-weighted cross section over the process-dependent parameters
described by the probability distribution functions introduced in equation
(4.1). In the previous equation, the quantities $V_{\rm seg},\,v_{\rm
rel},\,dn_{i,j}/dM_{i,j}$ depend on $M_{\rm b}^{\rm ini}$, and they are
estimated by relying on the limepy [110] public code, as described in section
3.2.

####  4.3.1 Direct capture

The cross section in the direct capture case is [128, 129]

$\sigma^{\mathrm{(DC)}}=2\pi\left(\frac{85\pi}{6\sqrt{2}}\right)^{2/7}\frac{G^{2}M^{12/7}\mu^{2/7}}{c^{10/7}v_{\rm
rel}^{18/7}}\,,$ (4.6)

where the total and reduced masses read as $M=M_{i}+M_{j}$ and
$\mu=M_{i}M_{j}/M$, respectively. In this case the velocity-weighted cross
section has to account only for the relative velocity probability distribution
function, i.e.,

$\left\langle\sigma^{(\mathrm{DC})}v_{\rm rel}\right\rangle=\int
dv_{\mathrm{rel}}p(v_{\mathrm{rel}})\sigma^{(\mathrm{DC})}v_{\rm rel}\,.$
(4.7)

We provide an analytical estimate of $p(v_{\mathrm{rel}})$ in appendix B.

####  4.3.2 Three-body encounters

Three-body encounters can result in the formation of bounded binary systems,
however only hard binaries are expected to survive [130]. Following the
current literature, we consider only binaries with hardness ratio $\eta$,
i.e., the ratio between the average potential and kinetic energy of the
environment [131, 111], that is larger than the typical minimum value
$\eta_{\mathrm{min}}\simeq 5$ [132]. Adapting the results of refs. [131, 132,
133] to our formalism, we define the average velocity-weighted cross section
for three-body encounters as

$\displaystyle\left\langle\sigma^{(\mathrm{3B})}v_{\rm rel}\right\rangle$
$\displaystyle=\sqrt{2}\pi^{2}\frac{G^{5}(M_{i}+M_{j})^{5}}{\left\langle
v_{\mathrm{rel}}\right\rangle^{9}}\frac{1+2\eta_{\rm min}}{\eta_{\min}^{5.5}}$
(4.8) $\displaystyle\qquad\qquad\sum_{k=\mathrm{ABH,PBH}}s^{-1}_{ij,k}\int
dM_{k}\frac{dn_{k}}{dM_{k}}\frac{V_{k,\mathrm{seg}}}{V_{\mathrm{seg}}}\left[1+2\eta_{\rm
min}\left(\frac{M_{i}+M_{j}+M_{k}}{M_{i}+M_{j}}\right)\right],$

where $\left\langle
v_{\mathrm{rel}}\right\rangle=2\sqrt{2/\pi}\sigma_{\mathrm{rel}}$ and
$\sigma_{\rm rel}^{2}=\sigma_{i}^{2}+\sigma_{j}^{2}$ are the average relative
velocity and the relative velocity dispersion between bodies $i$ and $j$,
respectively, $s_{ij,k}$ is a numerical factor that avoids double counting,777
We have $s_{ij,k}=3$ if $i=j=k$, $s_{ij,k}=2$ if $i\neq j$ and $s_{ij,k}=1$ if
$i=j\neq k$. and the three-body segregation volume is defined as $V_{k,\rm
seg}={\rm min}[r_{0i},r_{0j},r_{0k}]$.

####  4.3.3 Total binary contributions from globular clusters

Independently on the formation channel, all binaries emit the same kind of GW
signal during their inspiral, merger and ringdown phases. In this work we use
the templates provided in refs. [134, 135, 136] to describe the GW strain
$h^{\mathrm{(DC,3B)}}_{ij}$ across those three phases of the coalescence
process. Each phase has its own unique dependency on the compact object
masses; however, since in the frequency band of interest most of the emission
is concentrated during the inspiral phase, we use the inspiralling mass
dependence to compute the ABH effective mass. Therefore, for this purposes,

$\left|h^{\mathrm{(DC,3B)}}_{ij}\right|^{2}\propto\mu M^{2/3}.$ (4.9)

Figure 2: Binary formation rate density at redshift $z=0$ for direct capture
(left panel) and three-body (right panel) processes. In each panel we show
ABH-ABH (blue lines), ABH-PBH (green lines) and PBH-PBH (red lines) binaries,
both for $f_{\rm PBH}^{\rm cl}=1$ (solid lines) and $f_{\rm PBH}^{\rm
cl}=10^{-1}$ (dashed lines). The PBH mass is $M_{\rm PBH}=100\,M_{\odot}$ in
all cases reported.

In figure 2 we show the local, i.e., at redshift $z=0$, binary formation rate
density per cluster

$\mathcal{R}^{(\beta)}_{ij}(z,M^{\mathrm{ini}}_{\mathrm{b}})=\frac{dn_{\rm
cl}}{dM^{\mathrm{ini}}_{\rm b}}\Gamma_{ij,\rm tot}^{(\beta)}(M_{\rm b}^{\rm
ini})\,,$ (4.10)

for the direct capture and three-body channels. As showed in section 3.2, the
addition of a PBH component induces a non-trivial change in the GC structure,
both in terms of core densities and velocities, which translates into a non-
trivial behaviour in the binary formation rate density. First of all, we note
the presence of a hierarchy between the binary formation rates per cluster:
since PBHs form dense cores and dilute lighter ABHs, we have that the binary
formation rates of the former population is larger than that of the latter one
in both formation channels. Due to (the disappearance of) this dilution
effect, the binary formation rate of ABH-ABH systems increases when PBH
abundance in GC decreases, and vice versa, consistently with the picture
portrayed above. Moreover, the direct capture channel appears to be the
dominant contribution to the total binary formation rate in almost all cases
considered, except for GCs with low mass and low PBH abundance. While for
direct capture the binary formation rate per cluster decreases when
$f^{\mathrm{cl}}_{\mathrm{PBH}}$ does so, in the case of three-body processes
we observe the opposite behaviour for low mass GCs. This seemingly
counterintuitive effect is due to the fact that a lower PBH density is
effectively compensated by a lower relative velocity dispersion of PBHs due to
a shallower GC total potential well.

After they form, binaries merge in a characteristic time, which depends on the
initial conditions of the compact object system. Since those initial
conditions are stochastic, each formation channel is characterized by its own
time delay probability distribution function $p(t_{\rm d})$, see appendix A.
Therefore, the merger rate density is defined as

$R_{ij}^{(\beta)}(z)=\int dt_{\rm d}p(t_{\rm d})\int dM^{\mathrm{ini}}_{\rm
b}\mathcal{R}^{(\beta)}_{ij}(z_{\mathrm{bf}},M^{\mathrm{ini}}_{\mathrm{b}})\,,$
(4.11)

where the binary formation redshift $z_{\mathrm{bf}}$ is defined implicitly
for each time delay by $z=z_{\rm bf}+t_{\rm d}$, and

$t_{\rm d}=-\int_{z_{\mathrm{bf}}}^{z}\frac{dy}{(1+y)H(y)}\,.$ (4.12)

As we show in appendix A, time delays for our processes of interest are
negligible compared to the other timescales of the dynamical process,
therefore we approximate $p(t_{\rm d})\simeq\delta^{D}(t_{\rm d})$.

In table 1, we present the typical values of the local, i.e., at redshift
$z=0$, merger rate densities for our benchmark models.

$R_{ij}^{\rm(DC+3B)}(z=0)\ [{\rm Gpc}^{-3}\,{\rm yr}^{-1}]$ | ABH-ABH | ABH-PBH | PBH-PBH
---|---|---|---
$f^{\mathrm{cl}}_{\mathrm{PBH}}=1$ | $M_{\mathrm{PBH}}=1\ M_{\odot}$ | $2.8\times 10^{-3}$ | $1.2\times 10^{-3}$ | $4.4\times 10^{-3}$
$M_{\mathrm{PBH}}=10\ M_{\odot}$ | $5.7\times 10^{-3}$ | $2.4\times 10^{-2}$ | $4.5\times 10^{-2}$
$M_{\mathrm{PBH}}=100\ M_{\odot}$ | $1.3\times 10^{-5}$ | $7.8\times 10^{-4}$ | $1.2\times 10^{-1}$
$f^{\mathrm{cl}}_{\mathrm{PBH}}=10^{-1}$ | $M_{\mathrm{PBH}}=1\ M_{\odot}$ | $5.0\times 10^{-3}$ | $4.6\times 10^{-4}$ | $1.1\times 10^{-4}$
$M_{\mathrm{PBH}}=10\ M_{\odot}$ | $5.6\times 10^{-3}$ | $3.5\times 10^{-3}$ | $9.9\times 10^{-4}$
$M_{\mathrm{PBH}}=100\ M_{\odot}$ | $9.9\times 10^{-5}$ | $4.5\times 10^{-4}$ | $1.2\times 10^{-2}$
Table 1: Local binary merger rate summed over direct capture and three-body
processes for GCs with different PBH abundances and masses.

We note that in most of the cases where we have PBHs with
$M_{\mathrm{PBH}}\lesssim 10\ M_{\odot}$, the total BH (i.e., summing over
ABH-ABH, ABH-PBH, and PBH-PBH contributions) local merger rate is comparable
with the ABH-ABH merger rate $R^{\rm(DC+3B)}_{\rm ABH\ only}\simeq 5\times
10^{-3}\ \mathrm{Gpc^{-3}yr^{-1}}$ obtained considering GCs without PBHs.
While absolute numbers always depend on the model implemented and they should
be taken with caution, this trend seems to suggest that there might exist
scenarios where the knowledge of merger rates alone is not sufficient at all
in discriminating between radically different alternatives.

These systems emits in the Hz-kHz frequency band, which is targeted by current
and future ground-based interferometers. On the observational side, there are
already experimental upper limits on the GWB amplitude [137, 138, 139] in this
band, and a detection is within reach during next observational runs [140,
141]. However, in the following figures we include only the effective
“sensitivity” curves [142, 143] for a GW detector network given by Einstein
Telescope and two Cosmic Explorers (ET2CE). The sensitivity is determined by
the expected noise level for the experiments, and it accounts for one year of
observation time. These effective sensitivity curves
$\Omega^{\mathrm{eff}}_{n}(f)$ are a useful benchmark for comparison: any GWB
with a power-law-like frequency dependence (as in most of our cases) that
intersects these curves is expected to be detected with a signal-to-noise
ratio larger than unity after a year of observation. In figure 3, we show the
GWB generated by ABH-ABH, ABH-PBH, and PBH-PBH systems in GCs with different
PBH masses and abundances. We do not expect ABH-ABH and ABH-PBH contributions
to be detectable. On the other hand, if PBHs have masses
$M_{\mathrm{PBH}}\gtrsim 10\ M_{\odot}$, we observe that detection prospects
for the PBH-PBH component are in fact quite good, even for abundances as low
as $f^{\mathrm{cl}}_{\mathrm{PBH}}\simeq 10^{-1}$ in the most massive case
considered here. We stress again that the comparison between the benchmark
effective sensitivity curves and our GWB estimates is provided only to
visualize the relative importance of the different components, and to
understand which ones can provide a detectable signal. A complete study on the
detectability of these contributions, as well as the possibility of separate
them, is left for future work.

Figure 3: GWB generated by ABH-ABH (upper left panel), ABH-PBH (upper right
panel), and PBH-PBH (lower panel) binaries in GCs. In each panel solid and
dashed lines indicate $f_{\rm PBH}^{\rm cl}=1$ and $f_{\rm PBH}^{\rm
cl}=10^{-1}$ PBH abundances, respectively. Darker shades of red, green and
blue lines indicate heavier PBH masses in GCs. The gray dotted-dashed line
represents the forecasted sensitivity of the ET2CE GW detector network.

###  4.4 Signatures of primordial black hole binaries

We showed in the previous section that in many circumstances PBH binaries
formed in dense environments leave a detectable contribution in the GWB. It is
thus worth investigating whether this imprint is unique or not. Moreover, PBHs
form binaries via many different formation mechanisms, hence we necessarily
have to consider all of them for our analysis to be comprehensive. In
particular, although each source contributing to the GWB evolves independently
from the others, the joint set of channels has to be compatible with the local
GW merger rate $R^{\mathrm{LVK}}_{0}\simeq 20\,\mathrm{Gpc^{-3}\,yr^{-1}}$
measured by the LIGO-Virgo-KAGRA collaboration [144], i.e.,

$R^{\mathrm{LVK}}_{0}=R^{\mathrm{EPBH}}_{0}+R^{\mathrm{LPBH}}_{0}+R^{\mathrm{ABH}}_{0}+R^{\mathrm{GC}}_{0}\,,$
(4.13)

where the sum includes the early-time primordial, late-time primordial,
astrophysical, and dense environments contributions, respectively. The dense
environment contribution contains PBH-PBH, ABH-PBH, and ABH-ABH binaries.

The rates for early- and late-time PBH binary formation are computed as
discussed in sections 2.1 and 2.2. On the one hand, the late-time PBH binary
merger rate is independent of PBH mass (for the mass ranges considered here),
and it is given by $R_{0}^{\rm LPBH}(f_{\mathrm{PBH}}=1)=1.9\,{\rm
Gpc}^{-3}\,{\rm yr}^{-1}$ [52, 48, 27]. As for early binaries, even though
some initial theoretical estimates suggested large merger rates, recent
numerical results suggest that such values overestimate the real one [145,
146]. In a similar fashion to ref. [52], the parameter $\mathcal{A}_{m}$ is
defined for each fixed value of $M_{\rm PBH}$ in such a way that, when $f_{\rm
PBH}=1$, the total contribution from PBH binaries in the field constitutes a
certain fraction of the local GW merger rate [147]. As illustrative examples,
we consider the cases in which
$R^{\mathrm{EPBH}}_{0}+R^{\mathrm{LPBH}}_{0}=0.1\times R^{\mathrm{LVK}}_{0}$
and $R^{\mathrm{EPBH}}_{0}+R^{\mathrm{LPBH}}_{0}=0.3\times
R^{\mathrm{LVK}}_{0}$.888This value was suggested by some previous studies
[147], and we consider it an upper limit for the fraction of LVK events that
could be primordial. With these requirements, for $M_{\rm
PBH}=100\,M_{\odot}$, we set the early-time amplitude parameter to
$\mathcal{A}_{m}(f_{\mathrm{PBH}}=1)=12\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}$ for a
$30\%$ contribution, or $\mathcal{A}_{m}(f_{\mathrm{PBH}}=1)=0.37\,{\rm
Gpc}^{-3}\,{\rm yr}^{-1}$ for a $10\%$ contribution.

As we discussed in section 2, while the early- and late-time PBH binaries
merger rate depends on the $f_{\mathrm{PBH}}$ abundance parameter, the dense
environment contribution depends on the cluster parameter
$f^{\mathrm{cl}}_{\mathrm{PBH}}$. As we showed in section 3.3, as long as
$f_{\mathrm{PBH}}\gtrsim 2\times 10^{-5}$ the value of the parameter $f_{\rm
PBH}^{\rm cl}$ can be chosen independently from it. In other words, if PBHs
constitute a small fraction of the dark matter, the contribution of the late
and early formation channels to the GWB decreases; on the other hand, even if
the contribution of the binaries in the cluster is small, it remains constant
independently from the global $f_{\rm PBH}$, once that the abundance inside
the cluster is fixed through $f_{\rm PBH}^{\rm cl}$. Therefore, there are
regimes in which the relative contribution to the GWB of PBHs in dense
environments is more relevant than the one of the early and late formation
channels.

Figure 4: Ratio of the GWB from GC to the total GWB expected from early- and
late-time PBH binaries. The ratio is computed at the reference frequency
$f=50\,{\rm Hz}$ where the PBH backgrund is expected to peak for $M_{\rm
PBH}=100\,M_{\odot}$. The two panels showcase different scenarios, where the
total signal from early and late PBH binaries with $f_{\rm PBH}=1$ contributes
a $30\%$ (left panel) and $10\%$ (right panel) of the total merger rate. In
the left panel, $\mathcal{A}_{m}=12\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}$, and the
resulting early-time binaries local merger rate is $R_{0}^{\rm EPBH}(f_{\rm
PBH}=1)=4\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}$. In the right panel,
$\mathcal{A}_{m}=0.37\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}$, and $R_{0}^{\rm
EPBH}(f_{\rm PBH}=1)=0.1\,{\rm Gpc}^{-3}\,{\rm yr}^{-1}$. In both cases, the
late-binaries local merger rate is $R_{0}^{\rm LPBH}(f_{\rm PBH}=1)=2\,{\rm
Gpc}^{-3}\,{\rm yr}^{-1}$.

Figure 4 shows the relative importance of the GC contribution with respect to
the total PBH contribution from the field, with varying $f_{\rm PBH}$ and
$f_{\rm PBH}^{\rm cl}$. The figure shows the ratio of the GWB from GC to the
sum of the early- and late-time contributions, i.e., $\Omega^{\rm
GC}/\left(\Omega^{\rm EPBH}+\Omega^{\rm LPBH}\right)$, computed at the
reference frequency $f=50\,{\rm Hz}$ where the PBH background is expected to
peak for $M_{\rm PBH}=100\,M_{\odot}$. Depending on the value of both $f_{\rm
PBH}$ and $f_{\rm PBH}^{\rm cl}$, the GC contribution may be larger than the
one that is usually estimated in the literature from early- and late-time PBH
binaries.

Figure 5: GWB generated by early-time PBH (red lines), late-time PBH (green
lines), ABH (blue lines), and total GC (purple lines) binaries when the PBH
abundance is set to $f_{\rm PBH}=1$ (upper panel), $f_{\rm PBH}=10^{-1}$
(lower left panel), and $f_{\rm PBH}=10^{-2}$ (lower right panel). In all
cases PBHs have $M_{\mathrm{PBH}}=100\ M_{\odot}$, and we set
$\mathcal{A}_{m}$ following the procedure outlined in the main text. The
contribution from GCs assumes $f^{\mathrm{cl}}_{\mathrm{PBH}}=1$. The sum of
all contribution is represented by a black dashed line. The ET2CE effective
sensitivity curve for one year of observation time is showed as a dotted-
dashed gray line.

This relative importance must be considered in a wider scenario, that also
accounts for the other contributions to the total GWB. We show in figure 5 the
GWB generated by early-time, late-time, ABH, and GC binaries when
$M_{\mathrm{PBH}}=100\ M_{\odot}$ and $f_{\mathrm{PBH}}=1,10^{-1},10^{-2}$. In
the first case, the contribution from early- and late-time PBHs constitute a
$30\%$ of the total local binary merger rate – the other cases are rescaled
according to (2.2) and (2.3). While in the $f_{\rm PBH}=1$ case the sum of the
early- and late- binary contributions outnumbers the contribution PBHs inside
clusters, when $f_{\rm PBH}=10^{-1}$ the contributions are comparable, and for
$f_{\rm PBH}=10^{-2}$, PBH binaries inside dense environment contribute to the
GWB about one order of magnitude more than the other channels. These scenarios
showcase how even if $f_{\mathrm{PBH}}\ll 1$, the GC contribution can leave a
order tens-percent feature on the total GWB, especially in those cases where
PBHs have masses larger than those of ABHs999This conclusion relies on a
single effective mass to describe the ABH population. The implications of an
extended mass function will be investigated in a future work., since the
characteristic frequency of the merger scales as $f\propto M^{-1}$. In other
words, the presence of GWB features, which persist even in the case PBHs
contribute only to a fraction of the dark matter, can potentially provide an
indicator of more exotic binary formation channels, and deserves additional
scrutiny.

###  4.5 Hyperbolic encounters

Binary formation is not the only class of dynamical processes that take place
in dense environments; in fact, we are guaranteed to have also a second class
of events in GCs, i.e., GW bursts. When two unbound compact objects pass close
to each other, they emit gravitational radiation due to the reciprocal
deflection of their orbit, analogously to the bremsstrahlung emission in the
electromagnetic case [148]. In this section, we are interested in those cases
where the energy dissipated through GWs in the encounter is not sufficient to
create a bound pair, i.e., both initial and final eccentricities are larger
than unity. These events are called hyperbolic encounters ($\beta={\rm HE}$
hereafter), and they have been investigated for both ABHs [149, 150, 151, 152]
and PBHs [153, 154, 155]. Since multiple authors have provided slightly
different descriptions of the hyperbolic encounter dynamics and GW emission,
in appendix C we summarize the equations implemented in our code to model
these interactions.

The hyperbolic encounter cross section reads as [150]

$\sigma^{\rm(HE)}=\pi r^{2}_{\mathrm{ini}}\sin^{2}\theta_{\mathrm{ini}},$
(4.14)

where $r_{\mathrm{ini}}$ is the distance between BHs at the beginning of the
process, and
$\cos\theta_{\mathrm{ini}}=\hat{\mathbf{r}}_{\mathrm{ini}}\cdot\hat{\mathbf{v}}_{\mathrm{rel}}$
is the angle between the initial separation and relative velocity vectors. In
this case the velocity-averaged cross section reads as

$\left\langle\sigma^{(\mathrm{HE})}v_{\rm rel}\right\rangle=\int
dr_{\mathrm{ini}}d\theta_{\mathrm{ini}}dv_{\mathrm{rel}}p(r_{\mathrm{ini}},\theta_{\mathrm{ini}},v_{\mathrm{rel}})\sigma^{\rm(HE)}v_{\mathrm{rel}}\,,$
(4.15)

where we derive for the first time the probability distribution functions for
the initial distance, angle and relative velocity probability distribution
function. These are presented in appendix D.

In this scenario, it is more convenient to consider a single population of BHs
instead of two different ones, see the discussion in appendix D. Therefore,
instead of computing the effective mass for ABHs, we compute a global
effective mass $M^{\mathrm{eff}}_{\mathrm{HE}}$ for the ensemble of ABHs and
PBHs by solving the equation

$\Omega^{(\mathrm{HE})}\left(\frac{dn}{dM}=\frac{\rho_{0\rm ABH}+\rho_{0\rm
PBH}}{M_{\mathrm{HE}}^{\mathrm{eff}}}\delta^{D}(M-M_{\mathrm{HE}}^{\mathrm{eff}})\right)=\Omega^{(\mathrm{HE})}\left(\frac{dn_{\mathrm{tot}}}{dM}\right)\,,$
(4.16)

where

$\frac{dn_{\mathrm{tot}}}{dM}=\frac{dn_{\rm ABH}}{dM_{\rm ABH}}+\frac{dn_{\rm
PBH}}{dM_{\rm PBH}}=\frac{(2-\alpha)M^{-\alpha}\rho_{0\rm
ABH}}{(M^{\mathrm{max}}_{\mathrm{ABH}})^{2-\alpha}-(M^{\mathrm{min}}_{\mathrm{ABH}})^{2-\alpha}}+\frac{\rho_{0\rm
PBH}}{M_{\rm PBH}}\delta^{\rm D}(M-M_{\rm PBH})\,,$ (4.17)

and we use the fact that the GW strain scales as
$|h^{\mathrm{HE}}_{ij}|^{2}\propto\mu^{2}M^{4}$, as showed in equation (C.5),
for hyperbolic encounters.

Figure 6: Hyperbolic encounter GWB for PBH cluster abundances of $f_{\rm
PBH}^{\rm cl}=1$ (solid line) and $f_{\rm PBH}^{\rm cl}=10^{-1}$ (dashed
line). In all cases reported we assume $M_{\rm PBH}=100\,M_{\odot}$. The LISA
effective sensitivity curve for one year of observation time is showed as a
dotted-dashed gray line.

In this case, there is no time delay between the dynamical encounter and GW
emission, hence in the computation of the GWB we only need to implement the
interaction rate. This class of dynamical interactions is extremely common,
with interaction rates per cluster of order $\Gamma^{(\rm
HE)}_{\mathrm{tot}}(M_{\rm b}^{\rm ini}=10^{5}\,M_{\odot})\simeq 10^{-5}\,{\rm
yr}^{-1}$, approximately $10^{6}$ times larger than the direct capture
interaction rate per cluster. However, as showed in figure 6, the GWB produced
by hyperbolic encounters turns out to be many orders of magnitude below LISA
expected sensitivity. The GWB generated by lighter PBHs has a similar
magnitude to the ones reported in the figure. Strong GW bursts are generated
only by encounters with initial eccentricity very close to unity.
Unfortunately, that occurrence is very unlikely, as showed by the initial
eccentricity probability distribution functions in appendix D, explaining the
extremely weak GWB global signal. Since the effective sensitivity scales with
the square root of the observation time, we cannot foresee a detection of this
class of events even in the case of very massive PBHs.

##  5 Conclusions

In this work, we characterized the PBH contribution to the GWB, showing how
dense environments provide the perfect ground for dynamical interaction that
can ultimately be detected by current and next-generation GW observatories. In
particular, we showed that PBHs can leave a very characteristic imprint on the
GWB even if their abundance is well below the current upper bounds. Even in
the case where the well-known GWBs from early- and late-time PBH binaries are
below the sensitivity threshold of planned interferometers, PBHs in dense
environments can still provide a contribution to the background that is
comparable to the field-level ABH one, and therefore should not be neglected.
A more accurate analysis of the detectability of each individual background is
left for a future work currently in preparation.

In this work we modeled the internal structure of GCs, following what appears
to be the general consensus in the literature. However, we are well aware that
many properties of GCs are still an object of debate in the literature, and a
change in any of them could significantly boost/reduce the dynamical
interaction rate, and thus the GW signal. As an example, a more accurate
modeling of the signal, relying on numerical simulations, should also include
GC evaporation and the subsequent evolution of the mass function [156, 157,
158], that we neglected at this stage, as well as an extended mass function
for the ABHs. Nevertheless, we believe that semi-analytical models, as the
ones we propose, are still a valuable tool to explore the impact of
astrophysical parameters and assess which regions should be further
investigated numerically.

We showed that, if PBHs are present in the core of GCs, they have a
significant impact on their dynamical properties and internal structure,
causing the other populations to segregate within more or less diluted cores,
depending on whether they are lighter or heavier than PBHs. Therefore,
properly characterizing the GWB contribution from PBHs in dense environments
can also deepen our understanding of the internal structure of these
astrophysical systems, and on the amount of dark matter within their cores.

Moreover, even if in this paper we restricted our attention to GCs, we
acknowledge that our formalism can be readily applied to other kinds of
gravitationally bound system, such as nuclear clusters or even clusters made
exclusively of PBHs [159].

In an upcoming work in preparation, we will study the possibility of
disentangling the different channels and cases, so to investigate their
detectability and the cosmological tests that this observable will enable.

## Acknowledgments

The authors thank Konstantinos Kritos, Chiara Mingarelli, Emanuele Berti for
useful discussion.

SL thanks the John Hopkins University for the hospitality during the final
stages of this work. The visit was supported by a Balzan awarded scholarship,
funded by the Balzan Centre for Cosmological Studies. SL is supported by an
Azrieli International Postdoctoral Fellowship. NB is supported by PRD/ARPE
2022 “Cosmology with Gravitational waves and Large Scale Structure -
CosmoGraLSS” and acknowledges partial support from the National Science
Foundation (NSF) under Grant No. PHY-2112884. AR acknowledges funding from the
Italian Ministry of University and Research (MIUR) through the “Dipartimenti
di eccellenza” project “Science of the Universe”. L. V. acknowledges financial
support from the Supporting TAlent in Re-Search@University of Padova
(STARS@UNIPD) for the project “Constraining Cosmology and Astrophysics with
Gravitational Waves, Cosmic Microwave Background and Large-Scale Structure
cross-correlations”.

##  Appendix A Timescales

In this appendix we investigate different timescales of the systems of
interest, i.e., GCs and binaries, and we discuss their relevance for the
purpose of estimating the GWB generated by compact objects in dense
environments.

##### Stellar evolution timescale.

Even if at the time of formation we assume the GC to have exclusively a
stellar population, massive stars are expected to rapidly evolve and form a
population of ABHs. In other words, the massive star and ABH mass functions
are both evolving in time, with the former depleted in favor of the latter. We
estimate the ABH formation time as the minimum time required for a massive
star to burn Hydrogen, Helium, Carbon and all the other elements in its core,
i.e., $t_{\mathrm{ABH}}=t_{\mathrm{H}}+t_{\mathrm{He,C,...}}$, where the
Hydrogen burning time is given by [85]

$t_{\mathrm{H}}\simeq
35.3\left(\frac{8\,M_{\odot}}{M_{\star}}\right)^{\alpha_{\mathrm{H}}-1}\
\mathrm{Myr}\,,$ (A.1)

and the scaling exponent is $\alpha_{\mathrm{H}}=3.8$. The burning time of
heavier elements is of order $t_{\mathrm{He,C,...}}\simeq 0.2\times
t_{\mathrm{H}}$ [85]. Therefore, we find that the minimum ABH mass inside in a
GC is

$M^{\mathrm{min}}_{\mathrm{ABH}}=\left\\{\begin{aligned} &40\ M_{\odot}&\quad
t\leq t_{\mathrm{ABH}}(40\ M_{\odot})\,,\\\ &8\left[\frac{t_{\mathrm{ABH}}(8\
M_{\odot})}{t}\right]^{\frac{1}{\alpha_{\mathrm{H}}-1}}\ M_{\odot}&\quad
t_{\mathrm{ABH}}(40\ M_{\odot})<t<t_{\mathrm{ABH}}(8\ M_{\odot})\,,\\\ &8\
M_{\odot}&\quad t\geq t_{\mathrm{ABH}}(8\ M_{\odot})\,.\\\
\end{aligned}\right.$ (A.2)

Assuming that the GC and all stars inside it form at the same time, we have
that after a typical maximum timescale of order $t_{\mathrm{ABH}}\simeq
42\,\mathrm{Myr}$ the ABH population is fully formed. Given that this
timescale is only a fraction of the typical segregation timescale (e.g.,
$t_{\mathrm{ABH}}/t_{\mathrm{rel}}\simeq 10^{-1}$, see following paragraphs),
we can safely neglect the modelling of this transition time in our estimates.

##### Time delay between binary formation and merger.

Mergers of BH binaries typically do not happen immediately after binary
formation. Instead, binaries present a typical time delay distribution
determined by the initial properties of the system and of the environment.
Regarding direct capture processes, we follow ref. [160] and describe the
initial state of a binary using its initial semi-major axis and eccentricity

$a_{\mathrm{ini}}=-\frac{GM\mu}{2E_{\rm ini}},\qquad
e_{\mathrm{ini}}=\sqrt{1+2\frac{E_{\rm ini}b^{2}v_{\rm
rel}^{2}}{G^{2}M^{2}\mu}}\,,$ (A.3)

respectively, where $b$ is the impact parameter, $E_{\mathrm{ini}}$ is the
energy of the binary at the beginning of the inspiral phase and it reads as

$E_{\rm ini}=\frac{\mu v_{\rm rel}^{2}}{2}-\frac{85\pi
G^{7/2}}{12\sqrt{2}c^{5}}\frac{\mu^{2}M^{5/2}}{r_{p}^{7/2}}\,,$ (A.4)

and $r_{p}$ being the distance of closest approach, given by

$r_{p}\simeq\frac{b^{2}v_{\rm rel}^{2}}{2MG}\left(1-\frac{b^{2}v_{\rm
rel}^{4}}{4M^{2}G^{2}}\right)\,.$ (A.5)

Therefore, the delay time between binary formation and merger reads as [161]

$t_{\rm
d}=\frac{3}{85}\frac{a_{\mathrm{ini}}^{4}c^{5}}{G^{3}M^{2}\mu}\left(1-e_{\mathrm{ini}}^{2}\right)^{7/2}\,.$
(A.6)

We create a mock catalogue of binaries formed through direct capture process
and we derive their time delay distribution assuming that the impact parameter
is uniformly distributed between its minimum and maximum values

$\displaystyle b_{\mathrm{min}}$
$\displaystyle=\sqrt{12}\frac{GM}{cv_{\mathrm{rel}}}\,,\qquad
b_{\mathrm{max}}=\left(\frac{340\pi}{3}\right)^{1/7}\frac{GM^{6/7}\mu^{1/7}}{c^{2}}\left(\frac{v_{\mathrm{rel}}}{c}\right)^{-9/7}\,,$
(A.7)

respectively.

Figure 7: Time delay probability distribution for ABH-ABH (red), ABH-PBH
(green) and PBH-PBH (blue) binaries for a GC with
$M^{\mathrm{ini}}_{\mathrm{b}}=10^{5}\,M_{\odot}$,
$M_{\mathrm{PBH}}=10\,M_{\odot}$, and $f^{\mathrm{cl}}_{\mathrm{PBH}}=1$ (left
panel) or $f^{\mathrm{cl}}_{\mathrm{PBH}}=10^{-1}$ (right panel). In the case
of $f^{\mathrm{cl}}_{\mathrm{PBH}}=1$ all the binaries have time delay at most
of $1\,\mathrm{Myr}$, whereas in the $f^{\mathrm{cl}}_{\mathrm{PBH}}=10^{-1}$
case a negligible fraction of events has also larger time delays ($0.03\%$,
$0.01\%$ and $0.05\%$ for ABH-ABH, ABH-PBH and PBH-PBH, respectively).

We show the recovered time delay distributions in figure 7, and we observe
that for all classes of compact objects (i) the time delay distribution scales
approximately as $p(t_{\rm d})\propto t_{\rm d}^{-1}$, and (ii) the maximum
time delay is $\mathcal{O}(1)\,\mathrm{Myr}$; therefore, since the average
time delay is $\left\langle t_{\rm d}\right\rangle\ll 1\,\mathrm{Myr}$, we
conclude that we can safely ignore this time delay in our modelling of the
merger rate. On the other hand, regarding the three-body encounters binary
formation channel, we are not aware of any estimate on their expected time
delay distribution. However, we do not foresee the three-body time delay
distribution to be significantly different from the direct capture one, since
compact object participating in both kinds of interaction effectively populate
the same environment. Thus, also in the case of three-body encounters, we
neglect the time delay between binary formation and merger.

##### Mass segregation timescale.

The initial state of GCs are currently very uncertain, in particular little is
known about whether they form already mass-segregated or not, see, e.g., refs.
[162, 163] for a discussion on young globular clusters. If they do [164, 165],
density and velocity profiles are expected to be well described by the phase
space distribution in equation (3.6) during almost all the GC lifetime. In the
opposite case, the picture presented in section 3.2 is valid only after the GC
relaxed into its “equilibrium” configuration, hence after a typical timescale
from formation of the order of the relaxation time $t_{\mathrm{rel}}$ defined
in equation (3.10). In this work, we consider GCs that are not mass segregated
at formation time. However, since the typical relaxation time for GCs in the
mass range of our interest is of order $\mathcal{O}(300)\,\mathrm{Myr}$, and
it is significantly lower than the typical timescale over which binaries form
(see following paragraph), we neglect this effect in the modelling of the GW
emission. Note that even if very massive GC have a relaxation timescale of
Gyrs, the steepness of the GC mass function render their contribution to the
GWB subdominant. Therefore, we do not expect this approximation to
significantly alter our estimations.

##### Merger timescale.

The mass evolution of the GC component might be affected by the very same
merger events we are interested into. If mergers happen frequently, i.e., if

$H/\Gamma^{(\beta)}_{ij,\rm tot}\ll 1\,,$ (A.8)

the low-mass tail of the ABH and/or PBH mass function will get depleted of
compact objects in favor of the high-mass tail. Even if our initial assumption
was that the mass function is constant in time, we check a posteriori that
$H_{0}/\Gamma^{(\beta)}_{ij,\rm tot}\gtrsim 1$ in all binary formation
channels and for all combinations of sources. Only in the case of GCs with
$M_{\mathrm{PBH}}=100\,M_{\odot}$ and $f^{\mathrm{cl}}_{\mathrm{PBH}}=1$ we
observe $H_{0}/\Gamma^{(\mathrm{DC})}_{\rm PBH-
PBH,tot}\simeq\mathcal{O}(\mathrm{few})$ when
$M^{\mathrm{ini}}_{\mathrm{b}}\gtrsim 10^{6}\,M_{\odot}$, i.e., very massive
GCs might experience few PBH-PBH binary mergers per Hubble time. Since the GC
mass distribution is dominated by GCs with
$M^{\mathrm{ini}}_{\mathrm{b}}\simeq 10^{5}\,M_{\odot}$, we can safely
conclude that mergers do not significantly alter the compact object mass
function.

##  Appendix B Relative velocity probability distribution function

Cross sections of dynamical interactions described in section 4 depend on the
relative velocity of the compact objects involved in the process. In this
appendix we derive an approximate analytical expression for the relative
velocity probability distribution function starting from the velocity
distribution of BHs in the GC cores.

Figure 8: King and truncated Maxwell-Boltzmann (MB) ABH and PBH velocity
probability distribution functions in the core of GCs with
$M^{\mathrm{ini}}_{\mathrm{b}}=10^{5}\,M_{\odot}$,
$M_{\mathrm{PBH}}=10\,M_{\odot}$, and $f^{\mathrm{cl}}_{\mathrm{PBH}}=1$ (left
panel) or $f^{\mathrm{cl}}_{\mathrm{PBH}}=10^{-1}$ (right panel). Velocities
on the x-axis are rescaled by each respective population velocity scale.

We compare in figure 8 the velocity probability distribution function of
individual BH populations given by equation (3.6) to a truncated Maxwell-
Boltzmann distribution

$p(\boldsymbol{v}_{i})=\frac{\mathcal{N}_{i}}{(2\pi\sigma^{2}_{i})^{3/2}}e^{-v_{i}^{2}/2\sigma^{2}_{i}}\,,$
(B.1)

where the normalization factor reads as

$\mathcal{N}^{-1}_{i}=\frac{1}{(2\pi)^{3/2}\sigma^{3}_{i}}\int
d\Omega_{i}\int_{0}^{v_{\mathrm{esc}}}dv_{i}\
v^{2}_{i}e^{-v_{i}^{2}/2\sigma^{2}_{i}}=\mathrm{erf}\left(\frac{v_{\mathrm{esc}}}{\sqrt{2}\sigma_{i}}\right)-\sqrt{\frac{2}{\pi}}\frac{v_{\mathrm{esc}}}{\sigma_{i}}e^{-v^{2}_{\mathrm{esc}}/2\sigma^{2}_{i}}\,.$
(B.2)

As expected, since the typical velocity dispersion is smaller than the escape
velocity, i.e., $\sigma_{j}\ll v_{\mathrm{esc}}=\sqrt{2\Psi_{0}}$ for each
population of interest, we find that the truncated Maxwell-Boltzmann nicely
approximates the real phase space distributions. Since the difference between
phase space distributions appears to be relatively small, in the following we
adopt the truncated Maxwell-Boltzmann for our analytical estimates.

Under the change of variables defined by

$\boldsymbol{v}_{\mathrm{rel}}=\boldsymbol{v}_{2}-\boldsymbol{v}_{1}\,,\qquad\boldsymbol{v}_{c}=\frac{\sigma^{2}_{2}\boldsymbol{v}_{1}+\sigma^{2}_{1}\boldsymbol{v}_{2}}{\sigma^{2}_{\mathrm{rel}}}\,,$
(B.3)

we have that the original velocities read as

$\boldsymbol{v}_{1}=\boldsymbol{v}_{c}-\frac{\sigma^{2}_{1}}{\sigma^{2}_{\mathrm{rel}}}\boldsymbol{v}_{\mathrm{rel}}\,,\qquad\boldsymbol{v}_{2}=\boldsymbol{v}_{c}+\frac{\sigma^{2}_{2}}{\sigma^{2}_{\mathrm{rel}}}\boldsymbol{v}_{\mathrm{rel}}\,,\\\
$ (B.4)

with $\sigma^{2}_{\mathrm{rel}}=\sigma^{2}_{1}+\sigma^{2}_{2}$ and
$\sigma^{2}_{c}=\sigma^{2}_{1}\sigma^{2}_{2}/(\sigma^{2}_{1}+\sigma^{2}_{2})$.
Since the velocities of the two objects are bound to be lower than the escape
velocity $v_{\rm esc}$, the domain of the relative velocity is
$v_{\mathrm{rel}}\in[0,2v_{\mathrm{esc}}]$ while $v_{\rm c}$ has a “circular”
domain, i.e.,

$v^{2}_{c}\leq
v^{2}_{c,\mathrm{max}}(v_{\mathrm{rel}})=\frac{\sigma^{2}_{c}}{\sigma^{2}_{\mathrm{rel}}}\left[v^{2}_{\mathrm{esc}}\left(2+\frac{\sigma^{2}_{2}}{\sigma^{2}_{1}}+\frac{\sigma^{2}_{1}}{\sigma^{2}_{2}}\right)-v^{2}_{\mathrm{rel}}\right]\,.$
(B.5)

The integration domain of $\cos\theta_{c}=v_{c,z}/v_{c}$ is also restricted.
In particular, starting from the definitions of $\cos\theta_{c}$ and
$\cos\theta_{\mathrm{rel}}=v_{\mathrm{rel},z}/v_{\mathrm{rel}}$, we have

$\sigma^{2}_{\mathrm{rel}}v_{c}\cos\theta_{c}=\sigma^{2}_{1}v_{2}\cos\theta_{2}+\sigma^{2}_{2}v_{1}\cos\theta_{1}\,,\qquad\sigma^{2}_{2}v_{\mathrm{rel}}\cos\theta_{\mathrm{rel}}=\sigma^{2}_{2}v_{2}\cos\theta_{2}-\sigma^{2}_{2}v_{1}\cos\theta_{1}\,.$
(B.6)

By summing these two equations and later squaring them, we find

$\cos^{2}\theta_{c}+2\frac{\sigma^{2}_{2}}{\sigma^{2}_{\mathrm{rel}}}\frac{v_{\mathrm{rel}}}{v_{c}}\cos\theta_{\mathrm{rel}}\cos\theta_{c}+\frac{\sigma^{4}_{2}}{\sigma^{4}_{\mathrm{rel}}}\frac{v^{2}_{\mathrm{rel}}}{v^{2}_{c}}\cos^{2}\theta_{\mathrm{rel}}=\frac{v^{2}_{2}}{v^{2}_{c}}\cos^{2}\theta_{2}\leq\frac{v^{2}_{\mathrm{esc}}}{v^{2}_{c}}\,,$
(B.7)

which admits solutions for $\cos\theta_{c}$ in the range

$-\frac{\sigma^{2}_{2}}{\sigma^{2}_{\mathrm{rel}}}\frac{v_{\mathrm{rel}}}{v_{c}}\cos\theta_{\mathrm{rel}}-\frac{v_{\mathrm{esc}}}{v_{c}}\leq\cos\theta_{c}\leq-\frac{\sigma^{2}_{2}}{\sigma^{2}_{\mathrm{rel}}}\frac{v_{\mathrm{rel}}}{v_{c}}\cos\theta_{\mathrm{rel}}+\frac{v_{\mathrm{esc}}}{v_{c}}\,.$
(B.8)

When $v_{\mathrm{esc}}\gg\sigma_{1},\sigma_{2}$, as in our cases of interest,
we have $-1\lesssim\cos\theta_{c}\lesssim 1$, and the integration domain
becomes approximately unconstrained.

Figure 9: Relative velocity probability distribution function obtain from a
catalog (histograms) compared to the analytic prediction of equation (B.10)
(solid lines) for ABH-ABH (red), ABH-PBH (green) and PBH-PBH (blue) pairs in a
GC with $M^{\mathrm{ini}}_{\mathrm{b}}=10^{5}\,M_{\odot}$,
$M_{\mathrm{PBH}}=10\,M_{\odot}$, and $f^{\mathrm{cl}}_{\mathrm{PBH}}=1$ (left
panel) or $f^{\mathrm{cl}}_{\mathrm{PBH}}=10^{-1}$ (right panel).

The probability distribution function transform under such change of variables
as

$p(\boldsymbol{v}_{1})p(\boldsymbol{v}_{2})d^{3}v_{1}d^{3}v_{2}=p(\boldsymbol{v}_{c},\boldsymbol{v}_{\mathrm{rel}})d^{3}v_{c}d^{3}v_{\mathrm{rel}}=\frac{\mathcal{N}_{1}\mathcal{N}_{2}}{(2\pi)^{3}\sigma^{3}_{c}\sigma^{3}_{\mathrm{rel}}}e^{-v^{2}_{c}/2\sigma^{2}_{c}}e^{-v^{2}_{\mathrm{rel}}/2\sigma^{2}_{\mathrm{rel}}}d^{3}v_{c}d^{3}v_{\mathrm{rel}}\,,$
(B.9)

therefore, when $v_{\mathrm{esc}}\gg\sigma_{1},\sigma_{2}$, we find that the
relative velocity probability distribution function $p(v_{\mathrm{rel}})$
reads as

$\displaystyle p(v_{\mathrm{rel}})$
$\displaystyle=\mathcal{N}_{1}\mathcal{N}_{2}v^{2}_{\mathrm{rel}}\int
d\Omega_{\mathrm{rel}}\int_{0}^{v_{c,\mathrm{max}}(v_{\mathrm{rel}})}d\Omega_{c}dv_{c}v^{2}_{c}\frac{1}{(2\pi)^{3}\sigma^{3}_{c}\sigma^{3}_{\mathrm{rel}}}e^{-v^{2}_{c}/2\sigma^{2}_{c}}e^{-v^{2}_{\mathrm{rel}}/2\sigma^{2}_{\mathrm{rel}}}$
(B.10)
$\displaystyle=\mathcal{N}_{1}\mathcal{N}_{2}\frac{(4\pi)^{2}v^{2}_{\mathrm{rel}}e^{-v^{2}_{\mathrm{rel}}/2\sigma^{2}_{\mathrm{rel}}}}{(2\pi)^{3}\sigma^{3}_{c}\sigma^{3}_{\mathrm{rel}}}\int_{0}^{v_{c,\mathrm{max}}(v_{\mathrm{rel}})}dv_{c}v^{2}_{c}e^{-v^{2}_{c}/2\sigma^{2}_{c}}$
$\displaystyle=\mathcal{N}_{1}\mathcal{N}_{2}\frac{(4\pi)^{2}v^{2}_{\mathrm{rel}}e^{-v^{2}_{\mathrm{rel}}/2\sigma^{2}_{\mathrm{rel}}}}{(2\pi)^{3}\sigma^{3}_{\mathrm{rel}}}\left[\sqrt{\frac{\pi}{2}}\mathrm{erf}\left(\frac{v_{c,\mathrm{max}}(v_{\mathrm{rel}})}{\sqrt{2}\sigma_{c}}\right)-\frac{v_{c,\mathrm{max}}(v_{\mathrm{rel}})}{\sigma_{c}}e^{-v^{2}_{c,\mathrm{max}}(v_{\mathrm{rel}})/2\sigma^{2}_{c}}\right]\,.$

In figure 9 we show the comparison between the relative velocity distribution
given by the theoretical prediction of equation (B.10), and by generating a
catalog of events in which the velocities of the two bodies are separately
sampled from the individual Maxwell-Boltzmann distributions. For both
benchmark cases, we find that our analytic formula is accurate at percent
level, with deviations appearing only in the high relative velocity tail of
the probability distribution.

##  Appendix C Hyperbolic encounter dynamics

In this appendix we briefly summarize the dynamics of the hyperbolic encounter
between two bodies, mainly following refs. [166, 167]. In the rest frame of
the second body, the coordinates of the first one are parametrized as

$x(u)=a(e-\cosh u)\,,\qquad y(u)=a(e^{2}-1)^{1/2}\sinh u\,,$ (C.1)

where $a=GM/v^{2}_{\mathrm{rel}}$ is the semi-major axis, $e$ is the
eccentricity and $u$ is the anomaly. The latter implicitly parametrizes also
the time evolution of the system as

$t(u)=\sqrt{\frac{a^{3}}{GM}}(e\sinh u-u)\,,$ (C.2)

with $t_{0}=GM/v^{3}_{\mathrm{rel}}$ being a characteristic time of the
gravitational interaction. At time $t=0$ (corresponding to $u=0$) the two
bodies are at the minimal distance $r_{\mathrm{min}}=a(e-1)$.

The power emitted in GWs reads as

$P(t)=\frac{G}{45c^{5}}\partial^{3}_{t}D_{ij}\partial^{3}_{t}D_{ij}\,,$ (C.3)

where the quadrupole moment tensor is given by
$D_{ij}=3M_{ij}-\delta_{ij}M_{kk}$ and the second mass moment is defined as

$M_{ij}=\frac{1}{c^{2}}\int d^{3}x\ T^{00}(\boldsymbol{x})x_{i}x_{j}\,.$ (C.4)

In the simple case at hand we have $T^{00}(\boldsymbol{x})=\mu
c^{2}\delta^{D}\left(\boldsymbol{x}-\boldsymbol{x}(t)\right)$, therefore the
energy emitted though GWs during the encounter reads as

$\Delta
E=\int_{0}^{\infty}df\frac{dE}{df}=\frac{128\pi^{6}}{45}\frac{G^{5}\mu^{2}M^{4}}{c^{5}v^{8}_{\mathrm{rel}}}\int_{0}^{\infty}dff^{4}\mathcal{F}(2\pi
t_{0}f)\,,$ (C.5)

where

$\displaystyle\mathcal{F}(x)$
$\displaystyle=\left|\frac{3(e^{2}-1)}{e}H_{ix}^{(1)^{\prime}}(ixe)+\frac{e^{2}-3}{e^{2}}\frac{i}{x}H_{ix}^{(1)}(ixe)\right|^{2}$
(C.6)
$\displaystyle\qquad+\left|\frac{3(e^{2}-1)}{e}H_{ix}^{(1)^{\prime}}(ixe)+\frac{2e^{2}-3}{e^{2}}\frac{i}{x}H_{ix}^{(1)}(ixe)\right|^{2}$
$\displaystyle\qquad\qquad+\frac{18(e^{2}-1)}{e^{2}}\left|\frac{1}{x}H_{ix}^{(1)^{\prime}}(ixe)+\frac{e^{2}-1}{e}iH_{ix}^{(1)}(ixe)\right|^{2}\,,$

and $H^{(1)}_{\alpha}(\gamma)$ and
$H^{(1)^{\prime}}_{\alpha}(\gamma)=\frac{1}{2}\left[H^{(1)}_{\alpha-1}(\gamma)-H^{(1)}_{\alpha+1}(\gamma)\right]$
are Hankel functions of the first kind and their derivatives with respect to
the argument, respectively.

Depending on the value of $x$, both the order and the argument of the Hankel
function can be considerably large. While for small orders ($x\lesssim 100$)
it is still feasible to use the exact numerical result, for higher orders the
computation time becomes exceedingly long and we resort to analytical
approximations. First of all we introduce a new parameter $\beta$ via
$e=\cos^{-1}\beta$. If $x\tan^{3}\beta\gtrsim 1$ we use the standard
approximation for Hankel functions with argument greater than the order

$H^{(1)}_{ix}\left(\frac{ix}{\cos\beta}\right)\simeq-i\left(\frac{2}{\pi
x\tan\beta}\right)^{1/2}e^{-x(\tan\beta-\beta)}\sum_{n=0}^{\infty}\frac{\Gamma\left(n+\frac{1}{2}\right)}{\Gamma\left(\frac{1}{2}\right)}A_{n}(\beta)\left(\frac{-2}{x\tan\beta}\right)^{n}\,,$
(C.7)

where we checked that a good accuracy can be achieved stopping at $n=3$.
Coefficients in this approximation scheme read as

$\displaystyle A_{0}(\beta)$ $\displaystyle=1\,,\qquad
A_{1}(\beta)=\frac{1}{8}+\frac{5}{24\tan^{2}\beta}\,,$ (C.8) $\displaystyle
A_{2}(\beta)$
$\displaystyle=\frac{3}{128}+\frac{77}{576\tan^{2}\beta}+\frac{385}{3456\tan^{4}\beta}\,,$
$\displaystyle A_{3}(\beta)$
$\displaystyle=\frac{5}{1024}+\frac{1521}{25600\tan^{2}\beta}+\frac{17017}{138240\tan^{4}\beta}+\frac{17017}{248832\tan^{6}\beta}\,.$

Instead, if $x\tan^{3}\beta\lesssim 1$ we use the approximation for when order
and argument are nearly equal. More specifically, if
$x^{2/3}(1-\cos\beta)/(\cos\beta)^{2/3}\leq 1$ we use

$\displaystyle H^{(1)}_{ix}\left(\frac{ix}{\cos\beta}\right)$
$\displaystyle\simeq-\frac{2}{3\pi}\sum_{n=0}^{\infty}e^{\frac{2}{3}(n+1)\pi
i}\sin\left(\frac{n+1}{3}\pi\right)\Gamma\left(\frac{n+1}{3}\right)$ (C.9)
$\displaystyle\qquad\qquad\qquad\times\left(\frac{6\cos\beta}{ix}\right)^{(n+1)/3}B_{n}\left(\frac{ix(1-\cos\beta)}{\cos\beta}\right)\,,$

where in this case we stop at $n=5$ and the $B_{n}$ coefficients read as

$\displaystyle B_{0}(y)$ $\displaystyle=1\,,\quad B_{1}(y)=y\,,\quad
B_{2}(y)=\frac{y^{2}}{2}-\frac{1}{20}\,,\quad
B_{3}(y)=\frac{y^{3}}{6}-\frac{y}{15}\,,$ (C.10) $\displaystyle\quad B_{4}(y)$
$\displaystyle=\frac{y^{4}}{24}-\frac{y^{2}}{24}+\frac{1}{280}\,,\quad
B_{5}(y)=\frac{y^{5}}{120}-\frac{y^{3}}{60}+\frac{43y}{8400}\,.$

Otherwise we use the approximation

$H^{(1)}_{ix}\left(\frac{ix}{\cos\beta}\right)\simeq-i\frac{\tan\beta}{\sqrt{3}}e^{-x\left(\tan\beta-\beta-\frac{1}{3}\tan^{3}\beta\right)+\frac{2}{3}\pi
i}H^{(1)}_{1/3}\left(-\frac{x\tan^{3}\beta}{3}\right)\,.$ (C.11)

##  Appendix D Probability distribution functions for hyperbolic encounters

In this appendix we provide a semi-analytical estimate of
$p(r_{\mathrm{ini}},\theta_{\mathrm{ini}},v_{\mathrm{rel}})$, the probability
distribution function required to model the hyperbolic encounters in section
4.5. As a first approximation, we assume that the three stochastic variables
are statistically independent, thus the total distribution factorizes as
$p(r_{\mathrm{ini}},\theta_{\mathrm{ini}},v_{\mathrm{rel}})=p(r_{\mathrm{ini}})p(\theta_{\mathrm{ini}})p_{\mathrm{tot}}(v_{\mathrm{rel}})$.

Regarding the BH separation probability distribution function, instead of
considering the ABH and PBH populations separately, we just consider a BH
population with density

$n_{\mathrm{tot}}(r)\simeq\left\\{\begin{aligned}
&n_{\mathrm{ABH}}+n_{\mathrm{PBH}}\,,\qquad\qquad&r\leq
r_{\mathrm{0,PBH}}\,,\\\
&n_{\mathrm{ABH}}\,,\qquad\qquad&r_{\mathrm{0,PBH}}<r\leq
r_{\mathrm{0,ABH}}\,,\\\
&0\qquad\qquad&r_{\mathrm{0,ABH}}<r\,.\end{aligned}\right.$ (D.1)

Following the standard approach for an ideal gas of particles [168], we derive
the nearest neighbor distribution of BHs in the cluster core, i.e., the
initial distance distribution, which reads as

$p(r_{\mathrm{ini}})=\frac{3r^{2}_{\mathrm{ini}}\mathcal{N}_{r_{\mathrm{ini}}}}{a^{3}}e^{-(r_{\mathrm{ini}}/a)^{3}}\,,$
(D.2)

where $a=\left[3/(4\pi)(n_{\mathrm{ABH}}+n_{\mathrm{PBH}})\right]^{1/3}$ is
the Wigner-Seitz radius, and $\mathcal{N}_{r_{\mathrm{ini}}}$ is a
normalization constant obtained by restricting our region of interest to the
inner part of the cluster core, i.e., allowing only
$r_{\mathrm{ini}}\in[r_{\mathrm{min}},r_{\mathrm{max}}]$. The minimum radius
distance is defined as the one at which the dynamics of the two bodies is not
dominated anymore by the cluster, i.e., when the potential energies satisfy
$G(M_{1}+M_{2})/r_{\mathrm{min}}=\Psi_{0}$. On the other hand, we arbitrarily
define as maximum distance $r_{\mathrm{max}}=\bar{r}/\alpha_{\mathrm{cut-
off}}$, where $\bar{r}=a\times\mathrm{Gamma}(4/3)$ is the average distance
between BHs, $\mathrm{Gamma}(x)$ is the Gamma function and
$\alpha_{\mathrm{cut-off}}$ is a numerical factor that implicitly defines a
cut-off radius. We consider as threshold value $\alpha_{\mathrm{cut-off}}=10$
to ensure that there is less than $0.1\%$ probability of having a third BHs in
the same volume.

The total relative velocity probability distribution function can be obtained
by the applying Bayes theorem to the results derived in appendix B. Accounting
for a mixed population of ABHs and PBHs with two different relative abundances
and velocity dispersion, we find

$p_{\mathrm{tot}}(v_{\mathrm{rel}})=\pi^{2}_{\mathrm{ABH}}\,p_{\mathrm{ABH-
ABH}}(v_{\mathrm{rel}})+2\pi_{\mathrm{ABH}}\,\pi_{\mathrm{PBH}}\,p_{\mathrm{ABH-
PBH}}(v_{\mathrm{rel}})+\pi^{2}_{\mathrm{PBH}}\,p_{\mathrm{PBH-
PBH}}(v_{\mathrm{rel}})\,,$ (D.3)

where $\pi_{j}=N_{j}/N_{\mathrm{tot}}=n_{j}/n_{\mathrm{tot}}$ are relative
abundance probabilities.

Finally, since the vectors $\boldsymbol{r}_{\mathrm{ini}}$,
$\boldsymbol{v}_{1}$ and $\boldsymbol{v}_{2}$ are expected to be uniformly
distributed, it is straightforward to compute the initial angle
$\theta_{\mathrm{ini}}$ distribution. However, we have to include the
additional restriction that
$\cos\theta_{\mathrm{ini}}=\hat{\boldsymbol{v}}_{\mathrm{rel}}\cdot\hat{\boldsymbol{r}}<0$,
since in the opposite case the two bodies have already passed each other. As
for the case of the total relative velocity, the total initial angle
distribution reads as

$p_{\mathrm{tot}}(\theta_{\mathrm{ini}})=\pi^{2}_{\mathrm{ABH}}p_{\mathrm{ABH-
ABH}}(\theta_{\mathrm{ini}})+2\pi_{\mathrm{ABH}}\pi_{\mathrm{PBH}}p_{\mathrm{ABH-
PBH}}(\cos\theta_{\mathrm{ini}})+\pi^{2}_{\mathrm{PBH}}p_{\mathrm{PBH-
PBH}}(\theta_{\mathrm{ini}})\,.$ (D.4)

In order to ensure that the speed of the first BH remains subluminal we have
to impose the constrant

$r_{\mathrm{ini}}\sin\theta_{\mathrm{ini}}\gtrsim\frac{(e_{\mathrm{ini}}+1)^{3/2}}{2(e_{\mathrm{ini}}-1)^{1/2}}r_{S}\,,$
(D.5)

with $r_{S}=2GM/c^{2}$ being the sum of the Schwarzschild radii of the two
BHs, and $e_{\mathrm{ini}}$ is the initial eccentricity given by

$e_{\mathrm{ini}}=\left[1+\frac{v^{4}_{\mathrm{rel}}r^{2}_{\mathrm{ini}}\sin^{2}\theta_{\mathrm{ini}}}{G^{2}M^{2}}\right]^{1/2}\,.$
(D.6)

Since $r_{S}/r_{\mathrm{ini}}\ll 1$ is almost given by construction, we have
that the bound is non-trivial only for values of the initial eccentricity very
close to unity. In that case, by using the explicit expression for the initial
eccentricity, we obtain the bound

$\sin\theta_{\mathrm{ini}}\gtrsim\frac{r_{S}/r_{\mathrm{ini}}}{v_{\mathrm{rel}}/c}\,.$
(D.7)

Moreover, we have to ensure that the system remains unbound. The energy loss
due to GWs is

$\Delta
E=-\frac{8\mu^{2}v^{7}_{\mathrm{rel}}}{15Mc^{5}}g(e_{\mathrm{ini}})\,,$ (D.8)

where

$g(e_{\mathrm{ini}})=\frac{24\arccos\left(-\frac{1}{e_{\mathrm{ini}}}\right)\left(1+\frac{73}{24}e^{2}_{\mathrm{ini}}+\frac{37}{96}e^{4}_{\mathrm{ini}}\right)+\sqrt{e^{2}_{\mathrm{ini}}-1}\left(\frac{301}{6}+\frac{673}{12}e^{2}_{\mathrm{ini}}\right)}{(e^{2}_{\mathrm{ini}}-1)^{7/2}}\,.$
(D.9)

Since the energy of the system at initial time is
$E_{\mathrm{ini}}=\frac{1}{2}\mu v^{2}_{\mathrm{rel}}$, the condition to
impose to have $E_{\mathrm{fin}}>0$ is

$1-\frac{16}{15}\frac{\mu}{M}\left(\frac{v_{\mathrm{rel}}}{c}\right)^{5}g(e_{\mathrm{ini}})\gtrsim
0\,.$ (D.10)

Away from the $e_{\mathrm{ini}}\to 1$ limit, the condition is very easily
satisfied because of the relative velocity suppression. However, when the
eccentricity value gets close to unity, we find the extra condition

$\sin\theta_{\mathrm{ini}}\gtrsim\left(\frac{340\pi}{384}\right)^{1/7}\left(\frac{\mu}{M}\right)^{1/7}\frac{r_{S}/r_{\mathrm{ini}}}{(v_{\mathrm{rel}}/c)^{9/7}}\,,$
(D.11)

which turns out to be stronger than the previous condition, i.e., if this
condition is satisfied, the previous one is also satisfied.

Figure 10: Initial eccentricity distribution for GCs with
$M_{\mathrm{PBH}}=100\,M_{\odot}$ and different PBH abundance, i.e.,
$f^{\mathrm{cl}}_{\mathrm{PBH}}=1$ (left panel) and
$f^{\mathrm{cl}}_{\mathrm{PBH}}=10^{-1}$ (right panel).

Finally, we show in figure 10 the distribution of initial eccentricity for GCs
with different masses and PBH abundance. Also in this case, the eccentricity
probability distribution function is computed from a catalog of hyperbolic
encounters with parameters sampled from the distribution discussed above. The
emission of GW due to hyperbolic encounters is peaked at values of
eccentricity close to unity; however, as we can see from the figure, such
values are probable only in low-mass GCs. Therefore, it is reasonable to
expect that the overall GWB produced by this class of dynamical interactions
in this environment will not be comparable to other, more effective, emission
channels.

## References

  * [1] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “Observation of Gravitational Waves from a Binary Black Hole Merger”, Phys. Rev. Lett. 116 (Feb, 2016) 061102, arXiv:1602.03837.
  * [2] B. P. Abbott et al., “Multi-messenger Observations of a Binary Neutron Star Merger”, The Astrophysical Journal 848 no. 2, (Oct, 2017) L12, arXiv:1710.05833.
  * [3] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs”, Phys. Rev. X 9 (Sep, 2019) 031040, arXiv:1811.12907.
  * [4] The LIGO Scientific Collaboration and Virgo Collaboration, R. Abbott et al., “GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run”, Phys. Rev. X 11 (Jun, 2021) 021053, arXiv:2010.14527.
  * [5] The LIGO Scientific Collaboration, the Virgo Collaboration and the KAGRA Collaboration, R. Abbott et al., “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run”, Phys. Rev. X 13 (Dec, 2023) 041039, arXiv:2111.03606.
  * [6] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo”, The Astrophysical Journal 882 no. 2, (Sep, 2019) L24, arXiv:1811.12940.
  * [7] The LIGO Scientific Collaboration and Virgo Collaboration, R. Abbott et al., “Population Properties of Compact Objects from the Second LIGO–Virgo Gravitational-Wave Transient Catalog”, The Astrophysical Journal Letters 913 no. 1, (May, 2021) L7, arXiv:2010.14533.
  * [8] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “A gravitational-wave standard siren measurement of the Hubble constant”, Nature 551 (Nov, 2017) 85–88, arXiv:1710.05835.
  * [9] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “A Gravitational-wave Measurement of the Hubble Constant Following the Second Observing Run of Advanced LIGO and Virgo”, The Astrophysical Journal 909 no. 2, (Mar, 2021) 218, arXiv:1908.06060.
  * [10] The LIGO Scientific Collaboration and Virgo Collaboration, R. Abbott et al., “Constraints on the cosmic expansion history from GWTC-3”, The Astrophysical Journal 949 no. 2, (Jun, 2023) 76, arXiv:2111.03604.
  * [11] P. Creminelli and F. Vernizzi, “Dark Energy after GW170817 and GRB170817A”, Phys. Rev. Lett. 119 (Dec, 2017) 251302, arXiv:1710.05877.
  * [12] J. M. Ezquiaga and M. Zumalacárregui, “Dark Energy After GW170817: Dead Ends and the Road Ahead”, Phys. Rev. Lett. 119 (Dec, 2017) 251304, arXiv:1710.05901.
  * [13] T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, “Strong Constraints on Cosmological Gravity from GW170817 and GRB 170817A”, Phys. Rev. Lett. 119 (Dec, 2017) 251301, arXiv:1710.06394.
  * [14] G. Bertone, D. Croon, M. A. Amin, K. K. Boddy, B. J. Kavanagh, K. J. Mack, P. Natarajan, T. Opferkuch, K. Schutz, V. Takhistov, C. Weniger, and T.-T. Yu, “Gravitational wave probes of dark matter: challenges and opportunities”, SciPost Phys. Core 3 (2020) 007, arXiv:1907.10610.
  * [15] G. Bertone, D. Hooper, and J. Silk, “Particle dark matter: evidence, candidates and constraints”, Physics Reports 405 no. 5, (2005) 279–390, arXiv:hep-ph/0404175.
  * [16] G. Bertone and T. M. P. Tait, “A New Era in the Quest for Dark Matter”, Nature 562 (2018) 51–56, arXiv:1810.01668.
  * [17] G. Jungman, M. Kamionkowski, and K. Griest, “Supersymmetric dark matter”, Physics Reports 267 no. 5-6, (Mar, 1996) 195–373.
  * [18] L. Bergström, “Non-baryonic dark matter: observational evidence and detection methods”, Reports on Progress in Physics 63 no. 5, (Apr, 2000) 793–841, arXiv:hep-ph/0002126.
  * [19] M. S. Turner, “Windows on the Axion”, Phys. Rept. 197 (1990) 67–97.
  * [20] G. Carosi, A. Friedland, M. Giannotti, M. J. Pivovaroff, J. Ruz, and J. K. Vogel, “Probing the axion-photon coupling: phenomenological and experimental perspectives. A snowmass white paper”, arXiv:1309.7035.
  * [21] D. J. Marsh, “Axion cosmology”, Physics Reports 643 (Jul, 2016) 1–79, arXiv:1510.07633.
  * [22] A. Kusenko, “Sterile neutrinos: The dark side of the light fermions”, Physics Reports 481 no. 1-2, (Sep, 2009) 1–28, arXiv:0906.2968.
  * [23] Y. B. Zel’dovich and I. D. Novikov, “The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model”, Soviet Astronomy 10 (Feb, 1967) 602.
  * [24] S. Hawking, “Gravitationally collapsed objects of very low mass”, Monthly Notices of the Royal Astronomical Society 152 (1971) 75.
  * [25] B. J. Carr and S. W. Hawking, “Black Holes in the Early Universe”, Monthly Notices of the Royal Astronomical Society 168 no. 2, (Aug, 1974) 399–415.
  * [26] G. F. Chapline, “Cosmological effects of primordial black holes”, Nature 253 (Jan, 1975) 251.
  * [27] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, “Did LIGO Detect Dark Matter?”, Physical Review Letters 116 no. 20, (May, 2016) 201301, arXiv:1603.00464.
  * [28] S. Clesse and J. García-Bellido, “The clustering of massive Primordial Black Holes as Dark Matter: Measuring their mass distribution with advanced LIGO”, Physics of the Dark Universe 15 (2017) 142–147, arXiv:1603.05234.
  * [29] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, “Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914”, Phys. Rev. Lett. 117 (Aug, 2016) 061101, arXiv:1603.08338.
  * [30] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “Constraints on primordial black holes”, Rept. Prog. Phys. 84 no. 11, (2021) 116902, arXiv:2002.12778.
  * [31] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, “Primordial black holes - perspectives in gravitational wave astronomy”, Class. Quant. Grav. 35 no. 6, (2018) 063001, arXiv:1801.05235.
  * [32] The LISA Cosmology Working Group, E. Bagui et al., “Primordial black holes and their gravitational-wave signatures”, arXiv:2310.19857.
  * [33] C. Caprini and D. G. Figueroa, “Cosmological Backgrounds of Gravitational Waves”, Class. Quant. Grav. 35 no. 16, (2018) 163001, arXiv:1801.04268 [astro-ph.CO].
  * [34] The NANOGrav Collaboration, G. Agazie et al., “The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background”, Astrophys. J. Lett. 951 no. 1, (2023) L8, arXiv:2306.16213.
  * [35] The EPTA Collaboration and InPTA Collaboration, J. Antoniadis et al., “The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals”, Astron. Astrophys. 678 (2023) A50, arXiv:2306.16214.
  * [36] D. J. Reardon et al., “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array”, Astrophys. J. Lett. 951 no. 1, (2023) L6, arXiv:2306.16215.
  * [37] H. Xu et al., “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I”, Res. Astron. Astrophys. 23 no. 7, (2023) 075024, arXiv:2306.16216.
  * [38] The International Pulsar Timing Array Collaboration, G. Agazie et al., “Comparing recent PTA results on the nanohertz stochastic gravitational wave background”, arXiv:2309.00693.
  * [39] P. Amaro-Seoane et al., “Laser Interferometer Space Antenna”, arXiv:1702.00786.
  * [40] S. Phinney et al., “The Big Bang Observer: direct detection of gravitational waves from the birth of the Universe to the Present”, NASA Mission Concept Study (2004) .
  * [41] S. Kawamura et al., “The Japanese space gravitational wave antenna: DECIGO”, Classical and Quantum Gravity 28 no. 9, (Apr, 2011) 094011.
  * [42] S. Kawamura et al., “Current status of space gravitational wave antenna DECIGO and B-DECIGO”, arXiv:2006.13545.
  * [43] M. Punturo et al., “The Einstein Telescope: a third-generation gravitational wave observatory”, Classical and Quantum Gravity 27 no. 19, (Sep, 2010) 194002.
  * [44] M. Maggiore et al., “Science case for the Einstein telescope”, Journal of Cosmology and Astroparticle Physics 2020 no. 03, (Mar, 2020) 050, arXiv:1912.02622.
  * [45] D. Reitze et al., “Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO”, arXiv:1907.04833.
  * [46] M. Evans et al., “A Horizon Study for Cosmic Explorer: Science, Observatories, and Community”, arXiv:2109.09882.
  * [47] A. Raccanelli, E. D. Kovetz, S. Bird, I. Cholis, and J. B. Muñoz, “Determining the progenitors of merging black-hole binaries”, Phys. Rev. D 94 (Jul, 2016) 023516, arXiv:1605.01405.
  * [48] G. Scelfo, N. Bellomo, A. Raccanelli, S. Matarrese, and L. Verde, “GW$\times$LSS: chasing the progenitors of merging binary black holes”, Journal of Cosmology and Astroparticle Physics 2018 no. 09, (Sep, 2018) 039, arXiv:1809.03528.
  * [49] G. Scelfo, M. Spinelli, A. Raccanelli, L. Boco, A. Lapi, and M. Viel, “Gravitational waves × HI intensity mapping: cosmological and astrophysical applications”, Journal of Cosmology and Astroparticle Physics 2022 no. 01, (Jan, 2022) 004, arXiv:2106.09786.
  * [50] S. Mukherjee and J. Silk, “Can we distinguish astrophysical from primordial black holes via the stochastic gravitational wave background?”, Mon. Not. Roy. Astron. Soc. 506 no. 3, (2021) 3977–3985, arXiv:2105.11139 [gr-qc].
  * [51] S. Mukherjee, M. S. P. Meinema, and J. Silk, “Prospects of discovering subsolar primordial black holes using the stochastic gravitational wave background from third-generation detectors”, Mon. Not. Roy. Astron. Soc. 510 no. 4, (2022) 6218–6224, arXiv:2107.02181 [astro-ph.CO].
  * [52] M. Bosi, N. Bellomo, and A. Raccanelli, “Constraining extended cosmologies with GW×LSS cross-correlations”, Journal of Cosmology and Astroparticle Physics 11 (2023) 086, arXiv:2306.03031.
  * [53] S. Libanore, M. Liguori, and A. Raccanelli, “Signatures of primordial black holes in gravitational wave clustering”, Journal of Cosmology and Astroparticle Physics 2023 no. 08, (Aug, 2023) 055, arXiv:2306.03087.
  * [54] M. Raidal, V. Vaskonen, and H. Veermäe, “Gravitational waves from primordial black hole mergers”, Journal of Cosmology and Astroparticle Physics 2017 no. 09, (Sep, 2017) 037, arXiv:1707.01480.
  * [55] Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, “Merger rate of primordial black-hole binaries”, Phys. Rev. D 96 (Dec, 2017) 123523, arXiv:1709.06576.
  * [56] K. Kritos et al., “The Astro-Primordial Black Hole Merger Rates: a Reappraisal”, JCAP 05 (2021) 039, arXiv:2012.03585.
  * [57] X.-X. Kou, G. Fragione, and V. Mandic, “On the Stochastic Gravitational Wave Background from Binary Black Hole Mergers Dynamically Assembled in Dense Star Clusters”, arXiv:2401.04347.
  * [58] M. J. Benacquista and J. M. B. Downing, “Relativistic Binaries in Globular Clusters”, Living Reviews in Relativity 16 no. 1, (Mar, 2013) , arXiv:1110.4423.
  * [59] C. L. Rodriguez, M. Morscher, B. Pattabiraman, S. Chatterjee, C.-J. Haster, and F. A. Rasio, “Binary Black Hole Mergers from Globular Clusters: Implications for Advanced LIGO”, Phys. Rev. Lett. 115 (Jul, 2015) 051101, arXiv:1505.00792.
  * [60] G. Fragione and B. Kocsis, “Black Hole Mergers from an Evolving Population of Globular Clusters”, Phys. Rev. Lett. 121 (Oct, 2018) 161103, arXiv:1806.02351.
  * [61] R. M. O’Leary, B. Kocsis, and A. Loeb, “Gravitational waves from scattering of stellar-mass black holes in galactic nuclei”, Monthly Notices of the Royal Astronomical Society 395 no. 4, (05, 2009) 2127–2146, arXiv:0807.2638.
  * [62] M. C. Miller and V. M. Lauburg, “Mergers of stellar-mass black holes in nuclear star clusters”, The Astrophysical Journal 692 no. 1, (Feb, 2009) 917–923, arXiv:0804.2783.
  * [63] M. Zevin, C. Pankow, C. L. Rodriguez, L. Sampson, E. Chase, V. Kalogera, and F. A. Rasio, “Constraining Formation Models of Binary Black Holes with Gravitational-wave Observations”, The Astrophysical Journal 846 no. 1, (Sep, 2017) 82, arXiv:1704.07379.
  * [64] N. Neumayer, A. Seth, and T. Böker, “Nuclear star clusters”, The Astronomy and Astrophysics Review 28 no. 1, (Jul, 2020) , arXiv:2001.03626.
  * [65] M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez, F. Antonini, D. E. Holz, and C. Pankow, “One Channel to Rule Them All? Constraining the Origins of Binary Black Holes Using Multiple Formation Pathways”, The Astrophysical Journal 910 no. 2, (Apr, 2021) 152, arXiv:2011.10057.
  * [66] S. Banerjee, H. Baumgardt, and P. Kroupa, “Stellar-mass black holes in star clusters: implications for gravitational wave radiation”, Monthly Notices of the Royal Astronomical Society 402 no. 1, (02, 2010) 371–380, arXiv:0910.3954.
  * [67] S. Banerjee, “Stellar-mass black holes in young massive and open stellar clusters – IV. Updated stellar-evolutionary and black hole spin models and comparisons with the LIGO-Virgo O1/O2 merger-event data”, Monthly Notices of the Royal Astronomical Society 500 no. 3, (09, 2020) 3002–3026, arXiv:2004.07382.
  * [68] S. Rastello, M. Mapelli, U. N. Di Carlo, G. Iorio, A. Ballone, N. Giacobbo, F. Santoliquido, and S. Torniamenti, “Dynamics of binary black holes in low-mass young star clusters”, Monthly Notices of the Royal Astronomical Society 507 no. 3, (08, 2021) 3612–3625, arXiv:2105.01669.
  * [69] D. Chattopadhyay, J. Hurley, S. Stevenson, and A. Raidani, “Dynamical double black holes and their host cluster properties”, Monthly Notices of the Royal Astronomical Society 513 no. 3, (Jul, 2022) 4527–4555, arXiv:2202.08924.
  * [70] S. Torniamenti, S. Rastello, M. Mapelli, U. N. Di Carlo, A. Ballone, and M. Pasquato, “Dynamics of binary black holes in young star clusters: the impact of cluster mass and long-term evolution”, Monthly Notices of the Royal Astronomical Society 517 no. 2, (10, 2022) 2953–2965, arXiv:2203.08163.
  * [71] The Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters”, Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [72] N. Bellomo, J. L. Bernal, A. Raccanelli, and L. Verde, “Primordial black holes as dark matter: converting constraints from monochromatic to extended mass distributions”, Journal of Cosmology and Astroparticle Physics 2018 no. 01, (Jan, 2018) 004, arXiv:1709.07467.
  * [73] V. Bosch-Ramon and N. Bellomo, “Mechanical feedback effects on primordial black hole accretion”, Astron. Astrophys. 638 (2020) A132, arXiv:2004.11224.
  * [74] V. Bosch-Ramon, “3D hydrodynamical simulations of the impact of mechanical feedback on accretion in supersonic stellar-mass black holes”, Astron. Astrophys. 660 (Apr, 2022) A5, arXiv:2201.09601.
  * [75] L. Piga, M. Lucca, N. Bellomo, V. Bosch-Ramon, S. Matarrese, A. Raccanelli, and L. Verde, “The effect of outflows on CMB bounds from Primordial Black Hole accretion”, Journal of Cosmology and Astroparticle Physics 12 (2022) 016, arXiv:2210.14934.
  * [76] T. Nakamura, M. Sasaki, T. Tanaka, and K. S. Thorne, “Gravitational waves from coalescing black hole MACHO binaries”, Astrophys. J. Lett. 487 (1997) L139–L142, arXiv:astro-ph/9708060.
  * [77] K. Ioka, T. Chiba, T. Tanaka, and T. Nakamura, “Black hole binary formation in the expanding universe: Three body problem approximation”, Phys. Rev. D 58 (1998) 063003, arXiv:astro-ph/9807018.
  * [78] M. Raidal, C. Spethmann, V. Vaskonen, and H. Veermäe, “Formation and Evolution of Primordial Black Hole Binaries in the Early Universe”, Journal of Cosmology and Astroparticle Physics 02 (2019) 018, arXiv:1812.01930.
  * [79] G. Ballesteros, P. D. Serpico, and M. Taoso, “On the merger rate of primordial black holes: effects of nearest neighbours distribution and clustering”, Journal of Cosmology and Astroparticle Physics 10 (2018) 043, arXiv:1807.02084.
  * [80] V. De Luca, V. Desjacques, G. Franciolini, and A. Riotto, “The clustering evolution of primordial black holes”, Journal of Cosmology and Astroparticle Physics 11 (2020) 028, arXiv:2009.04731.
  * [81] J. Garriga and N. Triantafyllou, “Enhanced cosmological perturbations and the merger rate of PBH binaries”, Journal of Cosmology and Astroparticle Physics 09 (2019) 043, arXiv:1907.01455.
  * [82] K. Hayasaki, K. Takahashi, Y. Sendouda, and S. Nagataki, “Rapid merger of binary primordial black holes: An implication for GW150914”, Publ. Astron. Soc. Jap. 68 no. 4, (2016) 66, arXiv:0909.1738.
  * [83] S. Mashchenko and A. Sills, “Globular Clusters with Dark Matter Halos. I. Initial Relaxation”, The Astrophysical Journal 619 no. 1, (Jan, 2005) 243–257, arXiv:astro-ph/0409605.
  * [84] H. Baumgardt and S. Mieske, “High mass-to-light ratios of ultra-compact dwarf galaxies - evidence for dark matter?”, Monthly Notices of the Royal Astronomical Society 391 no. 2, (Dec, 2008) 942–948, arXiv:0809.2783.
  * [85] R. Kippenhahn, A. Weigert, and A. Weiss, Stellar structure and evolution.  Springer, 2012.
  * [86] K. El-Badry, E. Quataert, D. R. Weisz, N. Choksi, and M. Boylan-Kolchin, “The formation and hierarchical assembly of globular cluster populations”, Monthly Notices of the Royal Astronomical Society 482 no. 4, (11, 2018) 4528–4552, arXiv:1805.03652.
  * [87] E. Ardi and H. Baumgardt, “Depletion of dark matter within globular clusters”, J. Phys. Conf. Ser. 1503 no. 1, (2020) 012023.
  * [88] H. Wirth, K. Bekki, and K. Hayashi, “Formation of massive globular clusters with dark matter and its implication on dark matter annihilation”, Monthly Notices of the Royal Astronomical Society: Letters 496 no. 1, (May, 2020) L70–L74, arXiv:2005.07128.
  * [89] P. Boldrini, R. Mohayaee, and J. Silk, “Embedding globular clusters in dark matter minihaloes solves the cusp–core and timing problems in the Fornax dwarf galaxy”, Monthly Notices of the Royal Astronomical Society 492 no. 3, (Jan, 2020) 3169–3178, arXiv:1909.07404.
  * [90] R. G. Carlberg and L. C. Keating, “Simulating Globular Clusters in Dark Matter Subhalos”, The Astrophysical Journal 924 no. 2, (Jan, 2022) 77, arXiv:2105.13900.
  * [91] X. Ma, M. Y. Grudić , E. Quataert, P. F. Hopkins, C.-A. Faucher-Giguère, M. Boylan-Kolchin, A. Wetzel, J. hoon Kim, N. Murray, and D. Kereš, “Self-consistent proto-globular cluster formation in cosmological simulations of high-redshift galaxies”, Monthly Notices of the Royal Astronomical Society 493 no. 3, (Feb, 2020) 4315–4332, arXiv:1906.11261.
  * [92] F. Capela, M. Pshirkov, and P. Tinyakov, “Constraints on primordial black holes as dark matter candidates from capture by neutron stars”, Phys. Rev. D 87 no. 12, (2013) 123524, arXiv:1301.4984.
  * [93] D. C. Heggie and P. Hut, “Dark matter in globular clusters”, IAU Symp. 174 (1996) 303, arXiv:astro-ph/9511115.
  * [94] E. Hassani and S. M. G. Mousavi, “Dark Matter Effects on Stellar Populations in Globular Clusters”, arXiv:2106.03398.
  * [95] P. Kroupa, “On the variation of the initial mass function”, Monthly Notices of the Royal Astronomical Society 322 no. 2, (Apr, 2001) 231–246, arXiv:astro-ph/0009005.
  * [96] M. Mapelli, M. Dall’Amico, Y. Bouffanais, N. Giacobbo, M. Arca Sedda, M. C. Artale, A. Ballone, U. N. Di Carlo, G. Iorio, F. Santoliquido, and S. Torniamenti, “Hierarchical black hole mergers in young, globular and nuclear star clusters: the effect of metallicity, spin and cluster properties”, Monthly Notices of the Royal Astronomical Society 505 no. 1, (05, 2021) 339–358, arXiv:2103.05016.
  * [97] V. Pavlík, T. Jerábková, P. Kroupa, and H. Baumgardt, “The black hole retention fraction in star clusters”, A&A 617 (2018) A69, arXiv:1806.05192.
  * [98] M. Morscher, S. Umbreit, W. M. Farr, and F. A. Rasio, “Retention of stellar-mass black holes in globular clusters”, The Astrophysical Journal 763 no. 1, (Jan, 2013) L15, arXiv:1211.3372.
  * [99] D. Heggie and P. Hut, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics.  Cambridge University Press, 2003.
  * [100] J. E. Gunn and R. F. Griffin, “Dynamical studies of globular clusters based on photoelectric radial velocities of individual stars. I. M3.”, The Astrophysical Journal 84 (Jun, 1979) 752–773.
  * [101] G. Meylan and D. Heggie, “Internal dynamics of globular clusters”, The Astronomy and Astrophysics Review 8 (Feb, 1997) 1–143, arXiv:astro-ph/9610076.
  * [102] I. R. King, “The structure of star clusters. III. Some simple dynamical models”, The Astrophysical Journal 71 (Feb, 1966) 64.
  * [103] G. S. Da Costa and K. C. Freeman, “The structure and mass function of the globular cluster M3.”, The Astrophysical Journal 206 (May, 1976) 128–137.
  * [104] R. W. Michie, “On the Distribution of High Energy Stars in Spherical Stellar Systems”, Monthly Notices of the Royal Astronomical Society 125 no. 2, (08, 1962) 127–139.
  * [105] P. Miocchi, “Central energy equipartition in multimass models of globular clusters”, Monthly Notices of the Royal Astronomical Society 366 no. 1, (02, 2006) 227–234, arXiv:astro-ph/0511263.
  * [106] E. Khalisi, P. Amaro-Seoane, and R. Spurzem, “A comprehensive nbody study of mass segregation in star clusters: energy equipartition and escape”, Monthly Notices of the Royal Astronomical Society 374 no. 2, (12, 2006) 703–720, arXiv:astro-ph/0602570.
  * [107] M. Trenti and R. van der Marel, “No energy equipartition in globular clusters”, Monthly Notices of the Royal Astronomical Society 435 no. 4, (09, 2013) 3272–3282, arXiv:1302.2152.
  * [108] S. Torniamenti, G. Bertin, and P. Bianchini, “A simple two-component description of energy equipartition and mass segregation for anisotropic globular clusters”, A&A 632 (2019) A67, arXiv:1909.07221.
  * [109] S. Chandrasekhar, “Dynamical Friction. I. General Considerations: the Coefficient of Dynamical Friction.”, The Astrophysical Journal 97 (Mar, 1943) 255.
  * [110] M. Gieles and A. Zocchi, “A family of lowered isothermal models”, Monthly Notices of the Royal Astronomical Society 454 no. 1, (09, 2015) 576–592, arXiv:1508.02120.
  * [111] R. M. O’Leary, F. A. Rasio, J. M. Fregeau, N. Ivanova, and R. O’Shaughnessy, “Binary Mergers and Growth of Black Holes in Dense Star Clusters”, The Astrophysical Journal 637 no. 2, (Feb, 2006) 937, arXiv:astro-ph/0508224.
  * [112] V. Hénault-Brunet, M. Gieles, A. Sollima, L. L. Watkins, A. Zocchi, I. Claydon, E. Pancino, and H. Baumgardt, “Mass modelling globular clusters in the Gaiaera: a method comparison using mock data from an N-body simulation of M4”, Monthly Notices of the Royal Astronomical Society 483 no. 1, (11, 2018) 1400–1425, arXiv:1810.05167.
  * [113] J. Spitzer, Lyman, “Equipartition and the Formation of Compact Nuclei in Spherical Stellar Systems”, The Astrophysical Journal 158 (Dec, 1969) L139.
  * [114] J. Binney and S. Tremaine, Galactic dynamics.  Princeton University Press, 1987.
  * [115] M. Giersz and D. C. Heggie, “Statistics of N-body simulations – I. Equal masses before core collapse”, Monthly Notices of the Royal Astronomical Society 268 no. 1, (05, 1994) 257–275, arXiv:astro-ph/9305008.
  * [116] R. M. O’Leary, R. O’Shaughnessy, and F. A. Rasio, “Dynamical interactions and the black-hole merger rate of the Universe”, Phys. Rev. D 76 (Sep, 2007) 061504, arXiv:1606.07437.
  * [117] D. E. McLaughlin and S. M. Fall, “Shaping the Globular Cluster Mass Function by Stellar-Dynamical Evaporation”, The Astrophysical Journal 679 no. 2, (Jun, 2008) 1272, arXiv:0704.0080.
  * [118] S. S. Larsen, “Observational Constraints on Cluster Evolution”, in Mass Loss from Stars and the Evolution of Stellar Clusters, vol. 388, p. 279.  June, 2008.  arXiv:astro-ph/0609062.
  * [119] C. L. Rodriguez and A. Loeb, “Redshift Evolution of the Black Hole Merger Rate from Globular Clusters”, Astrophys. J. Lett. 866 no. 1, (2018) L5, arXiv:1809.01152.
  * [120] S. F. Portegies Zwart and S. McMillan, “Black hole mergers in the universe”, Astrophys. J. Lett. 528 (2000) L17, arXiv:astro-ph/9910061.
  * [121] F. Antonini and M. Gieles, “Merger rate of black hole binaries from globular clusters: Theoretical error bars and comparison to gravitational wave data from GWTC-2”, Phys. Rev. D 102 (Dec, 2020) 123016, arXiv:2009.01861.
  * [122] M. Boylan-Kolchin, “The globular cluster–dark matter halo connection”, Mon. Not. Roy. Astron. Soc. 472 no. 3, (2017) 3120–3130, arXiv:1705.01548.
  * [123] E. S. Phinney, “A Practical theorem on gravitational wave backgrounds”, arXiv:astro-ph/0108028.
  * [124] T. Regimbau, “The astrophysical gravitational wave stochastic background”, Res. Astron. Astrophys. 11 (2011) 369–390, arXiv:1101.2762.
  * [125] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford University Press, 2007.
  * [126] N. Bellomo, D. Bertacca, A. C. Jenkins, S. Matarrese, A. Raccanelli, T. Regimbau, A. Ricciardone, and M. Sakellariadou, “CLASS_GWB: robust modeling of the astrophysical gravitational wave background anisotropies”, Journal of Cosmology and Astroparticle Physics 2022 no. 06, (Jun, 2022) 030, arXiv:2110.15059.
  * [127] V. Baibhav, D. Gerosa, E. Berti, K. W. K. Wong, T. Helfer, and M. Mould, “The mass gap, the spin gap, and the origin of merging binary black holes”, Phys. Rev. D 102 no. 4, (2020) 043002, arXiv:2004.00650.
  * [128] G. D. Quinlan and S. L. Shapiro, “Dynamical Evolution of Dense Clusters of Compact Stars”, The Astrophysical Journal 343 (Aug, 989) 725.
  * [129] H. Mouri and Y. Taniguchi, “Runaway merging of black holes: analytical constraint on the timescale”, Astrophys. J. Lett. 566 (2002) L17–L20, arXiv:astro-ph/0201102.
  * [130] D. C. Heggie, “Binary Evolution in Stellar Dynamics”, Monthly Notices of the Royal Astronomical Society 173 no. 3, (12, 1975) 729–787.
  * [131] N. Ivanova, K. Belczynski, J. M. Fregeau, and F. A. Rasio, “The Evolution of binary fractions in globular clusters”, Mon. Not. Roy. Astron. Soc. 358 (2005) 572–584, arXiv:astro-ph/0501131.
  * [132] M. Morscher, S. Umbreit, W. M. Farr, and F. A. Rasio, “Retention of Stellar-Mass Black Holes in Globular Clusters”, Astrophys. J. Lett. 763 (2013) L15, arXiv:1211.3372.
  * [133] C. L. Rodriguez et al., “Modeling Dense Star Clusters in the Milky Way and beyond with the Cluster Monte Carlo Code”, The Astrophysical Journal Supplement Series 258 no. 2, (Jan, 2022) 22, arXiv:2106.02643.
  * [134] P. Ajith et al., “A Template bank for gravitational waveforms from coalescing binary black holes. I. Non-spinning binaries”, Phys. Rev. D 77 (2008) 104017, arXiv:0710.2335. [Erratum: Phys.Rev.D 79, 129901 (2009)].
  * [135] P. Ajith et al., “Inspiral-merger-ringdown waveforms for black-hole binaries with non-precessing spins”, Phys. Rev. Lett. 106 (2011) 241101, arXiv:0909.2867.
  * [136] P. Ajith, N. Fotopoulos, S. Privitera, A. Neunzert, and A. J. Weinstein, “Effectual template bank for the detection of gravitational waves from inspiralling compact binaries with generic spins”, Phys. Rev. D 89 no. 8, (2014) 084041, arXiv:1210.6666.
  * [137] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “Upper Limits on the Stochastic Gravitational-Wave Background from Advanced LIGO’s First Observing Run”, Phys. Rev. Lett. 118 (Mar, 2017) 121101, arXiv:1612.02029.
  * [138] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “Search for the isotropic stochastic background using data from Advanced LIGO’s second observing run”, Phys. Rev. D 100 (Sep, 2019) 061101, arXiv:1903.02886.
  * [139] The LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration, R. Abbott et al., “Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run”, Phys. Rev. D 104 (Jul, 2021) 022004, arXiv:2101.12130.
  * [140] T. Regimbau, M. Evans, N. Christensen, E. Katsavounidis, B. Sathyaprakash, and S. Vitale, “Digging Deeper: Observing Primordial Gravitational Waves below the Binary-Black-Hole-Produced Stochastic Background”, Phys. Rev. Lett. 118 (Apr, 2017) 151105, arXiv:1611.08943.
  * [141] The LIGO Scientific Collaboration and Virgo Collaboration, B. P. Abbott et al., “GW170817: Implications for the Stochastic Gravitational-Wave Background from Compact Binary Coalescences”, Phys. Rev. Lett. 120 (Feb, 2018) 091101, arXiv:1710.05837.
  * [142] E. Thrane and J. D. Romano, “Sensitivity curves for searches for gravitational-wave backgrounds”, Phys. Rev. D 88 (Dec, 2013) 124032, arXiv:1310.5300.
  * [143] The LISA Cosmology Working Group, N. Bartolo et al., “Probing anisotropies of the Stochastic Gravitational Wave Background with LISA”, Journal of Cosmology and Astroparticle Physics 2022 no. 11, (Nov, 2022) 009, arXiv:2201.08782.
  * [144] The LIGO Scientific Collaboration, the Virgo Collaboration and the KAGRA Collaboration, R. Abbott et al., “The population of merging compact binaries inferred using gravitational waves through GWTC-3”, arXiv:2111.03634.
  * [145] K. Jedamzik, “Primordial Black Hole Dark Matter and the LIGO/Virgo observations”, JCAP 09 (2020) 022, arXiv:2006.11172 [astro-ph.CO].
  * [146] K. Jedamzik, “Consistency of Primordial Black Hole Dark Matter with LIGO/Virgo Merger Rates”, Phys. Rev. Lett. 126 no. 5, (2021) 051302, arXiv:2007.03565 [astro-ph.CO].
  * [147] G. Franciolini, V. Baibhav, V. De Luca, K. K. Y. Ng, K. W. K. Wong, E. Berti, P. Pani, A. Riotto, and S. Vitale, “Searching for a subpopulation of primordial black holes in LIGO-Virgo gravitational-wave data”, Phys. Rev. D 105 no. 8, (2022) 083526, arXiv:2105.03349 [gr-qc].
  * [148] S. J. Kovacs and K. S. Thorne, “The generation of gravitational waves. III. Derivation of bremsstrahlung formulae.”, The Astrophysical Journal 217 (Oct, 1977) 252–280.
  * [149] B. Kocsis, M. E. Gaspar, and S. Marka, “Detection Rate Estimates of Gravity Waves Emitted during Parabolic Encounters of Stellar Black Holes in Globular Clusters”, The Astrophysical Journal 648 no. 1, (Sep, 2006) 411–429, arXiv:astro-ph/0603441.
  * [150] S. Capozziello, M. De Laurentis, F. De Paolis, G. Ingrosso, and A. Nucita, “Gravitational waves from hyperbolic encounters”, Mod. Phys. Lett. A 23 (2008) 99–107, arXiv:0801.0122.
  * [151] S. Mukherjee, S. Mitra, and S. Chatterjee, “Gravitational wave observatories may be able to detect hyperbolic encounters of black holes”, Mon. Not. Roy. Astron. Soc. 508 no. 4, (2021) 5064–5073, arXiv:2010.00916.
  * [152] S. Dandapat et al., “Gravitational waves from black-hole encounters: Prospects for ground and galaxy-based observatories”, Phys. Rev. D 108 no. 2, (2023) 024013, arXiv:2305.19318.
  * [153] J. García-Bellido and S. Nesseris, “Gravitational wave bursts from Primordial Black Hole hyperbolic encounters”, Physics of the Dark Universe 18 (2017) 123–126, arXiv:1706.02111.
  * [154] J. García-Bellido and S. Nesseris, “Gravitational wave energy emission and detection rates of Primordial Black Hole hyperbolic encounters”, Physics of the Dark Universe 21 (2018) 61–69, arXiv:1711.09702.
  * [155] G. Morrás, J. García-Bellido, and S. Nesseris, “Search for black hole hyperbolic encounters with gravitational wave detectors”, Physics of the Dark Universe 35 (2022) 100932, arXiv:2110.08000.
  * [156] J. P. Madrid, N. W. C. Leigh, J. R. Hurley, and M. Giersz, “Mass evaporation rate of globular clusters in a strong tidal field”, Monthly Notices of the Royal Astronomical Society 470 no. 2, (06, 2017) 1729–1737.
  * [157] K. Takahashi, “The Dynamical Evolution of Globular Clusters and the Stellar Initial Mass Function”, in Dynamics of Star Clusters and the Milky Way, S. Deiters, B. Fuchs, A. Just, R. Spurzem, and R. Wielen, eds., vol. 228 of Astronomical Society of the Pacific Conference Series, p. 69.  Jan., 2001.
  * [158] J. Shin, S. S. Kim, and K. Takahashi, “Dynamical evolution of the mass function and radial profile of the Galactic globular cluster system”, Mon. Not. Roy. Astron. Soc. 386 (2008) 67, arXiv:0805.0058 [astro-ph].
  * [159] K. M. Belotsky et al., “Clusters of primordial black holes”, Eur. Phys. J. C 79 no. 3, (2019) 246, arXiv:1807.06590.
  * [160] I. Cholis, E. D. Kovetz, Y. Ali-Haïmoud, S. Bird, M. Kamionkowski, J. B. Muñoz, and A. Raccanelli, “Orbital eccentricities in primordial black hole binaries”, Phys. Rev. D 94 no. 8, (2016) 084013, arXiv:1606.07437.
  * [161] P. C. Peters, “Gravitational Radiation and the Motion of Two Point Masses”, Phys. Rev. 136 (Nov, 1964) B1224–B1232.
  * [162] R. J. Allison, S. P. Goodwin, R. J. Parker, R. de Grijs, S. F. Portegies Zwart, and M. B. N. Kouwenhoven, “Dynamical mass segregation on a very short timescale”, Astrophys. J. Lett. 700 (2009) L99–L103, arXiv:0906.4806.
  * [163] S. Portegies Zwart, S. McMillan, and M. Gieles, “Young massive star clusters”, Ann. Rev. Astron. Astrophys. 48 (2010) 431–493, arXiv:1002.1961.
  * [164] R. S. Klessen, “The Formation of stellar clusters: Mass spectra from turbulent molecular cloud fragmentation”, Astrophys. J. 556 (2001) 837–846, arXiv:astro-ph/0104127.
  * [165] I. A. Bonnell and M. R. Bate, “Star formation through gravitational collapse and competitive accretion”, Mon. Not. Roy. Astron. Soc. 370 (2006) 488–494, arXiv:astro-ph/0604615.
  * [166] L. De Vittori, P. Jetzer, and A. Klein, “Gravitational wave energy spectrum of hyperbolic encounters”, Phys. Rev. D 86 (2012) 044017, arXiv:1207.5359.
  * [167] M. Gröbner, P. Jetzer, M. Haney, S. Tiwari, and W. Ishibashi, “A note on the gravitational wave energy spectrum of parabolic and hyperbolic encounters”, Classical and Quantum Gravity 37 no. 6, (Feb, 2020) 067002, arXiv:2001.05187.
  * [168] S. Chandrasekhar, “Stochastic problems in physics and astronomy”, Rev. Mod. Phys. 15 (1943) 1–89.