# Primordial Black Hole clusters, phenomenology & implications

José Francisco Nuño Siles<EMAIL_ADDRESS>Juan García-Bellido
<EMAIL_ADDRESS>Instituto de Física Teórica UAM-CSIC, Universidad
Autónoma de Madrid, Cantoblanco 28049 Madrid, Spain

###### Abstract

We present direct N-body simulations of black-hole-only clusters with up to
$2\cdot 10^{4}$ compact objects, zero natal spin and no primordial binaries as
predicted by various primordial black hole (PBH) Dark Matter models. The
clusters’ evolution is computed using NBODY6++GPU, including the effects of
the tidal field of the galaxy, the kicks of black hole mergers and orbit-
averaged energy loss by gravitational radiation of binaries. We investigate
clusters with four initial mass distributions, three of which attempt to model
a generic PBH scenario using a lognormal mass distribution and a fourth one
that can be directly linked to a monochromatic PBH scenario when accretion is
considered. More specifically, we dive into the clusters’ internal dynamics,
describing their expansion and evaporation, along with the resultant binary
black hole mergers. We also compare several simulations with and without black
hole merger kicks and find modelling implications for the probability of
hierarchical mergers.

††preprint: IFT-UAM/CSIC-24-61

##  I Introduction

Black holes (BHs) have long captured the attention of astrophysicists due to
their elusive and mysterious nature. Recent progress in observational
astronomy has revolutionized our ability to study them across various scales,
from detailed images of supermassive BHs at the centres of galaxies
Collaboration  _et al._ (2019) to the detection of microlensing events caused
by BHs passing in front of background sources Wyrzykowski, Ł.  _et al._
(2023). More relevant to this work, the detection of gravitational waves
emitted by merging binary black holes (BBHs) in Earth-based interferometers
Abbott  _et al._ (2023) has improved our understanding of their formation
mechanisms and population statistics while, at the same time, giving rise to
new fundamental questions Abbott  _et al._ (2020a, b, c); Clesse and García-
Bellido (2022), like the nature of Dark Matter Bird  _et al._ (2016); Clesse
and García-Bellido (2017); García-Bellido (2017), and the thermal evolution of
the early Universe Carr  _et al._ (2021).

In this context, the investigation of black hole-only clusters takes on
particular significance as a possible alternative origin capable of shedding
light onto some observations. Theoretical considerations, supported by
numerical simulations, have demonstrated the plausibility of forming such
clusters. For instance, BHs created during the radiation-dominated era are
naturally clustered Chisholm (2006, 2011) if they derive from large non-
Gaussian tails Ezquiaga  _et al._ (2020, 2023); Animali and Vennin (2024).
Other phenomena such as the appearance of closed domain walls Belotsky  _et
al._ (2019) and their collapse Rubin  _et al._ (2001) can lead to their
clustering too. We should also mention the possibility of the existence of
black-hole-only clusters originating via stellar evolution Gieles  _et al._
(2021) for which our analysis and conclusions remain valid with the
appropriate time delay.

Despite the progress in the field, our comprehension of the dynamics of black
hole-only clusters remains limited Trashorras  _et al._ (2021). The
gravitational interactions between thousands of individual BHs together with
the influence of the surrounding environment, present hard challenges for any
attempt at theoretical modelling given the vast range of scales needed to be
considered. However, these complexities also offer a unique opportunity for
N-body simulations to excel.

N-body simulations have proven instrumental in understanding a diverse range
of complex astrophysical phenomena. Using numerical methods to integrate the
gravitational equations of motion for all of the individual particles, we
obtain a detailed description of the dynamics governing the evolution of such
complex systems Inman and Ali-Haïmoud (2019). A complete introduction to the
theoretical foundations of the field of N-body simulations can be found in
Ref. Aarseth (2003) and a shorter, but less recent one, in Spurzem (1999).

In this paper, we present direct N-body simulations of black hole-only dense
clusters, focusing on systems with up to $2\cdot 10^{4}$ compact objects. Our
simulations incorporate physical effects such as the tidal field of a given
host galaxy, black hole relativistic merger kicks, and orbit-averaged energy
loss due to gravitational radiation from binary systems. By utilizing the
NBODY6++GPU code Wang, Long and Spurzem, Rainer and Aarseth, Sverre and
Nitadori, Keigo and Berczik, Peter and Kouwenhoven, M. B. N. and Naab,
Thorsten (2015), we accurately model the evolution of all the BHs for a Hubble
time. The text is structured as follows, in the first section II we describe
the initial conditions in detail, in the second section we analyze the cluster
dynamics, trying to understand why the clusters dissolve and expand. The
following section is focused on binary black hole mergers and the parameters
characterizing these events. In the last section, we summarize our findings
and conclude the work.

The main body of this work consists of the description of the physical
variables of interest regarding such hypothetical celestial bodies and their
behaviour in a Hubble time. We find promising hints of results that could be
linked with current observations but would need more in-depth research before
satisfying conclusions could be drawn. While there are possibly many more
findings in the data, we don’t embark on such extensive endeavour and treat
this work as a comprehensive proof of concept for future works to be carried
out.

##  II Methodology and initial mass function

To evolve our BHs-only clusters, we made use of the direct numerical
integrator NBODY6++GPU. The models we evolve consist of N black holes, where N
ranges from $\mathcal{O}(10^{3})$ to $\mathcal{O}(2\cdot 10^{4})$, drawn
randomly from 4 different initial mass functions plotted in Fig 1. The exact
numbers can be found in Table 1. The first one, a log-normal distribution with
$\mu=10$, $s=1.5$ and $\sigma=0.954$ is depicted in the right panel and
models, for instance, a monochromatic mass function when accretion is
considered and thus the masses of the BHs can increase. The other three
represent an approximation to the PBH Thermal model Carr  _et al._ (2021);
Bagui  _et al._ (2023) composed of a wide mass function with three widths for
the log-normal $\sigma=\\{0.5,1,1.5\\}$. See Carr  _et al._ (2024) for a
recent review of the model. We will refer to them as
{$M\&A$,$\sigma_{0.5}$,$\sigma_{1}$,$\sigma_{1.5}$} respectively. We set the
initial fraction of binaries as well as the individual spins of the BHs to
zero according to theoretical expectations. The BHs are then spatially
distributed in such a way that they follow a Plummer distribution H.C. (1911)
with a Plummer Radius $r_{p}=10~{}\mathrm{pc}$. Based on the assumption that
they could reside in the halo of a galaxy similar to our own, these clusters
are themselves immersed in a central gravitational potential with orbital
radius $R_{c}=34~{}\mathrm{kpc}$ and central mass $M=4.37\times
10^{10}~{}\mathrm{M}_{\odot}$ throughout the entire evolution. This is just a
point mass approximation which leads to a circular movement of period
$T=2.81\mathrm{Gyr}$. In Table 1 we explicitly write key statistics for the
individual clusters.

ID | $M_{\mathrm{total}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $R_{\mathrm{HM}}[pc]$ | $R_{\mathrm{HM}}[pc]$ | ID | $M_{\mathrm{total}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $R_{\mathrm{HM}}[pc]$ | $R_{\mathrm{HM}}[pc]$
---|---|---|---|---|---|---|---|---|---|---|---
M&A | $t=0$ | $t=0$ | $t=T_{U}$ | $t=0$ | $t=T_{U}$ | $\sigma_{0.5}$ | $t=0$ | $t=0$ | $t=T_{U}$ | $t=0$ | $t=T_{U}$
1295 | 16140 | 46.54 | 46.54 | 2.20 | 27.04 | 1520 | 1389 | 5.03 | 5.03 | 9.72 | 12.25
 2570 | 32166 | 66.10 | 97.80 | 2.38 | 27.17 | 3345 | 3007 | 4.50 | 5.08 | 3.23 | 15.55
4046 | 50088 | 63.64 | 91.98 | 4.27 | 26.74 | 5480 | 5009 | 6.67 | 6.67 | 0.77 | 16.21
 5199 | 64280 | 69.98 | 74.80 | 2.20 | 28.39 | 7678 | 7008 | 5.69 | 7.52 | 1.92 | 15.15
8077 | 100116 | 78.27 | 78.27 | 12.53 | 29.50 | 9937 | 9001 | 5.48 | 5.48 | 1.25 | 15.85
 8922 | 110196 | 76.15 | 81.52 | 13.31 | 32.01 | 12201 | 11001 | 5.00 | 7.77 | 2.90 | 15.29
10372 | 128337 | 104.26 | 104.26 | 6.20 | 27.54 | 14366 | 13005 | 5.42 | 8.22 | 2.65 | 16.04
 10535 | 130198 | 68.50 | 74.35 | 9.72 | 28.20 | 16428 | 14912 | 7.72 | 7.72 | 0.48 | 16.04
16159 | 200392 | 57.75 | 113.59 | 4.44 | 27.24 | 18776 | 17007 | 5.40 | 6.51 | 2.29 | 15.00
 20738 | 256346 | 134.10 | 145.70 | 1.37 | 26.65 | 20866 | 18918 | 5.31 | 7.36 | 3.91 | 15.35
ID | $M_{\mathrm{total}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $R_{\mathrm{HM}}[pc]$ | $R_{\mathrm{HM}}[pc]$ | ID | $M_{\mathrm{total}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $M_{\mathrm{max}}[M_{\odot}]$ | $R_{\mathrm{HM}}[pc]$ | $R_{\mathrm{HM}}[pc]$
$\sigma_{1}$ | $t=0$ | $t=0$ | $t=T_{U}$ | $t=0$ | $t=T_{U}$ | $\sigma_{1.5}$ | $t=0$ | $t=0$ | $t=T_{U}$ | $t=0$ | $t=T_{U}$
1505 | 2003 | 32.07 | 32.07 | 1.93 | 27.07 | 1220 | 3020 | 63.67 | 63.67 | 3.37 | 58.92
 3090 | 4004 | 18.56 | 18.56 | 7.52 | 22.21 | 3423 | 8010 | 105.74 | 105.74 | 7.15 | 41.00
5288 | 7007 | 34.14 | 34.14 | 7.60 | 21.97 | 5258 | 13014 | 541.18 | 541.18 | 1.95 | 40.15
 7507 | 10001 | 25.82 | 34.29 | 3.95 | 24.77 | 8025 | 20021 | 157.23 | 157.23 | 0.41 | 41.62
10663 | 14000 | 46.77 | 46.77 | 6.64 | 26.96 | 10011 | 24187 | 179.83 | 280.36 | 1.45 | 60.33
 12834 | 17004 | 32.81 | 58.06 | 4.41 | 23.91 | 12409 | 30007 | 114.51 | 114.51 | 0.07 | 55.98
15136 | 20000 | 58.31 | 58.31 | 0.47 | 25.47 | 14691 | 38027 | 731.69 | 1004.13 | 1.40 | 37.81
 17554 | 23008 | 45.68 | 64.90 | 0.82 | 25.62 | 17182 | 43013 | 554.27 | 554.27 | 5.59 | 46.36
20590 | 27008 | 51.60 | 76.41 | 4.70 | 26.64 | 20261 | 49021 | 657.47 | 784.02 | 3.52 | 41.41
Table 1: In this table we show the initial conditions for the clusters studied
in the paper. The ID also corresponds to the initial number of BHs. We also
provide the total mass of the cluster, the maximum mass of any single BH in
the cluster at t=0 as well as after we have evolved it for the age of the
Universe and the radius that encircles half of the total mass of the cluster
at both times. The reasons for the maximum mass differing at the beginning and
end of the evolution are due to mergers retained in the clusters or the escape
of the most massive body. Figure 1: Initial mass function of the clusters.
Those based on a typical PBH mass function ($\sigma_{i}$) are in the left
panel, while the monochromatic with accretion IMF ($M\&A$) is in the right
one. We also write the general formula for all the IMFs and the values for the
specific parameters to reproduce our results.

For these models to be of any real physical interest Primordial Black Holes
(PBHs) should have been created in large numbers at the beginning of the
Universe and close enough to each other so that they would have formed
clusters so dense as to overcome the Hubble expansion. Several mechanisms
could lead to such scenarios Rubin  _et al._ (2001); Khlopov  _et al._ (2002).
The idea of them making up a significant fraction of the Dark Matter would add
relevance to the study but is not essential. To understand the implications of
such a scenario, we can estimate the number of clusters needed to get the
observed DM mass. Given the clusters’ total masses range is
$\mathcal{O}(10^{3}-10^{5})\mathrm{M}_{\odot}$, we would need around
$\mathcal{O}(10^{6}-10^{9})$ clusters to conform the entirety of the dark
matter in our galaxy, where we have used a total DM mass of $1.2\times
10^{12}\mathrm{M}_{\odot}$ Fragione and Loeb (2017). Given the nature of the
clusters’ components, direct detection is intrinsically difficult; therefore,
our approach will focus on identifying potential footprints they may leave
behind.

##  III Dynamics

The dynamics governing self-gravitating systems are extremely non-linear,
thus, trying to infer general qualitative principles is the most we aim to do
in this paper. In this section, we will try to shed some light on the
understanding of the passage of time for these systems. We will do so by
following various global physical variables that would serve as a proxy for
the whole cluster.

From the point of view of the spatial extension of the clusters with time, we
can study global quantities such as the Lagrangian radii and the core radii as
defined in Aarseth (2003). We plot in the left panel of Fig 3 the Lagrangian
radii for 50% of the mass of the different clusters. This quantity is also
called the half-mass radius. The tendency is for the half-mass radius to grow
with time, with a steeper slope at the beginning of the simulation and
flattening towards the end, albeit always increasing. In the case of the core
radius, which is displayed for all the clusters in the right panel of Fig 3,
we observe that it remains constant, with a slight upward trend throughout the
evolution. This translates into a faster expansion of the outer layers,
converging towards an almost non-expanding core. Additionally, it is
interesting to note that most clusters exhibit a similar core radius
regardless of the initial number of BHs.

Now we analyse the stability of these clusters. When carrying out the
numerical evolution, it is evident from the beginning that BHs are constantly
escaping the combined gravitational influence of the cluster and the host
galaxy so that the total number being evolved by the code gets reduced with
time. In NBODY6++GPU, BHs are removed from the cluster evolution once they
have reached twice the tidal radius von Hoerner (1957) as calculated using the
mass of the cluster at the specific time. This will be our definition of an
escaper (single or binary).

These BHs acquire the needed escape velocity after one or several close
encounters with other BHs. This process happens so often that we can consider
all clusters to be metastable, that is, most of them don’t release enough BHs
to dissolve in the age of the Universe, but they will eventually do so. In the
left panel of Fig 2, we show the number of BHs left in the clusters as a
function of time, normalized to the initial number of BHs of each cluster so
that the various runs are comparable. We see that some of the lightest
clusters of the $\sigma_{i}$ type dissolve completely in the age of the
Universe and, in general, the larger the $\sigma$ for the IMF of the cluster,
the shorter the expected lifespan. The most stable set of initial conditions
seems to be that of $M\&A$. To explain this, we have to understand the main
mechanism behind the clusters’ dissolution and the difference in their
behaviour as a function of the IMF. We can imagine a simplistic scenario in
which we just consider 2-body interactions as the main driver of the
evaporation of the clusters. This is probably a good approximation as 3+ body
encounters are suppressed based on their scattering cross-section. Now, when
considering the different IMFs, we see that the main difference across the
various $\sigma$ is the mass ratio distribution for random pairs. The larger
the $\sigma$ is, the more extreme mass ratios exist, and the expected mass
ratio goes further away from 1. In the $M\&A$ case, most of the pairs
concentrate close to equal mass ratios. This led us to the conclusion that the
more extreme the mass ratios are, the more slingshots of the lighter BHs occur
which is translated into more BH escapers. This is also in agreement with the
fact that the lighter BHs escape earlier, leaving the heavier ones in the
cluster. This can directly be appreciated in the right panel of Fig 2 where we
plot the average mass of the BHs inside the half-mass radius. This quantity is
strictly increasing in the $\sigma_{i}$ cases with a slope that steepens in
the last stages of the cluster life due to the wider range of masses present
in those clusters. Given the almost monochromatic nature of the $M\&A$ IMFs,
the average mass tends to stay constant with a variance that decreases
inversely with the initial BH number.

Figure 2: In the left panel of this figure we plot the number of BHs for the
various clusters as a function of time normalized to the initial number of
BHs. The dark blue lines represent the lightest cluster for a given IMF
increasing upwards with no overlapping between the lines in the last $\sim 7$
Gyrs of evolution. In the right panel, we plot the average mass of BHs for the
various clusters inside the half-mass radius as a function of time.

Figure 3: In this figure we plot the Lagrangian radii for 50% of the mass and
the core radii for the various clusters as defined in Aarseth (2003). Most
lines overlap cluster-wise, so we don’t add a legend to distinguish them only
pointing out that the lighter clusters are also those showing the largest
variance. For the core radii, we include a zooming window to what looks like
possible oscillations for the heaviest of the $M\&A$ clusters, although it
could also be just shot noise in the calculation.

Figure 4: In this figure we plot the initial semimajor axis and eccentricity
($a_{0}$,$e_{0}$) of all the BHs that escape the gravitational pull of the
clusters within the age of the Universe. These are the binaries logged by the
Nbody code as escapers. In order from left to right and top to bottom, we are
showing the clusters {$M\&A$,$\sigma_{0.5}$,$\sigma_{1}$,$\sigma_{1.5}$}. The
dashed lines represent the values of ($a_{0}$,$e_{0}$) for which
$\tau_{\mathrm{merger}}(m_{1},m_{2},a_{0},e_{0})=\mathrm{T}_{Hubble}$ with
($m_{1}$,$m_{2}$)[$M_{\odot}$] being the pairs of number in the base of the
plot beside each line. These calculations, however, does not take into account
the time delay from the beginning of the simulation until the binary escapes
the cluster which can be quite substantial. The width of the dots represents
the total mass of the binary defined as $M_{\mathrm{Total}}=m_{1}+m_{2}$.

If a sufficient number of these clusters existed in our Universe, single black
hole escapers could potentially be detected from Earth through microlensing or
via stellar disruption events. Those that escape as binaries could also be
detected in the final stages of the inspiral, as well as the merger and the
ringdown in Earth-based gravitational waves interferometers. We will call
these binary escapers, off-cluster mergers if the initial conditions regarding
their orbits’ parameters once they have escaped from the gravitational pull of
the cluster are such that they coalesce in the age of the Universe. To
calculate the coalescence time, we assume orbit shrinking via gravitational
wave emission and use the formula as derived in Peters (1964) in post-
processing.

Starting with the single BH escapers, we first acknowledge the fact that the
distributions of masses are very similar to the IMF of the clusters with a
slight skew towards the lighter BHs due to usual gravitational mass
segregation. This means that enough observations of these rogue compact
objects would paint a good picture of the progenitor clusters’ distributions.

Regarding binary escapers, of which off-cluster mergers represent a subset, we
find a correlation in the distribution of semi-major axes with the initial BH
number. The larger the initial number of BHs in the cluster, the tighter the
binaries that escape (smaller semi-major axis) and vice-versa. This fact can
be visually spotted in figure 4. This is expected as the binaries that may
survive in a less dense cluster based on the rate of interactions, may not do
so in denser environments. In other words, the binding energy needed for a
binary to prevail and escape the cluster grows with the density. Regarding
their eccentricities ($e_{0}$), we see an excess close to maximum
eccentricity. This is just the result of the many interactions the binaries
need to endure before escaping the cluster.

Results also show that the absolute number of binary escapers increases with
the initial BH number within the same type of IMF. Going even further, it also
seems to correlate with the stability of the clusters surveyed. The more
stable the cluster type, the more binaries escape the gravitational influence
with the order from more to less stable being
{$M\&A$,$\sigma_{0.5}$,$\sigma_{1}$,$\sigma_{1.5}$}.

##  IV BBH mergers

During the numerical evolution of the clusters, the centres of two (or more)
BHs might get close enough due to the gravitational interactions that their
event horizons would merge and the code would start evolving them as a single
entity, thus declaring a coalescence. These events occur inside the clusters
Sigurdsson and Hernquist (1993); Zwart and McMillan (2002), mainly induced by
interactions where 3 or more compact objects are involved Heggie (1975), that
is, binary-single interactions or binary-binary interactions. We refer to them
as in-cluster mergers. The second type of merger, as described in Section III,
occurs when two BHs leave the cluster as a binary bound by their mutual
gravitational interaction. Due to gravitational radiation and under the
assumption that nothing else ever interacts with these pairs, they will
inevitably end up merging. If the time of merger is less than a Hubble time,
we refer to them as off-cluster mergers.

###  IV.1 Merger count

Understanding the difference between the various initial conditions can be
very useful in assessing the underlying mechanisms of the mergers. In the
right panel of Fig 5 we show four histograms with the number of mergers per
cluster. Discussing first the in-cluster mergers, we find that among the
$\sigma_{i}$ clusters there is a trend in the number of mergers inverse to the
width ($\sigma$) of the distribution. In the $\sigma_{0.5}$ and $\sigma_{1}$
cases, a maximum in the merger count seems to arise for
$N\sim\mathcal{O}(1.7\cdot 10^{4})$, although it could be just a statistical
anomaly. The total number count for comparable initial $N_{BH}$ lies below
that of the $M\&A$ cluster type. This can be explained based on the lower
stability of the clusters themselves due to the larger mass range present. The
large number of mergers for the $M\&A$ cluster type has as a result also the
existence of second-generation BHs, which we discuss in sec IV.2. For this
last case, there is a clear upward-sloping tendency in the number count both
for in and off-cluster mergers.

We observe a lower number count in the off-cluster merger statistics. This is
expected as the delay between the start of the simulations and the ejection of
the binary from the clusters can be significant, thus making the expected
merger time for the binary escapers way longer than Hubble time. This comes
together with the fact that gravitational wave emission is a slow mechanism
for shrinking the orbit of the binary in comparison with the interactions
governing the in-cluster mergers. Comparing cluster types, we find again fewer
counts for the $\sigma_{i}$ with no off-cluster mergers happening in any
$\sigma_{1,1.5}$ clusters. This means that no binary escaper ever merged in
the age of the Universe. One of the reasons for this is that the distribution
in the semi-major axis of these two cluster types lies at higher values
always, which is generally equivalent to larger merger times. Another reason
is that the average delay time before escaping is larger for $\sigma_{1,1.5}$
than for the other two clusters. It can also be explained by the fact that
there are fewer binary escapers for those two clusters, as can be seen in
figure 4.

Figure 5:  In the left panel we show a histogram of the number of in-cluster
mergers for the $M\&A$ cluster type. The bars with the blue border represent
the simulations used throughout the paper, where BH merger kicks are taken
into account. The bars with the black border are those same simulations when
the kicks are not considered. In the right panel, we plot a histogram of the
number of mergers per cluster type. We differentiate, in blue the in-cluster
and in red the off-cluster mergers.

###  IV.2 Hierarchical mergers

The product of any merger is another BH with the mass equal to the sum of the
masses and spins calculated using Jiménez-Forteza  _et al._ (2017). These
remnants may or may not stay in the cluster as the velocity imprinted on them
Lousto  _et al._ (2012, 2010); Campanelli  _et al._ (2007a) can be very high.
Those which not leave the cluster are subject to the possibility of being
involved in another merger. To understand how common this phenomenon is, we
show a histogram in the left panel of Fig 5 with the total number of BH
mergers identified in the lifetime of the clusters for the $M\&A$ cluster
type. There are two different bars per Cluster ID. The ones whose border is
black represent simulations in which BHs’ natal kicks are not taken into
account while bars with a blue border are the normal simulations that do
include such a physical phenomenon and thus represent a more complete
description of reality. The initial conditions are the same for both cases. We
find that when kicks are not present, the probability that a remnant stays in
the cluster is high as the velocity after the merger is obtained based on
linear momentum conservation. However, when kicks are included, relativistic
effects imprint a velocity in the remnant that is usually much larger than the
escape velocity of the cluster. This has as a consequence a sharp drop in the
probability of hierarchical mergers. We can see this reflected in the fact
that hierarchical mergers of up to $6^{\mathrm{th}}$ generation can happen in
most of the clusters if kicks are ignored but it’s only for N larger than
10535 BHs that we find hierarchical mergers when kicks are considered. What is
more, anything beyond $2^{\mathrm{nd}}$ generation is not on the cards, as the
kick velocity is even larger in the case where progenitor spin is non-zero
Campanelli  _et al._ (2007b). This has broader consequences in the
distribution of other binary parameters, such as the final spin. As can be
seen in Fig. 8 a limited number of BH remnants have $a_{f}>0.7$ and they are
all produced via second-generation mergers. A variation in the statistics of
such products would lead to critical differences in the detected population of
BHs via for example gravitational wave interferometers.

This result is of great importance as it highlights the significant impact
that an occasionally overlooked physical effect can have on observational
probes. Any simulation lacking the modelling of BH merger kicks is doomed to
overestimate the hierarchical merger rate and the skew to larger remnant
spins.

###  IV.3 Merger Rate

Cluster ID | Redshift | Rate z=0 | Cluster ID | Redshift | Rate z=0
---|---|---|---|---|---
$M\&A$ | last merger | (events/yr/Gpc3) | $\sigma_{0.5}$ | last merger | (events/yr/Gpc3)
2570 | 5.89 | 0.00 | 3345 | 1.71 | 0.00
 4046 | 2.79 | 0.00 | 5480 | 0.03 | 263.69
5199 | 0.18 | 1.36 | 7678 | 2.67 | 0.00
 8077 | 0.44 | 2.72 | 9937 | 0.83 | 80.93
8922 | 1.77 | 0.00 | 12201 | 1.65 | 0.00
 10372 | 0.43 | 1.36 | 14366 | 0.33 | 77.48
10535 | 0.04 | 5.45 | 16428 | 0.60 | 29.55
 16159 | 0.18 | 6.81 | 18776 | 0.34 | 45.33
20738 | 0.05 | 14.98 | 20866 | 0.44 | 18.34
Table 2:  We present the local merger rate and redshift of the last identified
in-cluster BH merger for every cluster with at least one merger within the age
of the Universe. We focus on the {$M\&A$,$\sigma_{0.5}$} cluster types due to
their larger merger counts. The rate is calculated using all the in-cluster
mergers with $z<1$ to approximate the local merger rate, assuming that all of
the Dark Matter is contained within such dense clusters.

Following the merger history of the various clusters throughout their
evolution could give us clues about possible probes of these models, as
currently, our best chance to infer the existence of these objects is via
gravitational waves. For that purpose, we plot in Fig 6 the in-cluster and
off-cluster merger rates as a function of time for the cluster types {$M\&A$,
$\sigma_{0.5}$} based on their larger merger statistics compared to the other
two. The first thing we observe is there is a pronounced peak for the in-
cluster merger rate at a time $T\sim\\{0.5,\,3\\}~{}\mathrm{Gyrs}$ or
$z\sim\\{9,\,2\\}$ respectively for the two cluster types and a distribution
for the off-cluster merger rate very broad with time. In Table 2 we include
the estimated local merger rate of the different clusters assuming they
individually comprise all of the dark matter in our Universe. This way of
presenting the rates for individual clusters rather than cumulatively for all
clusters within the same type facilitates comparison with the rate we could
estimate if we observed any one of these clusters in our Universe. Given the
limited statistics available in some cases, these estimates are susceptible to
shot noise for the single clusters.

From the in-cluster mergers perspective, current Earth-based interferometers
are on the verge of being able to detect events coming from as far as $z=2$,
with examples in O3b going as high as $z=1.18^{+0.73}_{-0.53}$ in the case of
GW190403$\\_$051519 Abbott  _et al._ (2023). Current star formation models
Madau and Dickinson (2014) predict a peak in the merger rate corresponding to
the peak in the star formation rate ($z\sim 1.8$ = Madau-Dickinson 3.5 Gyr
after BB) plus a time delay of about a Gyr due to the time of BH collapse,
binary formation and subsequent merger ($z\sim 1.1-0.9$). Finding an excess
beyond this redshift would possibly hint towards one of the models of dense
clusters of PBH analysed here.

In the case of off-cluster mergers, the rate is smaller, and instead of
increasing, their redshift dependence decreases towards larger redshifts. This
characteristic could be utilized to distinguish between their primordial and
astrophysical origins.

Figure 6: We show the in-cluster and binary escapers merger rate for the cases
of $M\&A$ and $\sigma_{0.5}$, which have the highest number of mergers,
allowing for more accurate kernel density plots of the underlying merger rate.
The binary escapers’ merger rates are represented as solid lines with similar
colours as the filled density plots depicting the in-cluster merger rates. We
observe a general trend where the in-cluster mergers peak around 0.54 Gyr
after the Big Bang for $M\&A$ clusters and 3.3 Gyr for $\sigma_{0.5}$
clusters, well before a Hubble time, while the binary escapers’ merger rate
extends significantly beyond it.

###  IV.4 Mass distribution

Another approach to analyzing BBHs is through their mass distribution, which
is highly dependent on the original distribution of masses due to the strong
suppression of hierarchical mergers in realistic simulations. In Fig 7 we
present a scatter plot of the mergers in the ($M_{1}$,$M_{2}$)[$M_{\odot}$]
plane. We can differentiate the two cases: in-cluster mergers directly
identified by the code and depicted as circles in the figure and off-cluster
mergers that we evolve in post-processing. A noticeable clustering around
equal mass ratios is observed, deviating from the expected binary distribution
one would obtain by randomly drawing pairs from the initial mass function.
This suggests a genuine preference for $q=1$ embedded in the mechanisms
responsible for these mergers. This preference implies that the binary
formation cross-section has a maximum for equal mass ratios, while unequal
ratios result in more hyperbolic encounters. The BBH merger with the minimum
mass ratio identified sits around $\sim 0.2$ in the $\sigma_{1.5}$, not shown
in the plot, which is above the minimum mass ratio detected in a confident BBH
which is $\sim 0.1$ for GW190814.

Figure 7: In this figure, we present the distribution of the masses of the BHs
that merge within the age of the Universe. They are categorized into two
groups: those merging inside the clusters, identified as collisions by the
code (in-cluster), and those accounted for as binary escapers and evolved
outside of the simulation (off-cluster). Some clusters don’t produce any
merger and thus, not every color in the legend is represented in the plots.
The left plot corresponds to the $M\&A$ type and the right plot represents the
$\sigma_{0.5}$ models.

###  IV.5 Spin distribution

Initially, the spins of all the BHs are set to 0 inspired by the isotropic
collapse of the primordial overdensities predicted by most PBH formation
theories. This assumption, nonetheless, may also apply to stellar BHs as the
efficiency of the angular momentum transport from the spin of the progenitor
star is still under debate (e.g. Cantiello  _et al._ (2014); Maeder and Meynet
(2000); Fuller  _et al._ (2019)). The various BHs can acquire spin via merging
with other BHs. The final spin of the remnant is calculated using Jiménez-
Forteza  _et al._ (2017), which depends also on the progenitors’ spins and
masses. It is also possible that due to the multiple close encounters in dense
clusters, the spin of the BH population is induced stochastically, with a
final spin distribution peaked at zero with a dispersion of about
$\sigma_{s}=0.2$ Jaraba and García-Bellido (2021). We did not consider in this
work the induced spin due to close hyperbolic encounters in the distribution
of spins of the BH population.

In Fig 8 we present the distribution for the absolute magnitude of the spin of
the remnant BHs. The left panel illustrates the distribution for in-cluster
mergers, while the right panel depicts the same for off-cluster mergers. We
observe that the bulk of the distribution is similar for both cases, as the
initial spins for all BHs are 0, and thus, the final spin depends entirely on
the initial masses. The most frequent value for the final spin is around
$a_{\mathrm{f}}=0.68$ as you need to reach very unequal mass ratios $q<0.4$ to
go below $a_{\mathrm{f}}=0.6$. There are however, some extreme spins, sitting
around $a_{\mathrm{f}}=0.9$ due to $2^{\mathrm{nd}}$ generation mergers where
one of the progenitors already had a non-negligible spin. They represent a
$4\%$ of all the BBHs that coalesce while representing a larger fraction of
the total number of off-cluster mergers. The smallest final spin sits around
$a_{\mathrm{f}}=0.4$, corresponding to an extreme mass ratio of $q=0.2$ in the
$\sigma_{1.5}$ type of clusters and the largest one lays very close to
$a_{\mathrm{f}}=0.9$, coming from a second generation merger and high total
mass. This distribution aligns with general observational results and our
expectations that most, if not all, BHs in our universe are Kerr BHs, as there
are many mechanisms to gain spin but very few to lose it.

Figure 8: In this plot we show the value for the final spins of the remnant
BHs as a function of the binary mass ratios of the progenitors. The
simulations assume zero spins as the initial condition. We distinguish between
mergers happening inside the cluster and outside and also between 1st and 2nd
generation mergers. Spins are calculated using Jiménez-Forteza  _et al._
(2017) either by the code itself or in post-processing.

###  IV.6 Residual eccentricity distribution

Another parameter crucial for understanding the nature of BBHs is the residual
eccentricity of the binaries. Looking at fig 6, we can differentiate the in-
cluster and the off-cluster mergers as two very different cases. The in-
cluster mergers are the result of complex dynamical processes where more than
two bodies are usually involved. The main process by which the BHs end up
merging is not gravitational wave emission, and thus, we expect to find highly
eccentric signals as eccentricity cannot be radiated away. This could be a
distinct feature if we were to detect an eccentric merger at a sufficiently
high redshift as we would not expect in-cluster mergers to happen before the
peak star formation rate era. This kind of highly eccentric merger, however,
presents challenges on its own to be measured (see Bustillo  _et al._
(2021a)). Given the lack of a significant inspiral phase in these events, in-
cluster mergers would likely manifest in Earth-based gravitational wave
interferometers as short bursts, of which one example might be GW190521 Abbott
_et al._ (2020c). This comes with various challenges as the very low number of
detectable cycles leaves the floor open to numerous alternative hypotheses
Bustillo  _et al._ (2021b); Gayathri  _et al._ (2022); Romero-Shaw  _et al._
(2020).

The other case, the off-cluster mergers, follow a completely different route
and the only mechanism by which they end up merging is via orbit shrinking
through GW emission. As a result, eccentricity is radiated away Peters (1964)
before the binary reaches the detectability threshold in current Earth-based
interferometers. To investigate the validity of such a hypothesis, we plot in
fig 9 the residual eccentricity for the off-cluster mergers, that is, the
eccentricity that, based on quadrupolar GW emission, the binaries would have
at the moment of the merger. Most of the cases that would concern us ($z<2$)
present a maximum residual eccentricity of $e\sim 10^{-4}$, which is hardly
detectable given current sensitivity and the fact that it might be confused
with spin effects Romero-Shaw  _et al._ (2023). The detectability threshold
lies at around $e\sim 10^{-3}$ for a high enough SNR event, justifying the use
of the quasi-circular approximation.

Figure 9: In this figure, we display the residual eccentricity of the off-
cluster mergers. We define residual eccentricity as the eccentricity the
binary would possess at the time of merger, after it has been radiated away
throughout the inspiral phase. This residual eccentricity could potentially be
detected from ground-based interferometers.

##  V Conclusions

In this paper, we have reported the results of our study of the phenomenology
of BHs-only clusters with different initial mass functions based on both
astrophysical and cosmological assumptions for their origin.

From a dynamics point of view, we have encountered differences in the
evaporation as well as the expansion rate of the clusters. The main
distinction between cluster types can be understood via their stability, which
seems to depend highly on the mass ratios of the typical encounters occurring
inside of them and driving the energy exchange between the BHs. In this
regard, we find that single BH escapers’ distributions resemble very closely
the original IMF except for a larger skew towards smaller masses and the total
number of binary escapers seems to be correlated again with the stability
criteria so that, the more stable the cluster is, the more binaries can fly
from its gravitational potential well.

Regarding BBH coalescences, we distinguished between in-cluster and off-
cluster mergers due to their significant differences in orbital parameters and
potential detectability from Earth. On the one hand, in-cluster mergers are
produced via 3+ body encounters inside the densest regions of the clusters,
with a peak in the merger rate at a relatively high redshift, generally above
the detectability horizon of current interferometers, although this is highly
model-dependent. These mergers are also characterized by very high residual
eccentricity and few detectable cycles from Earth. On the other hand, off-
cluster mergers occur outside of the cluster due to gravitational wave
emission. This mechanism leads to a very low residual eccentricity
($e<10^{-4}$), leading to quasi-circular BBH coalescence as observed from
Earth. The masses of the BBH mergers depend very much on the IMF although we
find a trend favouring mass ratios close to 1 independent of the IMF. Given
our zero-natal spin assumption, most of the final spins depend on the
progenitor masses with a very low probability of non-zero progenitor spin
($<4\%$).

We also conducted a comparison between the $M\&A$ clusters with and without
considering BH merger kicks, finding a strong dependence on the existence and
number of hierarchical mergers. This finding provides valuable insight into
the accuracy of numerical simulations. Predictions regarding the occurrence of
these types of mergers cannot be relied upon if merger kicks are not
incorporated.

In summary, this work aims to elucidate the phenomenology of complex non-
linear systems represented by gravitationally-bound BH-only clusters. We hint
at some results as promising candidates for unexplained observations and
expand the world of possibilities in numerical research.

## Acknowledgements

We would like to thank Rainer Spurzem and Manuel Arcca Sedda for their
valuable comments and suggestions during the drafting of this manuscript. The
authors acknowledge the use of the publicly available codes: NBODY6++GPU. They
acknowledge support from the research project PID2021-123012NB-C43 and the
Spanish Research Agency (Agencia Estatal de Investigación) through the Grant
IFT Centro de Excelencia Severo Ochoa No CEX2020-001007-S, funded by
MCIN/AEI/10.13039/501100011033. JFNS acknowledges support from MCIN through
Grant No. PRE2020-092571. The authors acknowledge the use of the Hydra cluster
at the IFT, on which some of the numerical computations for this paper took
place.

## References

  * Collaboration  _et al._ (2019) T. E. H. T. Collaboration, K. Akiyama, A. Alberdi, _et al._ , The Astrophysical Journal Letters 875, L1 (2019).
  * Wyrzykowski, Ł.  _et al._ (2023) Wyrzykowski, Ł., Kruszyńska, K., Rybicki, K. A., Holl, B., Lecœur-Taïbi, I., Mowlavi, N., Nienartowicz, K., Jevardat de Fombelle, G., Rimoldini, L., Audard, M., Garcia-Lario, P., Gavras, P., Evans, D. W., Hodgkin, S. T., and Eyer, L., Astronomy & Astrophysics 674, A23 (2023).
  * Abbott  _et al._ (2023) R. Abbott, T. D. Abbott, _et al._ (LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration), Phys. Rev. X 13, 041039 (2023).
  * Abbott  _et al._ (2020a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020a), arXiv:2001.01761 [astro-ph.HE] .
  * Abbott  _et al._ (2020b) B. Abbott _et al._ , The Astrophysical Journal 896, L44 (2020b).
  * Abbott  _et al._ (2020c) B. Abbott _et al._ , The Astrophysical Journal 900, L13 (2020c).
  * Clesse and García-Bellido (2022) S. Clesse and J. García-Bellido, Physics of the Dark Universe 38, 101111 (2022).
  * Bird  _et al._ (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haimoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Clesse and García-Bellido (2017) S. Clesse and J. García-Bellido, Phys. Dark Univ. 15, 142 (2017), arXiv:1603.05234 [astro-ph.CO] .
  * García-Bellido (2017) J. García-Bellido,  _Proceedings, 11th International LISA Symposium: Zurich, Switzerland, September 5-9, 2016_ , J. Phys. Conf. Ser. 840, 012032 (2017), arXiv:1702.08275 [astro-ph.CO] .
  * Carr  _et al._ (2021) B. Carr, S. Clesse, J. García-Bellido, and F. Kuhnel, Phys. Dark Univ. 31, 100755 (2021), arXiv:1906.08217 [astro-ph.CO] .
  * Chisholm (2006) J. R. Chisholm, Physical Review D 73 (2006), 10.1103/physrevd.73.083504.
  * Chisholm (2011) J. R. Chisholm, Physical Review D 84 (2011), 10.1103/physrevd.84.124031.
  * Ezquiaga  _et al._ (2020) J. M. Ezquiaga, J. García-Bellido, and V. Vennin, JCAP 03, 029 (2020), arXiv:1912.05399 [astro-ph.CO] .
  * Ezquiaga  _et al._ (2023) J. M. Ezquiaga, J. García-Bellido, and V. Vennin, Phys. Rev. Lett. 130, 121003 (2023), arXiv:2207.06317 [astro-ph.CO] .
  * Animali and Vennin (2024) C. Animali and V. Vennin, (2024), arXiv:2402.08642 [astro-ph.CO] .
  * Belotsky  _et al._ (2019) K. M. Belotsky, V. I. Dokuchaev, Y. N. Eroshenko, E. A. Esipova, M. Y. Khlopov, L. A. Khromykh, A. A. Kirillov, V. V. Nikulin, S. G. Rubin, and I. V. Svadkovsky, The European Physical Journal C 79 (2019), 10.1140/epjc/s10052-019-6741-4.
  * Rubin  _et al._ (2001) S. G. Rubin, A. S. Sakharov, and M. Y. Khlopov, Journal of Experimental and Theoretical Physics 92, 921 (2001).
  * Gieles  _et al._ (2021) M. Gieles, D. Erkal, F. Antonini, E. Balbinot, and J. Peñarrubia, Nature Astronomy 5, 957 (2021).
  * Trashorras  _et al._ (2021) M. Trashorras, J. García-Bellido, and S. Nesseris, Universe 7, 18 (2021), arXiv:2006.15018 [astro-ph.CO] .
  * Inman and Ali-Haïmoud (2019) D. Inman and Y. Ali-Haïmoud, Phys. Rev. D 100, 083528 (2019), arXiv:1907.08129 [astro-ph.CO] .
  * Aarseth (2003) S. J. Aarseth,  _Gravitational N-Body Simulations: Tools and Algorithms_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2003).
  * Spurzem (1999) R. Spurzem, Journal of Computational and Applied Mathematics 109, 407 (1999).
  * Wang, Long and Spurzem, Rainer and Aarseth, Sverre and Nitadori, Keigo and Berczik, Peter and Kouwenhoven, M. B. N. and Naab, Thorsten (2015) Wang, Long and Spurzem, Rainer and Aarseth, Sverre and Nitadori, Keigo and Berczik, Peter and Kouwenhoven, M. B. N. and Naab, Thorsten, Monthly Notices of the Royal Astronomical Society 450 (2015), 10.1093/mnras/stv817.
  * Bagui  _et al._ (2023) E. Bagui _et al._ (LISA Cosmology Working Group), (2023), arXiv:2310.19857 [astro-ph.CO] .
  * Carr  _et al._ (2024) B. Carr, S. Clesse, J. García-Bellido, M. Hawkins, and F. Kuhnel, Phys. Rept. 1054, 1 (2024), arXiv:2306.03903 [astro-ph.CO] .
  * H.C. (1911) P. H.C., Monthly Notices of the Royal Astronomical Society 71, 460 (1911).
  * Khlopov  _et al._ (2002) M. Y. Khlopov, S. G. Rubin, and A. S. Sakharov, “Strong primordial inhomogeneities and galaxy formation,”  (2002).
  * Fragione and Loeb (2017) G. Fragione and A. Loeb, New Astronomy 55, 32 (2017).
  * von Hoerner (1957) S. von Hoerner, Astrophys. J. 125, 451 (1957).
  * Peters (1964) P. C. Peters, Phys. Rev. 136, B1224 (1964).
  * Sigurdsson and Hernquist (1993) S. Sigurdsson and L. Hernquist, Nature (London) 364, 423 (1993).
  * Zwart and McMillan (2002) S. F. P. Zwart and S. L. W. McMillan, The Astrophysical Journal 576, 899 (2002).
  * Heggie (1975) D. C. Heggie, Monthly Notices of the Royal Astronomical Society 173, 729 (1975).
  * Jiménez-Forteza  _et al._ (2017) X. Jiménez-Forteza, D. Keitel, S. Husa, M. Hannam, S. Khan, and M. Pürrer, Physical Review D 95 (2017), 10.1103/physrevd.95.064024.
  * Lousto  _et al._ (2012) C. O. Lousto, Y. Zlochower, M. Dotti, and M. Volonteri, Phys. Rev. D 85, 084015 (2012).
  * Lousto  _et al._ (2010) C. O. Lousto, M. Campanelli, Y. Zlochower, and H. Nakano, Classical and Quantum Gravity 27, 114006 (2010).
  * Campanelli  _et al._ (2007a) M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Merritt, Phys. Rev. Lett. 98, 231102 (2007a).
  * Campanelli  _et al._ (2007b) M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Merritt, Physical Review Letters 98 (2007b), 10.1103/physrevlett.98.231102.
  * Madau and Dickinson (2014) P. Madau and M. Dickinson, Annual Review of Astronomy and Astrophysics 52, 415 (2014).
  * Cantiello  _et al._ (2014) M. Cantiello, C. Mankovich, L. Bildsten, J. Christensen-Dalsgaard, and B. Paxton, Astrophys. J. 788, 93 (2014), arXiv:1405.1419 [astro-ph.SR] .
  * Maeder and Meynet (2000) A. Maeder and G. Meynet, Annual Review of Astronomy and Astrophysics 38, 143 (2000), https://doi.org/10.1146/annurev.astro.38.1.143 .
  * Fuller  _et al._ (2019) J. Fuller, A. L. Piro, and A. S. Jermyn, Monthly Notices of the Royal Astronomical Society 485, 3661 (2019), https://academic.oup.com/mnras/article-pdf/485/3/3661/28230562/stz514.pdf .
  * Jaraba and García-Bellido (2021) S. Jaraba and J. García-Bellido, Phys. Dark Univ. 34, 100882 (2021), arXiv:2106.01436 [gr-qc] .
  * Bustillo  _et al._ (2021a) J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forné, and J. A. Font, Phys. Rev. Lett. 126, 201101 (2021a).
  * Bustillo  _et al._ (2021b) J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forné, J. A. Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, and S. H. W. Leong, Phys. Rev. Lett. 126, 081101 (2021b).
  * Gayathri  _et al._ (2022) V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczepanczyk, I. Bartos, M. Campanelli, S. Klimenko, C. Lousto, and R. O’Shaughnessy, “Eccentricity estimate for black hole mergers with numerical relativity simulations,” (2022), arXiv:2009.05461 [astro-ph.HE] .
  * Romero-Shaw  _et al._ (2020) I. Romero-Shaw, P. D. Lasky, E. Thrane, and J. C. Bustillo, The Astrophysical Journal Letters 903, L5 (2020).
  * Romero-Shaw  _et al._ (2023) I. M. Romero-Shaw, D. Gerosa, and N. Loutrel, Monthly Notices of the Royal Astronomical Society 519, 5352 (2023).