Dynamics and merger rate of primordial black holes in a cluster

The PBHs clusters can be source of gravitational waves, and the merger rate depends on the spatial distribution of PBHs in the cluster which changes over time. It is well known that gravitational collisional systems experience the core collapse, that leads to significantly increase of the central density and shrinking of the core. After core collapse, the cluster expands almost self-similarly (i.e. density profile extends in size without changing its shape). These dynamic processes affect at the merger rate of PBHs. In this paper the dynamics of the PBH cluster is considered using the Fokker-Planck equation. We calculate the merger rate of PBHs on cosmic time scales and show that its time dependence has a unique signature. Namely, it grows by about an order of magnitude at the moment of core collapse which depends on the characteristic of the cluster, and then decreases according to the dependence $\mathcal{R} \propto t^{-1.48}$. It was obtained for monochromatic and power-law PBH mass distributions with some fixed parameters. Obtained results can be used to test the model of the PBH clusters via observation of gravitational waves at high redshift.

The question of PBHs clustering might play a decisive role in observational effects. It is known that various models predict formation of PBH clusters both in the postinflation epoch (an initial clustering) [36][37][38][39][40][41][42] and in the early epoch of structure formation due to Poisson fluctuations [27,[43][44][45]. This study considers the model of cluster formation due to collapse of domain walls produced as a result of quantum fluctuations of scalar field(s) at the inflation stage [36][37][38]. Following this model, a falling power-law mass distribution is assumed here in a wide PBH mass range. The cluster decouples from the Hubble flow and forms a virialized system at the redshift z f ∼ 10 4 [11,46]. The parameters of the resulting cluster depend on the specific model of scalar field potential. Typically, it has characteristics close to globular star clusters with a wide range of PBH masses. However, the observed signatures (e.g., the merger rate) of the PBHs cluster depends on its structure at a specific redshift. Hence, the cluster evolution is important. The dynamics of a PBHs cluster was considered in the works [22,47,48] with help of the N-body simulation. However, such calculations are only suitable for clusters with a small number of black holes. In this work, one uses the Fokker-Planck equation to study dynamics of the PBHs cluster. This approach allows practically studying an arbitrary range of the possible cluster parameters in a short computational time.
Additionally, note that the picture can be changed by including other processes. It might be caused by taking into account baryons or/and the dark matter which change the dynamics of the cluster directly or through accretion processes. However, these effects are the question of a separate investigation. Here, the effect of the pure gravitational evolution of the PBH cluster is considered.
In this paper, we show that during the evolution, the spatial distribution of PBHs is changed due to two-body relaxation processes and the mass segregation. We estimate the time dependence of the PBH merger rate inside the cluster and show that it does not depend on the initial parameters of the cluster and has the unique signature. The obtained results might be used to test the model of the PBH cluster with the help of a future generation of gravitational waves detectors which will be able to detect black holes mergers at high redshifts. The effect under consideration is one of the rare possibly observable effects inherent in models of PBH clusters that may include other effects such as gravitational microlensing on the cluster [49]. In addition, effects giving constraints on PBHs as DM and constraints themselves in case of PBH clustering should be reconsidered [11].

The Fokker-Planck equation approach
The orbit-averaged Fokker-Planck (FP) equation is often used to study evolution of globular star clusters [50][51][52][53] and galactic nuclei [54][55][56][57]. In our case, this equation could be used to describe time evolution of the distribution function of PBHs f (E) as a result of diffusion in energy space. Note, the distribution function f (E) depends only on energy (per unit mass) E = v 2 /2 + φ(r). The FP equation for a multi-mass cluster is given by [55] where the index i refers to the i-th mass type of PBHs, the last loss-cone term νf i describes the absorption of PBHs by the massive central black hole (CBH) [55], and the coefficients D E and D EE are where the sum is over all the types of PBH. Γ = 4πG 2 ln Λ, ln Λ is the Coulomb logarithm, and G is the Newtonian gravitational constant. The expressions for q(E) and p(E) are where φ −1 (E) is the root of the equation E = φ(r). φ(r) is the gravitational potential which is given by Poisson equation solution The last term is the potential of the CBH with the mass M • , and ρ(r) is the density profile The evolution of cluster is described by both Fokker-Planck equation (1) and Poisson equation (6). The diffusion coefficients (2) and (3) are calculated via distribution function from a previous time step; therefore, (1) is solved as linear equation at a time step ∆τ . Then, (6) and (7) are solved by iterations, and so on. This calculation algorithm can be found in more detail in the papers [51,55].
The evolution of PBH cluster is similar to the dynamics of star clusters and is driven by twobody relaxation. This leads to the effect of a core collapse [51], as a result of which the radius and the mass of the core (the central part of the cluster) shrink to zero r c → 0, M c → 0, while the density of core goes to infinity ρ c → ∞. This phenomena is also known as "gravothermal catastrophe" [58]. The time dependence of r c (t) and M c (t) are given by [59] where t cc is the core collapse time Note, Eqs. (8) and (9) do not have a physically reasonable limit at t → t cc . In real stellar clusters, a core collapse stops due to the following reasons: they are formation of binary stars [60][61][62], or the presence of primordial binaries [63,64], or the influence of a massive CBH [55,65] which acts as a heating source. After the termination of the core collapse, cluster enters to a nearly self-similar expansion stage.

The merging of primordial black holes
The PBH cluster is the great place for intense black holes mergers. The cross section of binary PBH formation due to emission of gravitational waves is given by [52,66] where v rel is the relative velocity of PBHs, m and m are the masses of PBHs, and c is the speed of light. It is typically assumed that the forming binary is immediately merged [17,19,47,52,67]. Possible formation of binary systems due to triple PBH interactions are not taken into account here and commented in the end of this section. The merger rate of black holes per cluster will be where ρ is the density of PBHs. It is assumed for the sake of simplicity of calculations that all PBHs in the cluster have the same mass m. For estimation, one can write that inside the cluster core ρ(r) = ρ c and v rel = GM c /r c , where M c and r c are the mass and the radius of the core, respectively. Outside the core, the density falls as ρ ∝ r −β (β ≈ 2.2 is the typical value for a globular star cluster). Then, the merger rate can be written as Note, the estimation (13) depends only on the parameters of the central part of the cluster, even if the total mass of the cluster is much greater than the core mass.
Using the expressions (8) and (9), one can get the time evolution of the merger rate: Moreover, here we assume that the clusters contain central massive black hole. Then, an additional channel for black hole mergers is added, namely the capture of less massive PBHs by the CBH. The growth rate of the CBH mass can be obtained aṡ where σ c = 4πr 2 g (c/v rel ) 2 is the capture cross section of particles by the CBH [68] and r g is the gravitational radius of CBH. Using (8) and (9), one can geṫ Here, both (14) and (16) diverge at t → t cc , that is not a physical result. As noted above, when the core collapse is terminated in real systems, the merger rate decreases due to expansion of the cluster. However, one should keep in mind that the presented analysis is simplified and does not take into account the presence of a PBH mass spectrum in a cluster. The calculations performed for the cluster with the PBH mass spectrum are presented in this study. Nevertheless, this result can be useful for probing the model of PBHs clusters if an enhancement of the merger rate will be observed. Observation of an enhancement of a merger rate through gravitational waves with increasing redshift would support PBH cluster hypothesis. The calculation details of the merger rate for the PBH cluster with a spectrum mass is described in [69].
It is important to note that PBH binaries can be also formed through three-body interactions [60]. However, the evolution of such binary systems cannot be studied within the framework of the Fokker-Planck equation. These binaries can both be destroyed in a cluster and become tighter systems with subsequent merging through emission of gravitational waves. The study of these processes requires N-body simulation and is outside the scope of this work.

Evolution of PBHs clusters 4.1 PBHs clusters with monochromatic mass spectra
In order to show main features of the cluster dynamics, we first consider the simplified model. In this section, the masses of all PBHs are equal to M , the CBH mass is 100 M , and the total mass of the cluster is 10 5 M . We take the initial density profile of PBHs in the cluster in the following form where ρ 0 is the normalization constant. Note, that initial cusp ρ ∝ r −γ is inevitable due to presence of CBH. The cusp must be steeper than r −1/2 , otherwise the initial distribution will not be physical [70]. The evolution of the PBHs density profile is presented in Fig. 1. It should be noted, that the figure shows differential mass distribution dM (r)/dr ∝ r 2 ρ(r). One can see, that the slope of the density profile ρ ∝ r −7/4 appears in the central region of the cluster after ∼ 0.4 Myr. This distribution is typical for gravitationally interacting bodies around a massive black hole and is known as the cusp of Bahcall-Wolf. Then, the core is formed (yellow curve) outside the region of influence of the CBH. The further evolution follows the path of the classical gravothermal catastrophe: the inner part of the cluster is shrunk, and the central density increases. The green curve shows the final moment of the core collapse. After ∼ 600 Myr, the stage of near self-similar expansion begins (see Fig. 1b). The evolution of mass-shells are presented in Fig. 2. One can see that after ∼ 600 Myr all mass shells expand according to the law r ∝ t 2/3 which indicates the self-similar expansion character.  Figure 1: The evolution of the PBHs mass distribution in the cluster before (Fig. 1a) and after (Fig. 1b) the self-similar expansion phase is shown.
The merger rate of PBHs is presented in Fig. 3. One can see that the modern mergers are almost independent of the initial structure of the cluster (for the considered cluster mass). It is interesting to note that after the core collapse, all curves merge into one. It can be seen that the time evolution of the merger rate has a characteristic maximum at the moment of core collapse, see Fig. 3a. However, for cluster with a "steeper" density profile at the center, this maximum is not so pronounced (the case γ = 2). In addition, it should be noted that the capture rate of CBH has two peaks, see Fig. 3b. The second peak is similar to the peak shown in Fig. 3a, while the first one is caused by formation of the Bahcall-Wolf cusp. The asymptotic behavior of the merger rate is different for various masses of PBHs. Namely, the merger rate of approximately equal PBHs masses behaves as R ∝ t −1.48 . This dependence is very easy to understand due to the fact that R ∝ r  Thus, the study of the PBHs merger rate on cosmological time scales can be used to test the hypothesis of PBHs cluster existence. If future generations of gravitational waves experiments will detect the decreasing merger rate with decreasing redshift according to the obtained dependencies, it will be an indirect evidence of the existence of PBHs clusters. Moreover, if the dependence R ∝ t −1.48 is observed, it will indicate that PBHs clusters do not have a massive central BH.
In this case, the expansion is associated with formation of binary systems, similar to globular star clusters. Note, that the clusters can have different parameters and, respectively, positions of maxima of merger rate at the redshift axes. Integration over all the clusters should change the expected merger rate dependence from z. However, the asymptotic behavior at large times will follow the obtained dependence.

PBHs clusters with wide mass spectra
Cluster may contain PBHs distributed by masses in a wide range as it is predicted, e.g., in the model of the domain walls collapse [11]. This fact leads to a shift in time scales due to mass segregation and, as a consequence, to acceleration of the core collapse. In this section, we present calculation of PBHs cluster evolution with the mass spectrum dN/dm ∝ m −2 [11] and the mass range from 10 −2 M to 10 M with the CBH mass 100 M and the total cluster mass 10 5 M . The parameters of the density profile (17) are γ = 1, β = 5, r 0 = 1 pc. The evolution of the mass distribution of PBHs is shown in Fig. 4. One can see, that the dynamics are similar to a cluster with initially monochromatic PBH mass distribution, but the density profile is near to ρ ∝ r −2 . The presence of the mass spectrum of the PBHs leads to the mass segregation as result of which more massive BHs settle to the center cluster, see Fig. 5. In addition, the wide range of PBHs masses shifts relaxation timescales. The cluster enters the expansion phase after ∼ 30 Myr while for the cluster with monochromatic PBH mass distribution this phase occurs after ∼ 600 Myr. This fact is caused by the formation of a subcluster from heavy PBHs (which sink in the central part of the cluster due to dynamical friction) in the central region of the "host cluster" which evolves faster. The expansion of the cluster is also self-similar with the law r ∝ t 2/3 . In this work, we do not take into account the external conditions for cluster. At the stage of the structure formation, clusters become a part of the DM halo. This leads to the fact that the cluster is influenced by tidal forces as a result of which it cannot expand indefinitely. The outer layers will be captured by the gravitational potential of the host halo. Therefore, the cluster may be partially stripped. However, as noted earlier, the obtained results are valid in order of magnitude for the central part of a cluster. Hence, even if cluster has lost most of its mass, the obtained merger rate remains correct.
One can use the obtained merger rate in order to restrict the abundance of PBHs clusters. The modern merger rate of BHs evaluated from the LIGO/Virgo data is R ∼ 10 Gpc −3 y −1 [71] while, for instance, the merger rate for 10 M PBHs per the cluster is R cl ∼ 10 −12 y −1 . From the condition ρ crit Ω cl R cl /M cl R, one can obtain the density fraction of the PBHs cluster with the mass 10 5 M composed of PBHs with masses from 10 −2 M to 10M distributed as dN/dm ∝ m −2 with CBH mass 100M as Ω cl 0.01. The influence of binaries formed as a result of three-body interactions or primordial binaries formed before the formation of a cluster can have their contribution to the limit. Note that the limit on Ω cl from GW observation can be escaped by the choice of other parameters.

Conclusion
In this paper, the dynamics of the primordial black holes cluster was considered. We estimated the time dependence of the merger rate of PBHs in the cluster with the different parameters and showed that the cluster evolution significantly affects the merger rate. In addition, it was shown that the merger rate changes with time as R ∝ t −1.48 . Hence, the observation of GWs at high redshifts may shed light on the spatial distribution of PBHs in the Universe. The merger rate was considered for the PBHs cluster with both the narrow (monochromatic) and the wide mass distribution of BHs. In the first case, the merger rate has a peak at some redshift due to the core collapse. In the second case, the presence of mass distribution of PBHs inside a cluster leads to the settling of heavy BHs in the cluster center, while the lighter components tend to form outer layers of a cluster. Nevertheless, the main features of the evolution are the same as for the cluster model with the monochromatic mass distribution. We also showed that gravitational waves data constrain clusters with the total mass 10 5 M composed of PBHs with masses from 10 −2 M to 10M , distributed as dN/dm ∝ m −2 with CBH mass 100M as Ω cl 0.01. The last limit can be escaped by choosing other PBH cluster mass parameters. The search for similar effects as considered here or effect of gravitational microlensing on the cluster can allow testing the cluster model.

Funding
The work of V.D.S. was supported by a grant of Russian Science Foundation No 18-12-00213-P https://rscf.ru/project/18-12-00213/ and performed in Southern Federal University (SFEDU). The work of A.A.K. was supported by the Ministry of Science and Higher Education of the Russian Federation, Project "Fundamental properties of elementary particles and cosmology" № 0723-2020-0041.