Formation and clustering of primordial black holes in Brans-Dicke theory

The formation of primordial black holes in the early universe in Brans-Dicke scalar-tensor theory of gravity is investigated. Corrections to the threshold value of density perturbations are found. Above the threshold the gravitational collapse occurs after the cosmological horizon crossing. The corrections depend in a certain way on the evolving scalar field. They affect the probability of primordial black holes formation, and can lead to their clustering at large scales, if the scalar field is inhomogeneous. The formation of the clusters, in turn, increases the probability of black holes merge and the corresponding rate of gravitational wave bursts. The clusters can provide a significant contribution to the LIGO/Virgo gravitational wave events, if part of the observed events are associated with primordial black holes.


I. INTRODUCTION
The principal possibility of the primordial black holes (PBHs) formation in the early universe was stated in the works [1,2]. Several mechanisms were proposed, among them: collapses of adiabatic density perturbations [3,4], collapses of perturbations at early dust-like stages [5,6], collapses of domain walls [7][8][9], collapses of baryon charge fluctuations [10][11][12]. It is possible that some of the LIGO/Virgo events are explained by the merge of PBHs [13][14][15][16][17]. The value of density perturbations necessary for the formation of a PBH, can be achieved due to the presence of features in the inflationary potential [18,19], as well as in inflationary models with several scalar fields [20].
The process of PBHs formation depends on the underlying theory of gravity. In addition to numerous analytical and numerical calculations performed in the framework of General relativity, there are studies of the PBHs formation in modified gravity theories. The paper [21] considered the effect of the gravitational constant changes in the Brans-Dicke theory on the number of PBHs formed. The evolution of PBHs population in the Brans-Dicke cosmology was also studied in [22,23]. A PBH is formed only from density perturbation that is greater than certain threshold value. The number of PBHs is very sensitive to the threshold of their formation, and in General relativity, the threshold is ∼ 1/3 − 0.5. For some principal issues related to the threshold, see the paper [24]. The modification of the threshold in the Eddington-inspired-Born-Infeld gravity was investigated in [25], and it was shown that the threshold in this theory of gravity is not constant, but depends on the epoch of PBH formation.
In this paper, we also calculate corrections to the threshold of PBHs formation, but we do it within the framework of the Brans-Dicke theory. In contrast to [21][22][23], we write down equations for the evolution of spherical-symmetric perturbations in the Brans-Dicke theory and consider the gravitational collapse. The corrections depend on the scalar field of the Brans-Dicke theory. This dependence leads to the fact that PBHs in the regions with different scalar field are formed with different probability. The scalar field modulates the large-scale distribution of PBHs. The modulation effect is similar to the well-known biasing effect for galaxies. The galaxies are more likely to form in regions with higher density than in regions with lower density [26].
In [27] and [28], the attempts were made to explain the long-wave spectrum of perturbations (on the scales of galaxies and clusters of galaxies) by statistical fluctuations in the PBHs number density. As was shown in [29], this mechanism is ineffective. In [27], the initial (deterministic) long-wave fluctuations in a mixture of PBHs and radiation were considered. It was assumed that the perturbations in the PBHs have the same amplitude as the initial perturbations in the radiation. A more precise approach requires accounting for the biasing effect that accompanies the formation of PBHs at the background of long-wave perturbations. The possibility of the biasing for PBHs within the General relativity was proposed in [30]. However, later, by more precise calculations, it was shown that this effect is very small [31,32]. The effect, in a sense similar to the biasing effect, can also appear due to the non-Gaussian nature of initial perturbations [33][34][35].
A working mechanism for PBHs clustering was developed in [9], where it was shown that for a certain inflationary potential, spherical domain walls are formed, which break up into PBHs. The resulting PBH clusters may have a number of observational consequences for gravitational-wave astronomy and for the formation of early galaxies and quasars [36]. The influence of PBH clustering on the rate of their merge within the General relativity was discussed in [37].
In this paper, we investigate the new biasing effect for PBH in the Brans-Dicke theory. We show that the biasing (clustering) of the PBHs actually occurs. It can affect the rate of PBH merge, which may be one of the sources of gravitational waves detected by LIGO/Virgo. The Brans-Dicke theory [38,39] is one of the most well-known extensions of General relativity, reducing to the General relativity in the limiting case when one of the parameters of the theory ω → ∞. Therefore, as long as the experimental data do not exclude General relativity, the Brans-Dicke theory will also remain viable with the appropriate choice of parameters. Currently, the observational data give the following restrictions on the Brans-Dicke parameter ω > 40000 from the Cassini-Huygens measurements.
Perturbations of curvature and perturbations of scalar field that lead to the formation and clustering of PBHs are generated at the inflation stage. In this paper, we do not fix any particular theory of the inflation in scalar-tensor theory. We note only that, according to the theoretical calculations of [40], [41], in Brans-Dicke inflation theory, curvature perturbations and perturbations in a scalar field other than inflation can be independent. This is enough to demonstrate the new mechanism for PBHs clustering. Although the specific type of inflationary perturbations and quantitative characteristics of the clustering depend, of course, on the specific inflation theory.

A. Basic equations
We consider the Brans-Dicke theory with the Lagrangian where ω = const, and L m is the Lagrangian of matter. In calculations, we assume that the speed of light c = 1. Following approach of [4], we model a perturbed region of space collapsing to the PBH as a part of a closed universe with a metric where S = S(τ ) is the scale factor. And the metric outside the perturbed region is the metric of the flat universe R t S Τ Figure 1. The collapsing region is modelled by a part of a closed universe with the scale-factor S(τ ). Outside the perturbed region the flat universe is assumed with the scale-factor R(t). And some transition layers are possible.
where R = R(t), see Fig. 1. Between the regions some transition intermediate layers are possible, see discussion in [25]. In contrast to [4], we assume that the evolution is governed not by the General relativity, but the Brans-Dicke theory. In the Brans-Dicke theory, the cosmological equations for a flat and closed cosmological model have the form [42], where in the equation (4) a point means a derivative over t, and in (5) a point means a derivative over τ (and further similarly). We set the initial data at some early time t i (τ i ) after the inflation stage, but when all density perturbations are still small and described by linear theory. Values at the initial moment are marked with the index "i". The cosmological epoch of radiation dominance is considered (equation of state p = ε/3), therefore in the equations (4) and (5) one has The scalar field evolves asφR 3 = const,φS 3 = const [42], sȱ Note that for ωφ i t i /φ i ≪ 1, the last term on the right hand side (5) can be ignored throughout cosmological evolution, but in general the last term may be important. And the last closing equation gives a relation between τ and t [43]. This equation follows from the general expressions for the comoving coordinate system [44] and is valid both in General relativity and in Brans-Dicke theory. The value of the density perturbation is defined as In (8), the perturbation is taken at the initial moment t i (τ i ), and the initial constant-time hypersurface is chosen so that when t = t i (τ = τ i ), the following conditions are met [43] Note that in other works [45], [46], where the formation of PBHs was considered in the framework of General relativity, the constant-time hypersurface is chosen differently (we use the choice that coincides with the choice in [4] and [43]), and it must be taken into account when comparing the results of the zero approximation of General relativity. In the following sections, we first find the solution of (4) and (5) in the zero approximation atφ = 0 (coinciding with the General relativity), and then we solve the problem in the first linear approximation in Brans-Dicke theory.

B. Zero approximation
In the zero approximation,φ =φ = 0 and φ =φ = 1/G, where G is the Newtonian gravitational constant. This approximation matches the results obtained in General relativity in [4]. We will present here the known results from [4] and write down some new relations for the evolution of δ and for the threshold of PBH formation.
In the zero approximation, we denote A ≡ S (0) . The condition dA/dτ = 0 in (5) gives the scale factor at the moment τ max of the maximum expansion of the disturbed region because δ i ≪ 1, and from (10) we have The evolution of the scale factor in the zero approximation is described by the equation Introducing the parameter ψ, one get the solution in the parametric form [43] where the relation of the parameter ψ to τ has the form and the relationship between ψ and t is as follows In the last expression, in contrast to [43], we took into account the next order corrections for δ i , which will be required later to calculate δ(t).
The maximum expansion of the disturbed region corresponds to the moment of time and at the same time When selecting a constant-time hypersurface, the ratio of densities at the time of maximum expansion is ε/ε = (π/2) 4 . As already mentioned above, in other works (for example, [45,46]) a different choice of the hypersurface was used, in which the time in an undisturbed flat region of the universe is t = τ , and this choice results in ε/ε = 4.
The value of the perturbation is expressed as follows where ψ is related to t by (16). For t ≪ t max one get It means that in this coordinate system, the falling perturbation mode evolves as ∝ t −1/2 , in contrast to [46], where the law ∝ t −1 was obtained. And the evolution of a growing mode on scales larger than the horizon occurs according to the law ∝ t known for a synchronous reference frame. The PBH formation criterion was obtained analytically in [4]. The region of space with δ > 0 was modelled as part of a closed universe (Friedman model). The gravitational collapse of such a region with the formation of a PBH occurs if the relative value of the perturbation at the moment of horizon crossing δ H satisfies the following conditions: In early studies criterion for the threshold δ c was used, according to which PBH is formed, if the radius of the disturbed region at the time of of its expansion stop exceeds Jeans radius. It corresponds to the left inequality in (21). And the right inequality corresponds to the formation of a black hole, not a separate universe. As noted in [47], this criterion contains ambiguity in the expression for the Jeans radius and can therefore be considered as only an estimate. Further, when considering the scalar tensor theory, we will look for small corrections to the threshold of PBH formation, so in this paper we use the refined PBH formation criterion proposed in [47]. According to the improved criterion, a PBH is formed if the sound wave does not have time to reach the center from the periphery of the disturbed region by the time the expansion stops. In this case, the Jeans radius is assumed to be equal to the sound horizon.
We write the metric of the perturbed region in the form then the equation of the sound wave has the form [47] Substituting the above expressions, we get Integrating up to the moment of maximum expansion ψ = π/2, we have which in the case of the equation of state p = ε/3 matches the result obtained in [47]. However, we use a different choice of hypersurface, so the further calculation will be different. The conditions for the Jeans radius (sound horizon) and the condition for cosmological horizon crossing make up the system of equations where χ a is the boundary of the disturbed region, and ψ H is the value of the parameter ψ at the moment of the horizon crossing. From these two equations one get In the limit δ i ≪ 1 the (28) is transformed into the following nonlinear equation where Numerical solution gives the root of the above equation Substituting t H = ξ 0 t i /δ i in (19), we get the threshold for the PBH formation in the zero approximation Note that in other coordinate systems, the threshold value is different.

C. Corrections to the black hole's formation threshold
In the following approximation we take into account the evolution of the gravitational constant (change of the field φ) in the Brans-Dicke theory. We assume that all the corrections to General relativity are small.
After the variable replacing, we rewrite the equation (7) as follows Its solution is Substituting expressions for the zero approximation on the right, we get at ψ = π/2 Here, as before, the index "max" denotes the values at the moment when the expansion of the disturbed region stops, i.e. dS/dτ = 0. We introduce small corrections to the zero approximation as follows where s i is the relative deviation of φ from G −1 at the initial moment, p and α are new small functions, of which the first was already calculated above in (34). Let us write the equation (5) linearised over α: The last term contains the small parameterφ i t i G squared, but because of the possibility of large values of ω, we keep it. SubstitutingȦ from the zero approximation and introducing the notation x = A 2 /A 2 max , we rewrite (37) as where C 1,2,3 are some combinations of constants. General solution of the above equation is Using the x variable, the solution (34) is written as follows It shows the structure of the function p.
From the zero approximation, we have C 1 ≈ 1, Using these values, after integration we get At the moment of maximum expansion x = 1 one get α ≃ −µ − s i /2 + ωµ 2 δ i /3 from (41). For verification, let us calculate the same value in a different way. From the condition dS/dτ = 0 in (5), we obtain an algebraic equation of the 4th order for the value S max /S i , and its solution is Substituting expressions for the zero approximation, we obtain the same value α ≃ −µ − s i /2 + ωµ 2 δ i /3. Instead of (23), for the sound wave equation, we have the equation Through the variable x it looks as follows Note that when integrating this equation, the size of the sound horizon for the first and third terms in (41) is accumulated at early times τ (small x), and for the second term, the size is accumulated at later times just before the PBH formation. By performing the integration, we get the solution with the correction where was denoted Instead of (29), we get the equation The corrected solution can be found as follows where the derivative dξ/dγ| γ=0 ≃ 0.6255 is calculated by differentiating the equation (48). By a similar decomposition, we get the final correction for the threshold of PBH formation: where λ ≃ 0.432, and γ is given above in (47). The expression (50) is the main result of this work. In the following chapters, we will consider its application for the effect of PBH clustering and, accordingly, for the rate of gravitational bursts producing by the PBH merge. It is nontrivial to choose the initial moment of time t i , since in (47) the term in parentheses at the linear stage is invariant with respect to the choice of the initial moment and is equal to (Gφ max − 1), which follows from (34) and (35), but the other terms depend on t i . To fix the initial conditions, we assume that the moment t i corresponds to the end of inflation. This model is sufficient to demonstrate the effect of PBH clustering.

A. Perturbations in the PBH number density
To illustrate the effect, let's consider a simple model in which PBHs have a monochromatic mass spectrum, i.e. all PBHs are formed with the same mass M PBH . We assume also that the perturbations are Gaussian, and the mean squared value of perturbations on the scale of this mass at the moment of horizon crossing is denoted by σ H . For simplicity, let the scalar field have perturbations on some larger scale, which currently contains a mass of cold dark matter M ≫ M PBH . At the stage of inflation, perturbations with distinguished scales are generated if the inflationary potential has features, such as local flat segments [18,19].
The PBH is formed under the condition (21). The probability of PBH formation, i.e. the fraction of radiation transformed into the PBHs at the time of their formation, is written as where the last expression is obtained by considering the tail of the Gaussian distribution δ ≥ δ c ≫ σ H , and the upper limit of integration is not important. Now-days cosmological parameter of PBHs At the radiation-dominated stage, the scale factor of the universe is a(t) ∝ t 1/2 . If one assumes that the fraction ∼ 1 of all LIGO/Virgo events is due to PBHs with masses M PBH ∼ 30M ⊙ , then the PBHs constitute the fraction Ω PBH /Ω m ∼ 10 −3 of all dark matter [15], where Ω m is the cosmological density parameter of dark matter. From (51) and (52) one obtains σ H ≃ 0.013 in this case. Now we will show that perturbations in the scalar field that give a correction in (50) are translated into large-scale perturbations in the PBH number density, i.e. biasing effect is in place. The PBH distribution is modulated by the inhomogeneities of the scalar field and its time derivative in combination (47). Moreover, due to the threshold of PBH formation at the tail of the Gaussian distribution, it turns out that perturbations in the PBHs can be significantly amplified in comparison with perturbations in the scalar field.
Substitute the value δ c from (50) to (51). Then the perturbation of the PBH number density is written as The minus sign is explained by the fact that a positive value of γ increases the threshold, and PBHs form in a smaller number, and in the case of a negative γ, the opposite situation occurs. The same is valid for the mean root square values Since (δ (0) c ) 2 /σ 2 H ≫ 1, one has δ PBH ≫ λγ, that is, the perturbations in the PBH number density are much larger than the initial perturbations in the scalar field. For the example above (M PBH ∼ 30M ⊙ , Ω PBH /Ω m ∼ 10 −3 ) we have (δ However, the PBHs constitute only a small part of all dark matter, so the total value of the density perturbation is where Ω PBH is given by (52), and the last estimate corresponds to the example above. Here, fluctuations in the scalar field are translated into fluctuations in the total density of dark matter, although with less efficiency than in the PBH component alone, as it was the case (53). Previously, the bias effect for PBHs was already considered in General relativity [30]. The density perturbation was taken as the sum of the horizon-scale perturbation and long-scale perturbation. Therefore, on the background of a positive long-scale perturbation, the total perturbation can easier overcome the δ c threshold than on average. The usual bias effect applied to galaxies in galaxy clusters also consist in the fact that in region with high density, galaxies are formed more likely than in regions with low density [26,48]. In [30], a conformal Newtonian coordinate system was used, in which density perturbations at the radiation-dominated stage do not evolve and are equal to perturbations at the moment of the horizon crossing. The sum of perturbations in conformal Newtonian coordinate system was taken at the single time moment. It was shown later that the effect was overestimated by several orders of magnitude, see [31]. The easiest way to understand the confusion is to use a synchronous reference frame. Indeed, in the synchronous reference frame δ ∝ t on the super-horizon scales, therefore at the moment of PBH formation the long-wave perturbation has a much smaller values than at the moment of their horizon crossing. Therefore, the contribution of the long-wave perturbation in the total perturbation is much smaller than it was obtained in the naive calculation in the conformal-Newtonian frame, and the biasing effect is therefore very small. A more accurate method of summation of perturbations at different scales is required for the correct calculation in the conformal-Newtonian frame [31].

B. Influence of inhomogeneities on the rate of gravitational bursts
The number density of PBHs is increased in regions where γ has negative sign, see (53). Let's first consider the case when the clustering of the PBHs is so small that the specified regions did not form virialized clusters up to now.
Accidentally close PBHs can create a connected pair at the cosmological stage of radiation dominance, which is shrinking in size due to the emission of gravitational waves and, as a result, the PBH merge [13]. The elongation of their orbit and the time of merge is determined mainly by tidal forces from the third nearest PBH. The effect of inhomogeneities in the distribution of PBHs is not simply increases the probability of pair formation due to the closer location of PBHs in the inhomogeneities. It also plays a role that the formed pair must merge at the present time in order for the signal to be registered by the detector. This makes the dependence of merge rate on the clustering less obvious. The rate of PBH merges at the present time t 0 was found in [13] as where ρ c = 9.3 × 10 −30 g cm −3 is the critical density, Ω m ≈ 0.27, f is the fraction of PBHs in the dark matter, and the probability that the lifetime of a of PBH pair is less than t has the form where the average distance between the PBH isx where ρ eq is the average density of dark matter at the moment of transition to the dust-like stage of the universe evolution.
Since t max ≫ t 0 is valid for LIGO/Virgo signal parameters, the first term in (57) gives the dominant contribution to the merge rate. Then for a fixed f we have R b ∝x −12/37 . Since for a non-uniform distribution of the PBH x ∝ (1 + δ PBH ) −1/3 , where δ PBH is set by (53), the final expression R b ∝ (1 + δ PBH ) −4/37 depends only very weakly on δ PBH . Since it seems that the extremely large values of λγ cannot be reached, we conclude that the effect of PBH clustering on their collisions in bound pairs formed at the radiation-dominated stage is very small. This may not be in the case when PBHs form virialized clusters, which we will discuss in the next section.

C. Clusters of primordial black holes
Perturbations of PBH number density (53), which arose due to inhomogeneities of the scalar field, and the total perturbations (55) belong to the class of entropy perturbations (perturbations with constant curvature). They are independent on the curvature perturbations associated with the perturbations in the relic radiation density. At the stage of radiation dominance, these entropic perturbations almost do not evolve. Their value increases only 2.5 times due to the Meszaros effect. And after the onset of the cosmological stage of matter dominance, these perturbations evolve in the same way as the usual perturbations in dark matter ∝ t 2/3 .
Direct collisions of PBHs that are not part of binary systems can occur if the PBHs are located in virialized clusters. Let's consider the formation of such clusters and determine the rate of PBH collisions in them.
The cross section of the gravitational capture and subsequent collision of two PBHs [49] where r g ≡ 2GM PBH /c 2 , and v is the relative speed of the two PBHs, and the light speed was saved in the formulas of this section. The rate of PBH collisions in a single cluster with a total mass of M (the sum of the PBH masses and the mass of the rest of dark matter) where, n is the number density of PBHs, N = M (Ω PBH /Ω m )/M PBH ≃ (4π/3)R 3 n is the number of PBHs in the cluster, v ∼ (GM/R) 1/2 , R is the radius of the cluster, where, according to the spherical model of perturbation evolution, ρ cl ≃ 18π 2ρ (t f ),ρ(t) is the average cosmological density of dark matter at time t, and t f is the moment when the perturbation evolving according to the law ∝ t 2/3 reaches the value ≃ 2.81. The value Ω PBH /Ω m ≤ 10 −3 , because otherwise the rate of gravitational wave bursts due to initially binary systems would be too high (see the previous section). The rate of gravitational bursts in a unit volume is where B is the gain factor due to dynamic effects, discussed below, and the fraction of clusters in the dark matter composition where ∆ = δ 2 tot 1/2 . Observations give R b ∼ 10 yr −1 Gpc −3 [15]. Numerical calculation using the above equations gives the necessary value of perturbation and f cl < 1. If we ignore the dynamic effects and put B ∼ 1, we get δ tot f 14/31 cl ∼ 2. This means that even if we take a sufficiently large ∆ ∼ 0.1, it is impossible to get clusters of PBH that give the observed rate of gravitational bursts in this case.

IV. CONCLUSION
In this paper, within the framework of the Brans-Dicke theory, the PBH formation was studied. A region of space collapsing into a black hole was modelled as a part of a closed universe on the flat background. Corrections were found to the threshold of PBHs formation, which was previously calculated within the General relativity [4]. It turned out that these corrections depend on the scalar field of the Brans-Dicke theory and on its derivative taken at some initial time.
PBHs are formed in regions with sufficiently large curvature perturbations. If the scalar field is distributed statistically independent with respect to the distribution of curvature perturbations, then perturbations of the scalar field can modulate the number density of PBHs in the universe, leading to inhomogeneities or clustering of the PBHs. In the case of strong perturbations of scalar field, we can even expect the formation of virialized clusters of PBHs. The collisions of PBHs in clusters occur more frequently than on average, and this may affect the predicted rate of events for LIGO/Virgo detectors if some of the recorded events are explained by the PBHs. We showed that the effect of PBHs clustering is indeed the case in the Brans-Dicke theory. The formation of PBHs clusters in a minihalo of dark matter can significantly affect the rate of PBH mergers and the observed rate of LIGO/Virgo gravitational wave events.
PBH clustering can also influence the formation of structures in the dark matter at sub-galactic scales. This effect may be important for the formation of the first stars and massive black holes in gas clouds in the early universe, since the regions of clustered PBHs create seeds for large-scale density perturbations. An example of the formation of such an object from the dark matter with a cluster of PBHs in the center was discussed in the Section III C. However, a detailed study of this effect on a large-scale structure is beyond the scope of this paper. Note only that the results of the Section III C are applicable to clusters of PBHs that could have been formed due to the other possible mechanisms, and not only in the Brans-Dicke theory. The structure and the dynamical evolution of such clusters will be the same. It will be a PBH subsystem immersed in the dark matter halo and compressed due to the dynamic friction. In addition, in the outer regions, the formation of dark matter halos will continue due to the secondary accretion, as well as the baryonic matter will be captured. In these objects, the first stars and supermassive black holes could have appeared.
Calculations of [40] show that in the Brans-Dicke theory, the contribution of the entropy (constant curvature) perturbations generated at the inflation stage, to the total value of perturbations is negligible. In this paper, we have shown that the contribution can be amplified by about three orders of magnitude. The presence of such a contribution of entropy perturbations will not contradict the observation of anisotropy of relic radiation, if the perturbation spectrum in the Brans-Dicke scalar field is restricted by subgalactic scales only. At the same time, at larger scales, observed in cosmic microwave background, the differences with the standard picture may be below the level available for the current observations.