The clustering dynamics of primordial black boles in $ N $-body simulations

We explore the possibility that Dark Matter (DM) may be explained by a non-uniform background of approximately stellar-mass clusters of Primordial Black Holes (PBHs), by simulating the evolution them from recombination to the present with over 5000 realisations using a Newtonian $ N $-body code. We compute the cluster rate of evaporation, and extract the binary and merged sub-populations along with their parent and merger tree histories, lifetimes and formation rates; the dynamical and orbital parameter profiles, the degree of mass segregation and dynamical friction, and power spectrum of close encounters. Overall, we find that PBHs can constitute a viable DM candidate, and that their clustering presents a rich phenomenology throughout the history of the Universe. We show that binary systems constitute about 9.5% of all PBHs at present, with mass ratios of $ \bar{q}_{\mathrm{B}} = 0.154 $, and total masses of $ \bar{m}_{\mathrm{T},\mathrm{B}} = 303 \mathrm{M_\odot}$. Merged PBHs are rare, about 0.0023% of all PBHs at present, with mass ratios of $ \bar{q}_{\mathrm{M}} = 0.965 $ with total and chirp masses of $ \bar{m}_{\mathrm{T},\mathrm{M}} = 1670 \mathrm{M_\odot} $ and $ \bar{m}_{c,\mathrm{M}} = 642 \mathrm{M_\odot} $ respectively. We find that cluster puffing up and evaporation leads to bubbles of these PBHs of order 1 kpc containing at present times about 36% of objects and mass, with hundred pc sized cores. We also find that these PBH sub-haloes are distributed in wider PBH haloes of order hundreds of kpc, containing about 63% of objects and mass, coinciding with the sizes of galactic halos. We find at last high rates of close encounters of massive Black Holes ($ M \sim 1000 \mathrm{M_\odot} $), with $ \Gamma^{\mathrm{S}} = (1.2 \substack {+5.9 \\ -0.9}) \times 10^{7} \mathrm{yr^{-1} Gpc^{-3}}$ and mergers with $\Gamma^{\mathrm{M}} = 1337 \pm 41 \mathrm{yr^{-1} Gpc^{-3}} $.


I. INTRODUCTION
One of the most pressing problems in cosmology is that of the nature of Dark Matter (DM), first suggested in the thirties by Zwicky from the observation of motion anomalies in the Coma galaxy cluster in Ref. [1], and by Rubin in the seventies from the observation of an excess velocity in the tails of galaxy rotation curves in Ref. [2]. Evidence for DM is now abundant, from structure formation, galaxy dynamics, lensing, and Cosmic Microwave Background (CMB) observations (see Refs. [3][4][5][6]), largely in agreement with a DM density of Ω c h 2 = 0.120 ± 0.001 (see Ref. [6]).
For four decades now there have been a myriad attempts to explain the nature of DM, with possible candidates spanning a huge difference in order of magnitude, with masses from O 10 −22 eV to O(100) M . From new weakly interacting particles (WIMPs) such as ultra-light axions (see Ref. [7]) emerging in some extensions of the Standard Model of particle physics, to hypothesised Massive Halo Compact Objects (MACHOs) such as massive Black Holes (BHs) (see Ref. [8]) populating galactic halos.
However, and in spite of the very numerous efforts, no viable DM particle candidate has been found in experiments like the Large Hadron Collider (LHC), nor a convincing signal has been detected either in ground or underground facilities, directly or indirectly, with cryogenic detectors such as the Cryogenic Dark Matter Search * manuel.trashorras@csic.es † juan.garciabellido@uam.es ‡ savvas.nesseris@csic.es (CDMS), or liquid Xenon detectors like LUX, nor any annihilation probes like Fermi have found an excess signal in the sky that may be attributed to it. So overall, no experiment has yet observed a statistically significant direct detection of a DM particle candidate [9,10]. However, the XENON1T detector recently reported an excess over known backgrounds consistent with the solar axion model, albeit only at a low statistical significance of ∼ 3.5σ [11]. Yet, ever since the first LIGO observation of Gravitational Waves (GWs) (see Refs. [12][13][14]) from the merger of O(10) M BHs and subsequent detections, and additional merges observed in LIGO runs O1 and O2, another hypothesis, that is, that DM may be partially or entirely made of Primordial Black Holes (PBHs), is gaining ground, and a number of models been proposed (see Refs. [15][16][17][18][19][20]).
Experiments like LIGO-VIRGO herald a new era of GW astronomy in which new observations previously unavailable in the electromagnetic channel become now visible in the gravitational one, thus opening a powerful and far reaching window to the Universe. In particular, it is expected that, should DM be made of PBHs, then this new observational probe would be open to explore DM in novel ways.
Observational evidence of a DM component of the Universe is abundant: from purely dynamical probes such as galaxy rotation curves and star velocity dispersions to weak and strong gravitational lensing of distant galaxies, cosmic microwave background spectrum of temperature anisotropies, structure formation history, type Ia supernovae distance measurements and galaxy, cluster and quasars surveys probes, such as the baryon acoustic oscillation scale, redshift space distortions and Lyman-α spectroscopy. All of these probes allow for the possibility that all of DM may consist of BHs in the stellar massrange, particularly if BH clustering and BH mass spectra are considered (see Refs. [30][31][32][33]).
If so, these PBHs may offer a rich phenomenology not only able to explain the nature of DM, but also to explain some questions regarding early structure formation [34], the radial profiles of galactic cores, ultra-faint dwarf galaxy dynamics, correlations between the X-ray and Cosmic Infrared Background, as well as BH merger rates, mass spectrum and low spins (see Ref. [35] for a review of these issues).
Our paper is organised as follows: In Section II we present the theoretical background of PBH formation, while in Section III we describe how our simulations were designed and performed, their properties, the choice of the Initial conditions (ICs) and the Evolution Snapshots (ESs) as well as their limitations. We also, for better intuition of our results, provide animations of our simulations in the public repository here.
Once performed, the analysis of the simulations can then be grouped in two broad categories, one regarding the separate populations that arise in the simulations, and another regarding the dynamics.
Concerning the former, Section IV is devoted to the identification of differentiated populations arising in the simulations, which are mainly those of clustered and ejected PBHs, bounding and mergers of PBHs, evaporation rates, parent trees and merger trees, while Section V deals with the close encounter dynamics where two PBHs passing sufficiently near to each other produce a merger pair or a binary pair.
Concerning the latter, Section VI is dedicated to the study of the mass, position, velocity and density profiles of PBHs, segregated by their population, as well as to the degree of mass segregation and dynamical friction, while in Section VII we examine the orbital distributions of PBHs and its repercussions for both either transient or stable binaries. In Section VIII we present estimates of the hyperbolic encounter and merger event rates for a medium-sized galactic halo and comoving cosmological volume.
Finally, in Section IX we summarize our findings and comment on the feasibility of PBHs as a DM candidate according to our results.

A. Formation mechanisms
First proposed in the sixties when Zel'dovich and Novikov (see Ref. [36]) realised that BHs could very well have been formed in the early Universe, catastrophically growing by accreting the surrounding plasma during the radiation era from the extremely high gravitational forces. It was later, in the seventies, that Hawking (see Ref. [37]) and Carr (see Refs. [38,39]) found that a sufficiently overdense region or inhomogeneity in the early Universe could undergo gravitational collapse providing the seeds for formation of these BHs, more commonly referred to as PBHs in this context.
Several mechanisms have been proposed in the literature to form these PBHs in the early Universe (see Refs. [40,41]), such as domain walls, cosmic strings, and vacuum bubbles. In other models, however, they are formed by the gravitational collapse of overdense regions in the Universe either in the radiation era or, less frequently, the matter era. The existence of such overdensities is a consequence of certain inflationary models, and could be, if discovered, used to probe dynamics during inflation.
A typical scenario in which PBHs may form is as follows. During inflation high peaks in the primordial curvature power spectrum produce overdensities that grow and collapse gravitationally in the radiation era [41]. These high peaks are produced naturally by a number of inflationary models, such as in hybrid inflation by having a phase transition at the symmetry-breaking field (see Ref. [15]), single-field inflation with an inflection point (see Refs. [42][43][44]) or by having the inflaton coupled to gauge fields (see Refs. [45,46]). The clustering and mass spectrum of these PBH is highly dependent on the shape of the primordial curvature power spectrum; a narrow peak will typically produce a nearly monochromatic spectrum of PBHs, while a broader peak will produce PBHs with a wider mass range of O(0.01) M − O(1000) M . This mass power spectrum could then evolve by hierarchical merging (see Ref. [47]) right after recombination, but more importantly, from PBH masses that would greatly increase by accreting baryons from the dense plasma in the early Universe. The largest of these PBHs would then become Super Massive Black Holes (SMBHs) and InterMediate Mass Black Holes (IMBHs), provide the seeds for galaxies and accelerate cosmological structure formation at early times. These PBHs would still have a very small cross section with matter and are an ideal candidate for DM.
More generally speaking, PBHs may form with different masses greater than the Planck mass M P = 1.094 × 10 −38 M (see Ref. [48]) and with very low spin at different epochs as the size of the cosmological horizon increases, forming a broad spectrum, but their most frequent masses must be smaller than O 10 4 M at formation, since otherwise they would have altered significantly the early Universe dynamics by injecting energy to the medium before recombination which would be apparent in the Cosmic Microwave Background.
However, available PBHs masses in the late Universe are shifted with respect to their masses at formation due to Hawking evaporation, merger history and accretion. At present, PBHs masses may range from the largest supermassive BHs at the cores of the largest galaxies of about O 10 11 M , up to the tiniest possible BHs with masses larger than O 10 -18 M , which is the smallest non-evaporated BH mass possible at present (see Ref. [49]).
Note that this mass range is much broader than that of the more familiar astrophysical BHs formed with large spins as the remnants of the most massive stars in corecollapse type II supernova, and whose masses are restricted to be heavier than about 5 M but lighter than O(30) M (see Ref. [50]). Therefore, an observation of very low spin BHs and BH masses out of this range would provide significant support to the PBH hypothesis. This is particularly true for the lower bound, as there is no known astrophysical mechanism that may create BHs in the sub-solar mass range, while, through merger, spinning BHs in the LIGO-VIRGO mass range are possible, even if very unlikely, although impossible without spin, specially in the high mass range.

B. Observational constraints
There are many observational constraints on uniformly distributed PBHs which restrict their mass distribution. These constraints can be grouped into four different categories: i) Evaporation constraints arising from the fact that all PBH with a mass of less that O 10 −19 M would have evaporated by the Hawking mechanism at present-day (see Ref. [51]) ii) Lensing constraints from femtolensing of gammaray bursts (see Ref. [52]) and millilensing of compact radio sources (see Ref. [53]), as well as microlensing constraints from the EROS collaboration (see Ref. [54]), the OGLE collaboration (see Ref. [55]), HSC Andromeda observations (see Ref. [56]), Kepler observations (see Ref. [57]) and caustic crossings (see Ref. [58]).
Moreover, while a monochromatic PBH mass power spectrum of PBHs is excluded by observations, this kind of spectrum is not very realistic since gas accretion and merger histories will both broaden and shift the mass power spectrum over time, not mentioning that PBHs can form in a broad spectrum of masses to begin with, as previously mentioned. In either case it has been found that a PBHs mass spectrum spanning a few orders of magnitude in the range of O(0.01 − 1000) which for any particular mass does not reach more than 40% of the DM energy density is still viable (see Refs. [15,20,35,69,70]).
Also, these constraints are calculated for a monochromatic mass distribution. It has been found that a recalculation of these constraints for a wide mass distributions tends to have the effect of relaxing some of these exclusion curves (see Refs. [32,33,71]).
In addition to this, the constraints act on different epochs in the history of the Universe. For instance, most lensing and dynamical constraints such as wide halo binaries disruption act on redshifts z ≈ 0, constraints from type IA supernovae act on redshifts z ≈ 1, constraints from X-rays and the Cosmic Infrared Background act of redshifts z ≈ 10, and accretion constraints from the Cosmic Microwave Background act on redshifts O 1000 − 10 4 , see Ref. [72].
Overall, this amounts to a picture where PBHs could survive these constraints and explain all of DM either if the mass spectrum evolves in redshift through merging and accretion, enough not to be excluded by these constraints from accretion at early times as well as those from lensing and dynamical constraints at late times, PBHs mean masses growing from O(0.01 − 10) M in the early Universe to O(10 − 1000) M in the late Universe.
Moreover, clustered PBHs in the stellar mass range can comprise the totality of DM. Present constraints indicate that a monochromatic distribution of PBHs can account for 100% of DM in the stellar mass range of O(10) M , once one takes into account that clustering relaxes lensing constraints significantly [32,71].
More importantly, since these PBHs may form over an extended mass spectrum, it can be found that PBH may comprise the totality of DM even if a fraction of these objects have masses greater than O(100) M , where the accretion, dynamical friction, wide binaries disruption and ultra faint dwarfs constraints start to be significant for monochromatic PBH distributions. The LIGO-VIRGO discovery of BHs in the stellar mass cannot for the moment discern if these objects are of primordial or astrophysical origin.
That however, might change in the near future as it has been proposed that bursts with millisecond durations could be explained by the GW emission from the hyperbolic PBHs encounters in dense clusters, such as the ones considered this analysis, see Ref. [73,74]. These encounters have a very characteristic spectrum that could allow detection by both AdvLIGO and, in the future, LISA depending on the duration and peak frequency of the emission. Studying in detail these hyperbolic encounters and simulating these dense PBHs clusters is one of the main motivations of our work and we will present our results in detail in Sections V and VIII.

III. METHODOLOGY
We now proceed to describe our simulations, the essential properties of which are outlined in Section III A, the particular choice of ICs in Section III B and the information available in the time slices in Section III C along with comments on the advantages, disadvantages and limitations of our methodology.
In order to compute the trajectories of the gravitating bodies in our simulations we have made used of ASTROGRAV-3.7 (see Ref. [75]), a Java numerical Nbody code that makes use of a Verlet integration algorithm: a time-reversible symplectic integrator method for solving Newton's equations of motion with good numerical stability. As we discuss also later, this code is purely Newtonian and does not include any Generla Relativity (GR) effect, such as GW emission.
In particular, we set a single configuration of N O = 1000 bodies, and let it evolve throughout the entire history of the Universe, with ICs determined by the bodies' masses, positions, velocities and radii. These parameters, amounting to 7 degrees of freedom, are chosen as follows. iii) For the 3D velocity vectors: 3 degrees of freedom, taken again from a random realisation of an isotropic, centred Multi-Normal distribution, whose parameters can be determined under some assumptions by making use of the Virial Theorem.
iv) The radius correspond to the Schwarzschild radius for the object mass, which is completely determined by the object mass, and does not add any extra degrees of freedom.
Then, we repeat the process for N R = 5000 realisations in the same IC configuration, in order to improve our statistics. Note that these choices and their approximate range of values agrees with some inflationary models that produce PBHs [15]. Note however that, as we will see in Section III B, it is more practical for both positions and velocities to use instead of isotropic Multi-Normal distributions, the equivalent Maxwell-Boltzmann (MB) distributions over an isotropic 3D vector field.
The code then computes the evolution of the gravitational system in a purely classical framework within a static background, unlike other codes such as GEVOLUTION-1.1 (see Ref. [76]), a full relativistic code, or GADGET-2.0.7 (see Ref. [77]), a cosmological hydrodynamical code which adapts classical dynamics to an expanding background. The choice of ASTROGRAV-3.7 was made as it has a series of advantages over the aforementioned two codes.
First, the spatial regions where these PBH clusters may form are dense enough that they quickly decouple from the expansion of the Universe and stay decoupled throughout the DM and dark energy eras, so that the relativistic effects in an expanding background would not add up to significant increase in accuracy.
Second, the code computes the system evolution much faster, since the orbital dynamics of the orbiting bodies are treated classically, the full GR treatment being computationally far more time consuming. In the face of the large number of bodies in each realisation N O = 1000, the large number of realisations per model N R = 5000, and the long evolution times of O 10 10 yr, makes the simulation time-efficiency an item of crucial importance.
Third, even if GR effects are not properly simulated in the code, such as precession of the periastron, framedragging, binary spin-coupling, etc, are relevant in the strong field regime, the PBH clusters in our simulations are disperse enough to begin with and only disperse further over time, that only a very low number of these close approaches do occur in between time snapshots. This is so since the mean distance in between objects is typically O(0.1) pc and the mean minimum distance over the whole simulation time is around O(0.01) pc, which for the range of masses simulated in our code, of order O 10 1 M , leaves the magnitude of such effects completely negligible for the vast majority of trajectories.
Last, body collision and merger are, unlike in the alternatives, already implemented in the code, a feature of particular importance in the densest PBH clusters even if events of this kind are very rare, having found that the probability of a particular PBH to merge during the full time interval covered by the simulations is of O 10 −5 .
However, the collisions are treated classically in the code, and do not exhibit inspiralling binaries and loss of kinetic energy that leads to the emission of gravitational energy. The latter increases the likelihood for mergers as it tends to draw the bodies closer to each other, meaning the total amount of mergers in our simulations is underestimated and must be understood as a lower bound of the actual merger event rate.

A. Simulation properties
We proceed now to comment on the time and space resolution of our simulations, their implications for the data extraction, and constrain the range of validity where we may trust our findings, focusing in particular on whether the evolution of PBHs is dynamically constrained to a range where the code behaves as desired.

Time-resolution
The code has an adaptive time resolution δt, which depends on the forces at i th -step and is not to be confused with the time interval at which data is regularly stored, ∆t. Dependent on this δt, that is, the time-step in between the iterations of the Verlet integrator the global error in position and velocity grows with (δt 2 ).
We have, however, performed tests prior to the simulations to ensure that the error during the whole evolution period does not accumulate enough to meaningfully alter the simulation output by re-running the simulations with time-steps smaller than in the final simulations by a factor of a hundredth, with no significant effect, with position and velocities typically differing by less than 0.01% from their original values in the final output, but at the expense of both CPU time and memory.

Space-resolution
The code is a pure N -body integrator and therefore there is no theoretical minimum resolution or spatial separation as in other hydrodynamical codes. However, we can can consider a minimum spatial resolution, approximately corresponding to the maximum impact parameter derived from the effective cross section of two orbiting PBHs to become gravitationally bound. This relativistic constraint is due to the power emission of the incoming body, which is absent in the codes classic framework. The threshold cross section σ (see Ref. [78]) for the PBH pair binding is where b is the limiting impact parameter and the Schwarzschild's radius is Should the incoming object approach hyperbolically another object with a smaller impact parameter than the threshold, the simulated hyperbolic trajectory would clearly not track the actual bounded trajectory that one would observe. Considering the most unfavourable case in the range of masses and velocities in the simulations (see Figure VI.13a and VI.13b), then Eq. (1) renders an effective spatial resolution of ∆x = 100 AU, less than the mean minimum distance in-between bodies, of O(1000) AU.

B. Simulation Initial Conditions
We proceed now to describe the generation of the simulation ICs, and the choice of distributions for all the fundamental parameters within it. The dynamics of each of the i gravitating bodies, where i = 1, 2, ..., N O , in the simulation box is fully determined by the simulation ICs. The number of bodies in the simulation box is N O = 1000, and 4N O free parameters are enough to completely determine the realisation IC and its subsequent evolution i) N O free parameters for the mass, m i , one per body in the simulations, chosen from a random sample of a global Log-Normal distribution.
ii) 3N O free the position coordinates, x i , chosen from a random sample of a global Maxwell-Boltzmann distribution for the position moduli with isotropic randomised directions.
These parameters completely determine i) N O dependent parameters for the radii of the objects, r i , corresponding to the Schwarzschild radius of the individual body, thus naturally following a Log-Normal distribution of mass.
ii) 3N O dependent velocity coordinates, v i , chosen from a single random sample of a shell-specific isotropic Maxwell-Boltzmann distribution, determined by the number of bodies, N O , the individual masses, m i , and the individual positions x i , or conversely, an isotropic shell-specific Multi-Normal distribution, where by shell-specific we mean radially dependent with respect to the PBH cluster barycentre.
The choice of mass, position, velocity and radius parameters is then repeated for each of the N R = 5000 different realisations that comprise the simulations. An example of the cluster's position and velocity distribution is shown in Figure III.1. This choice of IC leads to an initial density ρ 0 ≈ 3N O µ m /4πa 3 x of O(1000) M pc −3 (see Figure VI.12a). In the remaining of this subsection, we proceed to detail the particular choice of the underlying distributions from which this parameters are extracted. A summary of the input distributions and realisations is given Table III.1.

Mass distribution
The initial mass profile D m (m i ) 0 of the cluster is obtained from a random realisation of a Log-Normal distribution for each realisation       where ∆m i = (log m i − µ m ) is the PBH logarithmic mass deviation from the mean and N m = (σ m √ 2π) −1 is a normalisation factor. The total mass in each realisation   Table III.1 and an example realisation, from the logarithmic mean µm = 2.0 M and deviation σm = 1.5 M . The thick, normal and thin grey lines signal the distribution mode, median and mean respectively. The orange and red lines signal the realisation median and mean respectively. Note the large hierarchy between modal, median and mean masses due to the large m → ∞ tail of the distribution, skewing the mean to large values, which results in a large variance of the total mass per realisation. The total mass is mtot is then, on average, M T = N O µ m . In particular, for our simulations, µ m = 1.5 M and σ m = 2.0 M , as laid out in Table III.1 and shown in Figure III.3.
In the previous distribution, however, one should be reminded that µ m is the mean of the natural logarithm of the mass, and not the mean mass itself. Since the actual mean massm is be a parameter of great importance in the computation of the background number of PBHs in a typical galactic halo and it is in any case of greater physical meaning, it will be computed from the Log-Normal distribution parameters µ m and σ m with the expressionm while, similarly, the actual variance of the mass can be computed from the distribution parameters with where in both cases the over-bar indicates that the quantity is computed in the linear (physical) scale. In the linear scale one can similarly extract the median value of the mass,m as well as the the modal value of the mass, as a function of the Log-Normal distribution parameters. Note that, generally speaking, a Log-Normal distribution has wide tails and as suchm >m >m.

Position distribution
The initial position profile D x (x i ) 0 of the cluster is obtained from a random sample of a 3D Multi-Normal distribution for each realisation, with density where ∆ x i = ( x i − µ x ) is the PBH position displacement from the cluster center of mass and N x = ( 2π |Σ x |) −1 is the normalisation factor. In our simulations, we consider purely spherical Σ x = σ 2 x 1 3 and centrally aligned µ x = 0 clusters of PBHs.
In practice, this profile is equivalent to that of isotropic vector field with the norms following a Maxwell-Boltzmann distribution, with the Maxwell-Boltzmann scale parameter a x given by where µ x is the mean position of the distribution. The variance of the position distribution can be directly computed from the scale parameter to which it is directly proportional with while the median position of the distribution can be computed from the scale parameter with the similarly linear relationx   Table III.1 and an example realisation, from the mean µx = 1.6 pc and scale parameter ax = 1.0 pc. The thick, normal and thin grey lines signal the distribution mode, median and mean respectively. The orange and red lines signal the realisation median and mean respectively. To to obtain the radial position abundance one needs only to multiply the previous distribution by the surface element dS = 4πx 2 , which enhances the tail. Note the hierarchy between modal, median and mean masses is much reduced with respect to the mass case, due to the smallness of the x → ∞ tail of the distribution, which results in a small variance of the total simulated volume covered in the ICs. and the modal position of the distribution is then as a function of the scale parameter. Note that a Maxwell-Boltzmann distribution has as wellx >x >x. Finally, the probability distribution density D x (x i ), as a function of the mean position µ x , is then given by where µ x is the mean position of the distribution and N x = 32/(π 2 µ 3 x ) is a the normalisation factor. This equivalence is due to the fact that the χ distribution for k = 3 degrees of freedom, i.e. the dimensionality of the simulation box, of the vector modules reduces precisely to the Maxwell-Boltzmann distribution, For simplicity, we compute the IC with the latter Maxwell-Boltzmann distribution as it is more computationally economic given the symmetry of the problem. In particular, the simulated PBH clusters have been chosen such that, µ x = 1.6 pc and a x = 1.0 pc, (see Table III

Velocity distribution
The initial velocity profile D v (v i ) 0 of the cluster is obtained as well from a random sample of a Maxwell-Boltzmann distribution for each realisation, which, as a function of the mean velocity µ v is given by where N v = 32/(π 2 µ 3 v ) is a normalisation factor. Note that both the distribution scale parameter, the velocity variance and the modal can be similarly obtained as it was the case in the position distribution in Eqs. (11), (12) and (14) respectively with the substitution x → v, being the same underlying distribution.
Assuming that the position distribution is a nearly homogeneous, spherically-symmetric matter distribution, that the cluster is sampled well enough that the distribution of objects traces accurately enough the prior distribution and that the Virial Theorem holds, then the mean velocity is calculated from the rotation curve of the galaxy, shown in Figure III.2. Such is the case of our simulations, where the masses of the individual PBHs in the cluster do not vary wildly in between them and there are enough of objects to accurately sample the distribution, with the mass evenly distributed across the cluster in a monopolar distribution.
The mean rotational velocity of the distribution, v i , in shells of radius x i ∈ (0, x max ) from the cluster barycentre is where the average cumulative mass, m i , in spheres of radius x i ∈ (0, r max ) from the cluster barycentre is then  Table III.1 distributions, at 64 linearly spaced radial bins in the range 0 ≤ x ≤ 4σx. The latter is obtained by integration of the shell-specific Maxwell-Boltzmann velocity distribution with a density kernel corresponding to the bin surface area dS = 4πx 2 , from x = 0 to x = 8σ. Both the core and the outskirts of the PBH cluster roughly understood as the regions with 0 ≤ x ≤ ax and 2ax ≤ x ≤ 4ax respectively) are not heavily populated, the core exhibits very low velocities in the range v ∈    Table III.1 and an example realisation, from the mean µv = 2.0 km/s and scale parameter av = 1.2 km/s. The thick, normal and thin grey lines signal the distribution mode, median and mean respectively. The orange and red lines signal the realisation's median and mean velocities respectively. Note, again, the hierarchy between modal, median and mean masses is much reduced with respect to the mass case, due to the smallness of the v → ∞ tail of the distribution, which results in a small variance of the total simulated velocities in the ICs.
It is more frequent, however, to express the mean velocity in terms of the Maxwell-Boltzmann scale parameter, a v , which, from Eq. (18) and Eq. (19), is given by Then, the initial, shell-dependent velocity distribution function, can be obtained by replacing the in velocity distribution of Eq. (17) the shell-dependent scale param-eter using Eq. (20) It can be shown, once Eq. (20) is integrated out to x i → ∞ that the resulting global velocity distribution is not a Maxwell-Boltzmann distribution, but that it can be well approximated by one with µ v = 2.0 km/s and a v = 1.2 km/s, (see Table III.1, Figure III Additionally, it should be noted that the individual positions and velocities of each object were displaced by means of a suitable change of frame of reference such that the PBHs in the simulation box, in each realisation, have zero total linear momentum i P i = 0, (21) so that the barycentre of all particles remains static throughout the simulation and the PBH cluster in particular does drift minimally from its initial position as it expels PBHs on one-to-one encounters. This transformation is performed as well to solve one practical problem, and that is that it reduces the computational time and data storage needed for the analysis of the output, since tracking individually the cluster and ejected PBHs with respect to the separate cluster frame of reference is faster in our analysis pipeline if the simulation as a whole is in the rest frame. This is because, as the positions and velocities transformation from the simulation box frame of reference to the cluster one skips the intermediate step of computing the simulation box barycentre at each time slice and it is precisely operating with position and velocities what takes up most of the computational time and storage.
We do as well perform as well a suitable rotation of the frame of reference so that there is zero total and angular momentum for similar reasons, so that there is no bulk total rotation throughout the simulation and the cluster and ejecta PBHs develop a minimal intrinsic rotation themselves as PBHs expelled from the cluster in one-to-one encounters. Note that again, the particular choice of a frame of reference in which the total angular momentum is zero has no impact on the simulation observables or dynamics, as it should be the case.

Body radii choice
In the code, mergers are treated classically and thus a merger occurs when the distance between the bodies centre of mass is less than the sum of their radius, set in Section III B as their Schwarzschild radius. Therefore, if two bodies (m i , r S i ) and (m j , r S j ) with i, j ∈ (1, 2, ..., N O ) approach each other to less than the threshold below which, a merger is then registered by the N -body code.

C. Simulation evolution snapshots
Once the simulation starts, N t − 1 = 64 time slices are produced and the masses, positions, velocities and other orbital parameters are recorded for all objects, until simulation time t 64 = 1.38 × 10 10 yr is reached.
Note that the position and velocity coordinates will henceforth be shown with respect to to the cluster barycentre, in the hereby called Cluster Frame (CF), a non-inertial frame of reference w.r.t the original Simulation Frame (SF), and which is calculated on a snapshotby-snapshot basis considering only the cluster objects, with positions an velocities transforming as and where the position and velocity centres of mass are computed, time by time-slice, with where δ Zi,C = 1 if object i in the simulation remains bounded to the cluster or δ Zi,C = 0 if said object has been ejected and is no longer bounded to the cluster. Then, by definition the Cluster Frame origin does not drift away from the cluster over time, unlike the case with the Simulation Frame, because of the ejected PBHs carrying mass away from the cluster in a non-perfectly isotropic manner.
We proceed now to describe the criteria chosen for the output choice of the time slices, the simulation averaging of the different realisations in order to provide our predictions with reliable confidence bands, the distinct phases through which the simulations go, as well as the algorithm design to reconstruct the merger and parent trees of each realisation.

Simulation period
All realisations of the simulated PBH cluster start with the IC at t 0 = 0 yr and end with the last snapshot exactly at t 64 = 1.38 × 10 10 yr, where 64 is the total number of evolution time-slices, the last of which corresponding to the current age of the Universe t U = (13.800±0.024)×10 10 yr according to Planck 2018 (CMB TT, TE, EE+lowE spectra, see Ref. [6]).
In order to capture the fast dynamics of PBHs at early simulation times and the much slower dynamics at late times, we settle then for a data collection algorithm in which the data collection time is updated in every timerun being multiplied by a factor of ten. Following this scheme, the simulation evolution is arranged in the following manner: i) Each simulation is divided into N r = 7 consecutive intervals, or runs, the first taking as the starting Simulation time-runs:   ii) Each successive run lasts for ten times the duration of the preceding one and outputs data every ninth of the run time, so that each run consists of nine non-overlapping time slices linearly binned in time.
iii) It can be shown, then, that per each realisation, there are in total, N t = 65 time-slices: 64 timeslices corresponding to the time slices at t j ∈ 13.8× [1000, 10 10 ]yr, the first of which is the effective IC for the first run, and an extra time-slices at the beginning for the global IC at t 0 = 0 yr, as detailed in Table III.2.

Evolution tests
It has been checked for a small number of realisations that both total mass i M i (t), total energy i E i (t), linear momentum i P i (t) = 0 and angular momentum i L i (t) = 0 are conserved for all times in the simulation period, in accordance with the ICs, where both quantities were set to zero by choosing an appropriate frame of reference as shown in Eq. (21) and Eq. (22). In particular we find that energy conservation is upheld up to O 10 −6 and stable in time evolution, while total linear and angular momenta are conserved to within 0.01% with respect to their corresponding initial values.

Realisation averaging
In order to better quantify and constrain the evolution of the simulated PBH clusters and in particular of their masses m, radii r S , absolute positions x = (x 1 , x 2 , x 3 ) and velocities v = (v 1 , v 2 , v 3 ) with respect to the simulation box frame that have been set at t 0 in the previous Section III B, we make exactly N R = 5000 different realisations, each having for the IC a different random realisation of the same probabilistic distribution of the physical parameters that fully determine the ICs. As will be later described in Section V, two distinct phases can be easily differentiated during this time evolution. Borrowing from similarly looking features in Monte Carlo Markov Chain simulations we describe the first phase as a "burn-in" (BI) period, followed by a "quasistatic" (QS) phase. The algorithm by which the time of transition between the burn-in period and the quasistatic phase of the simulation is found will be detailed next.

Burn-in period
One can easily observe in Figure IV.7 that the very first few stages of the time evolution, usually covering a fraction no more than the first O 10 −4 of the simulation time, show an exceedingly fast dynamic with respect to later times, forcefully expelling a large number of PBHs from the cluster to the background.
This feature occurs because, unlike other seemingly similar structures to this PBH clusters, like ordinary stellar globular clusters, the density in the ICs of this PBHs clusters is very large, much more so than in the case of stellar globular clusters.
In particular, in the simple model where there is no realisation variance and the IC random sample perfectly traces the IC distribution moments, at t 0 the total sim-ulation mass is mostly contained in a spherical region with In such a case, the mean total density of the PBH cluster at the initial time is in line with the mass density of the central regions of globular clusters, of ρ ≈ 1000 − 10 5 M pc −3 (see Ref. [79]), and orders of magnitude above that of the stellar density in the solar neighbourhood, ρ ≈ 0.040 ± 0.002 M pc −3 (see Ref. [80]). This results in that while PBHs move through their trajectories in between the first few snapshots at 0 yr ≤ t < 10 6 yr they observe large local gradients in the PBH cluster gravitational potential, often enough to expel these objects with speeds greater than the escape velocity of the cluster within their respective radial shell, and helping them escape as ejecta to the background. For similar reasons, the majority of hyperbolic encounters found within the simulations do occur within this short period, due to the much smaller mean inter-object distance in between PBHs. This is ultimately due to the fact that long term cluster stability imposes a harsh boundary to the minimum spatial volume of a cluster for a given mass. It is therefore natural, as it is also seen in Section VI, that after this brief period of fast evolution and complex dynamics, the cluster puffs-up and settles in the next stage.
The particular time at which the burn-in period is considered to transition to the quasi-static phase is found as follows. We compute the cluster population N C (t) and that of ejected objects N E (t) at any given simulation snapshot, and average this quantity over the N R = 1000 realisations so that the resulting evolution is well behaved and constrained, quasi-monotonously decreasing and increasing respectively for the two populations as it is shown in Figure IV Prior to this transition time t BI , there is essentially no evolution and their numbers remain static; after this the evolution of the cluster and ejected objects fit well to a power-law with α C > 0 for a decreasing cluster population, α E < 0 for an increasing ejecta population, and β C > β E > 0 since the initial ejecta population is for all realisations exceedingly small, sometimes nil. Due to the discrete nature of the output, the resulting transition time t BI at which the change between these two different regimes cannot be found exactly, but it can be estimated. At each time slice, the evolution of the cluster and ejecta populations is fitted to the logarithmic law of Eq. (31), from that particular time slice time to the end time of the run, of which there are N r = 7, in which such time slice is found, for all time slices N t − 1 = 64 time-slices.
The motivation to perform this fits in a run-by-run manner is to avoid overfitting some regions to the detriment of others since the separation in between time-slices varies by seven orders of magnitude in between the first and the last fit, and time-slices are globally neither linearly spaced nor logarithmically spaced, but rather a combination between the two. After the fitting, the array of times and determination coefficients, (t, r 2 ), is interpolated with a cubic spline, again in a run-by-run manner.
The resulting interpolated curve r 2 (t) generally has low values for t → t 0 , and unit values for t → t 64 . Inside each run i ∈ [1, 2, ..., N I ], the interpolated curve r 2 (t) will tend to have lower values for t → t 0,i , and higher values for t → t 64,i . The burn-in time, t BI , is then, the time where there is a first run with a maximum value of the r 2 (t) curve within 1% from the maximum of the range in the subsequent runs j ≥ i + 1 the results of which for the different populations and their combinations are shown in Table III.3. Note that, in all cases, t BI ≤ 5.93 × 10 6 yr.

Quasi-static phase
After the cluster has undergone the puffing up phase and expanded by a factor that varies greatly on the ICs, but usually up to O(10) during its very first stages, then the density decreases significantly and becomes comparable to that of stellar globular clusters.
Starting at this point, the evolution is much slower than in the previous phase. The fraction of merged objects per realisation barely amounts to 117 mergers out of a possible total of N R N O = 10 6 so there are barely any mergers in each of the individual clusters in the period that we have considered, so the number of objects states almost constant in our simulations.
Also starting at this point, the core of the cluster attracts the most massive and slower objects, sending the least massive and fastest of this objects to the periphery of the cluster through the well-known processes of mass segregation and dynamical friction, in a manner qualitatively similar to that of stellar globular clusters, but enhanced with respect to these due to the larger masses involved.
In this stage and for the rest of the PBH cluster evolution, the vast majority of would-be sources of gravitational radiation would come from hyperbolic encounters between pairs of objects, often ending in the slingshot and ejection of one of the PBHs at high speeds out of the cluster, resulting in a slow cluster population decrease, while releasing high-frequency GWs. Be reminded that our code does not include GR effects and will source cluster evaporation only as a consequence of close Newtonian 2-body encounters resulting in the ejection of a PBH, again in a very similar manner to the case of stellar globular clusters, with no emission of GWs and the associated velocity dampening. This slow cluster depletion of objects increases the matter density of the background at the expense of the cluster's own, which results in a bimodal distribution in the DM power spectrum.
Also in this stage a number of PBH binary systems are formed inside the cluster, with their survival depending on their possible ejection. Cluster binaries are rapidly disrupted by third-object encounters, usually lasting less that τ CB ∼ O 10 7 yr. However, having escaped to the much depopulated background, ejected binaries are stable and long lived since the time to merger is well approximated in Ref. [81]. In this scenario, only relatively infrequent massive, close PBHs could slowly inspiral into each other within a Hubble time, while emitting lowfrequency gravitational radiation with regular peaks at the periastron crossing.

Parent trees computation
In each time-slice we locate the possible bounded object sub-systems and their coordinates in their Orbital Frame (OF) of reference. Note that the orbital coordinate frame of reference is one centred on the barycentre of all the inferior objects of the sub-system, that is, all objects excluding those outside the shell defined by each PBH and centred around the cluster. This method of choice for the frame of reference is best when, as this is the case, no object is significantly more massive than all of the others combined, and the number of ejected objects dripping from the cluster grows overtime.
This way we identify sub-systems of objects within the simulations, corresponding to binary and multiple systems, either within the cluster or the ejecta, and even in a few brief instances sub-sub-systems corresponding to binaries orbiting a binary within the simulations. As for how the parent object is identified, in each time snapshot we find the parent tree of the simulation objects, that is, for all objects we find the most massive object that they are orbiting, in order to identify multiple systems. Depending on the different population to which the object may belong, then i) Cluster isolated PBHs do orbit the cluster barycentre, alone, and not as a part of more compact bounded sub-systems. They will then be assigned as their main-parent the cluster most massive PBH, centrally located, unless they they show up in bounded sub-systems in which a sequence of subparents follows this tracing the hierarchical substructure of the subsystem.
ii) Ejected isolated PBHs will always be assigned as their main parent again the cluster most massive, centred object unless they show up in bounded subsystems in which then a list of sub-parents follows just as it was the case for cluster objects. This is so even if for all other effects ejecta PBHs are not part of the cluster and exert only very weak and ever decreasing gravitational influence over it.
iii) Only the remaining fraction of cluster and ejecta objects, typically of order less than 10% for each of these respectively, do show up in the form of binary pairs of PBHs, and much more rarely, in the form of binary pairs already within a sub-system, with their sub-parent automatically assigned to the subsystem most massive object, and their parent tree tracing back to the aforementioned main parent of the cluster.
When an object belongs to a sub-system within a subsystem of say, the cluster or ejecta, we label the former as a 2 nd generation system and the latter we label a 1 st generation sub-system, while the cluster or ejecta is 0 th generation. Note that the parent tree registers all the hierarchy up to the most massive, centred object in the cluster, thus the tree always ends up with the same 0 th generation primogenitor for all objects, tracing generations in between to the last.
We compute orbital the semi-major axis, semi-minor axis and eccentricity from the relative position and velocity vectors in the Orbital Frame for each bounded pair i, j, which are computed similarly as those for the Cluster Frame in Eq. (25), and where again the position and velocity centres of mass are computed, time by time-slice, with where δ Zi,j ,B = 1 if the pair i, j is bounded and δ Zi,j ,B = 0 if the pair i, j is unbounded within a subsystem. From Eq. (36) the semi-major axis, semi-minor axis and eccentricity can be extracted by first computing the specific angular momentum of the pair with is the gravitational parameter. Then, the semi-major axis and orbital eccentricity can be computed with which, for sign(γ) < 0 yields an elliptical orbit, for sign(γ) > 0 yields a hyperbolic orbit, and γ → 0 constituting the parabolic limit where e = 1, a → ∞.

Merger trees computation
Also in each time-slice we compute the merger history of the bodies, tracking every occurring PBH absorption and in each merger event adding the identity of the least massive PBHs to the history of the more massive PBHs, thus providing a reconstruction of which particular objects have been merged to produce more massive ones. This is done by means of tracing every mass change in between every conservative pair of the N t − 1 = 64 last ESs, then solving N t − 1 = 64 constrained system of the diophantic equations, one per each pair of ESs, in to find out the identity of the pairs of objects that have merged and store the data where (∆M ) +,i is the mass gained by object i, (∆M ) −,j is the mass lost by object j, for all i = 1, 2, ..., N + objects that gain mass and for all j = 1, 2, ..., N − objects that loose mass in between snapshots, so N + ≤ N − . C is the mixing sparse matrix with elements C i,j = 0, 1 relating the absorber i with the absorbed j, of unit value in the case of a merger and zero in the opposite case. The mixing sparse matrix C can be very slow to solve in the case where a single PBH, or worse, many PBHs, have absorbed many other PBHs each in the same interval within snapshots. In that case, a number of matrix with elements C i,j could be found first before solving the full diophantic system of Eq. (43) and thus easing the problem if any of the easier merger cases or lack of thereof are identified first, mainly those of i) Impossible mergers: ii) Necessary mergers: iii) Single one-to-one mergers: iv) Multiple k -to-one mergers: After solving the aforementioned, one may proceed to solve the simplified diophantic system of Eq. (43).

IV. POPULATION STATISTICS
The time evolution of this positions and velocities is shown in Figure IV.7. Note that the cluster puffs-up to a factor of O(1000) to a final size of x C (1.38 × 10 10 yr) ≈ 100−1000 pc, decreasing its density by a factor of O 10 9 , while the ejecta expands further to distances in the range of x E (1.38 × 10 10 yr) ≈ 10 5 pc during the simulation period.
It is apparent in Figure IV.7 that the cluster timetrajectories exhibit a far wigglier pattern in comparison to the ejecta time-trajectories, as expected since the PBHs in the former ones move in quasi-random trajectories within the cluster while the latter ones very quickly arrive to the asymptotic regime where they move in straight lines.
The ejecta velocity dependence with time, however, is small in most cases, as it corresponds to asymptotically uniform rectilinear time-trajectories, the exception to this rule being the a few instances where wiggles are manifest, corresponding to PBHs expelled from the cluster in binary pairs, so that the velocity in the Cluster Frame frame is composed of two modes: a time-dependent orbital velocity over which the constant velocity of the binary barycentre is summed.
The Power Law (PL) fits of the time evolution of positions and velocities are given in Table IV.4. Note that no fit is made before the burn-in time for the ejecta positions and velocities of the simulations, given that the data variability is too large for this subset of objects, as very few of these exist in the time snapshots prior to the burn-in time. In particular, a threshold of 20 ejecta objects present collectively in all the N R = 5000 realisation considered has been required as the minimum over which the fits are computed both for the cluster and ejecta objects, which is only the case starting in the third time-run at the burn-in time t BI = 5.72 × 10 5 yr in Figure IV. 7 Now we proceed to identify the different populations arising in the simulations and show their evolution, which are the clustered vs. ejected PBH populations of Section IV A, the isolated vs. bounded PBH populations of     : zP(t ∈ r) = βi,z,Pt α i,z,P where z = x, v and r = 3, ..., 7 refers to the separate run where the fit is performed, βi,z,P corresponds to the amplitude of the PL while αi,z,P is the PL index, and P = C, E denotes the cluster and ejecta populations respectively.
Section IV B, and primitive vs merged PBH populations of Section IV C. The PBHs in our simulations can merge onto each other, so the total remaining PBH population is a monotonically decreasing function, while the absorbed PBH population is, in contrast, monotonically increasing over time. Once made this distinction, PBHs in the simulations can be, then, classified according to three different, independent criteria i) Depending on the eccentricity of the PBHs with respect to the cluster barycentre, one object belongs to the Cluster (C) for less than unity eccentricities: corresponding to instantaneously-closed orbits about the cluster barycentre, and it belongs to the Ejecta (E) for greater than unity eccentricities: corresponding to instantaneously-parabolic and hyperbolic orbits, where super-index 0 th denotes that it is the eccentricity with respect to the zero level, that is, the complete system, the one considering, and not with respect to the barycentre of a subsystem should the i PBH be part of one.
ii) Depending on the dimensionality of the parent tree, that is, on whether there is more than the unique 0 th PBHs common parent all objects share, the ob-ject may be classified as being part of the Bounded (B) population or the Isolated (I) population.
iii) Depending on the dimensionality of the merger tree, that is, on whether there is at least a 1 st PBHs merger event, the object may be classified as being part of the Merged (M) population or the Primitive (P) population.
There are a number of caveats to this classification algorithm that need further clarification.
First, and concerning the classification of PBHs into either cluster or ejecta objects, note that by "instantaneously" we mean that, at any given time, orbits can be assimilated with their respective Keplerian equivalents from their body to barycentre interaction, even though trough time no orbit will obviously remain a static Keplerian when successive interactions with other PBH takes place, and indeed, a PBH may flip from the cluster to ejecta population a number of times in short intervals, unless it quickly escapes the gravitational potential well of the cluster and joins indeed the ejecta population for the remainder of the simulations.
Second, and concerning the classification of PBHs into either bounded or isolated objects, note that any PBH part of a bounded system will be classified as such independently of whether the PBH is the more massive in the pair, in which case it will be labelled as the parent in the binary, or the least massive, in which case it will reference its parent in the parent tree. Also note that any bounded system save the cluster itself irrespective of the number of PBHs in it will be considered a binary system. However, in the unlikely event that the system is ternary or larger, this will be noted by subsequent parent trees following the hierarchy of masses in such multiple system, that is, by adding further levels in the parent trees.
Third, and concerning the classification of PBHs into either merged or primitive objects, note that, unlike in the other two cases there a PBH could flip-flop from one population to another and vice-versa, in this case the PBH can only transit either from the primitive to the merged or absorbed population it it is the more massive PBH in the merger or from the merged to the absorbed population it it is the less massive PBH in the merger, with no going back to the primitive population in any case, which is a decreasing function.
Note at last that all PBHs necessarily fall into one of the two options in the three separate categories, and therefore all PBHs must have exactly three tags: (C or E, B or I, M or P), Additionally, the tags R and A stand for the remaining and absorbed objects, and it is clear for all possible populations X = (C, E), (B, I), (M, P), where the populations in-between parenthesis are mutually exclusive.
We have quantified the statistical behaviour and evolution over time of all these populations by computing, snapshot-by-snapshot, their mean, median and 1σ and 2σ C.B., as well as fits to the PL expression of Eq. (31). Note two things, however, in relation with such statistical quantities i) First, that an outlier-exclusion algorithm has been applied to all data removed by more than 3σ from the best-fit of the calculated empirical distributions, which in effect removes at any given timesnapshot O(10) realisations from the computation of such statistical quantities, precisely those that exhibit very rare behaviour.
ii) Second, that a weak softening kernel has been applied to these quantities, in order to erase the shortlived, that is, single-snapshot, features that less rare but still atypical realisations have in the computation of such quantities, mainly in the form of irregularities of small amplitude and short wavelength that may render the corresponding timeevolved lines too wiggly.
iii) Third, that the kernel is weakest when applied to the total populations shown in the larger panels of Figure IV.8, since prior to softening the curves are already soft enough, and stronger when applied to the differential of populations shown in the smaller panels of same Figure IV.8, as those curves are, unsurprisingly, more wiggly.
The effects of both the outlier-exclusion algorithm and the softening kernel are, in fact, very minor, to the point of being barely visible in the case of the total population plots. They are, however, more relevant in the case of the differential population plots as they portray better the right levels of these statistical quantities, albeit at the expense of the small cyclic jumps visible in the 1σ and 2σ C.B. of Figure IV.8, a boundary effect due to the ten-fold increase in the time slices time-step ∆t r across different runs r = 0, 1, 2, ..., 7.

A. Cluster & ejecta populations
The simulated evolution of the cluster N C (t) and ejecta N E (t) population (see Figure IV.8a and Figure IV.8b respectively) shows, as it is expected, the gradual depopulation of the cluster and its loss of mass to the background, inter-cluster medium.
Note that the means, medians and confidence bands in said Figure IV.8 are cleaned with both the outlierexclusion and the softening kernel algorithms mentioned previously in Section IV, as otherwise they would show excessive variance in between time-slices due the very (a) Remaining cluster population of objects, N R,C (t).
(b) Remaining ejecta population of objects, N R,E (t).
Figure IV.8. Top: Total population of remaining cluster and ejecta objects in blue and teal respectively, NR,X (t) where X = C, E. The thick black line mark the mean of the discrete distribution of evolution trajectories per snapshot, while the thin black line in the same panel represents the median. Bottom: Differential of the total population of remaining cluster/ejecta objects, or loss/gain counts in between snapshots, rescaled by the simulation time, tδNR,X 2 (t) where X = C, E. The thick black line in the lower panels mark the mean differential counts, while the thin black line represents the median differential counts and the dotted black line corresponds the Bi coefficient of the logarithmic decrease/increase, or alternatively, the mean differential counts themselves. All: The darker and lighter regions represent the 68% and 95% C.B. of each individual evolution of the cluster and ejecta populations for all NR = 5000 realisations. The dotted vertical grey lines separate the Nr = 7 consecutive runs, and the dotted vertical grey line signals the burn-in time.
tBI = 2.64 × 10 5 yr Cluster population fits: Ejecta population fits:  31): NP(t ∈ r) = βi,Pt α i,P . αi,P is the PL index, and corresponds roughly to the constant value of the differential counts, while βi,P in turn, corresponds to the amplitude of the PL, r = 3, ..., 7 refers to the separate run where the fit is performed and P = (R, C), (R, E) denotes the remaining cluster and ejecta populations respectively.
Remaining cluster and ejecta population fraction and median:  Table IV.6. Mean remaining in-cluster and in-ejecta object population fractionsfi,R,P and medianfi,R,P where P = (C, E). The fraction is centred around the mean and shown are the 95% C.B. and is computed along with the median by averaging over the NR = 5000 realisations at selected times.
rare outlier O(0.001) realisation in which the cluster be-comes much more rapidly depopulated than usual, gen-erally because of the presence of an initially very massive PBH in the IC that slingshots many of its less massive companions, which is an artefact produced by the relatively large tails of the Log-Normal distribution that is used to sample the initial masses in the IC, and which will occasionally produce this very large PBHs with masses of O(0.001) M in a few of the N R = 5000 realisations. The population's evolution, showing an increase for the ejecta PBHs or decrease for the cluster objects has been fitted to the PL model of Eq. (31) and its results are shown in Table IV.5. The goodness of the fit is shown to be very good with less that O(0.001) variance unexplained by the PL model as shown by the goodness of the fit coefficient.
Note that the fits are done separately in each of the runs, and are computed only for the five lasts runs starting the burn-in time of the simulation t BI = 2.64 × 10 5 yr for the first of such runs. This is due to the fact that prior to that time the time-step of the simulation is too slow to capture the evolution within the time-runs as the characteristic time is during that period, as will be later computed and shown in Table VI.14 much greater than the time-step.
We provide as well the actual mean (95% C.L.) and median values of the fraction of cluster and ejecta objects within the simulations in Table IV.6, at the eight time-slices that bound the seven time-runs. Note that the cluster and the ejecta populations are almost, but not quite, complementary to each other. This is because some PBHs merge in the process, and is the case indeed here withṄ R (t) < 0, as will be seen in Section VI A 1. It will be seen later in Section IV C that the merging of PBH is so rare (thus the need for the high number of realisations to meaningfully capture it) that both fractions add up to unity.
In particular, we find in Figure IV.8 a high degree of cluster evaporation as roughly two thirds of the initial cluster objects are dispersed out of the cluster. The actual fraction and median number of PBH that remain in the cluster or join the ejecta population given in Table IV.6, show this trend indeed.
It will be later seen in Figure VI A 1 that the mass profiles of both the cluster and the ejected PBHs do not differ significantly throughout the simulations, meaning that both populations can be described at all times by the initial mass distribution, and so the mass fraction that is dispersed by the cluster is coinciding with that of the number of objects.

Cluster evaporation & virialisation
The phenomenon by which a group of gravitating bodies loses mass to the background, be it in stellar globular clusters or PBH clusters, is known as cluster "evaporation". In our particular case, starting at the burn-in time onwards , the evaporation rate is quite stable throughout the simulations. Our simulated PBH clusters are consis-tent with that fact, and indeed show up an evaporation rate roughly constant throughout the whole quasi-static phase, as it follows by the fact that the estimated confidence intervals of the β i coefficient of the PL expression of Eq. (31), shown in Table IV.5, are all quite close to each other, even if they quite not overlap due to the narrowness of the interval.
We find that, for the cluster, while the totality of PBHs are collectively bounded to the cluster at the initial times t < t BI,R = 4.92 × 10 5 yr, already in the third run the cluster starts to loose mass, doing so at a practically constant rate in proportion to the remaining mass, as it is shown in Figure IV.8a.
As for the ejecta, we find that the reverse is true: at the initial times t < t BI,E = 3.00 × 10 5 yr no object has acquired a speed high enough to escape the gravitational pull of the cluster in order to escape the cluster potential well. Only in the third run the background starts to be populated with ejected PBHs, doing so again at a practically constant rate in proportion to the remaining mass.
By the end of our simulations, the PBH cluster population has been reduced to fraction of objects Conversely, the PBH ejecta population has gained throughout the simulations a fraction of objects Note that both quantities being almost complementary to each other due to the lack of significant merging. Also note that, for the reasons shown in Section VI A 1, these fractions are good approximations to the fraction of mass, since the cluster and ejecta population don't develop differentiated mass profiles.
Cluster evaporation is apparent in Figure IV.8 starting from the burn-in time of t BI = 2.64 × 10 5 yr, a time at which the characteristic dynamical time of the cluster is comparable to the time-step in between the snapshots, and when the evolution is fast enough to be captured by the time-step. From that time onwards, then, the cluster proceeds to expel PBHs at a declining rate δN R,C ≈ −N R,E < 0, roughly given by The cluster and ejecta populations as a function of time are then well approximated by PLs, declining in the former case and growing in the latter case, whose fits are shown, time by time interval, in Table IV.5. It can be checked there that the fits are indeed very accurate with measures of the goodness of these with r 2 > 0.9 even right after the burn-in time.
Note that these PBH clusters are originated in the radiation era, and while they form this lasts, they do not, then, evaporate at all, their characteristic dynamical time being comparable to the total radiation domination era duration, and only right after cluster evaporation begins to be significant, therefore rendering this process irrelevant to cluster formation at least when the clusters are originated with parsec-sized characteristic length. Later evolution is then consistent with increasingly a majority of the PBH mass contained in diffuse haloes surrounding the proper PBH clusters.
It is shown in Figure IV.7a that ejecta haloes reach at present sizes of while cluster haloes reach at present sizes of Also shown in Figure IV.7a is the characteristic PBH ejecta halo typical velocities in the present to be in the rangev while the later cluster haloes reach typical velocities of v C,64 ∈ (0.030, 0.51) km/s.
This are indeed very interesting results, as one finds that ejecta haloes are comparable in size to typical galactic haloes, while cluster haloes reach sizes are comparable in size with dwarf galaxies and stellar globular cluster sizes, perhaps suggesting a link between the two.

B. Isolated & Bounded populations
A small portion of the PBHs in the simulations will end up in bounded systems, in which a larger object captures a smaller object and form a binary, as shown in Figure IV.9a and Figure IV.9b. These can be classified into two categories: transient, when the binaries last less that the simulation time, and stable, when they do last for more that the simulation time.
Note again that the means, medians and confidence bands in Figure IV.9 are cleaned with both the outlierexclusion and the softening kernel algorithms mentioned previously in Section IV in order to remove the noisier data, particularly in this case at the burn-in time crossing where the number of in-cluster binaries spikes to only to later be quickly reduced to very low levels as the IC is erased.
The population's evolution is fitted to PLs of Eq. (31) whose results are shown in Table IV.7. The fit is performed run by run, and for similar reasons as we had with the cluster and ejecta population, no fit is made before the burn-in time of the simulation t BI = 2.64×10 5 yr as the evolution is too slow to be captured by the timestep in the time-steps prior to it, or, in other words, the characteristic dynamical time the simulation prior to the burn-in time is, as will be later computed and shown in Table VI.14 much greater than the time-step at these time-slices. Figure IV.9a shows that while a respectably high fraction of PBHs appear in bounded systems initially in the cluster, the fraction quickly drops and stabilises at a low level of O(0.001) relative to the total number of PBHs.
As for the in-ejecta bounded PBHs, Figure IV.9b shows as it would be expected that their number consistently grows from almost none at all, as there is no ejecta population whatsoever in the initial stages of the simulations, to O(0.01) at present relative to the total number of PBHs. The actual fraction and median number of bounded PBH that remain in the cluster or join the ejecta population are given in Table IV.8.
Initially, the vast majority of the binaries are bound to be contained within the cluster, which contains a fraction of bounded objects of as there has not yet elapsed enough time to evaporate a significant fraction of objects to the ejecta, which contains a fraction of bounded objects of: both at the initial simulation time of t 0 = 0 yr, and with respect to the total population. Note that the large fraction of bounded cluster objects contained is an artefact of the IC, and quickly disappears after the burn-in time, which is natural considering that around this time the cluster undergoes a transformation, moving from the initial Multi-Normal position distribution of Figure III.1 to a cuspier one as will later be seen in VI.12a.
As a result from this, the density is increased at the PBH cluster core at the same time that the cluster slightly contracts about 5%-10% as seen in Figure IV.7a. As the PBHs fall to this dense environment they quickly acquire speeds as seen in Figure IV.7b that are able to disrupt the majority of bounded systems and transition from a Gaussian to a PL density profile.
After the transition and while in the cluster region, such binaries are typically very short lived, as they interact with neighbours multiple times in each run and so they are still disrupted easily. As a consequence of this, the population of cluster binaries is small overall and remains roughly constant, as the generation rate of binaries is matched by their disruption rate at any timeslice t > t BI .
The fraction of cluster binaries with respect to the total population varies then from an absolute minimum with respect to the total population of at t = 1.38 × 10 6 yr to a present relative maximum with  The thick black line mark to the mean of the discrete distribution of evolution trajectories per snapshot, while the thin black line represents the median. Bottom: Differential of the total population of bounded in-cluster/in-ejecta objects, or loss/gain counts in between snapshots, rescaled by the simulation time, tδNB,X (t) where X = C, E. The thick black line marks the mean differential counts, while the thin black line represents the median differential counts. The dotted black line corresponds the Bi coefficient of the logarithmic decrease/increase, or alternatively, the mean differential counts themselves. All: The darker and lighter regions represent the 68% and 95% C.B. of each individual evolution of the cluster and ejecta populations for all NR = 5000 realisations. The dotted vertical grey lines separate the Nr = 7 consecutive runs, and the dotted vertical grey line signals the burn-in time. The population's evolution is fitted to PLs in Table IV.7. tBI = 2.64 × 10 5 yr Bounded in-cluster population fits: Bounded in-ejecta population fits: at t = 1.38 × 10 10 yr, as seen in Figure IV.9a. Note that, as the cluster itself looses population, while the cluster bounded population remains nearly constant, then the fraction of PBHs contained in binaries within the cluster grows at a rate inverse to the evaporation rate, more that doubling from the burn-in time to present. However, when such a binary is ejected from the cluster then it very quickly ceases to be perturbed by nearest neighbour encounters and then, in the code classic treatment, its constituents remain in stable Keplerian orbits indefinitely, as seen in Figure IV.9b.
It follows from this fact that the population of binaries in the ejected background naturally increases overtime, eventually overcoming that of the population cluster binaries already at t ≈ O 10 8 yr and constituting at late times a non negligible fraction of all ejecta objects, which in our simulations is found to be at t = 1.38 × 10 6 yr to a present relative maximum of: at t = 1.38 × 10 10 yr, as seen in Figure IV.9b. Last, a more detailed analysis of binary pairs is given in Section V A 1, where we extract the mass profile of bounded PBHs and their population characteristics.

Simulation parent trees
While it is true that the majority of PBHs do not constitute a part of a binary system at any given time, a fraction of PBHs that form part of such systems has been found to lie in the interval after recombination at t = 1.38 × 10 6 yr and at the present time t = 1.38 × 10 10 yr, with f B = f B,C + f B,E , as previously explained in Section IV B. However, we have found as well a rich structure of subsequent binary sub-systems within multiple PBH systems. It has been found using the methods from Section III C 6 that a hierarchical substructure of binaries arises in the simulation IC and is maintained throughout the whole evolution period, the results of which can be seen in Table IV.9.
Generally speaking, we find that a majority of PBHs do not participate in any binary system. However, of those which take part in such systems, it has been found that i) A 99%-90% majority of those belong to binary (2component) systems in which two PBHs, one larger and one smaller, orbit each other, which we label on the binary hierarchy 0 th level and 1 st level respectively.
ii) Another 10%-1% of bounded objects are present in tertiary (3-component)  Thus, it is found that each subsequent level in the hierarchy is populated by roughly a hundredth the number of PBHs than the level immediately above. Note as well that, as shown in Table IV.9, the population in each level is less stable over time the deeper the level in the hierarchy.

C. Primitive & merged populations
Last, an even smaller portion of the PBHs in the simulations will, by the end of the simulation period, have merged with another and form larger PBHs. We have ii) (N R,C,B ) M = 4 such mergers occur in cluster, bounded objects leading to a total sub-fraction of merged objects of: iii) (N R,E,I ) M = 1 such mergers occur in ejecta, isolated objects leading to a total sub-fraction of merged objects of: iv) Last, (N R,E,B ) M = 62 such mergers occur in ejecta, bounded objects leading to a total sub-fraction of merged objects of:  Fraction of binaries at each of the eight time-slices dividing the seven runs. Each level lB indicates the level in the binary hierarchy, with lB = 2 indicating two bound PBHs in a binary system, lB = 3 indicating a binary system within another binary system, and lB = 4 indicating a binary system within another two binary systems. Note that there are Nr = 7 time-runs, for a total of eight bounding time-slices, corresponding to the times shown in Table III.2.
These results, summarised in Table IV.10, show that the large majority of merged PBHs belong to two populations of objects, mainly isolated, in-cluster PBHs that constitute about 43% of mergers and bounded, in-ejecta PBHs which constitute about 53% of mergers.
Indeed, the latter case of bounded, in-ejecta PBHs can be attributed to the fact that the PBHs that partake in the mergers are almost exclusively counted among the most massive in the simulations, as shown in Table V.12, and as the PBHs mutually approach each other within the cluster, they acquire large infall velocities as they fall in their absorber's potential well, that can overcome the cluster escape velocity and be expelled to the background. The former case of isolated, in-cluster PBHs, however, has a simpler explanation, and is merely due to the fact that cluster objects, by definition, live in a denser environment where collision is far more likely.
Even though the numbers are small, of the remaining cases, bounded, in-cluster PBHs merely constitute about 4% of PBH mergers, while isolated, in-ejecta objects which constitute a very minor 1% of PBH mergers. The cause of this can be attributed to the environment density where the mergers do take place.
In the former case, the merger takes place between bounded objects within the cluster sphere-of-influence, and their low number count can be attributed to the fact that the bounded in-cluster population is already quite small at any given time, even though it can be observed that events of this kind are enhanced with respect to the mergers arising from isolated in-cluster objects, as shown by the quotient of merger counts being larger than the population quotient at the last time-slice at present However, in the latter case, the merger takes place just within the periphery of the cluster with massive, high-velocity PBH, just evaporated and fast enough to overcome the cluster escape velocity so it is considered as part of the ejecta population even if its least massive merger pair still is part of the cluster, and so it cannot be consider a proper ejecta merger, consistent with the fact that, given the radial profile of ejecta velocities and the very low ejecta density, the chances of such mergers from objects incoming from the same single cluster are vanishingly small.
Merger tree hierarchy levels:  Be reminded that the number of PBH mergers that we find in our simulated cluster underestimates the amount of actual mergers that would be produced due to the non-relativistic behaviour of the code. Given the typically low object velocities shown in Figure IV.7b, it would expected first that many of the identified hyperbolic encounters and binary captures in our simulations would produce a binary capture or inspiral and merger in a relativistic code, and as such, the numbers provided in this section must be understood as lower bound to the actual merger amount.
A more detailed analysis of these merger events will given in Section V B 1, where we compute the mass profile of merging PBHs and their population characteristics.

Simulation merger trees
The identification of merger pairs in order to construct the merger trees has been done by the procedure described in Section III C 7, looking at mass differences of objects in between snapshots. Typically, the mergers happen at the last simulation time-runs, a time when average velocities are the smallest for all PBH populations at the core of the cluster, and where nearest neighbour interactions are most frequent, often arising from disruption of a transient binary.
Merger identification then straightforward, at least for most cases in our simulations, as there are not more that two mergers in between time snapshots in any realisation, so that, out of the (N R ) M = 117 merger events identified in our simulations, we find that i) (N R ) 1,1 M = 115 events are one-to-one mergers in which a single PBH absorbs another less massive PBH in between two particular time-slices belonging to the last simulation time-run.
ii) (N R ) 2,1 M = 2 events are two-to-one mergers, in which a single PBH subsequently absorbs two PBHs in between particular time-slices that may or may not be contiguous, within the last time-run.
Last, it should be noted that mergers are so rare in each of the simulations that, out of the methods devised in Section III C 7, only the first three of impossible, necessary and single-to-one mergers are relevant, as there has been no case found in all of the N R = 5000 realisations where, in between two consecutive time-slices, may have been many-to-one mergers.

V. CLOSE ENCOUNTERS
In this Section we take a closer look to the identified binary and merger pairs, in Sections V A and V B respectively. In particular, we study the collisional and shortdistance interaction dynamics of PBHs in these pairs, and extract the mass distributions of the PBHs participating in such systems. Also, we offer a lower bound prediction of the rates of occurrence of PBH mergers and of the PBH hyperbolic encounters per cluster. In particular, we find that, integrated across all the time-slices, binary pairs constitute about O(0.001) of all bodies in the simulations. About 85% of these remain in the cluster while the remaining 15% of these have been ejected in two-to-one encounters.
Note that this ratio is heavily influenced both by the IC, as the abundance of bounded in-cluster PBHs at early times is overwhelmingly larger than that of bounded inejecta PBHs, as well as the later evolution, as we had found in Section IV A that as the simulations approach the present time, the fraction of in-cluster PBHs remains roughly constant despite the gradual evaporation and depopulation of the cluster, while the corresponding number of in-ejecta PBHs grows monotonically, even if at an ever decreasing rate as shown in the PL-fits, also as a consequence of cluster evaporation.
If however, we restrict ourselves to time-slices posterior to the burn-in time, so that the IC has been already erased, then the numbers at present change significantly. About 65% of bounded PBHs do remain in the cluster while the remaining 35% of bounded PBHs will have been have been ejected in two-to-one encounters, as shown Section IV A.
In any case, the masses of the bodies participating in binary systems are naturally found to be substantially larger than those of the mean or median PBH in the simulations. As seen in Section III B 1, the PBH average and mean masses in the IC arem = 21.9 M and m = 7.55 M . Given both the lack of any meaningful number of mergers and that the mass profiles of the cluster and ejecta populations do not evolve differently throughout the simulations, the mean and median masses remain nearly constant for both the in-cluster and inejecta bounded populations up until the present time.
However, for binary pairs, the mean and median total masses in the pair arem T = 303 M andm T = 181 M , which is far more than twice the simulation mean or modal mass. This indicates that only very massive PBHs with sufficiently deep and far reaching gravitational potential wells are able to overcome the large, O(1) pc − O(1000) pc-sized and monotonically growing distances between PBHs in between the burn-in time and the present time as the cluster puffs up and become a parent in the binary pair.
Moreover, the binary pair mean and median mass ratios ofq T = 0.154 andq T = 0.0552 indicate that the less massive partner in the binary is typically between 6 and 20 times less massive than its parent, establishing a clear mass hierarchy between both, as more equally distribute binary pairs may be easier to disrupt in the long term.
In addition to this, and taking into account both the mean and median total mass and mass ratio, we find that just like its parent PBH, the secondary PBH in the pair is also typically larger than the typical PBH in the simulations, although to a lesser extent determined by the mass ratio. In particular, given that m 1 + m 2 = m T , m 2 /m 1 = q, and the LIGO convention m 1 ≥ m 2 , two PBHs in a binary pair, according to the previous results would have masses close to when computed from the mean or the median values of the total mass and mass ratio respectively, showing masses for the smaller partner that are bigger than the simulations' mean and median mass in both cases. This is maintained in the in-cluster and in-ejecta bounded populations of PBHs as well, although the degree to which the masses are larger than the typical mass or not varies a by a small amount. In particular, we find i) For the bounded, in-cluster PBHs,   Figure V.10, representing the binary profiles segregated by the population (cluster and ejecta) type integrated across all time-slices in the simulations starting at the burn-in time tBI = 2.64 × 10 5 yr up to the present, for a total of N R i=1 (NB)R = 2.32 × 10 6 binaries throughout all realisations. Note that we use our convention that C, E, P, R and B stand for the cluster, ejecta, primitive, remaining and bounded populations respectively. Note as well that we do not show the corresponding mean and median values of merged PBHs in bounded systems in this table since their number is negligible compared to the other cases, but we do show them however on Table V. 12. ii) For the bounded in-ejecta PBH, Overall, this indicates that the more one-sided distributions of mass in a binary pair as well as the more massive binary pairs are biased towards the ejecta population, as expected.

Hyperbolic encounter rates
We define a hyperbolic encounter between two bodies in the simulations as one in which the eccentricity of the approaching crosses the thresholds e = 1 in any direction in between snapshots t and t + ∆t or at any point during the close encounter, irrespective of the impact parameter b of such encounter or other orbital characteristics, and due to a close one-to-one body interaction. However, there are two issues with this approach that makes the interpretation of this quantity as the close encounter, GW-generating event a bit problematic.
First, that the number of hyperbolic encounters is by these means underestimated as one particular PBH with an orbital eccentricity crossing the parabolic threshold two times in opposite directions in between consecutive time slices will not register as a hyperbolic encounter at all despite that there have been, in fact two total encounters in such time interval.
Second, that such hyperbolic encounters most often happen at distances large enough that the amplitude of the emitted GWs will be too small to be registered on Earth detectors within the foreseeable future, meaning that this rates are not to be confused with the rate of detection of a gravitational wave event.
We can thus only provide a minimum threshold of the number of this encounters, given that an irreducible number of these hyperbolic encounters will be missing in any case from the fact that the time resolution is limited. This effect is not, in fact, very large, since, as seen in Table VI.14, the time-step in each runs is adapted to capture most of such interactions, being of an order of magnitude as the cluster crossing time, but being of the same order of magnitude the effect cannot be neglected either.
We now proceed to compute the minimum hyperbolic encounter rate per time-run and realisation as where N S r is the cumulative number of slingshots identified in the simulations with the aforementioned procedure at the end of time-run r and ∆t r = t i + −t i − the time-run total period, leading to a minimum hyperbolic encounter rate of: Γ S 4 (29 ≤ t i ≤ 37) = (0.9 +2.9 −0.6 ) × 10 4 Gyr −1 , Note that these slingshot rates have been extracted for non-relativistic simulations, adding up to the underestimation the underestimation of the rates as the emission of gravitational radiation results over time in a more compact clusters of PBHs, were one-to-one PBH encounters are more frequent.
However, given the large scale of the cluster, of O(1) pc already at the burn-in time, the effect is nearly negligible. Also, note that these rates have been computed for a single cluster. Latter, in Section VIII we will compute this quantity for a comoving cosmological volume of 1 Gpc 3 , and thus provide a more meaningful value for the slingshot rate of PBHs.

B. Merger pairs
Similarly to what we have done for the binary pairs of Section V A, we have computed the merger pair number and frequency in the simulations along with the total    Figure V.11, representing the merger profiles segregated by the population (cluster and ejecta) type integrated across all time-slices in the simulations starting at the time of first-merger (1M) tFM = 1.38 × 10 9 yr up to the present, for a total of N R i=1 (NM)R = 117 merger events throughout all realisations. Note that, as in the rest of the paper, we use our convention that C, E, I, B, R M stand for the cluster, ejecta, isolated, bounded, remaining and merged populations respectively.
In particular, we find that, integrated across all the time-slices, merged pairs constitute O 10 −5 of all bodies in the simulations, though note that, in contrast to what was the case with binary pairs, mergers occur late in the simulations, with most of them occurring in the last time-run, and none occurring during the first five runs.
Of these, about 46% remain in the cluster while the remaining 54% have been ejected in two-to-one encounters, occasionally because of the merger itself, as we had found in Section IV C. Unlike what happened with binary pairs, however, this ratio is not influenced at all by the IC as there are no mergers during that stage, but it may still be influenced by the later evolution as the cluster evaporates. Regrettably, with only 2.6% of mergers in the sixth run and the remaining 97.4% of mergers in the seventh run, we lack data to quantify if later evolution does indeed alter the balance between the in-cluster and in-ejecta merged population.
In any case, and for similar reasons as what was the case with binary pairs, the masses of the bodies participating in mergers are naturally found to be substantially larger than those of the mean or median PBH in the simulations, only even more so. We had extracted in Section III B 1 the PBH average and mean masses in the IC arem = 21.9 M andm = 7.55 M and found those to be practically invariant for the cluster and ejecta populations throughout time-slices up to the present.
For merged pairs, the mean and median total masses in the pair are indeed very large, far more so than for binary pairs, withm T = 1.67×10 3 M andm T = 1.51×10 3 M , which is 55−200 times more than the simulation mean or modal mass respectively. This is because only the truly massive PBHs, and most often the most massive PBH in the simulations, are sufficiently dominant over the other less massive PBHs that they may be able to overcome the order kpc-sized distances between PBHs in the later stages of evolution and merge still with a small number of them.
This typically happens for a merger pair of mean and median mass ratios ofq T = 0.965 andq T = 0.545 , indicating that the less massive partner in the merger is typically bigger than half the mass of the more massive part-ner, with no pronounced mass hierarchy between both, as it would be expected from the fact that in this classical computation the merger cross section is quadratic with the Schwarzschild radius and therefore quadratic with the mass, so the merger of two large PBHs is dramatically more likely to occur than the merger of a more uneven pair, especially so as in our non-relativistic simulations the bodies neither inspiral nor emit GWs.
Therefore, and taking into account the mean and median total mass and mass ratio, we find that, given that m 1 + m 2 = m T , m 2 /m 1 = q, and the LIGO convention when computed from the mean or the median values of the total mass and mass ratio respectively, showing masses for the smaller partner that are bigger than the simulation mean and median mass in both cases. We find as well that this feature is upheld in the incluster and in-ejecta bounded populations of PBHs as well, although the degree to which the masses are larger than the typical mass or not varies a by a small amount. In particular, we find that i) For the merged, isolated and in-cluster PBHs, This shows again that less one-sided distributions of mass than in the binary pairs as well as the fact that more massive binary pairs are biased towards the cluster population, which is to be expected from dynamical friction, significant in this mass range as will be shown in Section VI B. Also, we do not compute the cases of the merged, bounded, in-cluster PBHs and merged, isolated, in-ejecta PBHs as with only 4 and 1 cases respectively we are far from having enough sampling to give meaningful results.

Merger event rates
We have defined a merger between two bodies in the simulations in Section IV C 1, and Tables IV.10 and V. 12 show the frequency and population distribution of such events.
Note that while we on the one hand underestimate the number of PBH mergers in the simulations because of the fact that the N -body solver computes the evolution of the system in a non-relativistic manner, and as such there are no inspiral events that enhance the number of mergers from close binary systems, on the other hand we may be overestimating merging at later times since a rich merger history at early times will concentrate more mass at the core of the cluster and depopulate faster its outer layers, resulting in a diminished likelihood for merging at later times.
In this case, and considering the 117 merger events identified in total in the N R = 5000 realisations during the last two runs, we can roughly estimate the merger event rate per time-run and realisation with where P = (C, E), (I, B) is the population of interest and ∆t r = t i + − t i − the time-run total period. Then, replacing with the results of Table IV.10 we find that, We can also compute these results segregating by their origin, which can be of great importance when identifying one such merger as mergers arising from a hyperbolic encounter between two isolated PBHs or an inspiral encounter between two bounded PBHs are easy to discern as they leave very different signals in detectors. Then i) For all 51 mergers occurring from previously isolated PBHs, out of which 1 occurs in the sixth time-run and 50 do occur in the seventh time-run, we find that, Note that, in this last step we have taken into account the fact that most PBH mergers that are accounted in the ejecta do in fact occur within the cluster, only to be soon expelled to the ejecta as the simulations progress. This typically happens because of the large velocity kicks associated with the merger process, and so we have segregated the merged pair by whether they originate from single PBHs or bounded PBHs.
Also, these are the merger event rates for a nonrelativistic evolution for a single cluster. In Section VIII we will compute this quantity for a comoving cosmological volume of 1 Gpc 3 , and thus offer a more meaningful value for the merger event rate of PBHs, and one that can be interpreted as a minimum for the actual GR evolution.

VI. DYNAMICAL DISTRIBUTIONS
Now we proceed to extract in Section VI A the density, mass, position and velocity distributions of PBHs evolution, slice by slice in time and segregated by the sub-population type. We do as well quantify the cluster evaporation rate in Section VI B 1 and the degree of both mass segregation and dynamical friction within it in Section VI B 2.

A. Dynamical profiles
The PBH density profile, position profile and velocity profiles are shown in Figure VI.12 segregated by cluster and ejecta populations for selected times binding the runs, while the fits to the data are shown in Table VI.13.

Mass profile evolution
The evolved mass distribution of PBHs, ρ m,i | t , for i = C, E the cluster and ejecta populations respectively, has not been found to do develop differentiated characteristics or significant evolution throughout the whole simulation period and so no Figure is given showing it. Shown, however, in Table VI.13 are both the cluster and ejecta mass distribution parameters, which remain purely LN, and virtually indistinguishable from the masses IC (µ m , σ m ) = (2.0, 1.5) M . This is a somewhat surprising result, as it would be naively expected from dynamical arguments that the ejected object's mass profile would noticeably skew towards the lower end of the initial mass range, while the cluster object mass profile would in turn skew towards the higher end of the mass range, given the large mass differences in between PBHs, spanning O(4) orders of magnitude as seen in Figure III The underlying reasoning behind mass segregation is as follows. As more massive PBHs do concentrate in the inner core of the cluster, as shown in Figure VI  acquire in the process large kinetic energies, as shown in Figure VI.13b. Then, due to the large number of encoun-ters with massive PBHs in the core, the least massive objects are in turn expelled with increased likelihood from the core. In the end this results in an emerging massdepleted population of PBHs that would be segregated from the mass-enhanced population of PBHs in the cluster, and particularly, from the cluster core. Indeed, this is in a very restricted sense what happens. Despite the fact that the ejecta population comprises about 66% of PBHs at the last time-slice in the simulations, cases where the largest PBH ends up in the ejecta background are very rare, exceedingly so with respect to what would be expected if no correlation between mass and the likelihood of a PBH being slingshot away to the background is assumed.
In practice, this means that phenomena such as mass segregation and dynamical friction will be significant for about the top O(10) most massive objects in the simulations, which are well into the tail of the mass distribution. For the remaining O(1000) of objects, however, there is little evidence of significant mass segregation and dynamical friction, either in the cluster or in the ejecta populations. The cluster and ejecta-specific mass distributions do in fact skew towards larger and smaller masses respectively during the simulated period of time, but only so slightly, with their respective Log-Normal distribution parameters µ m C and σ m C increasing only by  The PL fits of the density profiles are shown in Table VI.13, once removed the outer 20% of the data points at which the cluster and ejecta densities suddenly fall off as there is less and less sampling of the density distribution as one transitions from the cluster sphere-ofinfluence to the ejecta sphere-of-influence in the first case, and from the ejecta sphere-of-influence to the void in the second place.
Note that we have tried pseudo-isothermal, Navarro-Frenk-White, and Einasto profiles for the fits, but it has been found that the both the cluster and ejecta density profiles were best described by simple PLs, with the PL index as a free parameter. Two different observations are noteworthy.
The first is that, for time-slices posterior to the reconfiguration of bodies that follows the burn-in time t BI when the IC is erased, the cluster volume-factored density profile exhibits a plateau up to the time-dependent cut-off radius x max C (t) where the cluster is sparse enough, that effectively can be considered as the cluster radius.
The second concerns the ejecta objects. The volumefactored distribution density profile of these, unlike in the previous case, does not exhibit any plateau, but rather a peak at a time-increasing length scale corresponding to the median distance travelled from the cluster barycenter. Before the peak, the actual ejecta density of objects gradually increases achieving a maximum at the cluster boundary. After the peak, the density PBHs falls abruptly, sampling only the more infrequent highvelocity PBHs. Last, the density profile arrives as well at each time-slice to a time-increasing cut-off length scale x max E (t) that corresponds to the maximum reach of PBHs expelled from the cluster.
Taking into account these two then we find that: i) The proper cluster density profile up to the cut-off distance roughly behaves as as inverse nearly-cubic PL: ii) The actual ejecta density profile up to the cut-off distance roughly behaves again as an inverse PL, with a slightly smaller power than in the cluster case: Note that the volume-factor x 3 that multiplies the actual density profiles ρ(x)| t masks the large range covered by the density within the cluster and ejecta spheres of influence, as distances extend to x max C → 10 4 pc and x max E → 10 4 pc respectively. In practice, we find in our simulations, after averaging all realisations, that densities are in the range of i) In the case of the initial time-slice of cluster objects: x ∈ (0.001, 10) pc, (112) t 0 = 0 yr.
iii) In the case of the first sufficiently populated timeslice of ejecta objects, with N E (t) > 100 : x ∈ (10, 10 7 ) pc, (121) It is worth noting, however, that cluster layers with densities smaller than ρ = 0.08 M pc −3 (see Ref. [80]) will already be disrupted at solar neighbourhood-like environments, and stripped of PBHs. Therefore, considering the more realistic case of clusters not in isolation but on top of a background mass density and in interaction with other clusters will put more stringent bounds both to the cluster and ejecta radial reach, coinciding with the local background matter density at their outer layers.

Position profile evolution
The PBH evolved position distributions, ρ x,i | t , with i = C, E being the cluster and ejecta populations, are given in Figure VI.12c and Figure VI.12b respectively at the five different time-slices that bound the last four time-runs and the IC. Log-Normal fits to these position profiles are shown in Table VI.13, when such fit is possible. There, a number of features become apparent.
The first is that right after the ICs, the cluster PBHs compress around 1.38 × 10 6 yr, only to later rebound and subsequently puffs up. At this time the PBHs do, on average, falling closest than they will ever do during the whole simulation period afterwards towards the cluster barycenter, approaching it to distances up to O(0.1) pc. This short compression-expansion phase is due to the choice of Maxwell-Boltzmann ICs in the process described in the previous section, which do not survive after this time.
Secondly, that as the cluster puffs-up, its characteristic length scale increases greatly as described in Eq. (164), with its radius expanding by three orders of magnitude, as seen for the last four time-slices in Figure VI.12c, after a brief phase of slight contraction on the second time-run after the IC. This transition happens right around the burn-in time t BI = 2.64×10 5 yr, and marks the transition from the Maxwell-Boltzmann IC that best describes the IC to the Log-Normal time slices later on, in the last four time-runs. This late-time cluster positions have been fitted to a variety of functions, including Maxwell-Boltzmann and Log-Normal distributions, finding out that it is the latter that best describes the data, in particular to accommodate the large tails in the later position distribution that develops as the cluster grows a subpopulation of far-off components with very elongated orbits with respect to the cluster barycenter, with periods comparable to the entire simulation time δt C | x x max C ≈ t 64 = 1.38 × 10 10 yr and nearly unit eccentricities e C | x x max C → 1 − . Cluster positions are constrained to be in the range x C (1.38 × 10 10 yr) ∈ (10, 2 × 10 3 ) pc, while the width of the constraint at a particular time is found to always be roughly of two orders of magnitude between the minimum and maximum values that enclose 99% of the variable's distribution. The cut-off distance, akin to the cluster radius, roughly behaves as Thirdly, that the stability of these PBHs in orbits bound to the cluster is almost guaranteed in our simulations, since the cluster lives in isolation, and the likelihood of an encounter that may transfer enough kinetic energy to overcome the cluster escape velocity is very low, as the ejected objects travel the cluster periphery where this far-off objects live very fast.
However, it will be seen in Section VIII A 1 that at late time-slices, in the very last run, the stability of cluster objects and particularly those on the periphery, far off the cluster barycenter, is far from assured in realistic scenarios where the cluster does not live in isolation, and inter-cluster interactions become likely.
A last interesting feature is that the ejecta reach, or alternatively, the ejecta sphere-of-influence, grows linearly in time, as it is natural, since the ejecta sphereof-influence is bounded by the fastest bodies that escape the cluster with isotropic and asymptotically constant velocities. As in the cluster case, the characteristic length scale increases linearly in time in the manner described in Eq. (174).
Ejecta positions are constrained to be in the range x E (1.38 × 10 6 yr) ∈ (0.8, 10) pc, x E (1.38 × 10 10 yr) ∈ (2, 100) × 10 3 pc, while the width of the previous constraint at any given time in the quasi-static phase is found to always be of one order of magnitude between the minimum and maximum bound that enclose 99% of the individual positions. The cut-off distance, which can be understood as the reach of the ejecta sphere-of-influence, scales as Late time ejecta positions have been as well been fitted to a large collection of distributions, including the MB, and Multi-Normal distributions.
However, and unlike the cluster case, no good single parametric fit has been found to accurately describe the at all post burn-in times simultaneously. The reason why the fits fail, particularly at late times, is that as PBHs are constantly emitted from the cluster at all times, the low end of the ejecta position distribution is constantly repopulated, which does not behave well with the exponentially-suppressed Log-Normal distribution in the logarithmic variable at that scale, while the long tails of the ejecta position distribution due to high speed objects are ill-described by the exponential suppression on the variable's higher end of the Maxwell-Boltzmann distribution, itself only good to describe the IC and the immediately following time slices.
Note that these ejecta PBHs remain in rectilinear asymptotically-uniform motion trajectories, gradually approaching their end-velocity in our simulations at infinity, and no longer taking part in the dynamics of the simulations, as seen in Figure IV.7b. This greatly increases the computational speed of the simulations at late times, since the characteristic discrete dynamical time-step in the N -body solver algorithm increases accordingly.
However, it will be seen in Section VIII A 2 that again in the last runs, the stability of such ejecta objects is no longer assured in realistic scenarios where the cluster does not live in isolation, as interactions between the ejecta objects and other clusters become likely at even earlier times.

Velocity profile evolution
The PBH evolved velocity distribution, ρ v,i | t , for i = C, E being the cluster and ejecta populations, is given in Figure VI.12e and Figure VI.12f respectively at the five different time-slices that bound the last four time-runs and the IC. Log-Normal fits to these velocity profiles are shown in Table VI.13, when such fits are appropriate, although, like it was the case with the position profiles computation of Section VI A 3, other fits to Maxwell-Boltzmann and Rayleigh distributions have been tried and proved to worse describe the data.
From Figure VI.12e and VI.12f it is apparent that the range covered by velocities for both the cluster and ejecta PBHs is much narrower than it was for masses or positions, as it would be expected from equipartition of energy arguments. In particular, cluster velocities are constrained to a narrow range of and the width of the constraint at any given time is estimated to be approximately of one order of magnitude between the minimum and maximum bound that enclose 99% of the individual velocities. A feature in Figure VI.12e stands out: right after the simulations starts, the cluster PBHs accelerate and reach the largest velocities they achieve on average during the whole simulated period of even O(10) km/s as the cluster contracts and then later rebounds due to the inward fall process described in Section VI A 3. After this phase, nonetheless, cluster velocities continue to decrease indefinitely as PBHs continue to expand and be diluted.
The time-dependent characteristic velocity scale is then Ejecta velocities are even more constrained than cluster velocities to a range of while the constraint width between the minimum and maximum bounds that enclose 99% of the distribution velocities is estimated to be approximately of half an order of magnitude at any given time after the burn-in period ends. The characteristic velocity scale is even less timedependent, and is given by The reason for this lack of time dependence in the ejecta velocity distribution is due to the balance between two competition effects. On the one hand, the dynamical relaxation of the cluster as it expands produces an ever increasingly slower motion of PBHs within the cluster. On the other hand, the cluster escape velocity, which is well approximated by given the symmetries of out clusters, and where X * C (t) is the time-dependent cluster radius and is the total mass contained in the sphere bounded by it, decreases over time. In the end, the difference between the decreasing kinetic energy kick to escape the cluster and the also decreasing typical cluster PBH kinetic energies is roughly constant and naturally produces in turn roughly constant ejecta PBH velocities at infinity.

B. Mass segregation & dynamical friction
In this section we comment on the degree to which mass segregation and dynamical friction are observed within the cluster and characterise both phenomena.

Cluster mass segregation
Mass segregation is the dynamical process by which the heavier bodies in a gravitationally bound system tend to fall to the central area of the system while the lighter bodies move farther away from the centre. It can be observed at a large range of scales is astrophysical systems, from star clusters to galaxy clusters and results in the negative cross-correlation of the mass and position distribution with heavier masses skewing towards smaller radial distances from the cluster barycentre.
The particular way this phenomenon works is as follows. During a close encounter of two bodies such as PBHs, the bodies exchange both kinetic energy and momentum. After a large number of encounters there is a tendency towards the bodies having similar energies from equipartition of energy arguments. Since energy is quadratic in velocity and linear in mass, the most massive objects move at slower speeds than the least massive objects in order to equipartition of energy be sustained.
This effect works for any range of body masses and radial positions, but is stronger the denser the medium and the more compressed the bodies' motions are, as shown by the time taken by the cluster bodies to achieve equipartition, called the relaxation time, which has been shown in Ref. [82] to approximately be where N R,O (t) is our time-dependent number of cluster PBHs, and t cross (t) is the time taken by PBHs to cross the cluster. In our simulations, Eq. (149) yields for the IC a relaxation time of t rel C,0 = (1.323 ± 0.071) × 10 7 yr, since at t 0 = 0 yr we have N R,C (t 0 ) = N O = 1000 ; and X C (t 0 ) = 1.710 ± 0.060 pc and V C (t 0 ) = 2.301 ± 0.091 pc from Table VIII.19 yields an initial crossing time of which is of the same order of magnitude than the burn-in time of t BI = 2.64 × 10 5 yr, an issue of great importance since, as Table VI. 15 shows, shortly after the burn-in time, the cluster starts to expand, which forces the crossing time to increase dramatically as cluster PBHs bulk speeds do not catch up with the increase in the characteristic length scale We illustrate cluster mass segregation in our simulations by extracting the ρ(m C , x C , t) distribution at the time-slices that divide the simulation time-runs. We show our results in Figure VI.13a and give the best-fit values in Table VI.16. There, a number of features are apparent i) First, that the phase space (m C , x C )| t occupied by the cluster PBHs undergoes an expansion of the position profile from a narrow Maxwell-Boltzmann distribution in the IC to a broader nearly Log-Normal distribution after the burn-in time at t BI = 2.64 × 10 5 yr, consistently with other observations.
ii) Second, that around the burn-in time, again the aforementioned phase in which the cluster contracts is apparent, in that the best-fit value of the distribution at t 28 = 1.34 × 10 6 yr has moved towards the lower end of the position range, the only time to do so throughout the simulations, again consistently with other observations throughout the simulations.
iii) Third, that for later times t > t 37 = 1.34 × 10 7 yr the cluster expands as the distribution's best-fit is displaced towards the higher end of the position range, while not changing significantly its mass value, in agreement with the finding in Section VI A 1 that neither the cluster nor the ejecta distributions develop a mass profile that is substantially different from the IC throughout the entire simulation period. iv) Last, that mass segregation is clearly apparent in Figure VI.13a in that both the 1σ and particularly the 2σ contours are increasingly tilted leftwards as the cluster evolves, suggesting the presence of a distinct population of massive PBHs in the cluster core.   Table VI.13. Note the mass distribution is not displayed since, first, the absence of a significant number of mergers implies that the evolution of the total mass profile of objects has an almost negligible evolution, and second, that as seen in said Table VI.13, the cluster and ejecta mass profiles do not segregate significantly over time.
This last point is evidenced too in that we have plotted on top of the distributions a random sample of 1000 PBHs, that is, 0.02% of available objects in each timeslice, with a PBH selection likelihood proportional to its   , where αi,P is the PL index, and corresponds roughly to the constant value of the differential counts, while βi,P in turn, corresponds to the amplitude of the PL. The mass, position and velocity profiles has been fitted to either a Log-Normal or a Maxwell-Boltzmann (MB) distribution. In the case of the Log-Normal fit, then ρz,i,P = (zσz,i, where µi,P is the expected value and σi,P corresponds to the standard deviation, both for the z variable's natural logarithm. In the case of a Maxwell-Boltzmann fit, however, In any case, z = x, v is the fitted variable, and P = C, E stands for the cluster and ejecta populations. The present table then shows the fit that better adjusts to the data at selected times. Overall, we have found that the Log-Normal fit beats the other choices, except for the initial time-slices before the burn-in time in position and velocities. mass to better capture the most massive of these. The largest of these PBHs are visible then in Figure VI.13a as black dots moving rightwards over time as the cluster puffs-up, but always to the left of the bulk of PBHs at their corresponding time-slice.

Cluster dynamical friction
Dynamical friction is the loss of momentum and kinetic energy of bodies by the gravitational interaction with their surrounding neighbours in space. Also called gravitational drag, this effect result on the negative cross-   Table VI.14. Shown are the crossing time, T cross P , and relaxation time, T rel P , of the cluster and ejecta bubbles that populate the galactic halo at selected times, and where P = C, E stands for the cluster and ejecta populations.  correlation of the mass and velocity distribution, with heavier masses skewing towards smaller cluster velocities. It has been observed in a variety of astrophysical systems, from galaxies and galaxy clusters at large scales to proto-planetary disks at smaller scales. In our case, where the PBHs cover a wide range of masses of O(0.1)M − O(1000)M , the lighter PBHs have an effect on the large PBHs being surrounded by them: the lighter, more abundant of these accelerate during slingshot encounters with their more massive partners gaining kinetic energy and momentum, slowing the more massive PBH of the pair in the process as energy is conserved.
Another way this mechanism is enhanced is conversely by the effect a large PBH has on the lighter surrounding PBHs on the denser cluster core: the larger of these PBHs wake up and accelerate the lighter partners towards himself, but by the time the lighter PBH arrives at its position the larger PBH has moved on already. This results in the tendency of large PBHs to be followed by an anisotropic trail of lighter PBHs that exert a gravitational pull opposite to the motion of the large PBH, slowing it down.
Generally speaking, both effects work the same for any range of masses and relative velocities between objects, but the effects are stronger the denser the medium and the slower the bodies' motions, since the gravitational pull is proportional to the square of the masses of the object and to the inverse square of the velocity. Moreover, dynamical friction is heavily suppressed at high velocities, since the larger the velocity, the less time available for the least massive objects to wake up and accelerate towards the most massive ones.
In our simulations, cluster dynamical friction is shown by extracting the ρ(m C , v C , t) distribution at time-slices that divide the simulation time-runs, shown in Figure VI.13b and whose best-fit values we give in Table VI.16, in which a number of features similar to those in the case for mass segregation of Section VI B 1 are apparent i) First is that the phase space (m C , v C )| t occupied by the cluster undergoes an expansion of the velocity profile from the IC's narrow Maxwell-Boltzmann distribution before the t BI = 2.64 × 10 5 yr to the later broader nearly Log-Normal distribution.
ii) Second, that around the burn-in time, just as the cluster contracts to later rebound, the velocities best-fit value of the distribution in the interval (t 28 , t 37 ) = 1.34 × (10 6 , 10 7 ) yr moves towards the higher end of the available range, in agreement with the accelerating infall velocities, the only time to do so throughout the simulations.
iii) Third, that for later times t > t 46 = 1.34 × 10 8 yr the distribution's velocity best-fit is displaced towards the higher end of the available range as the cluster expands, again while not changing significantly its mass value, consistently with the finding in Section VI A 1.
iv) Last, that Figure VI.13b shows velocities that are in fact, for the high mass objects, larger than the bulk velocities of the cluster as seen in that both the 1σ and particularly the 2σ contours are increasingly tilted rightwards as the cluster evolves to the later time intervals. This indicates that the population of core massive PBHs move at high speeds relative to each other.
Dynamical friction is then evidenced by the gradual kinetic energy loss of bulk cluster PBHs and the leftwardmoving velocity profile of Figure VI.13a, particularly so for the 1σ area, less so for the 2σ, but still clear if excluding the largest PBHs.
Note, however that, in contrast to the previous argument that gravitational drag is strongest for massive objects, Figure VI.13a shows that mass positively correlates with velocity. This is apparent too in that we have plotted on top of the distributions a random sample of 1000 PBHs or 0.02% of total objects in each time-slice, with a PBH selection likelihood proportional to its mass to better capture the most massive of these, showing that the most massive objects, with masses of O(1000) M are slingshot to high velocities up to O(100) km/s. The reason to this positive correlation stems from the fact that the cluster is not very dense, with a population of N O,0 = 1000 PBHs at most at the initial time-slice. Its quick expansion makes that the condition for gravitational drag from the trailing PBHs following massive PBHs is not realised in practise, since these large objects do not in fact move in a medium of lighter objects but in one of objects with comparable masses, falling into each other large potential wells and being accelerated to large velocities.

VII. ORBITAL DISTRIBUTIONS
Here in Sections VII A and VII B we find constraints to orbital parameter space by computing the 1σ and 2σ contours of the (a i , e i )| t distributions at the t logarithmically spaced time-slices that divide the simulation time-runs, for the i = C, E cluster and ejecta populations respectively, where e(t) and a(t) is the eccentricity and semimajor axis of PBHs. Also, in Section V A 1 we find the periods ∆t i at time t i , with the cluster and ejecta objects as well as the abundance of ejecta binaries and their orbital parameter distributions with respect to the binary barycentre.
We make a distinction between bounded and unbounded systems, splitting both the cluster and ejecta populations into isolated and bounded systems. as it has been found that they exhibit different properties from the earliest times. Note that eccentricities, semi-major axis and periods are computed with respect to the cluster barycenter in the case of cluster and ejecta isolated objects.
In the case of bounded cluster and ejecta objects, they are computed instead with respect to the bounded system barycenter, regardless of whether they may be part of binary systems, ternary, quaternary and quinary systems. Note yet that the orbital parameters of these bounded systems' barycenters are still included in the isolated population.

A. Cluster orbital properties
In Figure VII.14a we show the PBHs occupied orbital parameter space ρ(a C,I , e C,I , t), for the isolated cluster population, that, is, the PBH cluster population excluding binary, ternary, quaternary and quintic bounded systems. Also, the best-fit values of these distributions are given in Table VII.17. As this is the population of PBHs within the cluster, then it is clear that eccentricities lie in the elliptical range 0 ≤ e C,I , t) < 1 and semi-major axis track the typical cluster length scale at all times, a C,I,t ≈ X C .
A number of features do stand out. The distribution shape is largely uniform across all time-slices but the IC. It shows for t ≥ t 28 = 1.38 × 10 6 yr a peak at (a C,I , e C,I )| t ≈ (1, X C )| t , characterised by the distribution's best-fit value at (a C,I , e C,I )| BF t , and two tails i) The first tail is vertical in (a C,I , e C,I ) orbital phase space, of roughly normally distributed semi-major axes around the time-dependent semi-major axes best-fit value 0.865 pc ≤ a C,I | BF t ≤ 786 pc at each time-slice.
ii) The second tail is horizontal, roughly constant in eccentricities, truncated at the best fit value of e C,I | BF t = 1 and saturated at it. The distribution of eccentricities is then heavily skewed to the higher end of the available range at all times, with a distribution best-fit value at e C,I | t = 1.00 and PBHs consequently describing highly elliptical orbits around the cluster barycentre, and very rare circular orbits, with less than 5% of isolated PBHs with 0.00 ≤ e C,I , | t ≤ 0.15.

Custer binaries orbital properties
The complementary population of bounded cluster PBHs orbital parameter distributions is given in Figure VII.14c along with their best-fit values in Table VII.17, for t ≥ t 28 = 1.38 × 10 6 yr as there are not enough ejecta PBHs to extract reliable confidence regions before this time.
As it it shown there, the cluster, bounded PBHs distribution shape and position in (a C,B , e C,B ) of orbital parameter space does not differ substantially from that of cluster, unbounded PBHs, again exhibiting the aforementioned two tails, one horizontal of constant semi-major axis in the range of 0.128 pc ≤ a C,B | t ≤ 105 pc,    Table VI.16.
Cluster mass segregation: Cluster dynamical friction: asymptotically approaching the parabolic value over time, extending from the best-fit values a C,I | BF t . There is, however, a difference in one particular respect. As seen in Figure VII.14c and compared to Fig-ure VII.14a, binary PBH orbits tend to be more circular than those of cluster PBH orbits. This can be easily seen in the fact that, as opposed to the cluster, isolated case, 95% C.R. in orbital parameter space do reach indeed the (a) Cluster, isolated CF orbital parameters, ρ(e C , a I , t).
(d) Ejecta, bounded CF orbital parameters, ρ(e E , a B , t).   Figure VII.14, where i = C, E stand for the separate cluster and ejecta populations, and j = I, B stands for the separate isolated and bounded population distributions at selected times. e = 0 limit of circular orbits in each of the time-slices considered but that of the initial time, while the previous case eccentricities were, for all considered time-slices, above the threshold of e = 0.15 for 95% of these PBHs. Figure VII.15a shows the lifetimes of cluster bounded sub-systems, the majority of which coming as binary pairs as detailed in Section IV B.

Cluster transient binaries
We have located in total N C,B = 9.61 × 10 5 of such binary pairs by having, at each time-slice, identified all bounded in-cluster PBH sub-systems parent object and binary, ternary, quaternary and quintic pairs. Note that a large fraction of these sub-systems do exhibit continuity in the following time-slices, preserving the same primary and secondary objects identities, even if most do eventually disrupt, hence the label transient applies to them.
However, a non-negligible fraction of these do form by one time-slice only to disrupt by the following timeslice. This suggests again that the total number of binary pairs is underestimated by a non-negligible fraction of the aforementioned N C,B total since some binaries will form and disrupt between consecutive time-slices and will not show at all in our simulations.
We find that, prior to the burn-in time t BI = 2.64 × 10 5 yr, transients have lifetimes never exceeding the δt C,B = 3.0×10 7 yr, indicating that all of the sub-systems present in the IC are disrupted during the first four timeruns, none surviving to the present epoch.
Only later, for t > t BI do cluster binaries start forming with lifetimes typically ranging from a time-varying minimum value, coincident with the time-step at which the binary is formed δt C,B (t) = ∆t C,B (t), which is natural since, as explained, shorter-lived binaries cannot be captured in our simulations due to the discrete data extraction, up to timer greater than the simulation time of δt C,B = 1.38 × 10 10 yr.
Even after accounting for these missed binaries, it is clear from Figure VII.15a that these transient binary pairs do form with a lifetime distribution centred at a time value coinciding with the time-increasing characteristic time of evolution, which roughly equal to the simulation time. The tail of the lifetime distribution is large with lifetime values larger than the simulation time, as indeed some binaries formed as early as t = 9.66 × 10 4 yr, shortly-before the burn-in time, do indeed survive the whole simulation period as binary pairs in wide orbits around the PBH cluster, avoiding disruption by other cluster objects thanks to the low PBH density in such region.

B. Ejecta orbital properties
In Figure VII.14b we show the PBHs occupied orbital parameter space ρ(a E,I , e e,I , t) for the isolated ejecta population, considering only the PBH cluster population and excluding the binary, ternary, and quaternary bounded systems. Moreover, in Table VII.17 we give the best-fit values of these distributions. Since this is the population of ejecta PBHs, eccentricities are constrained to the hyperbolic range of 1 ≤ e C,I , t) ≤ ∞ and semi-major axis are inevitably larger than those of cluster PBHs, a C,I,t ≈ X E ≥ X C .
A number of issues are worth discussing. The distribution shape is again largely uniform across all time-slices, this time with less variation between the IC, and later evolution than it was the case for cluster PBHs. The shape is, naturally, very different in any case from that of cluster PBHs.
It does, however, exhibit some common characteristics. Just like in the cluster case, starting at t ≥ t 28 = 1.38 × 10 6 yr the distribution exhibits a peak at the bestfit value at (a E,I , e E,I )| BF t that asymptotically approaches the parabolic case e = 1, deviating from it by less that 5% by the time the simulations have entered the fourth run at t 37 = 1.38 × 10 7 yr.
Moreover, the distribution also exhibits two tails: i) A first slanted tail in (a E,I , e E,I ) orbital phase space towards smaller semi-major axis and higher eccentricities, starting at the best-fit value and including PBHs closer to the cluster with lower semi-majoraxes. In this case, it is found that eccentricity negatively correlate with semi-major axis and the best-fit value bound to the range 2.73 pc ≤ a E,I | t ≤ 200 pc, with the tails easily extending by a factor of O(100) of this value.
ii) A second tail that is horizontal in orbital phase space, roughly constant in eccentricities, truncated at the best fit value of e C,I | BF t = 1 and saturated at it just like it was the case for cluster objects. The distribution of eccentricities is this time heavily skewed to the lower end of the available range of eccentricities at all times, with a distribution best-fit value in the interval 1.00 ≤ e C,I | t ≤ 1.32 and PBHs consequently describing hyperbolic orbits, asymptotically approaching the inferior bound as time progresses.

Ejecta binaries orbital properties
The complementary orbital parameter distributions of for the population of bounded ejecta PBHs is given are Figure VII.14d along with their best-fit values in Table VII. 17. In this case, the ejecta, bounded PBHs distribution shape and position in (a E,B , e E,B ) of orbital parameter space does differ substantially from that of ejecta, unbounded PBHs, in contrast with cluster PBHs, where the distribution was roughly universal for both isolated and bounded objects. Now the shape does not exhibit the aforementioned two tails, but only the first slanted one, extending from the best-fit values a C,I | BF t . Orbital semi-major axis best-fits are then constrained to 0.217 pc ≤ a E,B | t ≤ 498 pc, while orbital eccentricities best fits are constrained to again asymptotically approaching the parabolic threshold as time progresses.
It is again the case that the orbital semi-major axis correlates negatively with eccentricities. Moreover, as seen in Figure VII.14d, and compared to Figure VII.14b, binary PBH orbits tend to approach the parabolic limit of hyperbolic orbits faster than those of cluster PBH orbits.
This can be easily seen in the fact that, as opposed to the ejecta, isolated case, the 95% C.R. in orbital parameter space do never exceed the e = 47 threshold for the first considered time-slice, decaying to the e = 20 level in the case of the last time-slice; while in the previous case eccentricities were below the threshold of e = 1.1 for the first considered time-slice, increasing to the e = 330 level in the case of the last time-slice, for 95% of these PBHs.

Ejecta stable binaries
Figure VII.15b shows the lifetimes of ejecta bounded sub-systems, the majority of which as simple binary pairs as detailed in Section IV B.
We have located in total N C,B = 1.87 × 10 5 of such binary pairs by having, at each time-slice, identified all bounded in-ejecta PBH sub-systems parent object and binary, ternary, quaternary and quintic pairs, in a similar way that we did with cluster objects. Note that, in this case, the vast majority of these sub-systems exhibit continuity in the following time-slices, preserving the same primary and secondary objects identities, given that only a negligible fraction of these do form by one time-slice only to be disrupted right after. This is due to the fact that once in the ejecta population, the chance of a PBH-PBH encounter that may disrupt the binary pair is extremely small and increasingly suppressed over time as our simulated cluster live in isolation and the PBH density quickly and monotonically decays as objects move away from the cluster core. Because of this, and unlike what was the case for cluster transients, this ejecta binary pairs are labelled stable, and are not underestimated by the discrete time-step of collecting information in the simulations.
We find that, prior the burn-in time t BI = 2.64×10 5 yr, transients have lifetimes never exceeding the δt C,B = 2.0×10 5 yr, indicating that all of the sub-systems present in the IC are disrupted during the first three time-runs, none surviving even to the burn-in time.
Only later, for t > t BI do ejecta binaries start forming with lifetimes typically ranging again, as it was the case for cluster binaries, from a time-varying minimum value coincident with the time-step at which the binary is formed δt C,B (t) = ∆t C,B (t), up to timer greater than the simulation time of δt C,B = 1.38 × 10 10 yr.
Note that, as Figure VII.15a suggests, this stable binary pairs do form with a lifetime distribution centred at a time coinciding with the monotonically-increasing characteristic time of evolution, which is roughly equal to the simulation time. The tail of the lifetime distribution is large with lifetime values larger than the simulation time, as indeed some binaries form as early as t = 4.14×10 5 yr.

C. Hyperbolic encounters power spectra
In this section we aim to extract the spectrum of hyperbolic encounters in our simulations. Ideally, we would extract the distribution by identifying the hyperbolic encounters themselves at the moment of the closest approach between the bodies, and computing the distribution of the relative orbital and dynamical parameters from there.
However, we had seen in Section V A 1 that a large majority of these hyperbolic encounters cannot be identified by this method as it is found that the time interval at which the simulation data is extracted is not short enough to capture accurately the trajectories of the participating PBHs. Two PBHs may be closest neighbours at one time-slice, yet at the next time-slice there may be tens of other PBHs closer to them even if they were close enough originally to have undergone a hyperbolic encounter themselves.
Instead, we aim to extract the spectrum of potential hyperbolic encounters in our simulations. We do this by computing the relative eccentricities and impact parameters of all cluster PBHs with respect to their closest 100 neighbours at selected time-slices, since those are assumed to be produce the majority of the actual encounters in the simulations. For each of the pairs, then we extract the masses of both PBHs in the pair m i , m j ), as well as the relative positions x rel i,j = x j − x i and relative velocities v rel i,j = v j − v i where i, j denote the j the closest PBH to the i the PBH in the cluster.
We compute the relative eccentricity of the pair then by making use of the expression where h rel i,j = x rel i,j × v rel i,j and G N is the gravitational constant in the appropriate units.
Then, we compute the pair semi-minor axis by making use of the expression which corresponds to the impact parameter of an hyperbolic encounter when the relative eccentricity exceeds the parabolic limit, for e rel i,j > 1. To ensure that our results are robust, we have repeated the power spectrum extraction algorithm for eccentricities and semi-minor axis where we had considered the 100 closest neighbours for the first, middle and last 1/3 of objects spanning the 1 st -to-33 rd , 34 th -to-67 th and 68 th -to-99 th closest neighbours and found that there power spectra does not change significantly in any of these cases, the only noticeable effect being a widening of the power spectrum to larger impact parameters if more than 100 neighbours are considered as it is be natural to expect in any case.  The result of this computation is presented next. We show the results of this computation, for a total of N rel i,j = 100N R N O = 10 8 pairs of PBHs per time-slice in all considered time-slices but the IC, in Figure VII.16a for the relative eccentricities and Figure VII.16b for the relative impact parameters. We do as well fits the resulting curves to Log-Normal distributions in Table VII.18, having tested prior to that other distributions and found that a Log-Normal distribution describes best the data.
We find in Figure VII.16a that, for all the considered time-slices, there is a significant and roughly timeindependent fraction of i, j pairs that may potentially produce hyperbolic encounters with e rel i,j > 1. Moreover, this encounters may be produced with eccentricities as large as large as e rel i,j = O(1000). Note that the fraction of hyperbolic encounters is small with a value of 34% at the IC at the initial time as one would expect from the smaller velocity dispersion of the Maxwell-Boltzmann distribution as compared to the velocity distribution at later stages of the evolution.
However, the fraction of hyperbolic encounters with respect to the total number of encounters does grow monotonically to a value of 64% at the beginning of the fourth time-run (t 28 = 1.38 × 10 6 yr) and reaching a final value of 71% at the present (t 64 = 1.38 × 10 10 yr).
Conversely, the fraction of potentially binding encounters decreases from 66% in the IC at the initial time to 36% at the beginning of the fourth time-run (t 28 = 1.38 × 10 6 yr) and finally arriving to 29% at the present (t 64 = 1.38 × 10 10 yr), in line with the overabundance of binaries in the IC and the monotonically decreasing rate of transients throughout as the cluster evolves. The bump in the power spectrum noticeable at eccentricities in the eccentricity range e rel i,j ≤ 0.1 is produced by real, and not potential, in-cluster binary pairs, and has, as expected, much less power than that the part of the power spectrum dominated by potential encounters in the ec-centricity range of e rel i,j ≤ O(0.01). We find as well in Figure VII.16b that, again for all the considered time-slices t i , i = 0, 28, 37, 46, 55, 64 that divide the last four runs and the IC, the range covered by the impact parameter for all potential encounters is constrained to the interval 0.001 pc ≤ b rel i,j ≤ 1000 pc. However, and unlike what was the case with eccentricities, the range varies greatly in time from the IC at t 0 = 0 yr where to the present time where tracking the characteristic length scale of the cluster in time.
It is clear, however, that the vast majority of all pairs here considered will not produce close encounters. Only a very small fraction of the total number of pairs will indeed be at distances of the order of 10 −5 pc ≈ 2 AU as for later times the power spectrum is noticeably suppressed at lengths scales of b rel i,j ≤ O(0.001) pc at the initial time, and of b rel i,j ≤ O(0.1) pc at the present time. Considering, however, that there are in total N rel i,j = 10 8 pairs potentially producing the encounters, then the likelihood of at least a few of these occurring during the whole evolution time remains large.
Note that the cluster infall and rebound stage described in Sections VI A 2-VI A 4 at around the burn-in time is apparent as well in the impact parameter profile evolution, which exhibits the same behaviour. The impact parameter profile at the IC is skewed towards values larger than the impact parameter distribution at the beginning of the fourth run (t 28 = 1.38 × 10 6 yr), indicative of the contraction of the cluster up to that time, and more importantly, of the concentration of PBHs in a cuspier core after the IC is erased after the burn-in time. Only after the beginning of the fifth run (t 37 = 1.38 × 10 7 yr) does the IC begin to skew towards values smaller than the impact parameter distribution at that time, increasingly so as the cluster puffs-up and the typical impact parameter grows accordingly.
Last, it is worth mentioning that no significant correlation has been found between eccentricity and semi-minor axis values for this number of pairs. Somewhat contrary to expectation, smaller impact parameters do not correlate with smaller eccentricities, as one would expect from dynamical arguments. The reason for this is that the number of neighbours considered for each PBH has been chosen, as explained in the beginning of the present section, to be small enough that the distributions are consistent with those extracted taking in only the nearest neighbours among those. The shell centred on each of the PBHs typically extends to a distance smaller than the cluster size by O(10), and is approximately homogeneous so it does not contain wildly different velocities or masses on average.

VIII. BACKGROUND DM IMPLICATIONS
In this section we aim to quantify the implications of the proposed PBHs clusters in the global DM distribution. In particular, Section VIII A will deal with predictions for the DM background of cluster and ejecta objects separately for a typically sized galactic halo, while in Section VIII B we will compute the distinct hyperbolic encounter and merger event rates of PBHs per comoving volume.

A. DM energy density
It was mentioned in Section VI A 3 that the trajectory stability of distant objects bound to the cluster or belonging to the ejecta sphere-of-influence is guaranteed in our simulations, since both the cluster and ejecta sphereof-influence lives in isolation and the likelihood of an encounter that may transfer enough kinetic energy to overcome the cluster escape velocity is very low, as the ejected objects travel the cluster periphery where this far-off objects live very fast, and also given that since the periods of these objects are very high.
However, in a realistic scenario this is not the case anymore. A single galactic halo may contain, as we will soon calculate, O 10 7 of these PBH clusters with their corresponding cluster and ejecta spheres-of-influence extending up to order kpc in the former case and reaching up to hundreds of kpc in the latter case. These numbers are large enough that eventually the volumes occupied by both the cluster and ejecta spheres-of-influence overlap, and so interaction between PBHs is no longer a remote possibility but a reality.
In the following, we aim to compute precisely the number of these PBH bubbles and the time at which they do interact or overlap within a typical galactic halo. For this computation, we take as our baseline a Milky Way-like halo, whose total mass Ref. [83] estimates to be at the 50% C.L. Moreover, we will assume a local DM fraction of f DM = Ω DM /Ω 0 | MW , and second, that these clusters constitute a particular DM fraction f PBH = Ω PBH /Ω DM .

Cluster DM energy density
Depending on the total mass profile of the PBH clusters, one can constrain their abundance and compute the mean inter-cluster distance. In our case, considering our clusters have, typically, a mean total mass of , (163) given that for our Log-Normal mass distribution the actual mass mean ism = exp µ m + σ 2 m /2 . The mass is itself contained in a sphere of a mean effective radius given by where the input logarithmic mean position, µ x and and deviation σ x can be used again to get the actual position mean,x = exp µ x + σ 2 x /2 , as the cluster profile was shown to be best described by a Log-Normal distribution in all time-slices but the very first ones before the IC is erased after the burn-in time.
The cluster mass, position and velocity distribution parameters µ i and σ i for i = m, x, v are given in Table VI.13, while the cluster characteristic scalesm,x and v are shown in Table VI. 15.
We now aim to estimate the number of simulated PBHs cluster bubbles in a typical galactic halo in order to infer what would be the population os these objects within a typical galaxy, their mass, mean inter-cluster distance and mean inter-PBH distance at the point where the cluster bubbles have expanded to the point of overlapping with each other. In order to do this, we will take in the following our own Milky Way halo as the standard of a typical galactic halo.
The number of these PBH clusters within the X < X * MW sphere centred at the Milky Way core is, considering the halo mass and effective radius from Ref. [83], found to be    Figure VII.16 distribution Log-Normal fits: ρz,i,P = (zσz,i,P √ 2π) −1 exp −(Log-Normalz − µz,i,P) 2 /2σ 2 z,i,P where µi,P is the expected value and σi,P corresponds to the standard deviation of the z variable's natural logarithm at selected times, and P = C, E stands for the cluster and ejecta populations respectively and z = e, b for the mass, position and velocity functions. Note that, fits to Rayleigh χ 2 2 and Maxwell-Boltzmann χ 2 3 distributions have also been tried for the position and velocity profiles, having found that, except for those time-slices before the burn-in time, a Log-Normal fit beats the former trials. Whenever none of the above fits meaningfully describe the data, r 2 < 0.8, no fit is provided.
where we have not corrected for the possible loss of ejecta PBHs out of the sphere since the fraction of M tot C (t) travelling such long distance is negligible, with about O(0.001) of PBH being slingshot away to distances larger than X < 100 kpc, being among the lightest of simulated PBHs on top of it. Correspondingly, the number of these in-cluster PBHs within the X < X * MW sphere centred at the Milky Way core is then Assuming now for simplicity a homogeneous distribution of PBH clusters within the halo sphere, then the mean inter-cluster distance within said sphere can be well approximated by X tot C (t) (V MW /N tot C (t)) 1/3 , where V MW stands for the Milky Way volume contained in the X < 100 kpc sphere, or after some substitutions and correspondingly, the mean inter-PBH distance within said sphere can be well approximated by assuming that the cluster have puffed-up enough to start overlapping with each other. Note that the mean inter-PBH distance X ind C (t) will be only of significance one clusters begin to overlap with each other, as prior to that point the PBH two-point correlation function will be bimodal function, and only when the clusters sufficiently overlap with each other will the PBH two-point correlation function be a unimodal function as PBH-to-PBH distances are all roughly the same. Now, each of these PBH cluster bubbles will contain a monotonically decreasing mass over given by the M tot C (t) of Eq. (163), while all the PBH clusters in the galactic halo will total an also monotonically decreasing joint mass of Last, the mean occupied PBH cluster volume fraction of said sphere is then F tot stands for the cluster volume contained in the cluster effective cluster radius X * C = 100 kpc, which leads to Assuming a DM fraction within the Milky Way halo sphere of f DM = 0.8 entirely due to PBHs, f PBH = 1 then Eq. (165), Eq. (167) and Eq. (170) yield a total cluster number per sphere, N tot C of separated on average by a distance, X tot C of that is, we have in the different realisations O 10 7 cluster bubbles with their own spheres of influence, separated by distances of O(100) pc, which evaporate and expand as the simulations evolve.
The estimated values of the mean individual PBH number N ind C (t) and mean individual PBH distance X ind C (t) are given in Table VIII. 19 for the cluster population, following the aforementioned procedure along with the mean mass per cluster bubble M tot C (t), mean in-cluster mass per halo M ind C (t), and cluster occupied halo volume fraction F tot C (t), at the te-slices corresponding to those that divide the last four time-runs plus the IC.
An interesting phenomenology arises from this computation. Table VIII. 19 shows that indeed the cluster objects are perfectly isolated initially, even though, at late times and soon after entering in the last time-run when t * ∈ (1.38, 10) Gyr, that is not longer the case, since F tot C,64 (t > t * ) ≈ 1 and these puffed-up clusters overlap. Then the typical inter-cluster PBHs and intra-cluster PBHs distances converge to a single value and the cluster distribution within the galactic halo approaches that of a homogeneous distribution where the outer layers of the clusters overlapping with each other. The cuspy cores of the original clusters remain distinct albeit less so as the simulations evolve approaches the present time. In particular, we find that the typical PBH-to-PBH distance regardless of whether the two originate in the same or separate cluster bubbles approximately is X ind C = 112 +20 −23 pc and of the same order of magnitude than the cluster-tocluster distance.

Ejecta DM energy density
We now aim to repeat the same analysis of Section VIII A 1 for ejecta objects, and constrain the abundance and compute the mean inter-ejecta spheres-ofinfluence distance, taking again into account that these PBH ejecta bubbles do not live in isolation.
Given the very low level of merging in our simulations, we can assume that the cluster and ejecta populations are complementary to each other and adding up to the initial cluster population since the number of ejecta objects is negligible at the time of the IC. Then, the ejecta will typically have a mean total mass of given that, as we had seen in the previous section, the PBH mean mass ism = exp µ m + σ 2 m /2 . The bulk of the mass is then contained in a sphere of a mean effective radius given by similarly to what was the case of the cluster bubble, and centred in it, but larger, by a growing factor of roughly O(100) at present times. Note that the mean position is x = exp µ x + σ 2 x /2 , as the ejecta position distribution is also best described by a Log-Normal distribution in all time-slices.
The ejecta mass, position and velocity distribution parameters µ i and σ i for i = m, x, v are given in Table VI.13, while the ejecta characteristic scalesm,x and v and are shown in Table VI.15. Since each ejecta bubble originates by the PBH expelled from the cluster at its centre, the number of these PBH ejecta bubbles within the X < X * MW sphere centred at the Milky Way core is, Correspondingly, the number of these in-ejecta PBHs within the X < X * MW sphere centred at the Milky Way core is then For a homogeneous distribution of PBH clusters within the halo sphere, then the mean inter-ejecta distance within said sphere need to track the mean inter-cluster distance, and so and correspondingly, the mean inter-PBH distance within said sphere can be well approximated by assuming that the ejecta have puffed-up enough to start overlapping with each other, a safe assumption given that it was shown that this was indeed the case for the cluster PBHs in the last run, and the fact that the ejecta bubbles reach to far larger distances from early on. Note that every PBH ejecta bubble will have a monotonically increasing mass over given by the M tot E (t) of Eq. (173), while all the PBH ejecta in the galactic halo will total as well an also monotonically increasing joint mass of The mean occupied PBH ejecta volume fraction of the ejecta sphere-of-influence will be then Last, assuming again a DM fraction within the Milky Way halo sphere of f DM = 0.8 entirely due to PBHs, f PBH = 1 then Eq. (175), Eq. (177) and Eq. (180) a total ejecta number per sphere, N tot C of N tot E = N tot C = (2.7 +1.7 −0.9 ) × 10 7 , separated on average by a distance, X tot C of X tot E = X tot C = 800 +190 −100 pc.
The estimated values of the mean individual PBH number N ind E (t) and mean individual PBH distance X ind E (t) are given in Table VIII. 19 for the ejecta population, along with the mean mass per ejecta bubble M tot E (t), mean in-ejecta mass per halo M ind E (t), and ejecta occupied halo volume fraction F tot E (t), at the time-slices corresponding to those that divide the last four time-runs plus the IC.
Many interesting features regarding the behaviour of the ejecta bubbles are already apparent in Table VIII.19, such that, while the ejecta bubbles sphere-of-influence are isolated initially, already by the fifth run they have started to overlap. By the present time the mixing between ejecta bubbles is complete as F tot E, 64 1 and it makes no longer sense to treat them as anything other than a uniform background, unlike the case of the cluster bubbles which still preserve a distinct, cuspy core.
In particular, we find that the typical PBH-to-PBH distance regardless of whether the two originate in the same or separate ejecta bubbles is approximately of O(100) pc, slowly decreasing as the now galactic halosized joint ejecta sphere-of-influence slow expansion is overcome by the gradual increase in population of the ejecta bubbles from continuing evaporation in the cluster, from at t 46 = 1.38 × 10 8 yr to a present value of X ind E = 93 +24 −13 pc.
Note, however, that in this computation we have neglected the fact that while our simulations have both the cluster an ejecta bubbles in isolation, at this scales that is clearly not the case anymore. Ejecta objects in particular will by now feel the potential wells of their neighbour bubbles and will at a point cease to expand when they finally populate the entire galactic DM halo. Individual ejected PBHs will not move in the asymptotically uniform rectilinear motion and their trajectories will bend to orbit the galaxy, most of them in the very little dense environment of the galactic halo outskirts, behaving as a virtually collisionless, cold DM component, given the extreme unlikelihood of a close encounter out if the cores of the cluster bubbles, let alone in the far less dense environment of the ejecta background.

B. Gravitational event rates
We aim now to compute a minimum threshold for the average hyperbolic encounter rate and merger event rate in time per comoving volume, starting from the computed rates for both per cluster in Sections V A 1 and V B 1.
In order to do so, we begin by computing the total DM mass contained in a comoving box of volume 1 Gpc 3 . The total DM energy density of the Universe is ρ DM = (3.965 ± 0.073) × 10 19 M Gpc −3 , according to Planck 2018 (TT, TE, EE+lowE, see Ref. [6]).
The initial cluster total mass is given by Eq (163), which reduces to M tot C,0 = N OmC (t 0 ) thus obtaining which is a large number indeed, albeit an expected one since it is known that a 1 Gpc 3 comoving box contains about O 10 7 Milky Way-sized galactic haloes.

Comoving hyperbolic encounter rates
We do now take the hyperbolic encounter rates of Section V A 1 and scale them to a 1 Gpc 3 -sized box with the transformation Γ S r → N tot C Γ S r . The final comoving slingshot rates are then These are indeed very large rates, however the same caveats that applied to the slingshot rates of Section V A 1 still apply here; mainly that while the number of hyperbolic encounters is underestimated because of the discrete time-step in each of the time-runs and the non-relativistic nature of the N-body code, we also find too many such events to be understood as gravitational wave generating events. Indeed most encounters most often happen at PBH-to-PBH distances way too large to detect the corresponding emission of GWs with current instruments.

Comoving merger event rates
We do at last take the merger event rates of Section V B 1 and scale them again to a 1 Gpc 3 -sized box with the transformation Γ M r → N tot C Γ M r . The final comoving merger event rates are then Γ M 6 (46 ≤ t i ≤ 55) = 852 ± 26 yr −1 Gpc −3 , The merger event rates segregated by their origin, would be, for all 51 mergers occurring from previously isolated PBHs Γ M,C,I 6 (46 ≤ t i ≤ 55) = 213 ± 11 yr −1 Gpc −3 , Note however, that the same caveats that applied to the merger event rates of Section V B 1 still apply here; mainly that we do underestimate mergers to some extent since the evolution is computed in a non-relativistic manner, and thus this rates may have to be interpreted in practice as lower bounds to the real PBH merger rates.
Indeed the lack of inspiral may be particularly relevant at the beginning of the simulations, as we have explained how an abundance of merging at early times has the effect of both increasing the mass density at the core of the cluster and decreasing the cluster total mass at later times by stripping the cluster of its outer layer due to in-falling PBHs acquiring large infall velocities only to be dispersed later on on a one-to-one encounter with another PBH at high relative speeds.
Note that, despite the large number of mergers produced per comoving volume, there is no tension between the total comoving merger rate from all population sources of PBHs found at present, Γ M 7 (56 ≤ t i ≤ 64) = 2409 ± 74 yr −1 Gpc −3 and the typically O(100) − O(1000) values estimated from LIGO's runs O1 and O2 (see Ref. [84]), since LIGO is not in any case sensitive merger events in which the participating BHs have total masses as large as those found in our simulation, typically of order 1000 M , though future experiments such as LISA will (see Ref. [85]).

IX. CONCLUSIONS
The dynamics of PBH in clusters is one of the cornerstones in our understanding of structure formation at small scales, as well as on the way classical PBH constraints can be resolved. Furthermore, PBHs have unique signatures that could be detected by current and future GW detectors, such as AdvLIGO and LISA. In particular, close hyperbolic PBH encounters in dense clusters, such as the ones we simulated here, emit bursts with millisecond durations which may be detected depending on the duration and peak frequency of the emission, but also the rates and distributions of the PBHs.
In this paper we found that PBHs can constitute a viable DM candidate, whose clustering offers a rich phenomenology. We have quantified the stability and rate of evaporation of PBH clusters, obtained the number of PBHs that remain in the cluster and those that are ejected from it to the background, also the fraction of these that show up in gravitationally-bound subsystems or merge with each other, as well as computing their parent and merger tree histories.
Furthermore, these PBH clusters are indeed stable and survive until present times retaining about a third of their total initial mass, while the remaining two thirds is transferred to the background ejecta, providing a more or less uniform spatial distribution of both free PBH and the cores of the original PBH clusters with a wide range of masses. In particular, we find that cluster puffing up and evaporation leads to PBH subhaloes of O(1) kpc in radius containing at present times about 36% of objects and mass, with O(100) pc cores. We also find that these PBH sub-haloes are distributed within larger PBH haloes of O(100) kpc, containing about 63% of objects and mass, coinciding with the sizes of galactic halos.
We have also extracted the component mass profiles and mass spectrum of binary and merger pairs in these clusters. In particular, we find that binary systems do appear frequently, constituting about 9.5% of all PBHs at present, with mean and median mass ratios ofq B = 0.154 andq B = 0.0552, and total masses ofm T,B = 303 M and m T,B = 183 M , and larger, less skewed binary systems for those PBHs in the ejecta rather than those in the cluster.
Alternatively, we find that mergers are very infrequent, with barely 0.0023% of all PBHs merged at present, having mean and median mass ratios ofq M = 0.965 and q M = 0.545. Merger total masses masses found arē m T,M = 1670 M andm T,M = 1510 M , while found chirp masses arem c,M = 642 M andm c,M = 567 M . Also, we find that the PBHs resulting from the largest mergers tend to segregate themselves to the cluster cores, and that merged PBHs escaping to the ejecta are often less massive and more skewed.
We have looked at the one-to-one close interactions of these PBHs and extracted lower bounds to the hyperbolic encounter and merger event rates produced. The simulations have a number of shortcomings that make a proper comparison with observations somewhat difficult, particularly so for the hyperbolic encounter rates. As for the merger event rates, we have found rates for massive black holes beyond the range of sensitivity of LIGO-VIRGO observations. We find close encounter rates to be Γ S = (1.2 +5.9 −0.9 ) × 10 7 yr −1 Gpc −3 and merger rates of Γ M = 1337 ± 41 yr −1 Gpc −3 . These values are thus far compatible with observations, and might be probed by future space-based observatories like LISA.
We have also computed the cluster and ejecta dynamical parameter space distributions, and found that the mass, position, velocity and density profiles and found them to be consistent with that of a cold DM component. Interestingly, at present times the clusters overlap just enough that their outer layers are indeed in contact with each other while retaining their separate cores, having length scales comparable to those of dwarf galaxies and globular clusters. Also interesting is the fact that the ejecta background once the many different ejecta bubbles merge has length scales comparable to medium sized galactic haloes. Moreover, we have studied the degree of cluster mass segregation and dynamical friction, and found those to be of no great significance.
We then determined the cluster and ejecta orbital parameter-space distributions of PBH clusters at selected times and the segregation by bounded and unbounded PBHs and studied their evolution. We have found that the profiles for both the orbital eccentricity and semimajor axis are compatible with current observations, that there is a large number of cluster transient binaries, disrupted by third encounters, and stable ejecta binaries that may be responsible for the events detected by LIGO. Last, we have computed the power spectrum of potential hyperbolic encounters between close cluster PBH pairs.
Overall, we have found that inhomogeneously distributed PBH clusters are a good candidate for a cold, nearly collisionless DM, consistent with astrophysical observations of the DM distribution and replicating the features that a good DM candidate must uphold, throughout the entire history of the Universe starting in the matter era.
It remains an open question how to extrapolate these results to a more general set of initial conditions in which the mass, position and velocity distribution of the PBHs may be different and the PBH cluster initially more concentrated. We leave the study of a more varied set of PBH clusters for future work.