Next Article in Journal
Machine Learning Using Rapidity-Mass Matrices for Event Classification Problems in HEP
Next Article in Special Issue
Quasi-Periodic Oscillatory Motion of Particles Orbiting a Distorted, Deformed Compact Object
Previous Article in Journal
A Superfluid Perspective on Neutron Star Dynamics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Clustering Dynamics of Primordial Black Holes in N-Body Simulations

by
Manuel Trashorras
,
Juan García-Bellido
* and
Savvas Nesseris
Instituto de Física Teórica UAM-CSIC, Universidad Autonóma de Madrid, Cantoblanco, 28049 Madrid, Spain
*
Author to whom correspondence should be addressed.
Submission received: 26 December 2020 / Revised: 4 January 2021 / Accepted: 7 January 2021 / Published: 15 January 2021
(This article belongs to the Special Issue Compact Astrophysical Objects)

Abstract

:
We explore the possibility that Dark Matter (DM) may be explained by a nonuniform background of approximately stellar mass clusters of Primordial Black Holes (PBHs) by simulating the evolution 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 q ¯ B = 0.154 , and total masses of m ¯ T , B = 303 M . Merged PBHs are rare, about 0.0023% of all PBHs at present, with mass ratios of q ¯ B = 0.965 with total and chirp masses of m ¯ T , B = 1670 M and m ¯ c , M = 642 M , 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 one-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 1000 M ), with Γ S = ( 1.2 + 5.9 0.9 ) × 10 7 yr 1 Gpc 3 and mergers with Γ M = 1337 ± 41 yr 1 Gpc 3 .

1. Introduction

One of the most pressing problems in cosmology is that of the nature of Dark Matter (DM), which was first suggested in the 1930s by Zwicky from the observation of motion anomalies in the Coma galaxy cluster in [1], and by Rubin in the 1970s from the observation of an excess velocity in the tails of galaxy rotation curves in [2]. Evidence for DM is now abundant, from structure formation, galaxy dynamics, lensing, and Cosmic Microwave Background (CMB) observations (see in [3,4,5,6]), largely in agreement with a DM density of Ω c h 2 = 0.120 ± 0.001 (see in [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 in [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 in [8]) populating galactic halos.
However, and in spite of numerous efforts, no viable DM particle candidate has been found in experiments like the Large Hadron Collider (LHC), nor has a convincing signal been detected either in ground or underground facilities, directly or indirectly, with cryogenic detectors such as the Cryogenic Dark Matter Search (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. Therefore, 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 in [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 in [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 [21].
Primordial Gravitational Waves (PGWs) arising from the interaction of these PBHs may also prove to be a useful tool to constrain the period of inflation in the early Universe, through their abundance, clustering, slingshot and merger event rates and dynamical profiles. Moreover, a Stochastic Gravitational Wave Background (SGWB) arising from the PBH–PBH interaction history from the Early to the late Universe may be detected in the future (see in [17,22]) by spaced-based interferometers, Fast Radio Bursts (FRB) [23,24], Pulsar Timing Arrays (PTA) [25,26] and 21 cm experiments [27,28,29,30].
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 mass range, particularly if BH clustering and BH mass spectra are considered (see in [31,32,33,34]).
If so, these PBHs may offer a rich phenomenology that can not only explain the nature of DM, but that also explains some questions regarding early structure formation [35]; 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 in [36] for a review of these issues).
Our paper is organised as follows. In Section 2, we present the theoretical background of PBH formation, while in Section 3 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 repository1. 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 4 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 5 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 6 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 7 we examine the orbital distributions of PBHs and their repercussions for either transient or stable binaries. In Section 8, we present estimates of the hyperbolic encounter and merger event rates for a medium-sized galactic halo and co-moving cosmological volume.
Finally, in Section 9 we present our conclusions, summarising our findings and comment on the feasibility of PBHs as a DM candidate according to our results.

2. Theoretical Background

2.1. Formation Mechanisms

First proposed in the 1960s when Zel’dovich and Novikov (see in [37]) 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 in [38]) and Carr (see in [39,40]) 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 in [41,42]), 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 [42]. 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 in [15]), single-field inflation with an inflection point (see in [43,44,45]) or by having the inflaton coupled to gauge fields (see in [46,47]). 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 in [48]) 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 in [49]) 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, as 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 in [50]).
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 core collapse type II supernova, and whose masses are restricted to be heavier than about 5 M but lighter than O ( 30 ) M (see in [51]). 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.

2.2. 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 in [52])
(ii)
Lensing constraints from femtolensing of gamma-ray bursts (see in [53]) and millilensing of compact radio sources (see in [54]), as well as microlensing constraints from the EROS collaboration (see in [55]), the OGLE collaboration (see in [56]), HSC Andromeda observations (see in [57]), Kepler observations (see in [58]) and caustic crossings (see in [59]).
(iii)
Dynamical constraints from white dwarf disruption (see in [60]), neutron star capture (see in [61]), wide halo binary disruption (see in [62]), dynamical friction on PBHs (see in [63]) and compact stellar systems in ultra-faint dwarf galaxies and Eridanus II (see in [64,65]).
(iv)
Accretion constraints from the energy injection by PBHs in the early Universe before photon–electron decoupling by spherical accretion (see in [66]) and accretion disks (see in [67]) on PBHs, as well as radio (see in [68]) and X-rays (see in [68,69]) emission from present PBHs in the late Universe.
Moreover, while a monochromatic PBH mass power spectrum of PBHs is excluded by observations, this kind of spectrum is not very realistic, as 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 in [15,20,36,70,71]).
Furthermore, 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 in [33,34,72]).
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 in [73].
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 [33,72].
More importantly, as 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.
This, 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 in [74,75]. 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 Section 5 and Section 8.

3. Methodology

We now proceed to describe our simulations, the essential properties of which are outlined in Section 3.1, the particular choice of ICs in Section 3.2 and the information available in the time slices in Section 3.3 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 in [76]), a Java numerical N-body 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 General 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.
(i)
For the masses: 1 degree of freedom, taken from a random realisation of a Log-Normal (LN) distribution with ( μ , σ ) m = ( 2.0 , 1.5 ) M .
(ii)
For the 3D position vectors: 3 degrees of freedom, taken from a random realisation of an isotropic, centred Multivariate-Normal (MN) distribution with M v = 0 pc and Σ v = 1.0 × 𝟙 3 pc 2 .
(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 3.2, 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 in [77]) a full relativistic code or GADGET-2.0.7 (see in [78]) a cosmological hydrodynamical code, which adapt 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, frame-dragging, 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.

3.1. 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.

3.1.1. Time-Resolution

The code has an adaptive time resolution δ t , which depends on the forces at ith-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.

3.1.2. 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 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 in [79]) for the PBH pair binding is
σ = 85 π 3 2 / 7 r S 2 v c 18 / 7 π Δ b 2 ,
where b is the limiting impact parameter and the Schwarzschild’s radius is
r S = 2 G m c 2 ,
where G is Newton’s constant, m the mass of the object and c the speed of light in vacuum. 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
0.1 M m 100 M ,
0.1 km / s v 1000 km / s ,
of masses and velocities in the simulations, then Equation (1) renders an effective spatial resolution of Δ x = 100 AU , less than the mean minimum distance in between bodies, of O ( 1000 ) AU .

3.2. 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 4 N 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)
3 N 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)
3 N 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 1. This choice of IC leads to an initial density ρ 0 3 N O μ m / 4 π a x 3 of O ( 1000 ) M pc 3 . 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 1.
Table 1. Mass, position and velocity IC distribution parameters. Distributions plots are shown in Figure 2a, Figure 3a and Figure 6a. A specific realisation of such distributions is shown in Figure 2b, Figure 3b and Figure 6b.
Table 1. Mass, position and velocity IC distribution parameters. Distributions plots are shown in Figure 2a, Figure 3a and Figure 6a. A specific realisation of such distributions is shown in Figure 2b, Figure 3b and Figure 6b.
Mass ICs Distribution:Position ICs Distribution:Velocity ICs Distribution:
D m m i 0     μ m     2.0 M D x x i 0     μ x     1.6 pc D v v i 0     μ v     2.0 km / s
[ LN ]     σ m     1.5 M [ MB ]     a x     1.0 pc [ MB ]     a v     1.2 km / s
Figure 1. Cluster IC positions and velocities. Both the position and velocities correspond to a sample of an isotropic Maxwell–Boltzmann distribution, further expanded in Figure 3. Note that there is no correlation at the IC time between the masses and positions, so more massive PBHs are not yet favoured to cluster closer to the cluster barycentre. Moreover, the velocity distribution computation is illustrated in Figure 5 and while Figure 6 shows the actual velocity profile.
Figure 1. Cluster IC positions and velocities. Both the position and velocities correspond to a sample of an isotropic Maxwell–Boltzmann distribution, further expanded in Figure 3. Note that there is no correlation at the IC time between the masses and positions, so more massive PBHs are not yet favoured to cluster closer to the cluster barycentre. Moreover, the velocity distribution computation is illustrated in Figure 5 and while Figure 6 shows the actual velocity profile.
Universe 07 00018 g001

3.2.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
D m m i = N m m i exp Δ m i 2 2 σ m 2 ,
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 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 1 and shown in Figure 2.
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. As the actual mean mass m ¯ 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 expression
m ¯ = exp μ m + 1 2 σ m 2 ,
while, similarly, the actual variance of the mass can be computed from the distribution parameters with
σ ¯ m 2 = exp σ m 2 1 exp 2 μ m + σ m 2 ,
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 ˜ = exp μ m ,
as well as the the modal value of the mass,
m ^ = exp μ m σ m 2 ,
as a function of the Log-Normal distribution parameters. Note that, generally speaking, a Log-Normal distribution has wide tails and as such m ¯ > m ˜ > m ^ .

3.2.2. 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
D x x i = N x exp Δ x i Σ x 1 Δ x i T 2 ,
where Δ x i = ( x i μ x ) is the PBH position displacement from the cluster centre of mass and N x = ( 2 π Σ x ) 1 is the normalisation factor. In our simulations, we consider purely spherical Σ x = σ x 2 𝟙 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
μ x = x ¯ = 2 a x 2 π ,
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
σ x 2 = 3 π 8 π a x 2 ,
while the median position of the distribution can be computed from the scale parameter with the similarly linear relation
x ˜ = 1.53 a x ,
and the modal position of the distribution is then
x ^ = 2 a x ,
as a function of the scale parameter. Note that a Maxwell–Boltzmann distribution has as well x ¯ > x ˜ > x ^ .
Finally, the probability distribution density D x x i , as a function of the mean position μ x , is then given by
D x x i = N x x i 2 exp 4 x i 2 π μ x 2 ,
where μ x is the mean position of the distribution and N x = 32 / ( π 2 μ x 3 ) 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,
D MN [ 0 3 , σ 2 𝟙 3 ] ( z ) D χ , 3 [ σ ] ( z ) = D MB [ σ ] ( z ) .
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 1 and Figure 3).

3.2.3. 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
D v v i = N v v i 2 exp 4 v i 2 π μ v 2 ,
where N v = 32 / ( π 2 μ v 3 ) is a normalisation factor. Note that the distribution scale parameter, velocity variance and modal can be similarly obtained as it was the case in the position distribution in Equations (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 4. 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
v i 2 = 12 π G μ M 5 x i 0 x i m ( x ) d x ,
where the average cumulative mass, m i , in spheres of radius x i ( 0 , r max ) from the cluster barycentre is then
m i ( x i ) = m i N O 0 x i D x x i 4 π x 2 d x .
It is more frequent, however, to express the mean velocity in terms of the Maxwell–Boltzmann scale parameter, a v , which, from Equations (18) and (19), is given by
μ v = 2 a v 2 π = 3 G m i ( x i ) 5 x i ,
Then, the initial, shell-dependent velocity distribution function can be obtained by replacing the in velocity distribution of Equation (17) the shell-dependent scale parameter using Equation (20)
It can be shown, once Equation (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 1, Figure 5 and Figure 6), which we denote by Quasi Maxwell–Boltzmann distribution (QMB).
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 ,
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; it does this by reducing the computational time and data storage needed for the analysis of the output, as individually tracking 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 also perform a suitable rotation of the frame of reference so that there is zero total and angular momentum
i L i = 0 ,
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.

3.2.4. 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 3.2 as their Schwarzschild radius. Therefore, if two bodies ( m i , r i S ) and ( m j , r j S ) with i , j ( 1 , 2 , , N O ) approach each other to less than the threshold
r min = r i S + r j S = 2 G m i + m j c 2 ,
below which, a merger is then registered by the N-body code.

3.3. 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 snapshot-by-snapshot basis considering only the cluster objects, with positions an velocities transforming as
x i CF = x i SF x i , C SF CM ,
v i CF = v i SF v i , C SF CM ,
and where the position and velocity centres of mass are computed, time by time-slice, with
x i , C SF CM = i N O δ Z i , C m i x i CF i N O δ Z i , C m i ,
v i , C SF CM = i N O δ Z i , C m i v i CF i N O δ Z i , C m i ,
where δ Z i , C = 1 if object i in the simulation remains bounded to the cluster or δ Z i , 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.

3.3.1. 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 in [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 time-run 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 point the IC, and each of the following runs taking as the starting point the last time slice of the preceding run.
(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 for each realisation, there are in total N t = 65 time slices: 64 time slices 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 2.

3.3.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 Equations (21) and (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.

3.3.3. 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 3.2, 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 5, 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 “quasi-static” (QS) phase. The algorithm by which the time of transition between the burn-in period and the quasi-static phase of the simulation is found will be detailed next.

3.3.4. Burn-in Period

One can easily observe in Figure 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 these 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 simulation mass is
m tot = N O m ¯ = N O exp ( μ m + σ m 2 / 2 ) 2.3 × 10 4 M ,
mostly contained in a spherical region with
x tot 2 x ¯ 3.2 pc .
In such a case, the mean total density of the PBH cluster at the initial time is
ρ = 3 m tot 4 π x tot 3 1.7 × 10 3 M pc 3 ,
in line with the mass density of the central regions of globular clusters, of ρ 1000 10 5 M pc 3 (see in [80]), and orders of magnitude above that of the stellar density in the solar neighbourhood, ρ 0.040 ± 0.002 M pc 3 (see in [81]).
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 6, 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 9.
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
N i = C , E ( t t BI ) β i t α i ,
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 Equation (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 fitting 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
t BI t : max [ r 2 ( t ) i ] t , j , t , j t , j 100 ,
t , j = min ( r 2 ( t ) j ) , j i + 1 ,
t , j = max ( r 2 ( t ) j ) , j i + 1 ,
the results of which for the different populations and their combinations are shown in Table 3. Note that, in all cases, t BI 5.93 × 10 6 yr .

3.3.5. 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.
Moreover, 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. Note 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.
Moreover, 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 [82]. In this scenario, only relatively infrequent massive, close PBHs could slowly inspiral into each other within a Hubble time, while emitting low-frequency gravitational radiation with regular peaks at the periastron crossing.

3.3.6. Parent Trees Computation

In each time slice we locate the possible bounded object subsystems 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 subsystem, 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 subsystems 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-subsystems 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 subsystems. They will then be assigned as their main parent the cluster most massive PBH, centrally located, unless they they show up in bounded subsystems in which a sequence of sub-parents 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 sub-systems 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 sub-system most massive object, and their parent tree tracing back to the aforementioned main parent of the cluster.
When an object belongs to a subsystem 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 subsystem, 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 Equation (25),
x i , j OF = x i CF x i , j , B CF CM ,
v i , j OF = v i CF v i , j , B CF CM ,
and where again the position and velocity centres of mass are computed, time by time slice, with
x i , j , B CF CM = i > j N O δ Z i , j , B m i x i CF i N O δ Z i , j , B m i ,
v i , j , B CF CM = i > j N O δ Z i , j , B m i v i CF i N O δ Z i , j , B m i ,
where δ Z i , j , B = 1 if the pair i , j is bounded and δ Z i , j , B = 0 if the pair i , j is unbounded within a subsystem.
From Equation (36), the semi-major axis, semi-minor axis and eccentricity can be extracted by first computing the specific angular momentum of the pair with
h i , j = μ i , j x i , j OF × v i , j OF ,
where γ = x i , j OF ( v i , j OF ) 2 2 μ i , j and μ i , j = G ( m i + m j ) is the gravitational parameter. Then, the semi-major axis and orbital eccentricity can be computed with
a = sign ( γ ) γ 1 μ i , j x i , j OF ,
e 2 = 1 + sign ( γ ) h i , j 2 a μ i , j ,
b 2 = a 2 sign ( γ ) ( 1 e 2 ) ,
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 .

3.3.7. Merger Trees Computation

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
( Δ M ) + , i = C i j ( Δ M ) , j ,
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 Equation (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:
( Δ M ) + , i < ( Δ M ) , j ,
C i , j = 0 .
(ii)
Necessary mergers:
( Δ M ) + , i > ( Δ M ) , j ,
i C i , j ( 2 , 3 , ) .
(iii)
Single one-to-one mergers:
( Δ M ) + , i = ( Δ M ) , j ,
j = j 1 ,
C i , j = 1 , C i , j j = C i i , j = 0 .
(iv)
Multiple k -to-one mergers:
( Δ M ) + , i = j ( Δ M ) , j ,
j = j 1 , , j k ,
C i , j = 1 , C i , j j = C i i , j = 0 .
After solving the aforementioned, one may proceed to solve the simplified diophantic system of Equation (43).

4. Population Statistics

The time evolution of this positions and velocities is shown in Figure 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 7 that the cluster time-trajectories exhibit a far wigglier pattern in comparison to the ejecta time-trajectories, which is as expected as the PBHs in the former move in quasi-random trajectories within the cluster while those in the latter 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 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 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 4.1, the isolated vs. bounded PBH populations of Section 4.2 and primitive vs. merged PBH populations of Section 4.3.
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 this distinction is made, 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:
N i C e i 0 < 1 ,
corresponding to instantaneously-closed orbits about the cluster barycentre, and it belongs to the Ejecta (E) for greater than unity eccentricities:
N i E e i 0 1 ,
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 object 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. Furthermore, 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 then that N R , X = N C and N A N M
N O = N R ( t ) + N A ( t ) ,
N R ( t ) = N C ( t ) + N E ( t ) , = N I ( t ) + N B ( t ) , = N P ( t ) + N M ( t ) ,
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 Equation (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 time-snapshot 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 short-lived, 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 time-evolved lines too wiggly.
(iii)
Third, that the kernel is weakest when applied to the total populations shown in the larger panels of Figure 8, as 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 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 8, a boundary effect due to the tenfold increase in the time slices time-step Δ t r across different runs r = 0 , 1 , 2 , , 7 .

4.1. Cluster and Ejecta Populations

The simulated evolution of the cluster N C ( t ) and ejecta N E ( t ) population (see Figure 8a,b 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 8 are cleaned with both the outlier-exclusion and the softening kernel algorithms mentioned previously in Section 4, as otherwise they would show excessive variance in between time-slices due the very rare outlier O ( 0.001 ) realisation in which the cluster becomes much more rapidly depopulated than usual, generally 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 Equation (31) and its results are shown in Table 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 13 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 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 N ˙ R ( t ) < 0 , as will be seen in Section 6.1.1. It will be seen later in Section 4.3 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 8a 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 6, show this trend indeed.
It will be later seen in Section 6.1.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 and 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 consistent 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 Equation (31), shown in Table 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 B I , 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 8a.
As for the ejecta, we find that the reverse is true: at the initial times t < t B I , 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
f C , 0 loss = | m C , 64 m C , 0 | m C , loss , 64 = 0.36 + 0.11 0.18 .
Conversely, the PBH ejecta population has gained throughout the simulations a fraction of objects
f E , 64 gain = | m E , 64 m C , 0 | m E , loss , 0 = 0.36 + 0.11 0.18 .
Note that both quantities being almost complementary to each other due to the lack of significant merging. Furthermore, note that, for the reasons shown in Section 6.1.1, these fractions are good approximations to the fraction of mass, as the cluster and ejecta population do not develop differentiated mass profiles.
Cluster evaporation is apparent in Figure 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, the cluster proceeds to expel PBHs at a declining rate δ N R , C N R , E < 0 , roughly given by
t δ N R , C ( t ) = t δ N R , E ( t ) = 60 + 65 145 .
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 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 7a that ejecta haloes reach at present sizes of
x ¯ E , 64 ( 7.0 × 10 3 , 3.5 × 10 5 ) pc ,
while cluster haloes reach at present sizes of
x ¯ C , 64 ( 80 , 2.3 × 10 3 ) pc .
Furthermore, shown in Figure 7a is the characteristic PBH ejecta halo typical velocities in the present to be in the range
v ¯ E , 64 ( 0.49 , 23 ) km / s ,
while the later cluster haloes reach typical velocities of
v ¯ C , 64 ( 0.030 , 0.51 ) km / s .
These 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.

4.2. 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 9a,b. 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 9 are cleaned with both the outlier-exclusion and the softening kernel algorithms mentioned previously in Section 4 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 Equation (31) whose results are shown in Table 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 time-step 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 13 much greater than the time-step at these time-slices. Figure 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 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 8.
Figure 9. Top: Total population of bounded in-cluster and in-ejecta objects in green and yellow, N B , X ( t ) where X = C , E . 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 δ N B , 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 B i 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 N R = 5000 realisations. The dotted vertical grey lines separate the N r = 7 consecutive runs, and the dotted vertical grey line signals the burn-in time. The population’s evolution is fitted to PLs in Table 7.
Figure 9. Top: Total population of bounded in-cluster and in-ejecta objects in green and yellow, N B , X ( t ) where X = C , E . 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 δ N B , 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 B i 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 N R = 5000 realisations. The dotted vertical grey lines separate the N r = 7 consecutive runs, and the dotted vertical grey line signals the burn-in time. The population’s evolution is fitted to PLs in Table 7.
Universe 07 00018 g009
Table 7. Mean bounded in-cluster and in-ejecta object population PL fits of Figure 9 from Equation (31): N P ( t r ) = β i , P t α i , P . α i , P is the PL index, again corresponding to the constant level of the differential counts, while β i , P , 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 bounded cluster and ejecta populations, respectively.
Table 7. Mean bounded in-cluster and in-ejecta object population PL fits of Figure 9 from Equation (31): N P ( t r ) = β i , P t α i , P . α i , P is the PL index, again corresponding to the constant level of the differential counts, while β i , P , 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 bounded cluster and ejecta populations, respectively.
    t BI = 2.64 × 10 5 yr Bounded in-Cluster Population Fits:Bounded in-Ejecta Population Fits:
r ( t i , t i + ) [yr]Fit: r i , B , C 2 β i , B , C α i , B , C Fit: r i , B , E 2 β i , B , E α i , B , E
3 ( t BI yr 1 , 1.30 × 10 6 ) [PL] 0.953 ( 4.4 ± 9.4 ) × 10 5 0.78 ± 0.16 [PL] 0.969 ( 1.6 ± 3.1 ) × 10 7 1.20 ± 0.15
4 1.38 × ( 10 6 , 10 7 ) [PL] 1.000 7.00 ± 0.47 0.0200 ± 0.0043 [PL] 1.000 ( 5.9 ± 1.5 ) × 10 4 0.550 ± 0.16
5 1.38 × ( 10 7 , 10 8 ) [PL] 1.000 9.80 ± 0.81 ( 1.0 ± 4.6 ) × 10 3 [PL] 0.999 0.044 ± 0.012 0.290 ± 0.015
6 1.38 × ( 10 8 , 10 9 ) [PL] 1.000 16.00 ± 0.84 0.0260 ± 0.0026 [PL] 1.000 0.520 ± 0.057 0.1600 ± 0.0053
7 1.38 × ( 10 9 , 10 10 ) [PL] 1.000 24.0 ± 1.4 0.0450 ± 0.0026 [PL] 1.000 2.20 ± 0.11 0.0920 ± 0.0023
Table 8. Mean bounded in-cluster and in-ejecta object population fractions f ˜ i , B , P and median f ˜ i , B , 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 N R = 5000 realisations at selected times.
Table 8. Mean bounded in-cluster and in-ejecta object population fractions f ˜ i , B , P and median f ˜ i , B , 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 N R = 5000 realisations at selected times.
Bounded in-Cluster and in-Ejecta Population Fraction and Median:
i t i [ yr ] f i , R , C f ˜ i , R , C f i , R , E f ˜ i , R , E i t i [ yr ] f i , R , C f ˜ i , R , C f i , R , E f ˜ i , R , E
00 0.086 + 0.032 0.030 0.085 0.0026 + 0.0040 0.0025 0.002 37 1.38 × 10 7 0.009 + 0.017 0.009 0.008 0.005 + 0.013 0.005 0.004
10 1.38 × 10 4 0.087 + 0.032 0.030 0.086 0.0027 + 0.004 0.0026 0.002 46 1.38 × 10 8 0.010 + 0.017 0.009 0.008 0.010 + 0.019 0.009 0.008
19 1.38 × 10 5 0.123 + 0.051 0.052 0.121 0.008 + 0.010 0.008 0.007 55 1.38 × 10 9 0.009 + 0.016 0.008 0.008 0.014 + 0.024 0.011 0.0012
28 1.38 × 10 6 0.008 + 0.010 0.007 0.007 0.0012 + 0.0044 0.0012 0.001 64 1.38 × 10 10 0.008 + 0.014 0.007 0.007 0.018 + 0.024 0.012 0.0016
Initially, the vast majority of the binaries are bound to be contained within the cluster, which contains a fraction of bounded objects of
f C , B , 0 = 0.086 + 0.0032 0.030 ,
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
f C , E , 0 = 0.0026 + 0.0040 0.0025 ,
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 1 to a cuspier one as will later be seen in Figure 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 7a. As the PBHs fall to this dense environment they quickly acquire speeds as seen in Figure 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 time slice 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
f C , B , 28 = N C , B , 2 N R , 0 = 0.008 + 0.010 0.007 ,
at t = 1.38 × 10 6 yr to a present relative maximum with respect to the total population of:
f C , B , 64 = N C , B , 64 N R , 0 = 0.008 + 0.014 0.007 ,
at t = 1.38 × 10 10 yr , as seen in Figure 9a. Note that, as the cluster itself loses 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 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
f E , B , 28 = N E , B , 28 N R , 0 = 0.0012 + 0.0044 0.0012 ,
at t = 1.38 × 10 6 yr to a present relative maximum of
f E , B , 64 = N E , B , 64 N R , 0 = 0.018 + 0.014 0.012 ,
at t = 1.38 × 10 10 yr , as seen in Figure 9b.
Last, a more detailed analysis of binary pairs is given in Section 5.1.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
f B , 28 = N B , 28 N R , 0 = 0.010 + 0.010 0.007 ,
after recombination at t = 1.38 × 10 6 yr and
f B , 64 = N B , 64 N R , 0 = 0.025 + 0.028 0.014 ,
at the present time t = 1.38 × 10 10 yr , with f B = f B , C + f B , E , as previously explained in Section 4.2.
However, we have found as well a rich structure of subsequent binary subsystems within multiple PBH systems. It has been found using the methods from Section 3.3.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 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 (2-component) 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) systems, in which the least massive PBH ( 2 nd -level) orbits the intermediate PBH ( 1 st -level), which itself orbits the most massive PBH ( 0 th -level), amounting to a second level in the hierarchy.
(iii)
Another 0.11–0.01% of bounded objects are present in quaternary (4-component) systems, adding an extra third level in the hierarchy.
(iv)
Last, 0.001–0.0001% of bounded objects are part of quinary 5-component) systems, adding a final forth level in the parent tree.
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 9, the population in each level is less stable over time the deeper the level in the hierarchy.

4.3. Primitive and 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 found in our simulation ( N R ) M = 117 such mergers in N R = 5000 realisations, leading to a total fraction of merged objects of f ( N R ) M = 2.34 × 10 5 . Out of these,
(i)
( N R , C , I ) M = 50 such mergers occur in cluster, isolated objects, leading to a subfraction of merged objects throughout all simulations at the last time slice of
f ( N R , C , I ) M = 1.00 × 10 5 .
(ii)
( N R , C , B ) M = 4 such mergers occur in cluster, bounded objects leading to a total subfraction of merged objects of
f ( N R , C , B ) M = 8.00 × 10 7 .
(iii)
( N R , E , I ) M = 1 such mergers occur in ejecta, isolated objects leading to a total subfraction of merged objects of
f ( N R , E , I ) M = 2.00 × 10 7 .
(iv)
Last, ( N R , E , B ) M = 62 such mergers occur in ejecta, bounded objects leading to a total subfraction of merged objects of
f ( N R , E , B ) M = 1.24 × 10 5 .
These results, summarised in Table 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 11, 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 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
( N R , C , B ) M ( N R , C , I ) M = 0.080 > 0.025 = f C , B f C , I f C , B 1 f C , B .
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; therefore, 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.
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 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 5.2.1, where we compute the mass profile of merging PBHs and their population characteristics.

4.3.1. Simulation Merger Trees

The identification of merger pairs in order to construct the merger trees has been done by the procedure described in Section 3.3.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 ) M 1 , 1 = 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 ) M 2 , 1 = 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 3.3.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.

5. Close Encounters

In this section, we take a closer look to the identified binary and merger pairs, in Section 5.1 and Section 5.2, respectively. In particular, we study the collisional and short-distance interaction dynamics of PBHs in these pairs, and extract the mass distributions of the PBHs participating in such systems. Moreover, we offer a lower bound prediction of the rates of occurrence of PBH mergers and of the PBH hyperbolic encounters per cluster.

5.1. Binary Pairs

Shown in Figure 10, we have computed the binary pair number and frequency in the simulations along with the total mass, chirp mass and mass ratio profile of the pair, as well as the component mass distributions. We have as well extracted the mean, median and modal values of the total mass, chirp mass and mass ratio profiles, which are shown in Table 12.
Table 11. Distribution mean and median values of the distributions in Figure 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) t FM = 1.38 × 10 9 yr up to the present, for a total of i = 1 N R ( N M ) 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 and M stand for the cluster, ejecta, isolated, bounded, remaining and merged populations, respectively.
Table 11. Distribution mean and median values of the distributions in Figure 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) t FM = 1.38 × 10 9 yr up to the present, for a total of i = 1 N R ( N M ) 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 and M stand for the cluster, ejecta, isolated, bounded, remaining and merged populations, respectively.
Merger Populations:Merger Distribution Parameters:
Pop.Frequency:Total Mass:Chirp Mass:Mass Ratio:
P i P j N R , M P i , P j f R , M i , j m ¯ T [ M ] m ˜ T [ M ] m ¯ C [ M ] m ˜ C [ M ] q ¯ q ˜
117 2.34 × 10 5 16701510642567 0.965 0.545
C I 50 1.00 × 10 5 21201690793667 0.517 0.517
C B 4 8.00 × 10 7 967925368368 0.634 0.641
E I 1 2.00 × 10 7 15901590689689 0.867 0.867
E B 62 1.24 × 10 5 13501230538477 0.594 0.561
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 in-ejecta PBHs, as well as the later evolution, as we had found in Section 4.1 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 4.1.
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 3.2.1, the PBH average and mean masses in the IC are m ¯ = 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 in-ejecta bounded populations up until the present time.
However, for binary pairs, the mean and median total masses in the pair are m ¯ T = 303 M and m ˜ 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 of q ¯ T = 0.154 and q ˜ 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
( m ¯ 1 , m ¯ 2 ) B = ( 261 , 40.9 ) M ,
( m ˜ 1 , m ˜ 2 ) B = ( 172 , 9.46 ) M ,
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,
( m ¯ 1 , m ¯ 2 ) B , C = ( 254 , 50.8 ) M ,
( m ˜ 1 , m ˜ 2 ) B , C = ( 177 , 7.90 ) M .
(ii)
For the bounded in-ejecta PBH,
( m ¯ 1 , m ¯ 2 ) B , E = ( 264 , 107 ) M ,
( m ˜ 1 , m ˜ 2 ) B , E = ( 186 , 13.5 ) M .
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.

5.1.1. 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 these 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, because, as seen in Table 13, the time-step in each run 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
Γ r S , P = N r S N r 1 S Δ t r ,
where N r S 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
Γ 1 S ( 0 t i 10 ) = ( 2.0 + 1.6 1.4 ) × 10 6 Gyr 1 ,
Γ 2 S ( 11 t i 19 ) = ( 2.2 + 1.2 1.8 ) × 10 6 Gyr 1 ,
Γ 3 S ( 20 t i 28 ) = ( 2.2 + 2.3 1.7 ) × 10 5 Gyr 1 ,
Γ 4 S ( 29 t i 37 ) = ( 0.9 + 2.9 0.6 ) × 10 4 Gyr 1 ,
Γ 5 S ( 38 t i 46 ) = ( 1.0 + 3.5 0.6 ) × 10 4 Gyr 1 ,
Γ 6 S ( 47 t i 55 ) = 98 + 400 63 Gyr 1 ,
Γ 7 S ( 56 t i 64 ) = 9 + 95 6 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. Furthermore, note that these rates have been computed for a single cluster. Latter, in Section 8 we will compute this quantity for a co-moving cosmological volume of 1 Gpc 3 , and thus provide a more meaningful value for the slingshot rate of PBHs.

5.2. Merger Pairs

Similar to what we have done for the binary pairs of Section 5.1, we have computed the merger pair number and frequency in the simulations along with the total mass, chirp mass and mass ratio profile of the pair, as well as the component mass distributions, shown in Figure 11. The mean, median and modal values of these distributions are shown in Figure 11.
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 4.3. 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 3.2.1 the PBH average and mean masses in the IC are m ¯ = 21.9 M and m ˜ = 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, with m ¯ T = 1.67 × 10 3 M and m ˜ T = 1.51 × 10 3 M , which is 55–200 times greater 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 of q ¯ T = 0.965 and q ˜ T = 0.545 , indicating that the less massive partner in the merger is typically bigger than half the mass of the more massive partner, 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 m 1 m 2 , two PBHs in a merged pair, according to the previous results would have masses close to
( m ¯ 1 , m ¯ 2 ) M = ( 850 , 820 ) M ,
( m ˜ 1 , m ˜ 2 ) M = ( 977 , 533 ) M ,
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 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 that
(i)
For the merged, isolated and in-cluster PBHs,
( m ¯ 1 , m ¯ 2 ) M , I , C = ( 1400 , 712 ) M ,
( m ˜ 1 , m ˜ 2 ) M , I , C = ( 1110 , 575 ) M .
(ii)
For the merged, bounded and in-ejecta PBH,
( m ¯ 1 , m ¯ 2 ) M , B , E = ( 847 , 503 ) M ,
( m ˜ 1 , m ˜ 2 ) M , B , E = ( 788 , 442 ) M .
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 6.2. Furthermore, we do not compute the cases of the merged, bounded, in-cluster PBHs and merged, isolated, in-ejecta PBHs as with only four and one cases, respectively, we are far from having enough sampling to give meaningful results.

5.2.1. Merger Event Rates

We have defined a merger between two bodies in the simulations in Section 4.3.1, and Table 10 and Table 11 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 as 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
Γ r M , P = ( l M 1 ) N r M , P N r 1 M , P Δ t r ,
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 10 we find that
Γ 6 M ( 46 t i 55 ) = 6.44 × 10 4 Gyr 1 ,
Γ 7 M ( 56 t i 64 ) = 1.82 × 10 3 Gyr 1 .
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
Γ 6 M , C , I ( 46 t i 55 ) = 1.61 × 10 4 Gyr 1 ,
Γ 7 M , C , I ( 56 t i 64 ) = 8.05 × 10 4 Gyr 1 .
(ii)
For all 66 mergers occurring originating in previously bounded systems, and out of which 3 occur in the sixth time-run and 63 do occur in the seventh time-run, we find that
Γ 6 M , E , B ( 46 t i 55 ) = 4.83 × 10 4 Gyr 1 ,
Γ 7 M , E , B ( 56 t i 64 ) = 1.01 × 10 3 Gyr 1 .
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.
Furthermore, these are the merger event rates for a non-relativistic evolution for a single cluster. In Section 8, we will compute this quantity for a co-moving 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.

6. Dynamical Distributions

Now, we proceed to extract in Section 6.1 the density, mass, position and velocity distributions of PBHs evolution, slice by slice in time and segregated by the subpopulation type. We do as well quantify the cluster evaporation rate in Section 6.2.1 and the degree of both mass segregation and dynamical friction within it in Section 6.2.2.

6.1. Dynamical Profiles

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

6.1.1. 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 14 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 2.
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 13a, and acquire in the process large kinetic energies, as shown in Figure 13b. Then, due to the large number of encounters 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 mass-depleted 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 μ C m and σ C m increasing only by
Δ μ C m = μ C , 64 m μ C , 0 m μ C , 0 m = 2.8 % ,
Δ σ C m = μ C , 64 m μ C , 0 m μ C , 0 m = 1.0 % ,
in the case of the in-cluster PBHs between the initial t 0 = 0 yr and final simulation times t 64 = 1.38 × 10 10 yr , and μ E m and σ E m conversely decreasing just by
Δ μ C m = μ C , 64 m μ C , 46 m μ C , 46 m = 0.20 % ,
Δ σ C m = μ C , 64 m μ C , 46 m μ C , 46 m = 0.13 % ,
in the case of the in-ejecta PBHs between the beginning of the t 46 = 1.38 × 10 8 yr and the final simulation times t 64 = 1.38 × 10 10 yr . Note that, the difference is very small for both parameters and populations, of O ( 0.01 ) in the case of mass expected values, and of O ( 0.001 ) in the case of the standard deviation of the logarithmic mass at the last time-slice, t 64 = 1.38 × 10 10 yr , consistent with a small degree of mass segregation between the cluster and ejecta populations.

6.1.2. Density Profile Evolution

Given in Figure 12a,b are the PBH volume-factored density profiles x 3 ρ i ( x ) | t , with i = C , E being the cluster and ejecta populations, computed at the five different time slices that bound the last four time-runs plus the initial time.
The PL fits of the density profiles are shown in Table 14, 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-of-influence 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 C max ( t ) where the cluster is sparse enough, that effectively can be considered as the cluster radius.
The second concerns the ejecta objects. The volume-factored 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 high-velocity PBHs. Last, the density profile arrives as well at each time-slice to a time-increasing cut-off length scale x E max ( 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:
ρ C ( x ) | t x 2.8 , t > t BI , x < x C max ( t ) .
(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:
ρ E ( x ) | t x 2.6 , t > t BI , x < x C max ( t ) .
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 C max 10 4 pc and x E max 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:
ρ C ( x ) | t 0 ( 10 6 , 10 6 ) M pc 3 ,
x ( 0.001 , 10 ) pc ,
t 0 = 0 yr .
(ii)
In the case of the last time-slice of cluster objects:
ρ C ( x ) | t 64 ( 10 6 , 10 22 ) M pc 3 ,
x ( 10 , 10 4 ) pc ,
t 64 = 1.38 × 10 10 yr .
(iii)
In the case of the first sufficiently populated time-slice of ejecta objects, with N E ( t ) > 100 :
ρ E ( x ) | t 28 ( 10 5 , 10 7 ) M pc 3 ,
x ( 0.001 , 10 ) pc ,
t 28 = 1.38 × 10 6 yr .
(iv)
In the case of the final time-slice of ejecta objects:
ρ E ( x ) | t 64 ( 10 8 , 10 32 ) M pc 3 ,
x ( 10 , 10 7 ) pc ,
t 64 = 1.38 × 10 10 yr .
It is worth noting, however, that cluster layers with densities smaller than ρ = 0.08 M pc 3 (see in [81]) 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.

6.1.3. 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 12b,c 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 14, 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.
Second, as the cluster puffs-up, its characteristic length scale increases greatly as described in Equation (164), with its radius expanding by three orders of magnitude, as seen for the last four time slices in Figure 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 C max t 64 = 1.38 × 10 10 yr and nearly unit eccentricities e C | x x C max 1 .
Cluster positions are constrained to be in the range
0.02 pc x C ( t ) 2 × 10 3 pc t > t 0 ,
x C ( 0 yr ) ( 0.02 , 20 ) pc ,
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
x C max ( t ) β x , C t α x , C pc , t > t BI ,
β x , C = ( 2.96 ± 0.56 ) × 10 4 pc ,
α x , C = 0.639 ± 0.010 .
Third, that the stability of these PBHs in orbits bound to the cluster is almost guaranteed in our simulations, as 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 8.1.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, as the ejecta sphere-of-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 Equation (174).
Ejecta positions are constrained to be in the range
0.4 pc x E ( t ) 10 5 pc t > t 0 ,
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
x E max ( t ) β x , E t α x , E pc , t > t BI ,
β x , E = ( 3.70 ± 0.80 ) × 10 5 pc ,
α x , E = 0.917 ± 0.012 .
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 7b. This greatly increases the computational speed of the simulations at late times, as the characteristic discrete dynamical time-step in the N-body solver algorithm increases accordingly.
However, it will be seen in Section 8.1.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.

6.1.4. 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 12e,f, 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 14, when such fits are appropriate, although, like it was the case with the position profiles computation of Section 6.1.3, other fits to Maxwell–Boltzmann and Rayleigh distributions have been tried and proved to worse describe the data.
From Figure 12e,f 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
0.01 km / s v E ( t ) 10 km / s t > t 0 ,
v C ( 0 yr ) ( 0.2 , 3 ) km / s ,
v C ( 1.38 × 10 10 yr ) ( 0.01 , 0.3 ) km / s ,
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 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 6.1.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
v C max ( t ) β v , C t α v , C pc , t > t BI ,
β v , C = ( 383 ± 259 ) × 10 4 km / s ,
α v , C = 0.276 ± 0.041 ,
Ejecta velocities are even more constrained than cluster velocities to a range of
0.08 km / s v E ( t ) 20 km / s t > t 0 ,
v E ( 1.38 × 10 6 yr ) ( 2 , 20 ) pc ,
v E ( 1.38 × 10 10 yr ) ( 0.08 , 8 ) × 10 3 pc ,
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 time-dependent and is given by
v E max ( t ) β v , E t α v , E pc , t > t BI ,
β v , E = ( 29.6 ± 4.8 ) × 10 4 km / s ,
α v , E = 0.0826 ± 0.0094 .
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
V C esc ( X C , t ) = 2 G N M ( x < X C ) X C 1 / 2 ,
given the symmetries of out clusters, and where X C ( t ) is the time-dependent cluster radius and
M ( x < X C ) = i = 1 N O m i θ ( x < X ) ,
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.

6.2. Mass Segregation and 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.

6.2.1. 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. As 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 [83] to approximately be
t rel ( t ) = N R , C ( t ) t cross ( t ) 8 log N R , C ( t ) ,
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, Equation (149) yields for the IC a relaxation time of
t C , 0 rel = ( 1.323 ± 0.071 ) × 10 7 yr ,
as 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 19 yields an initial crossing time of
t C , 0 cross = X C , 0 V C , 0 = ( 7.31 ± 0.41 ) × 10 5 yr ,
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 because, as Table 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
t C , 64 cross = X C , 64 V C , 64 = ( 1.72 ± 0.92 ) × 10 8 yr .
This increase, on top of the population loss from cluster evaporation, eventually increases the relaxation time by two orders of magnitude by the time of the last time-slice
t C rel ( t 64 ) = ( 1.32 + 2.00 2.10 ) × 10 9 yr .
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 13a and give the best-fit values in Table 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 6.1.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 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.
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 time-slice, with a PBH selection likelihood proportional to its mass to better capture the most massive of these. The largest of these PBHs are visible then in Figure 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.

6.2.2. 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-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, as 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, as 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 13b and whose best-fit values we give in Table 16, in which a number of features similar to those in the case for mass segregation of Section 6.2.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, 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 6.1.1.
(iv)
Last, that Figure 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 leftward-moving velocity profile of Figure 13a, particularly so for the 1 σ area and less so for the 2 σ , but still clear if excluding the largest PBHs.
Figure 13. Shown are the 1 σ and 2 σ contours for the cluster PBHs mass vs. position and velocity distributions at selected times. Fits for the mass, position, velocity and density probability distributions are given in Table 16. Black dots show the scatter a 1000 sample of individual PBHs in each time slice, randomly selected in the case of the bare distribution and selected with a probability proportional to their mass in the case of the mass weighted in order to better illustrate mass segregation and dynamical friction, as otherwise such very massive PBHs are easy to miss in the sample given that there are not more than O ( 1 ) m i > 1000 M massive objects per cluster. The Best-Fit (BF) values of the distributions, marked by the large coloured dots, are given in Table 16.
Figure 13. Shown are the 1 σ and 2 σ contours for the cluster PBHs mass vs. position and velocity distributions at selected times. Fits for the mass, position, velocity and density probability distributions are given in Table 16. Black dots show the scatter a 1000 sample of individual PBHs in each time slice, randomly selected in the case of the bare distribution and selected with a probability proportional to their mass in the case of the mass weighted in order to better illustrate mass segregation and dynamical friction, as otherwise such very massive PBHs are easy to miss in the sample given that there are not more than O ( 1 ) m i > 1000 M massive objects per cluster. The Best-Fit (BF) values of the distributions, marked by the large coloured dots, are given in Table 16.
Universe 07 00018 g013
Note, however, that, in contrast to the previous argument that gravitational drag is strongest for massive objects, Figure 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 for 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, as 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.
Table 16. Cluster mass segregation and dynamical friction from Figure 13. BF values for the mass segregation ρ i j ( x C , m C , t ) and dynamical friction ρ i j ( x C , m C , t ) plots, where j = w , b for the mass-weighed and bare (homogeneously weighted) distributions, respectively, from which we extract the BF values at selected times.
Table 16. Cluster mass segregation and dynamical friction from Figure 13. BF values for the mass segregation ρ i j ( x C , m C , t ) and dynamical friction ρ i j ( x C , m C , t ) plots, where j = w , b for the mass-weighed and bare (homogeneously weighted) distributions, respectively, from which we extract the BF values at selected times.
      Cluster Mass Segregation:Cluster Dynamical Friction:
i t i [ yr ] x i , C ( w )   [ pc ] m i , C ( w ) [ M ] x i , C ( b )   [ pc ] m i , C ( b ) [ M ] v i , C ( w )   [ km / s ] m i , C ( w ) [ M ] v i , C ( b ) [ km / s ] m i , C ( b ) [ M ]
00[BF] 1.35 125 1.82 6.62 [BF] 1.60 49.9 2.24 6.62
28 13.8 × 10 6 [BF] 1.07 49.3 1.35 6.62 [BF] 3.76 50.5 4.54 6.62
37 13.8 × 10 7 [BF] 4.93 77.3 6.33 6.62 [BF] 2.64 64.1 2.30 6.62
46 13.8 × 10 8 [BF] 27.2 103 24.4 6.62 [BF] 1.07 40.8 0.839 6.62
55 13.8 × 10 9 [BF] 86.1 57.2 107 6.87 [BF] 0.240 47.6 0.367 6.87
64 13.8 × 10 10 [BF]354 65.9 493 6.87 [BF] 0.101 41.3 0.166 6.87

7. Orbital Distributions

Here, in Section 7.1 and Section 7.2, 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 semi-major axis of PBHs. Furthermore, in Section 5.1.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.

7.1. Cluster Orbital Properties

In Figure 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. Furthermore, the best-fit values of these distributions are given in Table 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 ) | t BF , 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 | t BF 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 | t BF = 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 .

7.1.1. Custer Binaries Orbital Properties

The complementary population of bounded cluster PBHs orbital parameter distributions is given in Figure 14c along with their best-fit values in Table 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 do 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 ,
and one vertical of constant eccentricity with
0.916 e C , B | t BF 1.000 ,
asymptotically approaching the parabolic value over time, extending from the best-fit values a C , I | t BF .
There is, however, a difference in one particular respect. As seen in Figure 14c and compared to Figure 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 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.

7.1.2. Cluster Transient Binaries

Figure 15a shows the lifetimes of cluster bounded subsystems, the majority of which coming as binary pairs as detailed in Section 4.2.
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 subsystems parent object and binary, ternary, quaternary and quintic pairs. Note that a large fraction of these subsystems 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 time slice. This suggests again that the total number of binary pairs is underestimated by a non-negligible fraction of the aforementioned N C , B total as 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 subsystems present in the IC are disrupted during the first four time-runs, 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 because, 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 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 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 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.

7.2. Ejecta Orbital Properties

In Figure 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 17 we give the best-fit values of these distributions. As 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 best-fit value at ( a E , I , e E , I ) | t BF 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-major-axes. In this case, it is found that eccentricity negatively correlates 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 | t BF = 1 and saturated at 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.

7.2.1. Ejecta Binaries Orbital Properties

The complementary orbital parameter distributions of for the population of bounded ejecta PBHs is given are Figure 14d along with their best-fit values in Table 17.
In this case, the ejecta, bounded PBHs distribution shape and position in ( a E , B , e E , B ) of the orbital parameter space 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 | t BF . 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
1.00 e E , B | t BF 1.16 ,
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 14d, and compared to Figure 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 the orbital parameter space never exceeds 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.

7.2.2. Ejecta Stable Binaries

Figure 15b shows the lifetimes of ejecta bounded sub-systems, the majority of which as simple binary pairs as detailed in Section 4.2.
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 subsystems 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 subsystems 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 subsystems 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 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 .

7.3. 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 saw in Section 5.1.1 that a large majority of these hyperbolic encounters cannot be identified by this method as 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, as 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 i , j rel = x j x i and relative velocities v i , j rel = 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
e i , j rel = v i , j rel × h i , j rel G N ( m i + m j ) x i , j rel | x i , j rel | ,
where h i , j rel = x i , j rel × v i , j rel 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
e i , j rel = | h i , j rel | 2 G N ( m i + m j ) | 1 e i , j rel | 2 ,
which corresponds to the impact parameter of an hyperbolic encounter when the relative eccentricity exceeds the parabolic limit, for e i , j rel > 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 i , j rel = 100 N R N O = 10 8 pairs of PBHs per time slice in all considered time slices but the IC, in Figure 16a for the relative eccentricities and Figure 16b for the relative impact parameters. We also fit the resulting curves to Log-Normal distributions in Table 18, having tested prior to that other distributions and found that a Log-Normal distribution describes best the data.
We find in Figure 16a that, for all the considered time-slices, there is a significant and roughly time-independent fraction of i , j pairs that may potentially produce hyperbolic encounters with e i , j rel > 1 . Moreover, this encounters may be produced with eccentricities as large as large as e i , j rel = 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 grows 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 i , j rel 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 eccentricity range of e i , j rel O ( 0.01 ) .
We find as well in Figure 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 i , j rel 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
0.001 pc b i , j rel ( t 0 ) 1 pc ,
to the present time where
0.1 pc b i , j rel ( t 0 ) 1000 pc ,
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 i , j rel O ( 0.001 ) pc at the initial time, and of b i , j rel O ( 0.1 ) pc at the present time.
Considering, however, that there are in total N i , j rel = 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 Section 6.1.2, Section 6.1.3 and Section 6.1.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.

8. 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 8.1 will deal with predictions for the DM background of cluster and ejecta objects separately for a typically sized galactic halo, while in Section 8.2 we will compute the distinct hyperbolic encounter and merger event rates of PBHs per co-moving volume.

8.1. DM Energy Density

It was mentioned in Section 6.1.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, as both the cluster and ejecta sphere-of-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 as the periods of these objects are very high.
However, in a realistic scenario this is not the case any more. 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 [84] is estimated to be
M MW = M MW ( X < X MW = 100 kpc ) = ( 5.3 + 2.1 1.2 ) × 10 11 M ,
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 .

8.1.1. 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
M C tot ( t ) = i = 1 N O δ i , C ( t ) m i ( t ) N O m ¯ C ( t ) f C ( t ) N O exp μ m + σ m 2 / 2 ( 1 f C loss ( t ) ) ,
given that for our Log-Normal mass distribution the actual mass mean is m ¯ = exp μ m + σ m 2 / 2 . The mass is itself contained in a sphere of a mean effective radius given by
X C tot ( t ) = x ¯ C = exp μ x + σ x 2 / 2 ,
where the input logarithmic mean position, μ x and and deviation σ x can be used again to get the actual position mean, x ¯ = exp μ x + σ x 2 / 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 14, while the cluster characteristic scales m ¯ , x ¯ and v ¯ are shown in Table 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 the work in [84], found to be
N C tot = f DM f PBH M MW M C tot ( t 0 ) ,
where we have not corrected for the possible loss of ejecta PBHs out of the sphere as the fraction of M C tot ( 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
N C ind ( t ) = N O f DM f PBH M MW M C ind ( t ) .
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 C tot ( t ) ( V MW / N C tot ( t ) ) 1 / 3 , where V MW stands for the Milky Way volume contained in the X < 100 kpc sphere, or after some substitutions
X C tot X MW 3 N C tot 1 / 3 ,
and correspondingly, the mean inter-PBH distance within said sphere can be well approximated by
X C ind ( t ) X MW 3 N C ind ( t ) 1 / 3 ,
assuming that the cluster have puffed-up enough to start overlapping with each other. Note that the mean inter-PBH distance X C ind ( 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 C tot ( t ) of Equation (163), while all the PBH clusters in the galactic halo will total an also monotonically decreasing joint mass of
M C ind ( t ) N C tot M C tot ( t ) .
Last, the mean occupied PBH cluster volume fraction of said sphere is then F C tot ( t ) = N C tot ( t ) V C tot ( t ) / V MW , where V C tot ( t ) stands for the cluster volume contained in the cluster effective cluster radius X C = 100 kpc , which leads to
F C tot ( t ) N C tot X C ( t ) X MW 3 .
Assuming a DM fraction within the Milky Way halo sphere of f DM = 0.8 entirely due to PBHs, f PBH = 1 then Equations (165), (167) and (170) yield a total cluster number per sphere, N C tot of
N C tot = ( 2.7 + 1.7 0.9 ) × 10 7 ,
separated on average by a distance, X C tot of
X C tot = 800 + 190 100 pc ,
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 C ind ( t ) and mean individual PBH distance X C ind ( t ) are given in Table 19 for the cluster population, following the aforementioned procedure along with the mean mass per cluster bubble M C tot ( t ) , mean in-cluster mass per halo M C ind ( t ) , and cluster occupied halo volume fraction F C tot ( 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 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, as F C , 64 tot ( 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 C ind = 112 + 20 23 pc and of the same order of magnitude than the cluster-to-cluster distance.

8.1.2. Ejecta DM Energy Density

We now aim to repeat the same analysis of Section 8.1.1 for ejecta objects, and constrain the abundance and compute the mean inter-ejecta spheres-of-influence 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 as the number of ejecta objects is negligible at the time of the IC. Then, the ejecta will typically have a mean total mass of
M E tot ( t ) = i = 1 N O δ i , E ( t ) m i ( t ) N O m ¯ E ( t ) f E ( t ) N O ( exp μ m + σ m 2 / 2 ) f C gain ( t ) ,
given that, as we had seen in the previous section, the PBH mean mass is m ¯ = exp μ m + σ m 2 / 2 . The bulk of the mass is then contained in a sphere of a mean effective radius given by
X E tot ( t ) = x ¯ E = exp μ x + σ x 2 / 2 ,
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 + σ x 2 / 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 14, while the ejecta characteristic scales m ¯ , x ¯ and v ¯ and are shown in Table 15. As 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
N E tot = N C tot .
Correspondingly, the number of these in-ejecta PBHs within the X < X MW sphere centred at the Milky Way core is then
N E ind ( t ) = N O f DM f PBH M MW M E ind ( t ) .
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
X E tot = X C tot ,
and correspondingly, the mean inter-PBH distance within said sphere can be well approximated by
X E ind ( t ) X MW 3 N E ind ( t ) 1 / 3 ,
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 E tot ( t ) of Equation (173), while all the PBH ejecta in the galactic halo will total as well an also monotonically increasing joint mass of
M E ind ( t ) N E tot M E tot ( t ) .
The mean occupied PBH ejecta volume fraction of the ejecta sphere-of-influence will be then
F E tot ( t ) N E tot X E ( t ) X MW 3 .
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 Equations (175), (177) and (180) a total ejecta number per sphere, N C tot of
N E tot = N C tot = ( 2.7 + 1.7 0.9 ) × 10 7 ,
separated on average by a distance, X C tot of
X E tot = X C tot = 800 + 190 100 pc .
The estimated values of the mean individual PBH number N E ind ( t ) and mean individual PBH distance X E ind ( t ) are given in Table 19 for the ejecta population, along with the mean mass per ejecta bubble M E tot ( t ) , mean in-ejecta mass per halo M E ind ( t ) and ejecta occupied halo volume fraction F E tot ( 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 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 E , 64 tot 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 halo-sized 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
X E ind = 108 + 29 17 pc ,
at t 46 = 1.38 × 10 8 yr to a present value of
X E ind = 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 these scales that is clearly not the case any more. 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.

8.2. Gravitational Event Rates

We aim now to compute a minimum threshold for the average hyperbolic encounter rate and merger event rate in time per co-moving volume, starting from the computed rates for both per cluster in Section 5.1.1 and Section 5.2.1.
In order to do so, we begin by computing the total DM mass contained in a co-moving 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 Equation (163), which reduces to M C , 0 tot = N O m ¯ C ( t 0 ) thus obtaining
M C , 0 tot = ( 2.991 ± 0.074 ) × 10 4 M ,
for a total number of clusters per co-moving volume of N C tot = ρ DM / M C , 0 tot , which yields
N C , t = ( 1.325 ± 0.040 ) × 10 15 Gpc 3 ,
which is a large number indeed, albeit an expected one since it is known that a 1 Gpc 3 co-moving box contains about O ( 10 7 ) Milky Way-sized galactic haloes.

8.2.1. Co-Moving Hyperbolic Encounter Rates

We do now take the hyperbolic encounter rates of Section 5.1.1 and scale them to a 1 Gpc 3 -sized box with the transformation Γ r S N C tot Γ r S . The final co-moving slingshot rates are then
Γ 1 S ( 0 t i 10 ) = ( 2.7 + 7.6 2.1 ) × 10 12 yr 1 Gpc 3 ,
Γ 2 S ( 11 t i 19 ) = ( 3.0 + 2.5 1.6 ) × 10 12 yr 1 Gpc 3 ,
Γ 3 S ( 20 t i 28 ) = ( 2.9 + 3.9 2.6 ) × 10 11 yr 1 Gpc 3 ,
Γ 4 S ( 29 t i 37 ) = ( 1.2 + 4.8 0.9 ) × 10 10 yr 1 Gpc 3 ,
Γ 5 S ( 38 t i 46 ) = ( 1.4 + 5.2 1.1 ) × 10 9 yr 1 Gpc 3 ,
Γ 6 S ( 47 t i 55 ) = ( 1.3 + 5.5 1.0 ) × 10 8 yr 1 Gpc 3 ,
Γ 7 S ( 56 t i 64 ) = ( 1.2 + 5.9 0.9 ) × 10 7 yr 1 Gpc 3 .
These are indeed very large rates; however, the same caveats that applied to the slingshot rates of Section 5.1.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.

8.2.2. Co-Moving Merger Event Rates

We do at last take the merger event rates of Section 5.2.1 and scale them again to a 1 Gpc 3 -sized box with the transformation Γ r M N C tot Γ r M . The final co-moving merger event rates are then
Γ 6 M ( 46 t i 55 ) = 852 ± 26 yr 1 Gpc 3 ,
Γ 7 M ( 56 t i 64 ) = 2409 ± 74 yr 1 Gpc 3 .
The merger event rates segregated by their origin, would be, for all 51 mergers occurring from previously isolated PBHs
Γ 6 M , C , I ( 46 t i 55 ) = 213 ± 11 yr 1 Gpc 3 ,
Γ 7 M , C , I ( 56 t i 64 ) = 1066 ± 33 yr 1 Gpc 3 ,
while for all 66 mergers occurring originating in previously bounded systems, we find that
Γ 6 M , E , B ( 46 t i 55 ) = 639 ± 20 yr 1 Gpc 3 ,
Γ 7 M , E , B ( 56 t i 64 ) = 1337 ± 41 yr 1 Gpc 3 .
Note, however, that the same caveats that applied to the merger event rates of Section 5.2.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 co-moving volume, there is no tension between the total co-moving merger rate from all population sources of PBHs found at present, Γ 7 M ( 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 in [85]), as 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 in [86]).

9. Conclusions

Understanding the nature of Dark Matter remains one of the key questions in theoretical physics and cosmology and has resulted in a plethora of models ranging from weakly interacting massive particles (WIMPs), sterile neutrinos, axions, modifications of gravity like MOND but also PBHs, see in [87,88] for recent reviews. Furthermore, an interesting possibility is to have a mixed f ( R ) gravity model with axion, which could make it easier to be detected as it affects both the inflation and DE era [89,90].
Here, we focus on the PBHs as it also provides a model that connects the early time and late time physics and remains a plausible candidate for DM. In this case, PBHs not only play the role of DM, but they may also influence DE at late times, as by evaporating they inject energy and behave as a time-dependent DE equation of state, see in [91]. PBH may also induce non-adiabatic perturbations on DE, see in [92].
In this work, however, we focus on the dynamics of PBH in clusters, as it is one of the cornerstones in our understanding of structure formation at small scales, as well 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.
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 one-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 sub-haloes 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 of q ¯ B = 0.154 and q ˜ B = 0.0552 , and total masses of m ¯ 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 of q ¯ M = 0.965 and q ˜ M = 0.545 . Merger total masses masses found are m ¯ T , M = 1670 M and m ˜ T , M = 1510 M , while found chirp masses are m ¯ c , M = 642 M and m ˜ c , M = 567 M . Furthermore, 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. The fact that once the many different ejecta bubbles merge the ejecta background has length scales comparable to medium sized galactic haloes is also interesting. 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 semi-major 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 2.
Our work may also be extended in many other ways. For example, an important limitation in our current analysis is that we did not include GW emission from the PBHs as the AstroGrav code is purely Newtonian. This maybe resolved in the future by extending the equations of motion to higher order in the post-Newtonian expansion. Another possibility is to do simulations with a varying and higher number of PBHs in the simulation, so as to study more realistic clusters of different sizes and more varying in mass.

Author Contributions

Conceptualization, J.G.-B.; methodology, J.G.-B.; software, M.T.; validation, J.G.-B., S.N. and M.T.; formal analysis, M.T.; investigation, M.T.; resources, M.T.; data curation, M.T.; writing—original draft preparation, M.T.; writing—review and editing, J.G.-B. and S.N.; visualization, M.T.; supervision, J.G.-B. and S.N.; project administration, J.G.-B. and S.N.; funding acquisition, J.G.-B. and S.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Research Project PGC2018-094773-B-C32 [MINECO-FEDER]; Centro de Excelencia Severo Ochoa Program SEV-2016-0597; Ramón y Cajal program through Grant No. RYC-2014-15843; Grant No. BES-2016-077817 (MINECO-FPI).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors acknowledge support from the Research Project PGC2018-094773-B-C32 [MINECO-FEDER], the Centro de Excelencia Severo Ochoa Program SEV-2016-0597 and use of the Hydra HPC cluster at the Instituto de Física Teórica (UAM/CSIC). S.N. acknowledges support from the Ramón y Cajal program through Grant No. RYC-2014-15843. M.T. acknowledges support from the Grant No. BES-2016-077817 (MINECO-FPI).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zwicky, F. Die Rotverschiebung von extragalaktischen Nebeln. Helv. Phys. Acta 1933, 6, 110–127. [Google Scholar] [CrossRef]
  2. Rubin, V.C.; Ford, W.; Kent, J. Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions. Astrophys. J. 1970, 159, 379–403. [Google Scholar] [CrossRef]
  3. Massey, R.; Rhodes, J.; Ellis, R.; Scoville, N.; Leauthaud, A.; Finoguenov, A.; Capak, P.; Bacon, D.; Aussel, H.; Kneib, J.-P.; et al. Dark matter maps reveal cosmic scaffolding. Nature 2007, 445, 286. [Google Scholar] [CrossRef] [PubMed]
  4. Ade, P.; Aghanim, N.; Arnaud, M.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Banday, A.J.; Barreiro, R.B.; Bartolo, N.; Battaner, E.; et al. Planck 2015 results. XIV. Dark energy and modified gravity. Astron. Astrophys. 2016, 594, A14. [Google Scholar] [CrossRef] [Green Version]
  5. Aghanim, N.; Akrami, Y.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Ballardini, M.; Banday, A.J.; Barreiro, R.B.; Bartolo, N.; Basak, S.; et al. Planck 2018 results. V. CMB power spectra and likelihoods. arXiv 2019, arXiv:1907.12875. [Google Scholar]
  6. Aghanim, N.; Akrami, Y.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Ballardini, M.; Banday, A.J.; Barreiro, R.B.; Bartolo, N.; Basak, S.; et al. Planck 2018 results. VI. Cosmological parameters. arXiv 2018, arXiv:1807.06209. [Google Scholar]
  7. Hui, L.; Ostriker, J.P.; Tremaine, S.; Witten, E. Ultralight scalars as cosmological dark matter. Phys. Rev. D 2017, 95, 043541. [Google Scholar] [CrossRef] [Green Version]
  8. Silk, J.; Moore, B.; Diemand, J.; Bullock, J.; Kaplinghat, M.; Strigari, L.; Mellier, Y.; Merritt, D.; Bekenstein, J.; Gelmini, G.; et al. Particle Dark Matter: Observations, Models and Searches; Cambridge Univ. Press: Cambridge, UK, 2010. [Google Scholar] [CrossRef]
  9. Lin, T. Dark matter models and direct detection. PoS 2019, 333, 009. [Google Scholar] [CrossRef] [Green Version]
  10. Schumann, M. Direct Detection of WIMP Dark Matter: Concepts and Status. J. Phys. G 2019, 46, 103003. [Google Scholar] [CrossRef] [Green Version]
  11. Aprile, E.; Aalbers, J.; Agostini, F.; Alfonsi, M.; Althueser, L.; Amaro, F.D.; Antochi, V.C.; Angelino, E.; Angevaare, J.R.; Arneodo, F.; et al. Observation of Excess Electronic Recoil Events in XENON1T. arXiv 2020, arXiv:2006.09721. [Google Scholar]
  12. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abernathy, M.R.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X.; et al. Observation of Gravitational Waves from a Binary Black Hole Merger. Phys. Rev. Lett. 2016, 116, 061102. [Google Scholar] [CrossRef] [PubMed]
  13. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abernathy, M.R.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X.; et al. Astrophysical Implications of the Binary Black-Hole Merger GW150914. Astrophys. J. Lett. 2016, 818, L22. [Google Scholar] [CrossRef]
  14. Belczynski, K.; Holz, D.E.; Bulik, T.; O’Shaughnessy, R. The first gravitational-wave source from the isolated evolution of two 40-100 Msun stars. Nature 2016, 534, 512. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Clesse, S.; García-Bellido, J. Massive Primordial Black Holes from Hybrid Inflation as Dark Matter and the seeds of Galaxies. Phys. Rev. 2015, D92, 023524. [Google Scholar] [CrossRef] [Green Version]
  16. Bird, S.; Cholis, I.; Muñoz, J.B.; Ali-Haïmoud, Y.; Kamionkowski, M.; Kovetz, E.D.; Raccanelli, A.; Riess, A.G. Did LIGO detect dark matter? Phys. Rev. Lett. 2016, 116, 201301. [Google Scholar] [CrossRef]
  17. Clesse, S.; García-Bellido, J. Detecting the gravitational wave background from primordial black hole dark matter. Phys. Dark Univ. 2017, 18, 105–114. [Google Scholar] [CrossRef] [Green Version]
  18. Sasaki, M.; Suyama, T.; Tanaka, T.; Yokoyama, S. Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914. Phys. Rev. Lett. 2016, 117, 061101. [Google Scholar] [CrossRef] [Green Version]
  19. Carr, B.; Kuhnel, F.; Sandstad, M. Primordial Black Holes as Dark Matter. Phys. Rev. 2016, D94, 083504. [Google Scholar] [CrossRef] [Green Version]
  20. García-Bellido, J. Massive Primordial Black Holes as Dark Matter and their detection with Gravitational Waves. J. Phys. Conf. Ser. 2017, 840, 012032. [Google Scholar] [CrossRef]
  21. Clesse, S.; García-Bellido, J. GW190425, GW190521 and GW190814: Three candidate mergers of primordial black holes from the QCD epoch. arXiv 2020, arXiv:2007.06481. [Google Scholar]
  22. Clesse, S.; García-Bellido, J.; Orani, S. Detecting the Stochastic Gravitational Wave Background from Primordial Black Hole Formation. arXiv 2018, arXiv:1812.11011. [Google Scholar]
  23. Connor, L.; Lin, H.H.; Masui, K.; Oppermann, N.; Pen, U.L.; Peterson, J.B.; Roman, A.; Sievers, J. Constraints on the FRB rate at 700–900 MHz. Mon. Not. R. Astron. Soc. 2016, 460, 1054–1058. [Google Scholar] [CrossRef] [Green Version]
  24. Muñoz, J.B.; Kovetz, E.D.; Dai, L.; Kamionkowski, M. Lensing of Fast Radio Bursts as a Probe of Compact Dark Matter. Phys. Rev. Lett. 2016, 117, 091301. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Kramer, M.; Stappers, B. LOFAR, LEAP and beyond: Using next generation telescopes for pulsar astrophysics. PoS 2010, ISKAF2010, 034. [Google Scholar] [CrossRef] [Green Version]
  26. Schutz, K.; Liu, A. Pulsar timing can constrain primordial black holes in the LIGO mass window. Phys. Rev. D 2017, 95, 023002. [Google Scholar] [CrossRef] [Green Version]
  27. Madau, P.; Meiksin, A.; Rees, M.J. 21-CM tomography of the intergalactic medium at high redshift. Astrophys. J. 1997, 475, 429. [Google Scholar] [CrossRef] [Green Version]
  28. Tashiro, H.; Sugiyama, N. The effect of primordial black holes on 21 cm fluctuations. Mon. Not. R. Astron. Soc. 2013, 435, 3001. [Google Scholar] [CrossRef] [Green Version]
  29. Gong, J.O.; Kitajima, N. Small-scale structure and 21cm fluctuations by primordial black holes. JCAP 2017, 08, 017. [Google Scholar] [CrossRef] [Green Version]
  30. Hasinger, G. Illuminating the dark ages: Cosmic backgrounds from accretion onto primordial black hole dark matter. arXiv 2020, arXiv:2003.05150. [Google Scholar] [CrossRef]
  31. Belotsky, K.M.; Dokuchaev, V.I.; Eroshenko, Y.N.; Esipova, E.A.; Khlopov, M.Y.; Khromykh, L.A.; Kirillov, A.A.; Nikulin, V.V.; Rubin, S.G.; Svadkovsky, I.V. Clusters of primordial black holes. Eur. Phys. J. C 2019, 79, 246. [Google Scholar] [CrossRef] [Green Version]
  32. Desjacques, V.; Riotto, A. Spatial clustering of primordial black holes. Phys. Rev. D 2018, 98, 123533. [Google Scholar] [CrossRef] [Green Version]
  33. Calcino, J.; García-Bellido, J.; Davis, T.M. Updating the MACHO fraction of the Milky Way dark halo with improved mass models. Mon. Not. R. Astron. Soc. 2018, 479, 2889–2905. [Google Scholar] [CrossRef] [Green Version]
  34. García-Bellido, J.; Clesse, S.; Fleury, P. Primordial black holes survive SN lensing constraints. Phys. Dark Univ. 2018, 20, 95–100. [Google Scholar] [CrossRef] [Green Version]
  35. Inman, D.; Ali-Haïmoud, Y. Early structure formation in primordial black hole cosmologies. Phys. Rev. D 2019, 100, 083528. [Google Scholar] [CrossRef] [Green Version]
  36. Clesse, S.; García-Bellido, J. Seven Hints for Primordial Black Hole Dark Matter. Phys. Dark Univ. 2018, 22, 137–146. [Google Scholar] [CrossRef] [Green Version]
  37. Zel’dovich, Y.; Novikov, I. The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model. Sov. Astron. 1967, 10, 602. [Google Scholar]
  38. Hawking, S. Gravitationally collapsed objects of very low mass. Mon. Not. R. Astron. Soc. 1971, 152, 75. [Google Scholar] [CrossRef]
  39. Carr, B.J.; Hawking, S. Black holes in the early Universe. Mon. Not. R. Astron. Soc. 1974, 168, 399–415. [Google Scholar] [CrossRef] [Green Version]
  40. Carr, B.J. The Primordial black hole mass spectrum. Astrophys. J. 1975, 201, 1–19. [Google Scholar] [CrossRef]
  41. Carr, B.J. Primordial black holes and inflation. Current topics in astrofundamental physics. In Proceedings of the International School of Astrophysics *D. Chalonge*, 5th Course, Erice, Italy, 7–15 September 1996; pp. 306–321. [Google Scholar]
  42. García-Bellido, J.; Linde, A.D.; Wands, D. Density perturbations and black hole formation in hybrid inflation. Phys. Rev. 1996, D54, 6040–6058. [Google Scholar] [CrossRef] [Green Version]
  43. García-Bellido, J.; Ruiz Morales, E. Primordial black holes from single field models of inflation. Phys. Dark Univ. 2017, 18, 47–54. [Google Scholar] [CrossRef] [Green Version]
  44. Ezquiaga, J.M.; García-Bellido, J.; Ruiz Morales, E. Primordial Black Hole production in Critical Higgs Inflation. Phys. Lett. 2018, B776, 345–349. [Google Scholar] [CrossRef]
  45. Ezquiaga, J.M.; García-Bellido, J. Quantum diffusion beyond slow-roll: Implications for primordial black-hole production. JCAP 2018, 08, 018. [Google Scholar] [CrossRef] [Green Version]
  46. García-Bellido, J.; Peloso, M.; Unal, C. Gravitational waves at interferometer scales and primordial black holes in axion inflation. JCAP 2016, 12, 031. [Google Scholar] [CrossRef]
  47. García-Bellido, J.; Peloso, M.; Unal, C. Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter. JCAP 2017, 1709, 013. [Google Scholar] [CrossRef] [Green Version]
  48. Doctor, Z.; Wysocki, D.; O’Shaughnessy, R.; Holz, D.E.; Farr, B. Black Hole Coagulation: Modeling Hierarchical Mergers in Black Hole Populations. Astrophys. J. 2019. [Google Scholar] [CrossRef] [Green Version]
  49. Carr, B.J. Primordial black holes: Do they exist and are they useful? In Proceedings of the 59th Yamada Conference on Inflating Horizon of Particle Astrophysics and Cosmology, Tokyo, Japan, 20–24 June 2005. [Google Scholar]
  50. Page, D.N. Particle Emission Rates from a Black Hole: Massless Particles from an Uncharged, Nonrotating Hole. Phys. Rev. 1976, D13, 198–206. [Google Scholar] [CrossRef]
  51. Rhoades, C.E., Jr.; Ruffini, R. Maximum mass of a neutron star. Phys. Rev. Lett. 1974, 32, 324–327. [Google Scholar] [CrossRef] [Green Version]
  52. Carr, B.J.; Kohri, K.; Sendouda, Y.; Yokoyama, J. New cosmological constraints on primordial black holes. Phys. Rev. 2010, D81, 104019. [Google Scholar] [CrossRef] [Green Version]
  53. Barnacka, A.; Glicenstein, J.F.; Moderski, R. New constraints on primordial black holes abundance from femtolensing of gamma-ray bursts. Phys. Rev. 2012, D86, 043001. [Google Scholar] [CrossRef] [Green Version]
  54. Wilkinson, P.N.; Henstock, D.R.; Browne, I.W.A.; Polatidis, A.G.; Augusto, P.; Readhead, A.C.S.; Pearson, T.J.; Xu, W.; Taylor, G.B.; Vermeulen, R.C. Limits on the cosmological abundance of supermassive compact objects from a search for multiple imaging in compact radio sources. Phys. Rev. Lett. 2001, 86, 584–587. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Tisserand, P.; Le Guillou, L.; Afonso, C.; Albert, J.N.; Andersen, J.; Ansari, R.; Aubourg, E.; Bareyre, P.; Beaulieu, J.P.; Charlot, X.; et al. Limits on the Macho Content of the Galactic Halo from the EROS-2 Survey of the Magellanic Clouds. Astron. Astrophys. 2007, 469, 387–404. [Google Scholar] [CrossRef]
  56. Wyrzykowski, L.; Skowron, J.; Kozlowski, S.; Udalski, A.; Szymanski, M.K.; Kubiak, M.; Pietrzynski, G.; Soszynski, I.; Szewczyk, O.; Ulaczyk, K.; et al. The OGLE View of Microlensing towards the Magellanic Clouds. IV. OGLE-III SMC Data and Final Conclusions on MACHOs. Mon. Not. R. Astron. Soc. 2011, 416, 2949. [Google Scholar] [CrossRef] [Green Version]
  57. Niikura, H.; Takada, M.; Yasuda, N.; Lupton, R.H.; Sumi, T.; More, S.; Kurita, T.; Sugiyama, S.; More, A.; Oguri, M.; et al. Microlensing constraints on primordial black holes with the Subaru/HSC Andromeda observation. arXiv 2017, arXiv:1701.02151. [Google Scholar] [CrossRef]
  58. Griest, K.; Cieplak, A.M.; Lehner, M.J. Experimental Limits on Primordial Black Hole Dark Matter from the First 2 yr of Kepler Data. Astrophys. J. 2014, 786, 158. [Google Scholar] [CrossRef] [Green Version]
  59. Oguri, M.; Diego, J.M.; Kaiser, N.; Kelly, P.L.; Broadhurst, T. Understanding caustic crossings in giant arcs: Characteristic scales, event rates, and constraints on compact dark matter. Phys. Rev. 2018, D97, 023518. [Google Scholar] [CrossRef] [Green Version]
  60. Graham, P.W.; Rajendran, S.; Varela, J. Dark Matter Triggers of Supernovae. Phys. Rev. 2015, D92, 063007. [Google Scholar] [CrossRef] [Green Version]
  61. Capela, F.; Pshirkov, M.; Tinyakov, P. Constraints on primordial black holes as dark matter candidates from capture by neutron stars. Phys. Rev. 2013, D87, 123524. [Google Scholar] [CrossRef] [Green Version]
  62. Quinn, D.P.; Wilkinson, M.I.; Irwin, M.J.; Marshall, J.; Koch, A.; Belokurov, V. On the Reported Death of the MACHO Era. Mon. Not. R. Astron. Soc. 2009, 396, 11. [Google Scholar] [CrossRef]
  63. Carr, B.J.; Sakellariadou, M. Dynamical constraints on dark compact objects. Astrophys. J. 1999, 516, 195–220. [Google Scholar] [CrossRef]
  64. Brandt, T.D. Constraints on MACHO Dark Matter from Compact Stellar Systems in Ultra-Faint Dwarf Galaxies. Astrophys. J. 2016, 824, L31. [Google Scholar] [CrossRef]
  65. Li, T.S.; Simon, J.D.; Drlica-Wagner, A.; Bechtol, K.; Wang, M.Y.; García-Bellido, J.; Frieman, J.; Marshall, J.L.; James, D.J.; Strigari, L.; et al. Farthest Neighbor: The Distant Milky Way Satellite Eridanus II. Astrophys. J. 2017, 838, 8. [Google Scholar] [CrossRef]
  66. Ali-Haïmoud, Y.; Kamionkowski, M. Cosmic microwave background limits on accreting primordial black holes. Phys. Rev. 2017, D95, 043534. [Google Scholar] [CrossRef] [Green Version]
  67. Poulin, V.; Serpico, P.D.; Calore, F.; Clesse, S.; Kohri, K. CMB bounds on disk-accreting massive primordial black holes. Phys. Rev. 2017, D96, 083524. [Google Scholar] [CrossRef] [Green Version]
  68. Gaggero, D.; Bertone, G.; Calore, F.; Connors, R.M.T.; Lovell, M.; Markoff, S.; Storm, E. Searching for Primordial Black Holes in the radio and X-ray sky. Phys. Rev. Lett. 2017, 118, 241101. [Google Scholar] [CrossRef] [PubMed]
  69. Inoue, Y.; Kusenko, A. New X-ray bound on density of primordial black holes. JCAP 2017, 1710, 034. [Google Scholar] [CrossRef] [Green Version]
  70. Carr, B.; Clesse, S.; García-Bellido, J.; Kuhnel, F. Cosmic Conundra Explained by Thermal History and Primordial Black Holes. Phys. Dark Univ. 2021, 31, 100755. [Google Scholar] [CrossRef]
  71. García-Bellido, J. Primordial black holes and the origin of the matter–antimatter asymmetry. Philos. Trans. R. Soc. Lond. 2019, A377, 20190091. [Google Scholar] [CrossRef]
  72. García-Bellido, J.; Clesse, S. Constraints from microlensing experiments on clustered primordial black holes. Phys. Dark Univ. 2018, 19, 144–148. [Google Scholar] [CrossRef] [Green Version]
  73. García-Bellido, J. Primordial Black Holes. PoS 2018, EDSU2018, 042. [Google Scholar] [CrossRef]
  74. García-Bellido, J.; Nesseris, S. Gravitational wave bursts from Primordial Black Hole hyperbolic encounters. Phys. Dark Univ. 2017, 18, 123–126. [Google Scholar] [CrossRef] [Green Version]
  75. García-Bellido, J.; Nesseris, S. Gravitational wave energy emission and detection rates of Primordial Black Hole hyperbolic encounters. Phys. Dark Univ. 2018, 21, 61–69. [Google Scholar] [CrossRef] [Green Version]
  76. Calvert, R. AstroGrav. Available online: http://www.astrograv.co.uk (accessed on 12 January 2019).
  77. Adamek, J.; Daverio, D.; Durrer, R.; Kunz, M. GEVOLUTION: A cosmological N-body code based on General Relativity. JCAP 2016, 1607, 053. [Google Scholar] [CrossRef] [Green Version]
  78. Springel, V. The Cosmological simulation code GADGET-2. Mon. Not. R. Astron. Soc. 2005, 364, 1105–1134. [Google Scholar] [CrossRef] [Green Version]
  79. Mouri, H.; Taniguchi, Y. Runaway merging of black holes: Analytical constraint on the timescale. Astrophys. J. 2002, 566, L17–L20. [Google Scholar] [CrossRef] [Green Version]
  80. Chatterjee, S.; Umbreit, S.; Fregeau, J.M.; Rasio, F.A. Understanding the dynamical state of globular clusters: core-collapsed versus non-core-collapsed. Mon. Not. R. Astron. Soc. 2013, 429, 2881–2893. [Google Scholar] [CrossRef] [Green Version]
  81. Bovy, J. Stellar inventory of the solar neighbourhood using Gaia DR1. Mon. Not. R. Astron. Soc. 2017, 470, 1360–1387. [Google Scholar] [CrossRef]
  82. Peters, P.C. Gravitational Radiation and the Motion of Two Point Masses. Phys. Rev. 1964, 136, B1224–B1232. [Google Scholar] [CrossRef]
  83. Binney, J.; Tremaine, S. Galactic Dynamics, 2nd ed.; Princeton University Press: Princeton NJ, USA, 2008. [Google Scholar]
  84. Eadie, G.; Jurić, M. The Cumulative Mass Profile of the Milky Way as Determined by Globular Cluster Kinematics from GAIA DR2. Astrophys. J. 2019, 875, 159. [Google Scholar] [CrossRef]
  85. Gow, A.D.; Byrnes, C.T.; Hall, A.; Peacock, J.A. Primordial black hole merger rates: Distributions for multiple LIGO observables. JCAP 2020, 01, 031. [Google Scholar] [CrossRef] [Green Version]
  86. Bartolo, N.; Caprini, C.; Domcke, V.; Figueroa, D.G.; Garcia-Bellido, J.; Guzzetti, M.C.; Liguori, M.; Matarrese, S.; Peloso, M.; Petiteau, A.; et al. Science with the space-based interferometer LISA. IV: Probing inflation with gravitational waves. JCAP 2016, 1612, 026. [Google Scholar] [CrossRef]
  87. Arun, K.; Gudennavar, S.; Sivaram, C. Dark matter, dark energy, and alternate models: A review. Adv. Space Res. 2017, 60, 166–186. [Google Scholar] [CrossRef] [Green Version]
  88. Li, X.; Tang, L.; Lin, H.N. Comparing dark matter models, modified Newtonian dynamics and modified gravity in accounting for galaxy rotation curves. Chin. Phys. C 2017, 41, 055101. [Google Scholar] [CrossRef] [Green Version]
  89. Odintsov, S.; Oikonomou, V. f(R) Gravity Inflation with String-Corrected Axion Dark Matter. Phys. Rev. D 2019, 99, 064049. [Google Scholar] [CrossRef] [Green Version]
  90. Odintsov, S.; Oikonomou, V. Geometric Inflation and Dark Energy with Axion F(R) Gravity. Phys. Rev. D 2020, 101, 044009. [Google Scholar] [CrossRef] [Green Version]
  91. Nesseris, S.; Sapone, D.; Sypsas, S. Evaporating primordial black holes as varying dark energy. Phys. Dark Univ. 2020, 27, 100413. [Google Scholar] [CrossRef] [Green Version]
  92. Arjona, R.; García-Bellido, J.; Nesseris, S. Cosmological constraints on nonadiabatic dark energy perturbations. Phys. Rev. 2020, D102, 103526. [Google Scholar] [CrossRef]
  93. Jedamzik, K. Primordial Black Hole Dark Matter and the LIGO/Virgo observations. JCAP 2020, 09, 022. [Google Scholar] [CrossRef]
1.
2.
After our paper was submitted to the arXiv we received a message informing us that a related paper had appeared a few days before [93], addressing the problem of PBH clustering dynamics is a similar way.
Figure 2. Mass IC profile from Table 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 m tot = i m i = i ( m i ) = 2.19 × 10 4 M .
Figure 2. Mass IC profile from Table 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 m tot = i m i = i ( m i ) = 2.19 × 10 4 M .
Universe 07 00018 g002
Figure 3. Position IC profile from Table 1 and an example realisation, from the mean μ x = 1.6 pc and scale parameter a x = 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 d S = 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.
Figure 3. Position IC profile from Table 1 and an example realisation, from the mean μ x = 1.6 pc and scale parameter a x = 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 d S = 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.
Universe 07 00018 g003
Figure 4. Universal cluster cumulative mass and rotational velocity curves from Table 1, computed for a perfectly traced distribution of PBHs computed from Equation (19) where N O . Here, N O , the object mean logarithmic mass, μ m , and cluster mean displacement from the barycentre, σ x . The dark grey line marks the curve value at the radii containing 50% the cumulative mass, while the light grey lines do so for 68%, 95% and 99% the cumulative mass.
Figure 4. Universal cluster cumulative mass and rotational velocity curves from Table 1, computed for a perfectly traced distribution of PBHs computed from Equation (19) where N O . Here, N O , the object mean logarithmic mass, μ m , and cluster mean displacement from the barycentre, σ x . The dark grey line marks the curve value at the radii containing 50% the cumulative mass, while the light grey lines do so for 68%, 95% and 99% the cumulative mass.
Universe 07 00018 g004
Figure 5. Velocity IC shell-by shell (left) and integrated distributions (right) computed from from Table 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 d S = 4 π x 2 , from x = 0 to x = 8 σ . Both the core and the outskirts of the PBH cluster are roughly understood as the regions with 0 x a x and 2 a x x 4 a x , respectively, and are not heavily populated; the core exhibits very low velocities in the range v ( 0.0 , 1.0 ) km / s , while the outskirts exhibits intermediate velocities in the range v ( 1.0 , 3.0 ) km / s , whereas the populated intermediate bins do show the largest velocities in the range v ( 1.0 , 5.0 ) km / s .
Figure 5. Velocity IC shell-by shell (left) and integrated distributions (right) computed from from Table 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 d S = 4 π x 2 , from x = 0 to x = 8 σ . Both the core and the outskirts of the PBH cluster are roughly understood as the regions with 0 x a x and 2 a x x 4 a x , respectively, and are not heavily populated; the core exhibits very low velocities in the range v ( 0.0 , 1.0 ) km / s , while the outskirts exhibits intermediate velocities in the range v ( 1.0 , 3.0 ) km / s , whereas the populated intermediate bins do show the largest velocities in the range v ( 1.0 , 5.0 ) km / s .
Universe 07 00018 g005
Figure 6. Velocity IC distribution profile from computed from Table 1 and an example realisation, from the mean μ v = 2.0 km / s and scale parameter a v = 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 realisations 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.
Figure 6. Velocity IC distribution profile from computed from Table 1 and an example realisation, from the mean μ v = 2.0 km / s and scale parameter a v = 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 realisations 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.
Universe 07 00018 g006
Figure 7. Positions and velocities time evolution of a subset 10% of objects in one particular realisation. In blue and teal are the position and velocity curves for cluster and ejecta objects, respectively. For both positions and velocities, 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, and the thin black lines mark the 68% and 95% confidence regions. The mean fits obtained for the position and velocity profiles at each time slice are shown in.
Figure 7. Positions and velocities time evolution of a subset 10% of objects in one particular realisation. In blue and teal are the position and velocity curves for cluster and ejecta objects, respectively. For both positions and velocities, 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, and the thin black lines mark the 68% and 95% confidence regions. The mean fits obtained for the position and velocity profiles at each time slice are shown in.
Universe 07 00018 g007
Figure 8. Top: Total population of remaining cluster and ejecta objects in blue and teal, respectively; N R , 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 δ N R , 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 B i 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 N R = 5000 realisations. The dotted vertical grey lines separate the N r = 7 consecutive runs, and the dotted vertical grey line signals the burn-in time.
Figure 8. Top: Total population of remaining cluster and ejecta objects in blue and teal, respectively; N R , 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 δ N R , 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 B i 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 N R = 5000 realisations. The dotted vertical grey lines separate the N r = 7 consecutive runs, and the dotted vertical grey line signals the burn-in time.
Universe 07 00018 g008
Figure 10. Binary mass-related profiles segregated by their population (cluster and ejecta) type integrated across all time slices in the simulations from the burn-in time, t BI = 2.64 × 10 5 yr , to the present. Profiles are normalised by the total number of objects in the simulations at the initial time, N t 0 , and distribution mean and median values are given in Table 12. Note that panels (b,d,f) are arranged by choosing m 1 > m 2 to avoid double-plotting, and therefore all binary systems are represented by a single point at or below the first quadrant diagonal.
Figure 10. Binary mass-related profiles segregated by their population (cluster and ejecta) type integrated across all time slices in the simulations from the burn-in time, t BI = 2.64 × 10 5 yr , to the present. Profiles are normalised by the total number of objects in the simulations at the initial time, N t 0 , and distribution mean and median values are given in Table 12. Note that panels (b,d,f) are arranged by choosing m 1 > m 2 to avoid double-plotting, and therefore all binary systems are represented by a single point at or below the first quadrant diagonal.
Universe 07 00018 g010
Figure 11. Merger mass-related profiles segregated by the population (cluster and ejecta) type integrated across all time slices in the simulations staring at the time of first-merger t FM = 1.38 × 10 9 yr up to the present. Profiles are normalised by the total number of objects in the simulations at the initial time, N t 0 , and distribution mean and median values are given in Table 11. Note that panels (b,d,f) are arranged again as in Figure 10 by choosing m 1 > m 2 to avoid double-plotting and so again all binary systems are represented by a single point at or below the first quadrant diagonal.
Figure 11. Merger mass-related profiles segregated by the population (cluster and ejecta) type integrated across all time slices in the simulations staring at the time of first-merger t FM = 1.38 × 10 9 yr up to the present. Profiles are normalised by the total number of objects in the simulations at the initial time, N t 0 , and distribution mean and median values are given in Table 11. Note that panels (b,d,f) are arranged again as in Figure 10 by choosing m 1 > m 2 to avoid double-plotting and so again all binary systems are represented by a single point at or below the first quadrant diagonal.
Universe 07 00018 g011
Figure 12. Density, position and velocity quasi-static profiles segregated by population (cluster and ejecta) type at selected times. Fits for the mass, position, velocity and density probability distributions are given in Table 14. Note the mass distribution is not displayed as, 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, as seen in Table 14, the cluster and ejecta mass profiles do not segregate significantly over time.
Figure 12. Density, position and velocity quasi-static profiles segregated by population (cluster and ejecta) type at selected times. Fits for the mass, position, velocity and density probability distributions are given in Table 14. Note the mass distribution is not displayed as, 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, as seen in Table 14, the cluster and ejecta mass profiles do not segregate significantly over time.
Universe 07 00018 g012
Figure 14. Shown are the 1 σ and 2 σ contours PBHs orbital eccentricity vs. semi-major axis segregated by population (cluster and ejecta) type at selected times, t BI = 2.64 × 10 5 yr plus the initial time-slice t 0 . All panels: Black dots show the scatter a 100 sample of individual PBHs in each time-slice, randomly selected in the case of the bare distribution and selected with a uniform probability distribution. The Best-Fit (BF) values of the distributions, marked by the large coloured dots in the plots, are given in Table 17.
Figure 14. Shown are the 1 σ and 2 σ contours PBHs orbital eccentricity vs. semi-major axis segregated by population (cluster and ejecta) type at selected times, t BI = 2.64 × 10 5 yr plus the initial time-slice t 0 . All panels: Black dots show the scatter a 100 sample of individual PBHs in each time-slice, randomly selected in the case of the bare distribution and selected with a uniform probability distribution. The Best-Fit (BF) values of the distributions, marked by the large coloured dots in the plots, are given in Table 17.
Universe 07 00018 g014
Figure 15. All panels: Binary lifetime histogram for in-cluster and in-ejecta objects, for each time slice in the simulations. The total number of in-cluster binaries amounts to N C , B = 3.85 × 10 5 objects, while the total number of in-ejecta binaries amounts to N E , B = 7.46 × 10 4 objects, throughout all realisations and time-slices.
Figure 15. All panels: Binary lifetime histogram for in-cluster and in-ejecta objects, for each time slice in the simulations. The total number of in-cluster binaries amounts to N C , B = 3.85 × 10 5 objects, while the total number of in-ejecta binaries amounts to N E , B = 7.46 × 10 4 objects, throughout all realisations and time-slices.
Universe 07 00018 g015
Figure 16. Shown are the relative orbital eccentricity and impact distributions for the identified two-body encounters at selected times. These are computed by taking, for each realisation, the 100 PBHs closest to the cluster barycenter and calculating the array of 100 × 100 relative eccentricities among them, for the 6 time slices that split the last four runs and the IC. Fits for the distributions are given in Table 18.
Figure 16. Shown are the relative orbital eccentricity and impact distributions for the identified two-body encounters at selected times. These are computed by taking, for each realisation, the 100 PBHs closest to the cluster barycenter and calculating the array of 100 × 100 relative eccentricities among them, for the 6 time slices that split the last four runs and the IC. Fits for the distributions are given in Table 18.
Universe 07 00018 g016
Table 2. The N r = 7 time-runs, denoted by the index r = 1 , 2 , , 7 , into which each realisation is partitioned for a total of N t = 64 marked by the index i = i i = i = 1 , 2 , , 7 and the number of non-overlapping time slices in each run, as well as the initial time, ending time and time-step factored by the present age of the Universe, t 64 = 13.8 × 10 10 yr , which corresponds to the current age of the Universe t U = ( 13.800 ± 0.024 ) × 10 10 yr according to Planck 2018 (TT, TE, EE+lowE, see Ref. [6]). The time intervals and time-slices shown in Figure 7, displaying the PBHs positions and velocities at each of these.
Table 2. The N r = 7 time-runs, denoted by the index r = 1 , 2 , , 7 , into which each realisation is partitioned for a total of N t = 64 marked by the index i = i i = i = 1 , 2 , , 7 and the number of non-overlapping time slices in each run, as well as the initial time, ending time and time-step factored by the present age of the Universe, t 64 = 13.8 × 10 10 yr , which corresponds to the current age of the Universe t U = ( 13.800 ± 0.024 ) × 10 10 yr according to Planck 2018 (TT, TE, EE+lowE, see Ref. [6]). The time intervals and time-slices shown in Figure 7, displaying the PBHs positions and velocities at each of these.
Simulation Time-Runs:
r ( i , i + ) N t t i [yr] t i + [yr] Δ t r [ yr ] r ( i , i + ) N t t i [yr] t i + [yr] Δ t r [yr]
0 ( 0 , 1 ) 2 0.00 13.8 × 10 3 13.8 × 10 3 4 ( 29 , 37 ) 9 13.8 × 10 6 13.8 × 10 7 13.8 × 10 6
1 ( 2 , 10 ) 9 13.8 × 10 3 13.8 × 10 4 13.8 × 10 3 5 ( 38 , 46 ) 9 13.8 × 10 7 13.8 × 10 8 13.8 × 10 7
2 ( 11 , 19 ) 9 13.8 × 10 4 13.8 × 10 5 13.8 × 10 4 6 ( 47 , 55 ) 9 13.8 × 10 8 13.8 × 10 9 13.8 × 10 8
3 ( 20 , 28 ) 9 13.8 × 10 5 13.8 × 10 6 13.8 × 10 5 7 ( 56 , 64 ) 9 13.8 × 10 9 13.8 × 10 10 1.38 × 10 10
Table 3. Computed burn-in time of a number of populations of the remaining objects, mainly the cluster N C and ejecta N E populations, the population of objects in bounded N B systems, that is, those objects part of multiple, usually binary, system, or even binary within a binary system, and combinations of those. These burn-in times are shown in Figure 7 when computed from position and velocity evolution, along with Figure 8 and Figure 9 when computed from population evolution.
Table 3. Computed burn-in time of a number of populations of the remaining objects, mainly the cluster N C and ejecta N E populations, the population of objects in bounded N B systems, that is, those objects part of multiple, usually binary, system, or even binary within a binary system, and combinations of those. These burn-in times are shown in Figure 7 when computed from position and velocity evolution, along with Figure 8 and Figure 9 when computed from population evolution.
Snapshot Burn-in Times:
t BI [yr] N R N B N x N v
N C 4.92 × 10 5 3.50 × 10 5 2.95 × 10 5 2.95 × 10 5
N E 3.00 × 10 5 5.93 × 10 5 5.72 × 10 5 5.72 × 10 5
N T 0.00 3.59 × 10 5 5.72 × 10 5 5.72 × 10 5
Table 4. Mean cluster and ejecta position and velocity PL fits of Figure 7 from Equation (31): z P ( t r ) = β i , z , P t α 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, α i , z , P is the PL index and P = C , E denotes the cluster and ejecta populations, respectively.
Table 4. Mean cluster and ejecta position and velocity PL fits of Figure 7 from Equation (31): z P ( t r ) = β i , z , P t α 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, α i , z , P is the PL index and P = C , E denotes the cluster and ejecta populations, respectively.
    t BI = 5.72 × 10 5 yr Cluster Objects Positions:Cluster Objects Velocities:
r ( t i , t i + ) [yr]Fit: r i , x , C 2 β i , x , C   [ pc ] α i , x , C Fit: r i , v , C 2 β i , v , C   [ km / s ] α i , v , C
3 ( t BI yr 1 , 1.38 × 10 6 ) [PL] 1.000 ( 3.13 ± 0.90 ) × 10 5 0.788 ± 0.021 [PL] 1.000 ( 3.8 ± 1.1 ) × 10 3 0.468 ± 0.021
4 1.38 × ( 10 6 , 10 7 ) [PL] 1.000 ( 2.34 ± 0.21 ) × 10 4 0.6407 ± 0.0055 [PL] 1.000 577 ± 49 0.3311 ± 0.0055
5 1.38 × ( 10 7 , 10 8 ) [PL] 1.000 ( 2.374 ± 0.081 ) × 10 4 0.6410 ± 0.0019 [PL] 1.000 911 ± 23 0.3602 ± 0.0014
6 1.38 × ( 10 8 , 10 9 ) [PL] 1.000 ( 2.553 ± 0.061 ) × 10 4 0.6374 ± 0.0012 [PL] 1.000 912 ± 23 0.3599 ± 0.0013
7 1.38 × ( 10 9 , 10 10 ) [PL] 1.000 ( 2.525 ± 0.069 ) × 10 4 0.6374 ± 0.0012 [PL] 1.000 499 ± 31 0.3309 ± 0.0028
    t BI = 5.72 × 10 5 yr Ejecta Objects Positions:Ejecta Objects Velocities:
r ( t i , t i + ) [yr]Fit: r i , x , E 2 β i , x , E   [ pc ] α i , x , E Fit: r i , v , E 2 β i , v , E   [ km / s ] α i , v , E
3 ( t BI yr 1 , 1.38 × 10 6 ) [PL] 1.000 ( 2.29 ± 0.84 ) × 10 7 1.236 ± 0.026 [PL] 1.000 63 ± 22 0.128 ± 0.025
4 1.38 × ( 10 6 , 10 7 ) [PL] 1.000 ( 7.79 ± 0.83 ) × 10 6 0.9917 ± 0.0066 [PL] 1.000 48.6 ± 4.3 0.1072 ± 0.0057
5 1.38 × ( 10 7 , 10 8 ) [PL] 1.000 ( 1.956 ± 0.058 ) × 10 5 0.9362 ± 0.0016 [PL] 1.000 59.0 ± 2.3 0.1201 ± 0.0022
6 1.38 × ( 10 8 , 10 9 ) [PL] 1.000 ( 1.896 ± 0.072 ) × 10 5 0.9370 ± 0.0018 [PL] 1.000 30.0 ± 1.3 0.0841 ± 0.0022
7 1.38 × ( 10 9 , 10 10 ) [PL] 1.000 ( 1.555 ± 0.018 ) × 10 5 0.94624 ± 0.00049 [PL] 1.000 18.26 ± 0.39 0.06032 ± 0.00095
Table 5. Mean remaining in-cluster and in-ejecta object population PL fits of Figure 8 from Equation (31): N P ( t r ) = β i , P t α 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.
Table 5. Mean remaining in-cluster and in-ejecta object population PL fits of Figure 8 from Equation (31): N P ( t r ) = β i , P t α 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.
    t BI = 2.64 × 10 5 yr Cluster Population Fits:Ejecta Population Fits:
r ( t i , t i + ) [yr]Fit: r i , R , C 2 β i , R , C α i , R , C Fit: r i , R , E 2 β i , R , E α i , R , E
3 ( t BI yr 1 , 1.30 × 10 6 ) [PL] 1.000 ( 1.982 ± 0.049 ) × 10 3 0.0549 ± 0.0018 [PL] 0.991 ( 5.4 ± 8.7 ) × 10 5 1.00 ± 0.12
4 1.38 × ( 10 6 , 10 7 ) [PL] 1.000 ( 3.037 ± 0.085 ) × 10 3 0.0847 ± 0.0018 [PL] 0.999 0.320 ± 0.095 0.0410 ± 0.018
5 1.38 × ( 10 7 , 10 8 ) [PL] 1.000 ( 3.8684 ± 0.0088 ) × 10 3 0.09968 ± 0.00013 [PL] 1.000 9.9 ± 1.4 0.2000 ± 0.0076
6 1.38 × ( 10 8 , 10 9 ) [PL] 1.000 ( 4.210 ± 0.019 ) × 10 3 0.10418 ± 0.00023 [PL] 1.000 46.0 ± 2.9 0.1200 ± 0.0031
7 1.38 × ( 10 9 , 10 10 ) [PL] 1.000 ( 4.874 ± 0.082 ) × 10 3 0.11107 ± 0.00075 [PL] 1.000 100.0 ± 3.1 0.0780 ± 0.0013
Table 6. Mean remaining in-cluster and in-ejecta object population fractions f ˜ i , R , P and median f ˜ i , 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 N R = 5000 realisations at selected times.
Table 6. Mean remaining in-cluster and in-ejecta object population fractions f ˜ i , R , P and median f ˜ i , 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 N R = 5000 realisations at selected times.
Remaining Cluster and Ejecta Population Fraction and Median:
i t i [ yr ] f i , R , C f ˜ i , R , C f i , R , E f ˜ i , R , E i t i [ yr ] f i , R , C f ˜ i , R , C f i , R , E f ˜ i , R , E
00 0.9975 + 0.0002 0.0025 1.000 0.0025 + 0.0025 0.0002 0.000 37 1.38 × 10 7 0.76 + 0.08 0.13 0.763 0.25 + 0.12 0.09 0.237
10 1.38 × 10 4 0.9950 + 0.0015 0.0015 0.999 0.0050 + 0.0015 0.0015 0.001 46 1.38 × 10 8 0.59 + 0.11 0.16 0.595 0.40 + 0.16 1.10 0.505
19 1.38 × 10 5 0.990 + 0.002 0.010 0.991 0.009 + 0.011 0.004 0.009 55 1.38 × 10 9 0.47 + 0.11 0.17 0.476 0.53 + 0.17 0.11 0.524
28 1.38 × 10 6 0.912 + 0.051 0.072 0.910 0.092 + 0.072 0.050 0.090 64 1.38 × 10 10 0.36 + 0.11 0.18 0.368 0.63 + 0.17 0.11 0.632
Table 9. Fraction of binaries at each of the eight time-slices dividing the seven runs. Each level l B indicates the level in the binary hierarchy, with l B = 2 indicating two bound PBHs in a binary system, l B = 3 indicating a binary system within another binary system and l B = 4 indicating a binary system within another two binary systems. Note that there are N r = 7 time-runs, for a total of eight bounding time slices, corresponding to the times shown in Table 2.
Table 9. Fraction of binaries at each of the eight time-slices dividing the seven runs. Each level l B indicates the level in the binary hierarchy, with l B = 2 indicating two bound PBHs in a binary system, l B = 3 indicating a binary system within another binary system and l B = 4 indicating a binary system within another two binary systems. Note that there are N r = 7 time-runs, for a total of eight bounding time slices, corresponding to the times shown in Table 2.
Parent Tree Hierarchy Levels:
i t i [ yr ] l B = 2 l B = 3 l B = 4 i t i [ yr ] l B = 2 l B = 3 l B = 4
00 8.07 × 0.01 9.88 × 10 4 2.14 × 10 6 37 1.38 × 10 7 1.13 × 0.01 4.25 × 10 5 0.00
10 1.38 × 10 4 8.07 × 0.01 9.88 × 10 4 1.98 × 10 6 46 1.38 × 10 8 1.67 × 0.01 1.49 × 10 4 5.63 × 10 7
19 1.38 × 10 5 9.45 × 0.01 1.06 × 10 3 1.56 × 10 6 55 1.38 × 10 9 2.29 × 0.01 3.28 × 10 4 2.93 × 10 6
28 1.38 × 10 6 1.10 × 0.1 3.21 × 10 3 2.65 × 10 5 64 1.38 × 10 10 2.90 × 0.01 6.16 × 10 4 6.45 × 10 6
Table 10. Fraction of mergers grouped from the time slices dividing the first five time-runs. Each level l M indicates the level in the merger tree hierarchy, with l M = 2 indicating the merger two PBHs onto a final one, and l M = 3 indicating the merger of three PBHs onto single PBH, accumulated in the runs.
Table 10. Fraction of mergers grouped from the time slices dividing the first five time-runs. Each level l M indicates the level in the merger tree hierarchy, with l M = 2 indicating the merger two PBHs onto a final one, and l M = 3 indicating the merger of three PBHs onto single PBH, accumulated in the runs.
Merger Tree Hierarchy Levels:
r t i [ yr ] t i + [ yr ] l M = 2 l M = 3
1 5 1.38 × 10 3 1.38 × 10 9 0.00 0.00
6 1.38 × 10 8 1.38 × 10 9 6.00 × 10 7 0.00
7 1.38 × 10 9 1.38 × 10 10 2.24 × 10 5 4.00 × 10 7
Table 12. Distribution mean and median values of the distributions in Figure 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 t BI = 2.64 × 10 5 yr up to the present, for a total of i = 1 N R ( N B ) 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 as their number is negligible compared to the other cases, but we do show them, however, on Table 11.
Table 12. Distribution mean and median values of the distributions in Figure 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 t BI = 2.64 × 10 5 yr up to the present, for a total of i = 1 N R ( N B ) 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 as their number is negligible compared to the other cases, but we do show them, however, on Table 11.
Binary Populations:Binary Distribution Parameters:
Pop.Frequency:Total Mass:Chirp Mass:Mass Ratio:
P i P j N R , B P i , P j f R , B i , j m ¯ T [ M ] m ˜ T [ M ] m ¯ C [ M ] m ˜ C [ M ] q ¯ q ˜
1.27 × 10 6 2.55 × 10 3 303181 45.0 25.2 0.154 0.0552
C P 1.05 × 10 5 2.10 × 10 3 305185 39.9 26.2 0.200 0.0446
E P 2.27 × 10 4 4.54 × 10 4 371202 61.2 29.0 0.404 0.0715
Table 13. Shown are the crossing time, T P cross , and relaxation time, T P rel , of the cluster and ejecta bubbles that populate the galactic halo at selected times, where P = C , E stands for the cluster and ejecta populations.
Table 13. Shown are the crossing time, T P cross , and relaxation time, T P rel , of the cluster and ejecta bubbles that populate the galactic halo at selected times, where P = C , E stands for the cluster and ejecta populations.
      Cluster Time Scales:Ejecta Time Scales:
i t i [ yr ] T C cross [ yr ] T C rel [ yr ] T E cross [ yr ] T E rel [ yr ]
00 ( 4.2 ± 0.7 ) × 10 5 ( 2.3 + 0.7 2.3 ) × 10 5
28 1.38 × 10 6 ( 3.5 ± 0.5 ) × 10 5 ( 8.0 + 5.1 7.2 ) × 10 5 ( 1.72 ± 0.20 ) × 10 6 ( 2.93 + 0.41 0.43 ) × 10 7
37 1.38 × 10 7 ( 2.22 ± 0.33 ) × 10 6 ( 1.21 + 0.52 0.73 ) × 10 7 ( 1.65 ± 0.26 ) × 10 7 ( 2.42 + 0.52 0.53 ) × 10 8
46 1.38 × 10 8 ( 2.4 ± 0.5 ) × 10 7 ( 2.05 + 0.63 0.86 ) × 10 8 ( 1.61 ± 0.28 ) × 10 8 ( 1.89 + 0.51 0.59 ) × 10 9
55 1.38 × 10 9 ( 8.4 ± 2.7 ) × 10 7 ( 8.9 + 3.4 4.0 ) × 10 8 ( 1.64 ± 0.28 ) × 10 9 ( 1.57 + 0.50 0.56 ) × 10 10
64 1.38 × 10 10 ( 1.5 + 1.6 1.5 ) × 10 8 ( 1.9 + 2.0 1.9 ) × 10 9 ( 1.80 ± 0.40 ) × 10 10 ( 1.39 + 0.49 0.68 ) × 10 11
Table 14. Density, mass, position and velocity quasi-static profile fits segregated by population (cluster and ejecta) type of Figure 12. The density profile has been fitted to a PL [PL] in the range delimited by the 0 th and 80 th percentiles in position so as to ensure that the density is computed within the cluster and ejecta spheres. The fit is, following Equation (31), ρ z ( t i ) = β i , P t i α i , P , 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 , P 2 π ) 1 exp ( ( Log Normal z μ z , i , P ) 2 / 2 σ z , i , P 2 ) , 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, χ 3 2 = 2 / π z 2 a z , i , P 3 exp ( z 2 / 2 a z , i , P 2 ) . 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.
Table 14. Density, mass, position and velocity quasi-static profile fits segregated by population (cluster and ejecta) type of Figure 12. The density profile has been fitted to a PL [PL] in the range delimited by the 0 th and 80 th percentiles in position so as to ensure that the density is computed within the cluster and ejecta spheres. The fit is, following Equation (31), ρ z ( t i ) = β i , P t i α i , P , 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 , P 2 π ) 1 exp ( ( Log Normal z μ z , i , P ) 2 / 2 σ z , i , P 2 ) , 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, χ 3 2 = 2 / π z 2 a z , i , P 3 exp ( z 2 / 2 a z , i , P 2 ) . 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.
      Cluster Density Distributions:Ejecta Density Distributions:
i t i [ yr ] Fit: r i , C , ρ 2 β i , C ( ρ ) [ pc ] α i , C ( ρ ) Fit: r i , E , ρ 2 β i , E ( ρ ) [ pc ] α i , E ( ρ )
00[PL] 0.998 3.84 ± 0.75 2.879 ± 0.050
28 13.8 × 10 6 [PL] 1.000 5.85 ± 0.90 3.023 ± 0.031 [PL] 0.999 0.183 ± 0.027 2.898 ± 0.027
37 13.8 × 10 7 [PL] 1.000 1.93 ± 0.17 3.005 ± 0.030 [PL] 1.000 0.2487 ± 0.0019 2.7023 ± 0.0046
46 13.8 × 10 8 [PL] 1.000 0.113 ± 0.016 3.257 ± 0.022 [PL] 1.000 0.035018 ± 0.000062 2.6077 ± 0.0047
55 13.8 × 10 9 [PL] 1.000 0.1057 ± 0.0097 3.246 ± 0.019 [PL] 1.000 ( 3.166 ± 0.013 ) 0.001 2.5875 ± 0.0022
64 13.8 × 10 10 [PL] 1.000 0.0588 ± 0.0020 3.2288 ± 0.0074 [PL] 1.000 ( 2.592 ± 0.046 ) 10 4 2.6072 ± 0.0041
      Cluster Mass Distributions:Ejecta Mass Distributions:
i t i [ yr ] Fit: r i , C , m 2 μ i , C ( m ) [ pc ] σ i , C ( m ) Fit: r i , E , m 2 μ i , E ( m ) [ pc ] σ i , E ( m )
00[LN] 0.997 2.115 ± 0.019 1.602 ± 0.010
28 13.8 × 10 6 [LN] 0.997 2.126 ± 0.019 1.604 ± 0.011 [LN] 0.998 2.039 ± 0.020 1.584 ± 0.011
37 13.8 × 10 7 [LN] 0.997 2.128 ± 0.022 1.605 ± 0.012 [LN] 0.997 2.093 ± 0.020 1.598 ± 0.011
46 13.8 × 10 8 [LN] 0.996 2.134 ± 0.023 1.605 ± 0.013 [LN] 0.998 2.097 ± 0.019 1.600 ± 0.011
55 13.8 × 10 9 [LN] 0.995 2.152 ± 0.025 1.611 ± 0.014 [LN] 0.997 2.094 ± 0.019 1.598 ± 0.011
64 13.8 × 10 10 [LN] 0.994 2.175 ± 0.027 1.617 ± 0.015 [LN] 0.997 2.093 ± 0.020 1.598 ± 0.011
      Cluster Position Distributions:Ejecta Position Distributions:
i t i [ yr ] Fit: r i , C , x 2 a i , C ( x ) [ pc ]
00[MB] 0.988 0.813 ± 0.021
i t i [ yr ] Fit: r i , C , x 2 μ i , C ( x ) [ pc ] σ i , C ( x ) Fit: r i , E , x 2 μ i , E ( x ) [ pc ] σ i , E ( x )
28 13.8 × 10 6 [LN] 0.994 0.448 ± 0.078 1.043 ± 0.036 [LN] 0.994 2.368 ± 0.085 0.943 ± 0.044
37 13.8 × 10 7 [LN] 0.997 1.791 ± 0.086 1.014 ± 0.044 [LN] 0.996 4.37 ± 0.12 0.909 ± 0.062
46 13.8 × 10 8 [LN] 0.999 3.306 ± 0.088 1.052 ± 0.040 [LN] 0.998 6.50 ± 0.12 1.003 ± 0.063
55 13.8 × 10 9 [LN] 0.999 4.781 ± 0.087 1.066 ± 0.038 [LN] 0.999 8.54 ± 0.12 1.009 ± 0.058
64 13.8 × 10 10 [LN] 1.000 6.260 ± 0.095 1.083 ± 0.049 [LN] 0.999 10.75 ± 0.14 1.071 ± 0.076
      Cluster Velocity Distributions:Ejecta Velocity Distributions:
i t i [ yr ] Fit: r i , C , v 2 a i , C ( v ) [ pc ]
00[MB] 0.968 0.913 ± 0.023
i t i [ yr ] Fit: r i , C , v 2 μ i , C ( v ) [ pc ] σ i , C ( v ) Fit: r i , E , v 2 μ i , E ( v ) [ pc ] σ i , E ( v )
28 13.8 × 10 6 [LN] 0.989 1.638 ± 0.096 0.883 ± 0.047 [LN] 0.989 2.137 ± 0.063 0.478 ± 0.038
37 13.8 × 10 7 [LN] 0.985 1.00 ± 0.10 0.986 ± 0.051 [LN] 0.990 1.715 ± 0.079 0.693 ± 0.042
46 13.8 × 10 8 [LN] 0.965 0.17 ± 0.11 1.008 ± 0.054 [LN] 0.990 1.51 ± 0.10 0.881 ± 0.055
55 13.8 × 10 9 [LN] 0.825 0.00 + 0.25 0.00 1.34 ± 0.13 [LN] 0.989 1.21 ± 0.10 0.912 ± 0.051
64 13.8 × 10 10 [LN] 0.753 0.00 + 0.63 0.00 1.89 ± 0.44 [LN] 0.985 1.05 ± 0.12 0.949 ± 0.062
Table 15. Cluster and ejecta mean PBHs mean mass, position and velocity at selected time slices. These are extracted from the fitted Log-Normal and Maxwell–Boltzmann distribution parameters of Table 14 and by making use of Equations (6) and (11) of Section 3.2.
Table 15. Cluster and ejecta mean PBHs mean mass, position and velocity at selected time slices. These are extracted from the fitted Log-Normal and Maxwell–Boltzmann distribution parameters of Table 14 and by making use of Equations (6) and (11) of Section 3.2.
      Cluster Characteristic Profiles:Ejecta Characteristic Profiles:
i t i [ yr ] m ¯ i , C [ M ] x ¯ i , C [ pc ] v ¯ i , C [ km / s ] m ¯ i , E [ M ] x ¯ i , E [ pc ] v ¯ i , E [ km / s ]
00 29.92 ± 0.75 1.39 ± 0.16 1.462 ± 0.040
28 1.38 × 10 6 30.30 ± 0.84 2.70 ± 0.23 3.85 ± 0.13 26.90 ± 0.83 16.7 ± 1.6 9.50 ± 0.56
37 1.38 × 10 7 30.45 ± 0.94 10.0 ± 1.0 4.42 ± 0.38 29.14 ± 0.79 119 ± 16 7.12 ± 0.59
46 1.38 × 10 8 30.6 ± 1.0 47.0 ± 5.1 1.97 ± 0.24 29.33 ± 0.81 ( 1.100 ± 0.150 ) × 10 3 6.66 ± 0.75
55 1.38 × 10 9 31.5 ± 1.1 210 ± 20 2.51 ± 0.68 29.14 ± 0.84 ( 8.5 ± 1.1 ) × 10 3 4.15 ± 0.58
64 1.38 × 10 10 32.5 ± 1.2 940 ± 100 6.4 ± 5.5 29.10 ± 0.80 ( 8.3 ± 1.3 ) × 10 4 4.49 ± 0.58
Table 17. Shown are the cluster and ejecta BF values for the orbital eccentricity and semi-major axis distribution, ρ ( e i , j , a i , j , t ) , from Figure 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.
Table 17. Shown are the cluster and ejecta BF values for the orbital eccentricity and semi-major axis distribution, ρ ( e i , j , a i , j , t ) , from Figure 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.
      Cluster Objects Distributions:Ejecta Objects Distributions:
i t i [ yr ] a i , C , I   [ pc ] e i , C , I a i , C , B   [ pc ] e i , C , B a i , E , I   [ pc ] e i , E , I a i , E , B   [ pc ] e i , E , B
00[BF] 0.865 1.000 0.128 1.000 [BF] 0.226 1.49
28 13.8 × 10 6 [BF] 1.77 1.000 0.283 0.918 [BF] 2.73 1.00 0.0546 2.82
37 13.8 × 10 7 [BF] 7.14 1.000 1.16 0.916 [BF] 2.73 1.00 2.26 1.04
46 13.8 × 10 8 [BF] 30.3 1.000 8.43 1.000 [BF] 14.4 1.00 35.6 1.06
55 13.8 × 10 9 [BF]150 1.000 24.5 1.000 [BF] 81.8 1.32 150 1.00
64 13.8 × 10 10 [BF]786 1.000 105 1.000 [BF]200 1.32 744 1.07
Table 18. Nearby objects eccentricity and impact parameter distributions of Figure 16 distribution Log-Normal fits: ρ z , i , P = ( z σ z , i , P 2 π ) 1 exp ( ( Log Normal z μ z , i , P ) 2 / 2 σ z , i , P 2 ) 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 χ 3 2 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.
Table 18. Nearby objects eccentricity and impact parameter distributions of Figure 16 distribution Log-Normal fits: ρ z , i , P = ( z σ z , i , P 2 π ) 1 exp ( ( Log Normal z μ z , i , P ) 2 / 2 σ z , i , P 2 ) 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 χ 3 2 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.
      Relative Eccentricities Fits:Relative Impact Parameters Fits:
i t i [ yr ] Fit: r i , e 2 μ i , e σ i , e Fit: r i , b 2 μ i , b   [ pc ] σ i , b   [ pc ]
00[LN] 0.992 1.81 ± 0.14 1.559 ± 0.040 [LN] 0.856 0.0 ± 5.7 1.45 ± 0.14
28 13.8 × 10 6 [LN] 0.998 6.27 ± 0.31 2.342 ± 0.088 [LN] 0.509 0.0 + 3.0 0.0 1.83 ± 0.77
37 13.8 × 10 7 [LN] 0.998 6.57 ± 0.29 2.208 ± 0.083 [LN] 0.771 9.8 ± 3.8 4.0 ± 1.0
46 13.8 × 10 8 [LN] 0.999 6.74 ± 0.23 2.268 ± 0.068 [LN] 0.789 0.00 + 0.94 0.00 1.71 ± 0.24
55 13.8 × 10 9 [LN] 0.999 6.72 ± 0.23 2.166 ± 0.064 [LN] 0.975 3.43 ± 0.21 1.084 ± 0.054
64 13.8 × 10 10 [LN] 0.999 7.33 ± 0.25 2.330 ± 0.072 [LN] 0.987 5.55 ± 0.51 1.45 ± 0.15
Table 19. Shown are the number of bubbles, N P , segregated by the population (cluster and ejecta) type, in a Milky Way-like galactic halo if all DM is made of PBHs, as well as the volume fraction occupied by the spheres of influence of these bubbles in such galaxy, F P at selected times, and where P = C , E stands for the cluster and ejecta populations. Shown as well are the individual masses of such bubbles, M P ind [ M ] , and the total mass of all such bubbles in a milky Way-like galactic halo, M P tot [ M ] . Last, the average distance between individual PBHs in such bubbles, X P ind [ pc ] , are shown for the time slices in which the bubble’s spheres on influence have undergone sufficient expansion that they start to merge with each other, a point where the volume fraction holds F P 1 .
Table 19. Shown are the number of bubbles, N P , segregated by the population (cluster and ejecta) type, in a Milky Way-like galactic halo if all DM is made of PBHs, as well as the volume fraction occupied by the spheres of influence of these bubbles in such galaxy, F P at selected times, and where P = C , E stands for the cluster and ejecta populations. Shown as well are the individual masses of such bubbles, M P ind [ M ] , and the total mass of all such bubbles in a milky Way-like galactic halo, M P tot [ M ] . Last, the average distance between individual PBHs in such bubbles, X P ind [ pc ] , are shown for the time slices in which the bubble’s spheres on influence have undergone sufficient expansion that they start to merge with each other, a point where the volume fraction holds F P 1 .
Cluster Background Number:Cluster Background Properties:
i t i [ yr ] N C tot F C tot M C ind [ M ] M C tot [ M ] X C ind [ pc ]
00 ( 2.7 + 1.7 0.9 ) × 10 10 ( 5.3 + 4.0 2.7 ) × 10 9 ( 3.0 + 2.6 1.4 ) × 10 4 ( 8.0 + 5.0 2.6 ) × 10 11
28 1.38 × 10 6 ( 2.4 + 1.5 0.8 ) × 10 10 ( 3.8 + 2.9 1.8 ) × 10 8 ( 2.8 + 2.4 1.3 ) × 10 4 ( 7.4 + 5.0 2.5 ) × 10 11
37 1.38 × 10 7 ( 2.0 + 1.3 0.7 ) × 10 10 ( 2.0 + 1.5 1.0 ) × 10 6 ( 2.3 + 2.0 1.1 ) × 10 4 ( 6.1 + 4.0 2.3 ) × 10 11
46 1.38 × 10 8 ( 1.6 + 1.0 0.7 ) × [ 10 10 ( 2.1 + 1.6 1.0 ) × 10 4 ( 1.8 + 1.6 1.0 ) × 10 4 ( 4.9 + 3.2 2.0 ) × 10 11
55 1.38 × 10 9 ( 1.3 + 0.8 0.6 ) × 10 10 0.018 + 0.14 0.009 ( 1.5 + 1.3 0.9 ) × 10 4 ( 4.0 + 2.6 1.9 ) × 10 11
64 1.38 × 10 10 ( 9.7 + 7.2 6.0 ) × 10 9 1.6 + 1.3 0.8 ( 1.2 + 1.1 0.8 ) × 10 4 ( 3.2 + 2.2 1.8 ) × 10 11 112 + 29 23
Ejecta Background Number:Ejecta Background Properties:
i t i [ yr ] N E tot F E tot M E ind [ M ] M E tot [ M ] X E ind [ pc ]
28 1.38 × 10 6 ( 6.7 + 8.1 2.6 ) × 10 7 ( 9.0 + 7.1 3.8 ) × 10 6 67 + 89 34 ( 1.8 + 2.1 0.7 ) × 10 9
37 1.38 × 10 7 ( 6.7 + 5.2 3.2 ) × 10 9 ( 3.3 + 2.7 1.8 ) × 10 3 ( 7.0 + 7.1 4.2 ) × 10 3 ( 1.9 + 1.5 0.9 ) × 10 11
46 1.38 × 10 8 ( 1.12 + 0.81 0.39 ) × 10 10 2.6 + 2.1 1.5 ( 1.2 + 1.1 0.6 ) × 10 4 ( 3.2 + 2.3 1.3 ) × 10 11 108 + 29 17
55 1.38 × 10 9 ( 1.4 + 1.0 0.5 ) × 10 10 ( 1.21 + 0.99 0.67 ) × 10 3 ( 1.5 + 1.4 0.8 ) × 10 4 ( 4.1 + 2.9 1.6 ) × 10 11 99 + 26 14
64 1.38 × 10 10 ( 1.7 + 1.2 0.6 ) × 10 10 ( 1.1 + 1.0 0.7 ) × 10 6 ( 1.8 + 1.7 0.9 ) × 10 4 ( 4.9 + 3.4 1.8 ) × 10 11 93 + 24 13
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Trashorras, M.; García-Bellido, J.; Nesseris, S. The Clustering Dynamics of Primordial Black Holes in N-Body Simulations. Universe 2021, 7, 18. https://0-doi-org.brum.beds.ac.uk/10.3390/universe7010018

AMA Style

Trashorras M, García-Bellido J, Nesseris S. The Clustering Dynamics of Primordial Black Holes in N-Body Simulations. Universe. 2021; 7(1):18. https://0-doi-org.brum.beds.ac.uk/10.3390/universe7010018

Chicago/Turabian Style

Trashorras, Manuel, Juan García-Bellido, and Savvas Nesseris. 2021. "The Clustering Dynamics of Primordial Black Holes in N-Body Simulations" Universe 7, no. 1: 18. https://0-doi-org.brum.beds.ac.uk/10.3390/universe7010018

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop