Clusters of Primordial Black Holes

The Primordial Black Holes (PBHs) are a well-established probe for new physics in the very early Universe. We discuss here the possibility of PBH agglomeration into clusters that may have several prominent observable features. The clusters can form due to closed domain walls appearance in the natural and hybrid inflation models whose subsequent evolution leads to PBH formation. The dynamical evolution of such clusters discussed here is of crucial importance. Such a model inherits all the advantages of uniformly distributed PBHs, like possible explanation of supermassive black holes existence (origin of the early quasars), the binary black hole mergers registered by LIGO/Virgo through gravitational waves, which could provide ways to test the model in future, the contribution to reionization of the Universe. If PBHs form clusters, they could alleviate or completely avoid existing constraints on the abundance of uniformly distributed to PBH, thus allowing PBH to be a viable dark matter candidate. Most of the existing constraints on uniform PBH density should be re-considered to the case of PBH clustering. Furthermore, unidentified cosmic gamma-ray point-like sources could be (partially) accounted. We conclude that models leading to PBH clustering are favored compared to models predicting the uniform distribution of PBHs.


Introduction
The idea of the Primordial Black Holes (PBHs) formation was predicted five decades ago [1]. Then this scenario was grounded and developed in [2][3][4] and in many other works. The PBHs have not been identified in observations, although some astrophysical effects can be attributed to the PBHs, e.g., supermassive black holes in early quasars. Therefore till now, the PBHs give information about processes in the early Universe only in the form of restrictions on the primordial perturbations [5] and on physical conditions at different epochs. There are also hopes that the evidence about PBHs existence will be found soon. It is important now to describe and develop in detail models of PBHs formation and their possible effects in cosmology and astrophysics.
A new wave of interest to PBHs has been induced by the recent discovery of gravitational waves from BH mergers, raising discussions of PBHs with mass 10-100M as dark matter (DM) candidate (see, e.g. [6][7][8][9][10]). This is still a controversial question (e.g. [11][12][13][14][15][16][17][18]), to which we pay attention below. PBHs as DM at other mass values (M ∼ 10 17 g, M ∼ 10 20−24 g) are still topical (see e.g. [19] and references given below). The search for PBHs as the DM is one of the hot-topical issues which gets a new power in the context of PBH cluster models.
There are several models of PBHs formation. PBHs can be formed during the collapses of adiabatic (curvature) density perturbations in relativistic fluid [20]. They can be formed also at the early dust-like stages [21,22] and rather effectively on stages of a dominance of dissipative superheavy metastable particles owing to a rapid evolution of star-like objects that such particles form [23,24]. There is also an exciting model of PBHs formation from the baryon charge fluctuations [25,26].
Another set of models uses the mechanism of domain walls formation and evolution with the subsequent collapse [27][28][29]. As we show below, quantum fluctuations of a scalar field near a potential maximum or saddle point during inflation lead to a formation of closed domain walls. After the inflation is finished, the walls could collapse with black holes in the final state. There is a substantial amount of the inflationary models containing a potential of appropriate shape. The most known examples are the natural inflation [30,31] and the hybrid inflation [32] (and its supergravity realization [33]). The landscape string theory provides us with a wide class of the potentials with saddle points, see review [34] and references within.
In this review, we rely on the models that lead to closed wall formation. To figure out the idea, let us assume that the potential of a scalar field possesses at least two different vacuum states. In such a situation there are two qualitatively different possibilities to distribute these states in the early Universe. The first possibility implies that space contains approximately equal amounts of both vacuum states, which is typically the case at the usual thermal phase transition. It is supposed that after the phase transition, the two vacuum states are separated by a vacuum wall (kink). This leads to the wall dominated Universe that sure contradicts observational data. The alternative possibility corresponds to the case when the two vacuum states are populated with different probabilities. In this case, islands of a less probable vacuum state appear, surrounded by the sea of another, more probable vacuum state. The initial distribution of the scalar field which breaks a symmetry between the vacuum states is the essential condition for such an asymmetric final state.
In the framework of our model, the formation of such islands starts at the inflationary stage. Inflationary fluctuations permanently change the value of a scalar field in different space regions leading to the asymmetrical field distribution at the end of inflation. As a result, a distribution of the scalar field is already settled to the end of inflation leading in the future to the islands of one vacuum in the sea of another one. The phase transition accompanied by the walls formation takes place only after the end of inflation, deeply in the Friedman-Robertson-Walker (FRW) epoch. Therefore, the distribution of the scalar field which was formed during the inflation allows the formation of closed walls in the FRW stage. At some instant after crossing the horizon, such walls become causally connected as a whole and begin to contract because of the surface tension. As a result, the energy of a closed wall may be focused within a small volume inside the gravitational radius which is the necessary condition for a black hole creation.
The mechanism of the wall formation is quite general and can be applied to those inflationary models containing extrema of their potential. Some of them where developed in [35][36][37][38][39] where the potentials with saddle points were considered. It has been revealed there that the PBH mass spectrum and PBH ability depend strongly on the potential parameters and an initial field value, see discussion in Sect. 4.1. The PBH mass can vary in order of magnitudes. Depending on the chosen parameters, this mechanism can be applied to the description of the dark matter in the form of small BHs or the first quasar formation for massive BHs. Moreover, the cluster distribution of PBH is realized in a natural way.
The idea of PBH clusterization appears to be rather promising [40][41][42]. Indeed, such structures could explain the early quasars observations [43], contribute to the dark matter, produce gravitational waves bursts, be a reason for reionization and point-like cosmic gamma-ray sources.
Another model of the PBH cluster structure was also suggested in the framework of mechanism of black hole formation from strong primordial density inhomogeneities, see [44][45][46][47][48][49][50], and discussion in [46]. The idea of the research is based on the assumption that density peaks in the primordial density fluctuations are clustered more intensively than the ordinary density fluctuations [51]. The latter may imply a deviation from the Poisson distribution of the created PBHs in space. However, the Poisson effect is also supposed to hold for the PBH clusterization. Besides that, a seed effect, when PBH attracts surrounding matter along with other PBHs, is considered to be the possible independent reason for clusterization. The last two effects are used in the series of papers mentioned above for the large-scale structure formationgalaxies and their clusters (see also [52]).
Heating of the surrounding matter is the inherent property of the domain wall mechanism of PBH (cluster) formation. While collapsing, the domain wall partially transfers its kinetic energy to the ambient matter. It could allow discriminating between the models by observations.
The idea of PBH clustering gives a natural explanation of the origin of massive binary BH. Recent indication [53] of the possible existence of ∼ 10 4 black holes with a nearly solar mass within one parsec around galactic center could also be interpreted in the given cluster framework.
There are three aims of this review. Firstly we describe the mechanism that leads to the formation of the PBHs clusters. Secondly, we discuss the space and mass evolution of PBHs in such cluster provided that initial PBHs distribution is assigned. The third aim is to discuss observable effects caused by the PBH clusters. The discussion of cluster evolution and observational implications, conducted in Sects. 3, 4, 5, being rather general, is illustrated by a specific model described in Sect. 2.4.

PBH formation
One of the reasons for primordial black holes formation is phase transitions during and after inflation. This is threestep process. As discussed in the Introduction, closed walls connecting two vacuum states appear after the inflation is terminated. The shape of most of them is far from beeing spherical. At the second step the walls evolve in the plasma of relativistic particles that are produced by the reheating stage of the inflation. Internal tension leads to decreasing the surface square transmitting its potential energy into the kinetic energy of the wall. In its turn, the kinetic energy is also decreasing due to friction of the surrounding matter. As the result, the wall minimizes its surface energy by acquiring spherical shape and shrinking. This step is highly model dependent. To proceed, we suppose that the wall acquires the spherical shape soon after it crosses the horizon. The third step consists of wall shrinking due to its internal tension. It is finished by a black hole creation if the gravitational radius of the wall is larger than the wall's width.
One can conclude that production of closed walls lead to strong inhomogeneities able to collapse into PBH. Multiple quantum fluctuations of a scalar field result in PBHs clusters formation, with a scalar field not necessarily being an inflaton. In the review we will focus basically on this scenario of PBHs formation discussed in [35,37,54].
The mechanism of the PBH cluster formation is represented in sections A, B, C in general form so that it can be used in a wide class of potentials. In section D we perform numerical analysis on the basis of specific model. This model is also used in Sect. 3. It is shown there that the mass spectrum of cluster can be varied in wide range depending on Lagrangian parameters. Starting from Sect. 4 we discuss different applications for the PBH clusters characterised by different mass spectrum.

Multiple fluctuations on inflationary stage
Consider a scalar field φ with potential possessing at least one maximum or a saddle point. As a result of inflation the observ-able Universe is formed from a single causally connected region with a characteristic size of H −1 . Let our Universe be nucleated with an initial field value φ u close to a potential maximum or saddle point. Due to the classical motion, most of the space will be finally filled with a field lying in a specific minimum. Nevertheless, quantum fluctuations during the inflation gradually bring the field over the potential barrier in individual space regions. After inflation, classical field evolution in such space regions rolls the field down to another minimum, resulting in the formation of a wall between such domain and the outer Universe. For future discussion, let us denote as a set of those field values φ c , (i.e. φ c ∈ ) that lead to the formation of the space islands (domains) after the inflation is finished.
Multiple quantum fluctuations during inflation can be described as random walks [55]. One can solve the Fokker-Planck equation for the field φ [56] and get f (φ, t) -field φ values distribution by the moment of time t. Hereinafter we assume by definition the relation N ≡ Ht of the number of e-folds N to the time t passed from the beginning of inflation. Neglecting the shape of the potential (for more detailed calculations see [37,[56][57][58][59]) and considering φ u as an initial field value at the beginning of inflation one can get the probability density of finding field φ at an arbitrary space point The formula (1) is valid for the field φ, a mass of which can be neglected in comparison with the Hubble parameter during inflation, [60]. The corresponding formula for m = 0 was obtained in [61,62]. The probability of finding the field within the region by the moment of time t is: where φ cr is the boundary field value, separating the region in a case of one-field potential. The space of the size of the Hubble radius is divided into n u (t) = e 3Ht causally independent patches by the moment t. Hence, the number of critical fluctuations (or the number of critical regions) is: The size of fluctuations is limited by the size of the horizon H −1 . Hereinafter we assume the number of critical regions to be much less than the number of all causally independent regions. Otherwise walls abundance would contradict the observable rate of the Universe expansion and the CMB data.

Distribution of regions over the size
Let the critical fluctuation takes place at the moment of time t from the beginning of inflation. After inflation, its size will be: where N inf 60 is the total number of e-folds needed for the formation of the visible Universe. After the inflation is finished, the field tends to different minima in different areas.
Main part of the Universe contains the field in a specific minimum while the regions filled by the critical fluctuations form the domains of another vacuum separated from the surrounding space by walls. As has been shown in [63,64], the walls are not spherical just after formation. On the contrary, the fractal-like surface is much more natural for this mechanism. Here we suppose that the energy of the wall fluctuations is quickly transferred to matter that leads to spherical walls surrounded by heated matter. Closed wall starts to collapse after crossing its horizon which is expanding on the RD-stage as r hor (τ ) = 2τ , where τ is the period of time passed after the beginning of the RDstage. The domain wall itself expands as r (τ ) = r inf √ τ/t inf on this stage, where t inf = N inf /H . By defining the moment of time τ when the horizon reaches the size of the wall, and substituting it into r (τ ) we get the size of the wall at the beginning of the collapse: The expression (5) relates the moment of the domain formation during inflation t to the final size of the domain r . By substituting t ≡ t (r ) into (3) we get the distribution of walls over the size n c (r ) = P t (r ) e 3Ht (r ) .

Distribution of regions over the mass
The fluctuating mechanism discussed here leads to the PBH formation caused by the closed walls collapse. The PBHs mass spectrum is the result of complex processes that take place during inflation, just after the inflation and is strongly model dependent. Indeed, walls are created having the form far from the spherical one [35,36,65,66]. Their small inhomogeneities are damped transferring the energy to the surrounding media and heating it. The larger perturbations of wall surface survive by the moment t h of the horizon crossing. The nonspherical collapse with the BHs in final stage has been discussed in [65,66].
For the spherical walls, there are three regimes of walls evolution depending on their surface energy density μ. If the wall energy is too small, it cannot be concentrated to form the black hole (see formula (13) and discussion around it). If the wall is too massive ("supercritical" in terms of [64]), inequality holds and according to [64], such a wall leads to the creation of baby universes with wormholes inside them. For an exterior observer, a black hole is formed. Here t h is the time scale when the horizon reaches a wall size. A wall for which t h < t μ is called "subcritical". Its gravitational field can be neglected before horizon crossing the wall scale. If its form is near spherical, it shrinks to a size smaller than the corresponding Schwarzschild radius, forming a black hole with ordinary internal structure. The final masses of the black holes depending on the initial data and parameters of the model are also analyzed in [63,64,67].
For the qualitative analysis, we suggest the case of the subcritical spherical walls. Also, we consider the case of the field potential with the energy-degenerated minima.
Radius of a spherical wall with the mass m is r (m) = √ m/4πμ, where μ is the surface energy density. The latter depends on specific form of the potential. By substituting r ≡ r (m) into (6) we get the distribution over the mass: The expression (8) is the direct consequence of expression (3). It gives the total number of critical regions that have arisen by the time t. Therefore expression (8) denotes the total number of PBHs with masses greater than m -the integral distribution. Differential distribution can be obtained by differentiating expression (8).
2.4 Model with specific shape of the potential Let us illustrate the speculations above by a specific potential satisfying the necessary conditions -the presence of extremum. Consider a complex scalar field with the potential where φ = re iθ . This field coexists with an inflaton field which drives the Hubble constant H during the inflationary stage. Last term in (9) reflects the renormalization effects to the Lagrangian (see details and refs in [68]) and has been used in the model of the baryon asymmetry observed in the Universe [69]. This term is assumed to be negligible at the inflationary stage. The parameter appears as a result of the instanton effects and renormalization, so that its value cannot be large and we assume that H, f . The second term in (9) plays a significant role in the post-inflationary stage, when the Hubble parameter decreases with time (e.g. H = 1/2t during the radiation dominated stage, H = 2/3t during the matter dominated stage).
When the inflation is finished, the field is captured in the potential minima φ = f √ 2 e 2π ni , n = 0, 1, 2, . . .. If neighbour space areas are filled by different vacuum states, say θ = 0 and θ = 2π , the second term in (9) leads to the effective Lagrangian for the dynamical value χ = f θ . The one-dimensional kink solution for the cosine potential is well known [70]: It is assumed throughout the paper that the wall width d is much smaller than the wall size. Therefore the wall is almost plane for an observer near the wall and we can apply formula (11) for describing the distribution of phase across a plane border separating two areas with different vacuum states. Knowledge of solution (11) allows one to find the wall's surface energy density Energy of this field configuration is concentrated in the plane of width d = 2 f / 2 forming a wall with the center at z 0 . The domain size immediately after the end of inflation markedly exceeds the horizon size at the FRW expansion stage. The overall contraction of the closed wall may begin only when the horizon size will be equal to the domain size R w with total energy E w ∼ 4π R 2 w μ. Up to this moment, the characteristic domain size increases with the expanding Universe. Evidently, internal stresses developed in the wall after crossing the horizon and external friction leads to minimization of the wall's surface. This implies that, having entered the horizon, the wall tends, firstly, to acquire a spherical shape and, secondly, to contract toward the center. The wall contracts up to the minimal size of the order of the wall width d. A black hole can be formed if this scale is smaller than the gravitational radius of a wall. Following [71], the black holes are formed under condition Here the Planck mass is denoted as M Pl . Indeed, if the width of the wall exceeds its gravitational radius, the energy of the wall can not be concentrated within it. The wall size d depends only on the physical parameters of the Lagrangian while the total wall energy can be macroscopically large. The inequality (13) gives lower limit for the black hole mass. We choose the following values of the model parameters, which agree with the data on the observable CMB anisotropy: H = 10 13 GeV, N inf = 60, field parameter = 0.05GeV, and initial value θ u = 0.05π . The distribution is very sensitive to the parameter values.
It can be seen ( Fig. 1) that number of PBH depends crucially on the vacuum expectation value f . Moreover, the number of PBH is affected significantly by the initial field value θ U : the closer is the initial value to the critical value π -the more PBH of all masses will be in the Universe. The wall's density, which depends on a value of the parameter , also affects how massive PBH can be formed, but this affection is not as significant as the changing in θ U or f . Domains with low masses do not collapse into black holes because their gravitational radius is less than the thickness of their wall, see inequality (13). With the chosen model parameters values, the border is located at the mass of ∼ 10 −2 M as is shown in Fig. 1. The considered mechanism of PBH formation leads to nontrivial situation when massive PBH can exist along with total absence of lower masses PBH. Too large walls begin to dominate locally in the Universe at the RD-stage, which causes them to stop expanding with a subsequent collapse into the PBH as follows from (7). Therefore, there is an upper limit on the size of the walls and, as a consequence, on the mass of the PBH. With the chosen model parameters values the upper border is located at the mass of ∼ 10 4 M .

PBH cluster formation: internal structure of cluster
As has been shown in [35] black holes are created in a company with substantial amount of smaller masses BHs. These walls form their common gravitational well where they merge to create BHs of larger masses. The subsequent evolution of such PBH cluster depends on its space/mass distribution, Lagrangian parameters and initial conditions.
One has to conclude that the formation of the PBH cluster is a complicated multi-step process. In this review, we will base on the assumption that a part of the closed walls are able to get rid of the surface perturbations and acquire a spherical shape.
Additional remark is necessary. The mechanism of closed wall formation is based on the quantum fluctuations of a scalar field near maximum of the potential V (φ) where the Fig. 1 PBHs mass spectrum for the visible Universe: f = 1.4H line represents realistic PBH distribution. In this case, the number of PBHs with a mass ∼ 10 2 M corresponds to the number of galaxies in the observable Universe. The role of other lines is to demonstrate a strong dependence of the distribution on the parameter f potential derivative is small. This means that the classical motion of the field is also slow. At the same time, the energy density fluctuations δρ/ρ ∝ 1/|φ|. The immediate conclusion is that the energy density fluctuations increase rapidly in the area of interest that could lead to another mechanism of BH formation acting simultaneously with the main one.
In fact, the model [72] is based on this mechanism.
As it has been told earlier, critical region will be formed with higher probability due to multiple fluctuations instead of direct tunneling. Hence, the closer the field value inside the space region-predecessor is to the critical value, the higher is the probability of the critical region formation. Thus, the probability of new PBH formation around the existing one (the parent PBH) increases. PBHs form clusters with fractal structure. Subsequent evolution of these clusters lead to the merging of black holes into more massive ones.
PBH mass and distance distributions inside the cluster Fractal structure of clusters implies that around each PBH there is an aggregate of smaller PBHs, around each of which -an aggregate of even smaller PBHs and etc. To study space and mass structure of a cluster let us choose some «seed» (parent) PBH with the mass M 0 and consider a cluster of smaller PBHs around it. Let us fix a size r cl of this region by the end of inflation. Then according to (4), the moment of time when this region was formed is t cl = t (r cl ). We also need the moment of time when the main PBH was formed. It can be derived from (5): t 0 = t (r (M 0 )). As long as such a PBH should be the only one to form, during the period of time t 0 − t cl only one fluctuation should occur. One can write down the following equation by using (3): Condition (14) connects the size of a chosen region r cl to the initial field value φ inside it in the way that only a single PBH with the mass M 0 forms in this region after inflation. By solving it one can get the dependence φ(M 0 , r cl ). Now we can calculate the spatial-mass distribution of PBHs inside specific region. It is described by the expression derived from (3) by substituting the initial field value φ(M 0 , r cl ) and time t − t cl : Thus, by defining the mass of a «seed» PBH one can obtain the distribution of the rest PBHs in the region r cl over the mass. It is worth noting that similarly to (8) this distribution is integral, i.e. gives the number of PBHs with masses higher than m inside the region r cl . Considering (15) as a function of r = r cl one can obtain space-mass distribution.
Consider as the example specific PBH cluster with a «seed» black hole of the mass 10 2 M . The calculations were performed on the basis of model (9). Spatial and mass PBH distributions inside the cluster are shown in Figs The distribution of PBHs just after their detachment from the Universe expansion (at z ∼ 10 4 ). Characteristic cluster radius is ∼ 1 pc. After detachment, the central region appears to be more dense (though it is not clearly seen in the figure) while the total number of the PBHs inside the cluster remains almost the same sequential detachment of the spherical layers of the cluster from the Universe expansion. The procedure and theory of the calculation is discussed in the next section.

Internal dynamics of PBHs inside a cluster
Previous sections were devoted to the description of the mechanism of PBH formation as the result of walls collapse. The discussion was illustrated by specific model (9) with the parameters leading to small masses of PBH. The latter could be used in the discussion on the dark matter origin. As is shown in this section, the same mechanism could be applied for the problem of early quasars. To succeed, the PBH should be made more massive.
The clusters with massive PBHs can be obtained in the framework of not only the same mechanism but even with the same specific model (9). The only thing that should be altered is the model parameters. We demonstrate this in our first paper [29]. In this section, we alter significantly the parameter only -from = 0.05 in previous sections to = 5 in this section.
It will be shown that the virialized cluster contains a central supermassive black hole that can serve as the seed for an early quasar. The origin of the dark matter is not discussed in this section. The model for explaining both the dark matter and the early quasars on the basis of the same approach is elaborated.

Cluster detachment from the Hubble flow
Let us consider a detachment of large cluster from the mean Hubble flow. In this case the result is more clear and affects the dark matter distribution. We describe the gravitational dynamics of some particular spherical layer, after the moment when this layer has entered under the cosmological horizon, i.e. its radius r < ct. In addition to BHs, there is dark matter inside and around the cluster of BHs, provided that BHs do not constitute all the dark matter. The cluster represents the initial density perturbation (of the entropy or isocurvature type) which determines the subsequent evolution of this composite system, and this perturbation competes with ordinary inflationary perturbations which produce CMB anisotropy and most of the large scale structures as in the standard cosmological scenario. The initial mass profile of the BHs cluster M h (r i ) is shown in Fig. 4 in comparison with dark matter mass profiles M DM (r i ) inside and around the same cluster.
Let us consider the spherical shell with the total BHs mass M h inside the sphere of radius r , radiation density ρ r , dark matter density ρ DM , and Lambda-term density ρ . The latter quantity is important only at the nearest epoch z ∼ 1 and can be neglected at the early stages z 1 of the Universe evolution. The radiation inside the cluster is homogeneous with high accuracy, and its gravitation can be taken into account by the substitution ρ r → ρ r + 3 p r /c 2 = 4ρ r /3 into Newtonian equations, and analogously for the Lambla-term. The evolution equation for the shell has the form Following [73] we use the parametrization where ξ is the comoving radius, a(t) is the scale factor of the universe, and the function b(t) describes the inhomogeneous evolution. The DM mass inside the shell is By using the Friedmann equation the (16) can be rewritten as where and δ h = M h /M DM . At the limit ρ → 0 the Eq. (18) is analogous to the equation obtained in the work [73] for the case of axionic dark matter. But we consider much more higher densities (δ h > 10 4 ) where the approximate fitting solutions of [73] are not applicable, so one should solve Eq. (18) numerically. We start the numerical solution from the moment of the horizon crossing of each particular shell with the initial conditions shown in Fig. 4. Some particular shell stops expanding, thenṙ = 0 (db/dz = b/(1 + z)) at some radius r = r s . After the contraction from r s to r c = r s /2 the shell virializes and fixes finally at the radius r c . Therefore the mean density ρ of the formed cluster is 8 times its density at the moment of maximum expansion and the virial radius is equal to The shells detachments start from the center of the cluster. At the very early stage of evolution the central mass M c 1.1 × 10 4 M have occurred under the Schwarzschild radius. This is the mass of the central BH which can grow further due to the smaller BHs capture and baryons accretion.
The shells with δ h > 1 (M DM < M h ) detached at the radiation dominated stage, but the remaining shells detached later at the matter dominated stage. For the outer shells the dark matter has the principle influence on the gravitational dynamics.
The dark matter accompanies the BHs shells detachments with the formation of the density profile in analogy with secondary accretion mechanism. But the density profile does not follow the usual secondary accretion low ρ ∝ r −9/4 due to the non-compactness of the seed mass where M DM (r c ) is obtained from the solution of (18). This density profile is shown in Fig. 5. The virialized region is including more and more DM shells until some outer shell is distorted by the gravitation from the surrounding density perturbations of inflationary origin. The equality of the counteracting forces from the total inner mass and from the outer perturbations determines the final radius of the system. The accretion of the gas onto the central BH can explain the early activity. And at the present epoch these objects look like dense galaxies with dormant BHs at the centers.

Internal dynamics after detachment
The study of the gravitational dynamics of PBHs clusters after detachment reduces to the problem of N -body simulation that can be solved only numerically. For these purposes we have used the NBODY6++ code [74], aimed at the simulation of globular stars clusters. We have made some changes to the program: all bodies have been deprived of their stel- Fig. 6 The typical spatial black holes distribution within cluster, the result of numerical simulations. Black holes «sizes» show its mass distribution but do not correspond to real ones lar characteristics (luminosity, types of star and their evolution, metallicity, etc.), objects' «sizes» have been redefined according to their gravitational radius (r i = r g i = 2G M i /c 2 , where r g is the gravitational radius), merging mechanism has been added.
Here we present the first results of the N -body simulation of the PBHs cluster evolution. In the numerical simulation, we included only BHs with M > 0.1M due to limited computer resources. Such BHs define the main inhomogeneities in the density structure (see Figs. 3b, 6). Further, PBHs with other masses will be taken into account as far as possible since they can give noticeable dynamical friction effects and modify mass distribution due to changing merger process probability.
Let us consider a typical example of simulation with N = 10 4 initial black holes in the cluster and the mass distribution The spatial distribution is chosen in the form Here R = 1 pc is the cluster size. The initial velocities were chosen according to the Maxwell distribution with the virial velocity v 0 ∼ 1 km/s. We assume the most massive black hole to be in the center of the cluster with zero initial velocity.
The PBH distribution in Fig. 3b is approximated by expressions (23) (23) is typical of the considered model with simple potential (9). The mass distribution within the cluster slightly changes due to mergers of PBHs and their escape from the cluster.
The results are presented in Fig. 7. One can see that cluster as a whole does not change its structure significantly: the mass spectral index is not changed essentially (α f ≈ −2.4), the number of particles and the total cluster mass is reduced by ∼ 15% . The cluster is not destroyed and remainsa gravitationally bound system by the moment z ∼ 20 − 40, when, according to estimates, the accretion starts influencing considerably the growth of the most massive black holes (see [76]). The numerical calculations are to be further developed in this field. Existing N -body simulation codes like NBODY6++ should be adapted for these purposes. In particular, modeling the dynamics of massive black holes is still quite problematic.
Notice that the total mass of PBH cluster is enough to form the early quasar (at z > 6). The problem of the early quasar origin is widely discussed. The simulations show that early quasars formation is possible in the case if by the moment z ∼ 20 − 40 there had already existed their «seeds» -either «abnormal» baryonic objects that further collapse into black holes or black holes themselves with masses 10 3−5 M (such masses are required to get quasars with mass M ∼ 10 9 M , which are mostly observed at z ∼ 6 [76]). However, the formation of massive «seeds» at such high redshifts requires additional explanation. Our analysis shows that we can offer one more scenario for the natural formation of these «seeds».

Observable properties
As already mentioned, the PBH clusters could be a clue to the solution of different problems in cosmology and astrophysics, identified from different observations. These are dark matter (DM) problem, supermassive black holes (SMBH) origin, black hole mergers registered by LIGO/Virgo, also the origin of unidentified cosmic gamma-ray sources, reionization of the Universe and maybe others. Here we list them as possibly applied to the PBH cluster model.

Supermassive black holes
The fact that each galaxy contains a supermassive black hole is almost no longer in doubt (see, e.g., the latest results of Chandra experiment [77]). Standard mechanism of SMBH formation requires a quite long time to reach so large black hole mass and is hardly consistent with SMBH existence at the redshift z 4 − 5 [76,78]. However now a few dozens In PBH cluster mechanism, such quasars can be formed as a result of accretion of seed PBH [26,43,62,[80][81][82]. The mass of the seed PBH is defined by cluster's formation considered above and its subsequent evolution.
As was said just above, getting PBH of the mass M ∼ 10 3 M at the redshift z ∼ 20 − 40 could give rise to SMBH of M ∼ 10 9 M at z ∼ 7 due to accretion in critical regime (see, e.g. [76,82]).
A hypothesis on the primordial origin of SMBH is suggested to be probed in the measurement of 21 cm hydrogen line [83]. The observational effect should be induced by the accretion of surrounding matter onto PBH during the Dark Age epoch.
Existent correlation between X-ray and infrared backgrounds can be treated in favour of massive PBH as a seed of star formation at high redshift (z ∼ 20) [10].
If a cluster contains many PBHs of intermediate mass, it can give the opportunity to check or limit the model with the help of LIGO and LISA [84,85]. In the recent work [86] such attempt is undertaken concerning the clusters, which is commented below.
The intriguing possibility which would be worth considering in the framework of the cluster model [87] is the effect of a spin alignment of the large groups of quasars which was recently observed for two quasar groups (at z ≈ 1.3) [88].

PBH clusters and observational limits on PBH dark matter
If PBHs are formed in an appropriate amount, they can play the role of the dark matter (DM) in the Universe [89]. How- Fig. 8 Existing constraints on relative contribution of the uniform PBHs distribution into dark matter density with monochromatic mass distribution. The most robust constraints in their regions are basically depicted following the notation of [91]: from extragalactic gamma-ray background (EG) observation, femtolensing (F), neutron star destruction (NS), searches for gravitational microlensing events by Subaru Hyper Suprime-Cam (HSC), Kepler satellite (K), MACHO (ML), quasars millilensing (MLQ), wide binary (WB) and star cluster in Eridanus dwarf galaxy (E) destruction, CMB distortion due to accretion effects (WMAP), effects of dynamical friction in our Galaxy (DF) and large-scale structure (LSS). CMB constraint was taken for (more minimal) case of spherically symmetric accretion [17] in comparison with the less minimal ("CMB disk") one [91] (shown by dashed line). Double-headed arrow shows region where the constraints are to be reconsidered for PBH clusters ever, the «fine tuning» is required to avoid the PBH 1 or PBH 1 situations. As a dominant form of DM, PBHs have been proposed a long time ago; special activity started from the first results on the MACHO search [90]. By now there have appeared many constraints on such candidate in the very different mass ranges (see Fig. 8).
A comprehensive review on that can be found in [91]; there is one more appeared recently [92]. Since the first review (mid-2016), several new constraints had time to appear. The search for microlensing events of stars in M31 with HSC telescope [93] constrains PBHs with mass 10 20 /10 27 g (interval should be cut from below with 10 23 g because of violation of geometric optics approximation [94]). Also, additional analysis of ionization and thermal history of the Universe through the data on cosmic microwave background (CMB) anisotropy allows updating constraint on PBH abundance in different mass ranges. In the light mass range, the analysis provides constraints in a narrow range around M ∼ 10 15 g a little bit stronger than existed [95,96], which originates from CMB distortion by Hawking evaporation. Re-consideration of accretion effects taking into account that they may proceed in a disk regime (not spherically symmetric one) may lead to ruling out PBHs as DM for M M (plus-minus an order of magnitude depending on the velocity parameters) [97] from CMB data. Also higher mass interval of PBH density is claimed to be constrained by the data on Super Novae (SN) gravitational lensing (for M 10 −2 M ) [98], and by possible heat deposition in baryons limited with 21-cm absorption line observation by EDGES experiment at z ≈ 17 (for M ∼ 10 M ) [99]. However, SN constraint is criticized [100], and EDGES data are to be confirmed. There is also some dispute about microlensing constraints for M ∼ 10M [40,101]. However, they may be not so crucial in this mass range. Though some new constraint at M ∼ 10 −6 −10 −3 M appears from OGLE data on microlensing events [102], it allows interpretation as the indication of the Earth mass PBH with percent fraction in DM. There appear new constraints from BH merger events registered by LIGO. Work [103] puts some constraint on PBH amount at M ∼ 10M considering the contribution of gravitational waves from BH mergers to relativistic matter component in the Universe; we comment it together with the constraint from [12] for the cluster below.
Imposition of all the restrictions does not seem to leave a chance for PBH to be a dominant DM candidate. However, the current situation is more complicated. Many works were arguing that constraints coming from CMB, star cluster disruption in Eridanus dwarf galaxy and microlensing effects in the mass range ∼ 1−100M for uniformly distributed PBHs can be weaker. Therefore, both dark matter and gravitational waves from BH merger could be explained by PBHs (see, e.g. [10,52,104]). Also, PBHs of the mass within two gaps around 10 17 g and 10 20−24 g are so far being considered as DM candidates. Also, there are suggested possible constraints to be obtained by LISA in future coming from observation of gravitational waves signal from PBHs sinking to SMBH in our Galaxy [85,105] (in the ranges 10 20 − 10 24 g and 10 −2 − 10 −1 M ), and pulsar timing effects (delay of signal from pulsar) induced by PBH moving along the line of sight [106] (in the PBH mass interval 1 − 10 3 M ). EDGE experiment is suggested to put a constraint on PBH of masses 10 −5 − 10M [107] by investigating density inhomogeneity  (25)) through observation of 21 cm hydrogen line as it was done with Ly-α clouds [44] . These new possible constraints may overlap some still not yet (fully) forbidden mass ranges.
However, one should note that all the constraints mentioned above have been obtained basically for monochromatic PBH mass distribution. As it was shown above it is not the case of PBH cluster as well as of many other PBH models [91]. In the papers [92,[108][109][110] the constraints are tried to be generalized for extended mass spectrum (everywhere a lognormal distribution is considered and also power-law and exponential ones in the last reference). No weakening restriction was virtually obtained.
But not every type of constraint can be generalized «universally»: i.e., as if PBHs with different mass values contribute to the same observably constrained value (say, number of lensing events, or upper limit on the number of destroyed white dwarfs or neutron stars, etc.). For instance, constraint coming from gamma-ray background relevant for PBH with the mass M ∼ 10 15 − 10 17 g cannot be reduced to the «universal» generalization, since PBHs with different mass contribute to different parts of energy γ -spectrum. As was shown in [111], with power-law mass distribution one can gather all DM in the form of PBHs. The upper limit of M in (25) is defined by a mass value starting from which constraint from a femtolensing comes into force. We quote figure here 9, which shows a quite big parameter range (blue region in the coloured figure) where PBHs can substantially contribute to DM and simultaneously to reionization of the Universe (see Fig. 3 left of [111] and discussion below in Sect. 5.5). The PBH as DM with the mass around the same value is also considered in other works, e.g., [94,112]. Contribution to DM can also be maximized with several delta-functions (peaks) at different masses of PBH mass spectrum [111,113].
In the cluster model, a mass distribution can be very different depending on the model parameters (mostly of a scalar field potential). The case considered above gives nearly decreasing power-law from PBH mass (α ≈ −(2−3)). However, one can get growing power law as well as distribution with several peaks at different mass values. However, it is not the only opportunity to change or even avoid PBH abundance limitations in case of a cluster.
Almost all the existing constraints on PBH abundance with monochromatic mass distribution have to be reconsidered in the case of space and mass distributed PBH clusters. For instance, all constraints from femto-/micro-/milli-lensing can be irrelevant since it is a cluster as a whole that should give a lensing effect, while its size may exceed the Einstein radius. Besides, the probability of lensing for smaller PBHs can be suppressed because of their inhomogeneous (clustered) distribution in space with effectively diminished number density to the cluster's number density. However one should take into account that there can be an unclustered component. The latter could consist of those PBH that escaped by clusters, and the clusters themselves could be disrupted in Galaxy as well.
The investigation of the lensing effect for the cluster which size exceeds the Einstein radius is underway now by our group. The opposite case is considered in [40,41]. They take the cluster model [44,45] with the lognormal mass distribution of PBHs inside the cluster. Some increase of total PBH abundance, in this case, can be reached due to an increase of PBH amount inside the cluster which shifts and broadens the mass distribution of the clusters what may make limitations a little more flexible.
As to dynamical constraints, they should also be reconsidered. These are disruption of white dwarfs, neutron stars, wide binary systems, star clusters, etc. Constraint coming from star cluster disruption in dwarf galaxy looks like one of the strongest and reliable for the PBH clusters. The constraint could relate to the PBHs of masses around ∼ 10M inside the cluster.
Also in similar mass range, existing data on gravitational waves can constrain the cluster abundance, as claimed in the recent work [86]. However, the process of PBHs pairs formation considered there seems to be more relevant at early stages (before virialization and decoupling from Hubble flow) of cluster evolution. Also, one should note that these constraints, if they hold, can be ineffective for the PBH cluster model considered here. PBHs of the mass ∼ 10M are not expected to give a dominant contribution to DM density provided by falling power-law mass spectrum (31), (25). If we take α = −2.4, as obtained above, and the mass range the relative contribution of PBHs with M > M to their total density is less than 10 −7 . Note that the range (26) is different from the one above. Nonetheless, it could be got so by customizing initial parameters. However, it would strongly complicate the simulation of the cluster.
As to the constraint from CMB, based on consideration of accretion by PBH with mass M (the most robust constraint for the most massive PBHs), it can be weakened for clusters due to a motion of PBH inside them, changing the regime of accretion [41]. At the same time, if one takes the very massive PBHs, their amount is not supposed to be large to explain seeding SMBH in galaxy centers.
Finally, the majority of constraints considered and obtained in [91] for the monochromatic PBH mass distribution is worth being re-considered for PBH clusters with extended mass distribution both for PBHs inside the cluster and for clusters themselves. In Fig. 8 we point out the PBH mass interval of constraints obtained for monochromatic mass distribution which is to be re-considered for clusters. We suppose that this interval can be significant.
There is another advantage of considering the (clustered) PBHs as the dark matter. As was noticed in [10], their gravitational interaction allows circumventing the problems of standard non-interacting CDM scenario like cusp crisis, dwarf galaxies overproduction, etc.

Gravitational waves from PBH coalescence
By now ten BH merger events have been detected by LIGO/Virgo [114]; the first five announced [115][116][117][118][119] are given in the Table 1. It seems from viewpoint of standard astrophysical origin to be rather unnatural, that they are pairs of heavy black holes ( 10M ) with low initial (effective inspiral) spins. More arguments contra BH of star origin and pro primordial one for given events are collected in [10].
Although it is possible that binary BHs with masses of ∼ 30M can be formed after deaths of Population III stars with low metallicity. But there still exist large uncertainties in theoretical predictions of the merger rate by a factor of −0.06 +0.14 −0.14 O(100) due to unknown astrophysical parameters [127,128]. Such possible origin could be confirmed by observation of an electromagnetic signal, but no electromagnetic counterparts were observed for BH merger events [129]. It is worth mentioning once again counter-arguments for star origin given in [10]. The rate of PBHs collisions in clusters was calculated in [130] for the planned LISA (Laser Interferometer Space Antenna) instrument sensitivity curve in the frequencies range 10 −4 −1 Hz. Although the LIGO/Virgo interferometers fulfilled the first detection of gravitational waves, the primary approach of [130] is applicable to LIGO/Virgo too, but the preferred masses for detection move to lower values due to the smaller working frequencies of the ground-based detectors. We plan to perform the similar calculations for LIGO/Virgo in future works, but here we outline the general ideas and present the numerical results, obtained in [130] for LISA. These predictions can be used if the LISA project will finally be realized.
The initial mass spectrum of PBHs in the cluster is supposed to be the same as in the paper [131]. According to this work, the distribution just after the cluster virialization can be fitted as The total mass of the PBHs prevails that of the dark matter at the distances r < 1.6 pc from the center, and we restrict our analysis to this region. It should be taken into account that two-body relaxation in the cluster leads to the contraction of the inner region and the merging of the inner shells with a central black hole. For each moment, we find the radius of the full relaxation and exclude its interior from the calculations. Then we calculate the rate and the distribution of the PBHs collisions. The cross-section for two PBHs capture and subsequent merging is [132] σ mer = 2π 85π where M 0 , M are the masses of the PBHs, and v rel is their relative velocity which is assumed to be of the order of virial velocities in the cluster center. With this cross-section, the rate of the PBHs collisions in the cluster can be calculated as where dn is the PBHs number density with masses in the interval M/M + d M. It was found in [130] that the main signal (for LISA) comes from the collisions of PBHs in the cluster with the central most massive black hole. Note that in the case of LIGO/Virgo the mutual collisions between smaller black holes in the cluster would be more important.
The signal-to-noise ratio can be calculated as [133] ρ where S h ( f ) is the noise spectrum of the detector. The functionh( f ) represents the spectrum of the signal. It was calculated in [133] for the local population of sources (z 1) and was obviously generalized by [130] for arbitrary red-shifts z. For each PBH's mass M there is a distance threshold -maximum redshift z from which it can be detected. We divide the full mass interval into several smaller mass ranges and calculate the gravitational bursts rate for each of the intervals. We also suppose that the merging of the intermediate-mass black holes in the cluster finally produces the supermassive black holes in the typical galaxies. This assumption allows one to fix the number density of the clusters in the universe. Abolishing this assumption one obtain different normalization of the events rate.
The results of the calculations are presented in Fig. 10. One can see that the detection of the gravitational bursts by LISA would be possible during the reasonable time of observations. The similar study for LIGO/Virgo is desirable and is planned soon. The dependence of the rate on the red-shift z in our model differs sufficiently from the one obtained in the model [134], where astrophysical black holes, formed from the collapses of the gas clouds, were considered. This specific redshift dependence of the merger rate is determined by the early beginning of the two-body relaxation in the cluster centers. It led to the central density increase and the merger rate amplification. The difference in the redshift dependence could help to discriminate between models of black holes formation.
Notice that there is a contribution to the gravitational waves from the very light PBHs which can be experimentally accessible in future [12,135].

PBH clusters as point-like gamma-ray sources
Due to weak accretion of surrounding matter, observation of BHs with masses less than 10 4 M can be quite difficult. An alternative method for the BHs detection through the Hawking's radiation [136] is useful only for BHs with low masses that are situated near the Earth. From [137,138] one can conclude that situation can be different for PBH clusters with a large number of small BHs (with masses 10 15 g). The integral Hawking's radiation of such a cluster could then be sufficient for detection by modern gamma-ray telescopes on the Earth (in the near-earth space). Now there is a big database on unidentified point-like cosmic gamma-ray sources (PGRS) obtained by Fermi-LAT [139] which has more than 500 such sources.
The most of γ -ray sources (more than 60% of their total amount) are associated with active galactic nuclei, γ -ray pulsars and some other Galactic sources. However, there is an essential part of PGRS without counterparts found in other wavelengths and thus non-identified. They could have a conventional origin and can not be identified because of insufficient angular resolution but also could be of unknown physical origin, which attracts so far special interest (e.g. [140][141][142]).
We made estimations for specific case of cluster model [137,138] following [131]. Fractal-like PBH cluster structure is supposed to be formed around the most massive PBH being the progenitor of SMBH. It, in turn, is assumed to give eventually rise to a galaxy formation including ours.
In given case we had about 1500 clusters per our Galaxy, each with total mass about 10M and contained numerous small PBHs. Initial PBH mass distribution was obtained there in the form of a power-law just like we get here A total amount of photons from the cluster with E > 100 MeV is obtained to bė Then conditioṅ defines the maximum distance at which the cluster can be seen, where F min is the respective Fermi lower threshold. It allows us estimating roughly possible number of PBHs which could be observed as gamma-ray sources, where n PBH is the PBH cluster number density in Galaxy which was supposed to be homogeneous. So, PBH cluster at r max has a minimal flux I accessible to observation, which defines the lower boundary of the gray band in the Fig. 11. The closer the PBH clusters, the brighter they are, but the smaller their number is. At some distance only one cluster can be expected to be observed, the intensity, in this case, corresponds to the upper boundary of the band in Fig. 11.
In the future, one has to take into account the evaporation of light PBHs from a cluster which gives the main effect in gamma-radiation. Nonetheless, at the chosen parameters, it was found out that Fermi/LAT could observe the clusters as PGRS with power spectra as Fig. 11 shows. So (at least) part of the unidentified PGRS observed by Fermi/LAT could be explained by the PBH clusters.
In Fig. 11, sensitivities of X-ray telescopes (INTEGRAL [143] and OSSE [144]) are also shown, which turned out to be well above the expected flux which Fermi/LAT can see.
As one can see in Fig. 11, at the given parameters, the gamma-ray spectrum behaves from energy as E −3 at high energetic part (at low one it goes as ∝ E because of change of PBH mass spectrum (31) due to Hawking evaporation). Unfortunately, not for all PGRS observed a spectrum was Fig. 11 The expected gamma-ray spectra from PBH clusters (the shaded region). Below the shaded region the flux from PBH is lower than LAT integral sensitivity. Above it, too few clusters are expected to be observed as PGRS, namely less than one (see the text for more details). The differential sensitivity of detectors INTEGRAL and OSSE is also shown Fig. 12 Unidentified PGRS seen by Fermi LAT (blue crosses) on celestial sphere (in galactic coordinates b and l) are shown. The red filled markers are the sources with spectrum indices 3 within 1σ error, unfilled markers are the sources with spectrum indices 3 within 3σ error measured, but one can select those of them where it was done, and select from them those with a similar spectrum, Fig. 12. One can note that such sources are distributed over the sky more or less homogeneously as it would be expected for clusters.

Reionization
Reionization of the Universe happened at z ∼ 7 − 8 is still possibly one of the problems of cosmology and astrophysics. This fact was established by different observations, mainly through Gunn-Peterson effect and CMB data (see, e.g. [145,146]). There are attempts to explain it with dwarf galaxies, quasars, and stars of Population III [147][148][149]. abs and the total rate, for M = 5 × 10 16 g. Light gray lines relate to the respective losses obtained approximately (see [152,153]).˙ ev is shown to illustrate total evaporation rate into all species and e ± + γ only PBHs could be at least one of the contributors into reionization of the Universe, and both light and massive PBHs could do it. The massive ones -due to radiation induced by accretion, [17,97]. However, space correlations of the X-ray and Infrared backgrounds could instead favour this opportunity [10,150,151].
Light PBHs could also contribute, they can be divided into two components. As was shown above, the PBH cluster loses predominantly light PBHs during its evolution, which then may be supposed to spread more or less homogeneously in space. However, at the same time the total cluster mass changes insignificantly, so the lightest PBHs escaped the cluster make up small DM fraction and may avoid constraints coming from gamma-radiation as well as those inside it.
All the light PBHs, both inside and outside the cluster, can contribute to reionization of the Universe due to Hawking's radiation. Reionization effect due to homogeneously distributed PBHs have been considered by our group [111,[153][154][155]. PBH of the mass M ∼ 10 15 − 10 17 g would produce ionizing Hawking's radiation in the form of (sub)relativistic electrons/positrons and photons (along with non-ionizing neutrinos and gravitons). They lose their energy due to different processes: e ± -due to inverse Compton scattering off CMB (Cosmic Microwave Background) photons, ionization losses, redshift; γ -due to Compton scattering off electrons of the medium and redshift, and photons from e +annihilation have the same losses. All the energy loss rates are shown in Fig. 13 depending on the redshift for monochromatic PBH mass distribution with M = 5 × 10 16 g and their density as big as possible from data on the gamma-ray background. In Fig. 13  abs is the rate of energy deposition in baryonic matter per unit volume for i-th process, ρ crit is the critical density of the Universe.
As one can see from the figure, these are e ± -ionization losses which give the main contribution to energy transfer from PBH radiation to the matter. Re-ionization effect is reached in given approximation for PBH mass around (3..7) × 10 16 g.
The effect can be enhanced if one uses less trivial mass distributions, e.g. two-peak or power-law ones [111]. Reionization (enhancement of contribution in it) could be reached at a quite noticeable region of parameters for a power-law mass distribution (see Fig. 3 left of the reference above).
However, within our approximation, it was supposed that all energy deposits in matter via a single channel -heat, while it should be shared between three -heat, ionization (what goes to overcoming of ionization potential), atom excitation [95,96,[156][157][158]. The effect could be weaker if we took into account all the channels.
Also one should take into account the data on baryon temperature at z ∼ 17, deduced from 21 cm absorption line of hydrogen observation [159], which may give extra restriction for PBH in a given mass range [160].

Conclusion
The necessary components of the modern theory of the structure and evolution of the Universe are the inflation, the baryosynthesis, and the dark matter/energy. The physical basis of this currently Standard cosmological paradigm implies physics beyond the Standard Model of fundamental natural forces. The extension of the Standard Model of elementary particles is necessary to explain the mechanisms of inflation and baryosynthesis, as well as physics of dark matter and dark energy, inevitably contains additional cosmologically viable consequences. Primordial black holes are among the most important probes for such features of new physics in the very early Universe.
The progress in observations of distant object unambiguously points to the early formation of the massive and supermassive black holes. The origin of black holes of 10-30 Solar masses, gravitational wave signal of whose coalescence was detected by advanced LIGO and Virgo, or the origin of massive black hole seeds for AGNs may cause problems for their explanation based on stellar evolution in galaxies and can find their nature in primordial black holes.
A substantial amount of models are elaborated to explain the phenomena of primordial black hole formation. Some of them contain an opportunity of PBH cluster formation. The set of such models is discussed in this review. We focus mainly on one of them related to the phase transition with the formation of the closed walls. There are three steps of the PBH formation. Firstly, the quantum fluctuations during the inflation produce the space areas with scalar field values at both sides of a potential extremum. Secondly, after the end of inflation, the field settled in different vacua forming the closed walls. At the third stage, these walls move through the surrounding media, heating it and finally forming the black holes.
Field fluctuations at successive steps of inflation lead to the change of vacuum at smaller scales so that after the second phase transition the largest wall becomes inevitably accompanied by a cloud of smaller walls. Provided that the width of the wall does not exceed its gravitational radius, the walls collapse into smaller black holes, so that PBH cluster is formed.
Potential properties and initial conditions determine the range of PBH masses in clusters. The minimal mass is determined by the condition that the gravitational radius of the wall exceeds its width, while the condition that the wall does not start to dominate before it enters the horizon, puts an upper limit on the PBH mass.
Inevitable ingredient of the model is the internal evolution of the PBH cluster. After a cluster separates from the expansion, the formation of BH binaries inside it is strongly facilitated as compared with the case of stochastic PBH distribution. The PBHs merge and scatter each other, some of them evaporate from the cluster. The first numerical simulations presented in this review indicate that space and the mass distributions are not altered significantly. We show in the Sect. 4.2 that only ∼ 15% of black holes leave the cluster. Remaining BHs form the compact cluster with sizes ∼ 1 pc and the total mass 10 3 M at the moment z ∼ 20−40. Such clusters could serve as the «seeds» of the early quasars observed at z > 6.
The features of PBH clustering can shed new light on the observational data and provide us with new degrees of freedom in their analysis. PBH cluster model inherits all the advantages of those models predicting uniform PBH distribution like a possible explanation of the existence of supermassive black holes (origin of the early quasars), binary BH merges registered by LIGO/Virgo through gravitational waves (GW), contribution to reionization of the Universe, but also has additional benefits. The cluster could alleviate or avoid at all the existing constraints on the uniformly distributed monochromatic PBH abundance so possibly becoming a real DM candidate. Most of the existing constraints on PBH density should be re-considered if to apply to the clusters. The first attempts, mentioned in this review, have been already done to limit possible abundance due to observed BH merger rate, but PBHs of respective mass range could be much below that limit. Also unidentified cosmic gammaray point-like sources could be (partially) accounted for by them. One can conclude that it seems to be a much more viable model as compared with the models predicting uniform PBH distribution.
It is worth noting that the gravitational waves can become the specific prospective tool for the cluster model validation. The PBHs merges into pairs into the clusters could provide the sources of GW bursts in addition to the binary black holes of stellar origin. One of the most interesting effects is the possibility of recurrent (multiple) events of GW bursts from the same sufficiently dense cluster. On the contrary, the binary systems of stellar origin produce only single events from different galaxies. Therefore, future observations of the anticipated multiple events can provide the decisive test of the PBHs clustering discussed in this review.
The discussion made above is applicable to any compact cluster model. To date, several such models, different from the one we basically focused on here, have been proposed. Their formation mechanisms are based on strong primordial density fluctuations. According to these models, PBHs are clustered in the course of evolution due to their Poisson distribution or due to the seed effect. Such mechanisms are also applied to big clusters on the scale of large scale structure to explain its formation.
The mechanism of the PBH cluster formation discussed in this review can be tested by future observations. Indeed, the closed walls are its necessary attribute. They are nucleated having the shape far from the spherical one [35,36,65,66]. The wall inhomogeneities are damped transferring the energy to the surrounding media and heating it. According to our estimation, the heating could be extremely effective so that the temperature of matter around such PBH cluster could be high. This allows one to separate our mechanism from the others. We develop the theme, understanding its importance.
In the context of cosmoparticle physics, studying the fundamental relationship of cosmology and particle physics, search for PBH clusters and their study acquire the meaning of nontrivial probe of new physics of the very early Universe.