Formation and Clustering of Primordial Black Holes in Brans-Dicke Theory

Abstract

The formation of primordial black holes in the early universe in the 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.

1. 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], and 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 that are 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 that are 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 only formed from density perturbation that is greater than the 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 that are 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 it 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,28], 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. This mechanism is ineffective, as shown in [29]. 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 that is 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 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] 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 $\omega \to \infty$. 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 provide the following restrictions on the Brans-Dicke parameter $\omega >\mathrm{40,000}$ 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, in Brans-Dicke inflation theory, curvature perturbations and perturbations in a scalar field other than inflation can be independent, according to the theoretical calculations of [39,40]. 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.

2. Primordial Black Holes in the Brans-Dicke Theory

2.1. Basic Equations

We consider the Brans-Dicke theory with the Lagrangian

$$\mathcal{L}=\varphi R-\frac{\omega g^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi }{\varphi }+16\pi {\mathcal{L}}_m,$$

(1)

where $\omega =const$ and ${\mathcal{L}}_m$ is the Lagrangian of matter. In calculations, we assume that the speed of light $c=1$.

Following the approach of [4], we model a perturbed region of space collapsing to the PBH as a part of a closed universe with a metric

$$d{s}^2=d{\tau }^2-S^2\left(\frac{d{r}^2}{1-r^2}+r^{2}d{\Omega }^2\right),$$

(2)

where $S=S\left(\tau \right)$ is the scale factor. Additionally, the metric outside the perturbed region is the metric of the flat universe

$$d{s}^2=d{t}^2-R^2\left(d{r}^2+r^{2}d{\Omega }^2\right),$$

(3)

where $R=R\left(t\right)$, see Figure 1. Between the regions some transition intermediate layers are possible, see the 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 models have the form [41],

$$\frac{1}{R^2}{\left(\frac{dR}{dt}\right)}^2=\frac{8\pi }{3}\frac{1}{\overline{\varphi }}\overline{\epsilon }-\frac{\dot{\overline{\varphi }}}{\overline{\varphi }}\frac{\dot{R}}{R}+\frac{\omega }{6}\frac{\dot{\overline{{\varphi }^2}}}{\overline{{\varphi }^2}},$$

(4)

$$\frac{1}{S^2}{\left(\frac{dS}{d\tau }\right)}^2+\frac{1}{S^2}=\frac{8\pi }{3}\frac{1}{\varphi }\epsilon -\frac{\dot{\varphi }}{\varphi }\frac{\dot{S}}{S}+\frac{\omega }{6}\frac{\dot{{\varphi }^2}}{{\varphi }^2},$$

(5)

where, in the Equation (4), a point means a derivative over t, and in Equation (5) a point means a derivative over $\tau$ (and further similarly).

We set the initial data at some early time $t_i$ (${\tau }_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=\epsilon /3$); therefore, in the Equations (4) and (5), one has

$$\overline{\epsilon }={\overline{\epsilon }}_i\frac{R_i^4}{R^4},\epsilon ={\epsilon }_i\frac{S_i^4}{S^4}.$$

(6)

The scalar field evolves as $\dot{\overline{\varphi }}{R}^3=const$, $\dot{\varphi }{S}^3=const$ [41], so

$$\dot{\overline{\varphi }}={\dot{\overline{\varphi }}}_i\frac{R_i^3}{R^3},\dot{\varphi }={\dot{\varphi }}_i\frac{S_i^3}{S^3}.$$

(7)

Note that, for $\omega {\dot{\varphi }}_i{t}_i/{\varphi }_i\ll 1$, the last term on the right hand side Equation (5) can be ignored throughout cosmological evolution, but, in general, the last term may be important.

Additionally, the last closing equation

$$\frac{d\tau }{dt}=\frac{S}{R{(1+{\delta }_i)}^{1/4}}$$

(8)

gives a relation between $\tau$ and t [42]. This equation follows from the general expressions for the comoving coordinate system [43] and it is valid both in General relativity and in Brans-Dicke theory. The value of the density perturbation is defined as

$$\delta =\frac{\epsilon -\overline{\epsilon }}{\overline{\epsilon }}.$$

(9)

In Equation (8), the perturbation is taken at the initial moment $t_i$ (${\tau }_i$), and the initial constant-time hypersurface is chosen, so that, when $t=t_i$ ($\tau ={\tau }_i$), the following conditions are met [42]

$$S_i=R_i,{\left.dS/d\tau \right|}_{{\tau }_i}={\left.dR/dt\right|}_{t_i}.$$

(10)

Note that, in other works [44,45], 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,42]), 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 Equations (4) and (5) in the zero approximation at $\dot{\varphi }=0$ (coinciding with the General relativity), and we then solve the problem in the first linear approximation in Brans-Dicke theory.

2.2. Zero Approximation

In the zero approximation, $\dot{\varphi }=\dot{\overline{\varphi }}=0$ and $\varphi =\overline{\varphi }=1/G$, where G is the Newtonian gravitational constant. This approximation matches the results that were obtained in General relativity in [4]. Here, we will present the known results from [4] and write down some new relations for the evolution of $\delta$ and for the threshold of PBH formation.

In the zero approximation, we denote $A\equiv S^{\left(0\right)}$. The condition $dA/d\tau =0$ presented in Equation (5) gives the scale factor at the moment ${\tau }_{\max }$ of the maximum expansion of the disturbed region

$$A_{\max }=\frac{{(1+{\delta }_i)}^{1/2}{A}_i}{{\delta }_i^{1/2}}\approx \frac{A_i}{{\delta }_i^{1/2}},$$

(11)

because ${\delta }_i\ll 1$, and from Equation (10), we have

$${\delta }_i=4G{\varphi }_i\frac{t_i^2}{A_i^2}=4\frac{t_i^2}{A_i^2}.$$

(12)

The evolution of the scale factor in the zero approximation is described by the equation

$$\frac{dA}{d\tau }={\left(\frac{A_{\max }^2}{A^2}-1\right)}^{1/2}.$$

(13)

Introducing the parameter $\psi$, one obtains the solution in the parametric form [42]

$$A=A_{\max }sin\psi ,$$

(14)

where the relation of the parameter $\psi$ to $\tau$ has the form

$$\tau =A_{\max }(1-cos\psi ),$$

(15)

and the relationship between $\psi$ and t is as follows

$$\psi =arcsin\left(\frac{{\delta }_i^{1/2}}{{(1+{\delta }_i)}^{1/2}}\right)+\frac{{\delta }_i^{1/2}}{{(1+{\delta }_i)}^{1/4}}\left[{\left(\frac{t}{t_i}\right)}^{1/2}-1\right].$$

(16)

In the last expression, in contrast to [42], we took the next order corrections for ${\delta }_i$ into account, which will be required later in order to calculate $\delta \left(t\right)$.

The maximum expansion of the disturbed region corresponds to the moment of time

$$t_{\max }=t_i{\left(\frac{\pi }{2}\right)}^2\frac{{(1+{\delta }_i)}^{1/2}}{{\delta }_i},$$

(17)

and at the same time

$${\tau }_{\max }=t_{i}2\frac{{(1+{\delta }_i)}^{1/2}}{{\delta }_i}.$$

(18)

When selecting a constant-time hypersurface, the ratio of densities at the time of maximum expansion is $\epsilon /\overline{\epsilon }={(\pi /2)}^4$. In other works (for example, [44,45]), a different choice of the hypersurface was used, in which the time in an undisturbed flat region of the universe is $t=\tau$, and this choice results in $\epsilon /\overline{\epsilon }=4$, as already mentioned above.

The value of the perturbation is expressed, as follows

$$\delta =\frac{\epsilon -\overline{\epsilon }}{\overline{\epsilon }}={\left(\frac{t}{t_i}\right)}^2\frac{{\delta }_i^2}{(1+{\delta }_i)}\frac{1}{{sin}^4\psi }-1,$$

(19)

where $\psi$ is related to t by Equation (16). For $t\ll t_{\max }$, one obtains

$$\delta \approx \frac{1}{3}{\delta }_i{\left(\frac{t}{t_i}\right)}^{-1/2}+\frac{2}{3}{\delta }_i\left(\frac{t}{t_i}\right)$$

(20)

It means that, in this coordinate system, the falling perturbation mode evolves as $\propto t^{-1/2}$, in contrast to [45], where the law $\propto t^{-1}$ was obtained. Additionally, the evolution of a growing mode on scales larger than the horizon occurs according to the law $\propto t$ known for a synchronous reference frame.

The PBH formation criterion was analytically obtained in [4]. The region of space with $\delta >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 ${\delta }_{\mathrm{H}}$ satisfies the following conditions:

$${\delta }_c\le {\delta }_{\mathrm{H}}\le 1.,$$

(21)

In early studies, criterion for the threshold ${\delta }_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 Equation (21). Additionally, the right inequality corresponds to the formation of a black hole, not a separate universe. As noted in [46], this criterion contains ambiguity in the expression for the Jeans radius and can, therefore, only be considered as an estimation. 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 that was proposed in [46]. 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

$$d{s}^2=d{\tau }^2-A^2\left(d{\chi }^2+{sin}^2\chi d{\Omega }^2\right),$$

(22)

then the equation of the sound wave has the form [46]

$$A\frac{d\chi }{d\tau }=-\frac{1}{\sqrt{3}}.$$

(23)

Substituting the above expressions, we obtain

$$\frac{d\chi }{d\psi }=-\frac{1}{\sqrt{3}}.$$

(24)

Integrating up to the moment of maximum expansion $\psi =\pi /2$, we have

$${\chi }_{\max }=\frac{\pi }{2}\frac{1}{\sqrt{3}},$$

(25)

which, in the case of the equation of state $p=\epsilon /3$, matches the result obtained in [46]. However, we use a different choice of hypersurface, so the further calculation will be different. The conditions for the Jeans radius (sound horizon) and condition for cosmological horizon crossing make up the system of equations

$$R_{\mathrm{J}}=A_{\max }sin{\chi }_{\max }=A_{\max }sin{\chi }_a,$$

(26)

$$R_{\mathrm{H}}=2t_H=A_{\max }sin{\chi }_{a}sin{\psi }_{\mathrm{H}},$$

(27)

where ${\chi }_a$ is the boundary of the disturbed region and ${\psi }_{\mathrm{H}}$ is the value of the parameter $\psi$ at the moment of the horizon crossing. From these two equations, one obtains

$$2t_H=A_{\max }sin\left(\frac{\pi }{2}\frac{1}{\sqrt{3}}\right)sin{\psi }_{\mathrm{H}}.$$

(28)

In the limit ${\delta }_i\ll 1$, the Equation (28) is transformed into the following nonlinear equation

$$\xi =sin({\xi }^{1/2})sin\left(\frac{\pi }{2}\frac{1}{\sqrt{3}}\right),$$

(29)

where

$$\xi =\frac{{\delta }_i{t}_H}{t_i}.$$

(30)

Numerical solution gives the root of the above equation

$$\xi ={\xi }_0\simeq 0.52.$$

(31)

Substituting $t_H={\xi }_0{t}_i/{\delta }_i$ in Equation (19), we obtain the threshold for the PBH formation in the zero approximation

$${\delta }_c^{\left(0\right)}=0.42.$$

(32)

Note that, in other coordinate systems, the threshold value is different.

2.3. Corrections to the Black Hole’s Formation Threshold

In the following approximation, we take the evolution of the gravitational constant (change of the field $\varphi$) in the Brans-Dicke theory into account. We assume that all of the corrections to General relativity are small.

After the variable replacing, we rewrite the Equation (7), as follows

$$\frac{d\varphi }{dA}={\dot{\varphi }}_i\frac{A_i^3}{A^3}{\left(\frac{A_{\max }^2}{A^2}-1\right)}^{1/2}.$$

(33)

Its solution is

$$\varphi ={\varphi }_i+\frac{{\dot{\varphi }}_i{A}_i^3}{A_{\max }^2}\left(cot{\psi }_i-cot\psi \right).$$

(34)

Substituting expressions for the zero approximation on the right, we get at $\psi =\pi /2$

$${\varphi }_{max}={\varphi }_i+2{\dot{\varphi }}_i{t}_i.$$

(35)

Here, as before, the index “max” denotes the values at the moment when the expansion of the disturbed region stops, i.e., $dS/d\tau =0$.

We introduce small corrections to the zero approximation, as follows

$${\varphi }_i=G^{-1}(1+s_i),\varphi =G^{-1}(1+s_i+p),S=A(1+\alpha ),$$

(36)

where $s_i$ is the relative deviation of $\varphi$ from $G^{-1}$ at the initial moment, p and $\alpha$ are new small functions, of which the first was already calculated above in Equation (34). Let us write the Equation (5) that is linearised over $\alpha$:

$$2\frac{\dot{A}}{A}\dot{\alpha }-2\alpha \frac{1}{A^2}=-\frac{8\pi G}{3}{\epsilon }_i{S}_i^4\frac{1}{A^4}(s_i+p+4\alpha )-{\dot{\varphi }}_{i}G{S}_i^3\frac{\dot{A}}{A^4}+\frac{\omega }{6}{\dot{\varphi }}_i^2{G}^2{S}_i^6\frac{1}{A^6}.$$

(37)

The last term contains the small parameter ${\dot{\varphi }}_i{t}_{i}G$ squared, but, because of the possibility of large values of $\omega$, we keep it. Substituting $\dot{A}$ from the zero approximation and introducing the notation $x=A^2/A_{\max }^2$, we rewrite Equation (37), as

$$4\frac{d\alpha }{dx}(1-x)-2\alpha =-C_1\frac{s_i+p+4\alpha }{x}-C_2\frac{{(1-x)}^{1/2}}{x^{3/2}}+C_3\frac{1}{x^2},$$

(38)

where $C_{1,2,3}$ are some combinations of constants. General solution of the above equation is

$$\begin{array}{ccc}\hfill \alpha & =& {\alpha }_i{\left(\frac{x_i}{x}\right)}^{C_1}{\left(\frac{1-x}{1-x_i}\right)}^{(2C_1-1)/2}+\hfill \\ & +& \frac{1}{4}{x}^{-C_1}{\left(1-x\right)}^{(2C_1-1)/2}\underset{x_i}{\overset{x}{\int }}dy{y}^{C_1}{\left(1-y\right)}^{-(2C_1+1)/2}\left(-C_1\frac{s_i+p}{y}-C_2\frac{{(1-y)}^{1/2}}{y^{3/2}}+C_3\frac{1}{y^2}\right).\hfill \end{array}$$

(39)

Using the x variable, the solution Equation (34) is written, as follows

$$\varphi =G^{-1}\left[1+s_i+2\mu -2\mu {\delta }_i^{1/2}\frac{{(1-x)}^{1/2}}{x^{1/2}}\right].$$

(40)

It shows the structure of the function p.

From the zero approximation, we have $C_1\approx 1$, $C_2\approx 2\mu {\delta }_i^{1/2}$, $C_3\approx 2\omega {\mu }^2{\delta }_i$, where $\mu ={\dot{\varphi }}_i{t}_{i}G$. Using these values, after integration we obtain

$$\begin{array}{ccc}\hfill \alpha & \simeq & {\alpha }_i\frac{x_i}{x}{\left(\frac{1-x}{1-x_i}\right)}^{1/2}-(\mu +s_i/2-\omega {\mu }^2{\delta }_i/3)\left(x^{-1}-x^{-1}{(1-x)}^{1/2}{(1-x_i)}^{-1/2}\right)+\hfill \\ & +& \frac{1}{6}\omega {\mu }^2{\delta }_i{x}^{-1}{(1-x)}^{1/2}ln\left[\frac{1-{(1-x)}^{1/2}}{1-{(1-x_i)}^{1/2}}\frac{1+{(1-x_i)}^{1/2}}{1+{(1-x)}^{1/2}}\right].\hfill \end{array}$$

(41)

At the moment of maximum expansion $x=1$, one gets $\alpha \simeq -\mu -s_i/2+\omega {\mu }^2{\delta }_i/3$ from Equation (41). For verification, let us calculate the same value in a different way. From the condition $dS/d\tau =0$ in Equation (5), we obtain an algebraic equation of the fourth order for the value $S_{\max }/S_i$, and its solution is

$${\left(S_{\max }/S_i\right)}^2=\frac{4\pi S_i^2{\epsilon }_i}{3{\varphi }_{\max }}\left(1+\sqrt{1+\Sigma }\right),\Sigma =\frac{3\omega {\dot{\varphi }}_i^2}{32{\pi }^2{S}_i^2{\epsilon }_i^2}.$$

(42)

Substituting expressions for the zero approximation, we obtain the same value $\alpha \simeq -\mu -s_i/2+\omega {\mu }^2{\delta }_i/3$.

Instead of Equation (23), for the sound wave equation, we have the equation

$$A(1+\alpha )\frac{d\chi }{d\tau }=-\frac{1}{\sqrt{3}}.$$

(43)

Through the variable x it looks, as follows

$$d\chi \simeq -\frac{dx(1-\alpha )}{2\sqrt{3}{x}^{1/2}{(1+x)}^{1/2}}.$$

(44)

Note that, when integrating this equation, the size of the sound horizon for the first and third terms in Equation (41) is accumulated at early times $\tau$ (small x), and, for the second term, the size is accumulated at later times, just before the PBH formation. By performing the integration, we obtain the solution with the correction

$${\chi }_{\max }\simeq \frac{\pi }{2\sqrt{3}}\left(1+\gamma \right),$$

(45)

where was denoted

$$\gamma =\frac{1}{\pi }(s_i+2\mu )+\frac{2}{3\pi }\omega {\mu }^2{\delta }_i^{1/2}-\frac{2}{\pi }{\alpha }_i{\delta }_i^{1/2}.$$

(46)

Instead of Equation (29), we obtain the equation

$$\xi =sin({\xi }^{1/2}(1+\gamma ))sin\left(\frac{\pi }{2}\frac{1}{\sqrt{3}}(1+\gamma )\right).$$

(47)

The corrected solution can be found, as follows

$$\xi ={\xi }_0+\gamma {\left.\frac{d\xi }{d\gamma }\right|}_{\gamma =0},$$

(48)

where the derivative ${d\xi /d\gamma |}_{\gamma =0}\simeq 0.6255$ is calculated by differentiating the Equation (47). By a similar decomposition, we get the final correction for the threshold of PBH formation:

$${\delta }_c={\delta }_c^{\left(0\right)}(1+\lambda \gamma ),$$

(49)

where $\lambda \simeq 0.432$, and $\gamma$ is given above in Equation (46). The expression Equation (49) is the main result of this work. In the following sections, 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 Equation (46), the term in parentheses at the linear stage is invariant with respect to the choice of the initial moment and it is equal to $(G{\varphi }_{\max }-1)$, which follows from Equations (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 for demonstrating the effect of PBH clustering.

Note that the corrections to the threshold in the Eddington-inspired-Born-Infeld gravity [25] depend on the epoch of PBH formation. However, our result (Equations (46) and (49)) is independent of the epoch of collapse. In Brans-Dicke gravity, the threshold only depends on the initial conditions.

In this section, we have considered small corrections to General relativity from Brans-Dicke theory. It is interesting to qualitatively discuss what would happen in the case of larger corrections. The equations for the metric in the static spherically-symmetric case, if $\dot{\varphi }=0$, was obtained in [41] (Equations (9.9.24)–(9.9.26)). In this case, there is no new characteristic scale, but only power-law corrections to the post-Newtonian approximation. The situation is different when $\dot{\varphi }\ne 0$. If, for example, the term $1/S^2$ in Equation (5) is of the order of the last term, one has the characteristic scale

$$S_D\sim c{|\dot{\varphi }/\varphi |}^{-1}{(6/\omega )}^{1/2}.$$

(50)

One could formally relate this scale to the radius at which a flat rotation curve begins to appear in galaxies, which is attributed to dark matter in General relativity. According to astronomical constraints, $\dot{G}/G<{10}^{-11}$ yr${}^{-1}$. Therefore,

$$S_D>4\times {10}^2{\left(\omega /40000\right)}^{-1/2}\mathrm{Mpc},$$

(51)

where $G\sim {\varphi }^{-1}$. This scale is of the order of large galaxy superclusters scale. In modified gravity theories [47,48,49], other expressions for the new scale were obtained in terms of parameters of a particular theory and the new scales could be smaller.

3. Clustering of Black Holes

3.1. Perturbations in the PBH Number Density

To illustrate the effect, let us consider a simple model in which PBHs have a monochromatic mass spectrum, i.e., all PBHs are formed with the same mass $M_{\mathrm{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 ${\sigma }_{\mathrm{H}}$. For simplicity, let the scalar field have perturbations on a larger scale, which currently contains a mass of cold dark matter $M\gg M_{\mathrm{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 Equation (21). The probability of PBH formation, i.e., the fraction of radiation that transformed into the PBHs at the time of their formation, is written as

$$\beta \left({\delta }_c\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{H}}}\underset{{\delta }_c}{\overset{1}{\int }}d\delta exp\left(-\frac{{\delta }^2}{2{\sigma }_{\mathrm{H}}^2}\right)\simeq \frac{{\sigma }_{\mathrm{H}}}{{\delta }_c\sqrt{2\pi }}{e}^{-{\delta }_c^2/\left(2{\sigma }_{\mathrm{H}}^2\right)},$$

(52)

where the last expression is obtained by considering the tail of the Gaussian distribution $\delta \ge {\delta }_c\gg {\sigma }_{\mathrm{H}}$, and the upper limit of integration is not important. Nowadays, the cosmological parameter of PBHs

$${\Omega }_{\mathrm{PBH}}\simeq \beta \frac{a\left(t_{\mathrm{eq}}\right)}{a\left(t_{\max }\right)}.$$

(53)

At the radiation-dominated stage, the scale factor of the universe is $a\left(t\right)\propto t^{1/2}$.

If one assumes that the fraction ∼1 of all LIGO/Virgo events is due to PBHs with masses $M_{\mathrm{PBH}}\sim 30M_{\odot }$, then the PBHs constitute the fraction ${\Omega }_{\mathrm{PBH}}/{\Omega }_m\sim {10}^{-3}$ of all dark matter [15], where ${\Omega }_m$ is the cosmological density parameter of dark matter. From Equations (52) and (53), one obtains ${\sigma }_{\mathrm{H}}\simeq 0.013$ in this case.

Now, we will show that perturbations in the scalar field that give a correction in Equation (49) are translated into large-scale perturbations in the PBH number density, i.e., biasing effect is in place. The inhomogeneities of the scalar field and its time derivative in combination Equation (46) modulate the PBH distribution. Moreover, it turns out that perturbations in the PBHs can be significantly amplified in comparison with perturbations in the scalar field due to the threshold of PBH formation at the tail of the Gaussian distribution.

Substitute the value ${\delta }_c$ from Equation (49) to Equation (52). Subsequently, the perturbation of the PBH number density is written as

$${\delta }_{\mathrm{PBH}}\left[\gamma \right]=\frac{\beta \left({\delta }_c\left[\gamma \right]\right)}{\beta \left({\delta }_c[\gamma =0]\right)}-1\simeq -\frac{{({\delta }_c^{\left(0\right)})}^2}{{\sigma }_{\mathrm{H}}^2}\lambda \gamma .$$

(54)

The minus sign is explained by the fact that a positive value of $\gamma$ increases the threshold, and PBHs form in a smaller number, and, in the case of a negative $\gamma$, the opposite situation occurs. The same is valid for the mean root square values

$${\langle {\delta }_{\mathrm{PBH}}^2\rangle }^{1/2}\simeq \frac{{({\delta }_c^{\left(0\right)})}^2}{{\sigma }_{\mathrm{H}}^2}\lambda {\langle {\gamma }^2\rangle }^{1/2}.$$

(55)

Because ${({\delta }_c^{\left(0\right)})}^2/{\sigma }_{\mathrm{H}}^2\gg 1$, one has ${\delta }_{\mathrm{PBH}}\gg \lambda \gamma$, which is, the perturbations in the PBH number density are much larger than the initial perturbations in the scalar field. For the example above ($M_{\mathrm{PBH}}\sim 30M_{\odot }$, ${\Omega }_{\mathrm{PBH}}/{\Omega }_m\sim {10}^{-3}$), we have ${({\delta }_c^{\left(0\right)})}^2/{\sigma }_{\mathrm{H}}^2\sim {10}^3$, i.e., fluctuations of the scalar field are translated into fluctuations of the PBHs, amplified by three orders of magnitude, according to Equation (54). This is the biasing effect in application to the process of PBH formation in the Brans-Dicke theory.

However, the PBHs only constitute a small part of all dark matter, so the total value of the density perturbation is

$${\delta }_{\mathrm{tot}}\simeq -\frac{{\Omega }_{\mathrm{PBH}}}{{\Omega }_m}\frac{{({\delta }_c^{\left(0\right)})}^2}{{\sigma }_{\mathrm{H}}^2}\lambda \gamma \sim -\lambda \gamma ,$$

(56)

where ${\Omega }_{\mathrm{PBH}}$ is given by Equation (53) 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 was in the case Equation (54). The clustering effect, which we consider, may also arise in other theories of gravity with additional degrees of freedom, for example, with new scalar fields. The PBH distribution can be modulated by the large-scale inhomogeneities, if they exist through these degrees of freedom.

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 more easily overcome the ${\delta }_c$ threshold than on average. The usual bias effect applied to galaxies in galaxy clusters also consists in the fact that, in a region with high density, galaxies are formed more likely than in regions with low density [26,50]. In [30], a conformal Newtonian coordinate system was used, in which density perturbations at the radiation-dominated stage do not evolve and they 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 $\delta \propto 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 therefore, the biasing effect is 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].

3.2. Influence of Inhomogeneities on the Rate of Gravitational Bursts

The number density of PBHs is increased in regions where $\gamma$ has negative sign, see Equation (54). Let us 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 mainly determined 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

$$R_b={\left.\frac{{\rho }_c{\Omega }_{m}f}{M_{\mathrm{BH}}}\frac{dP(<t)}{dt}\right|}_{t=t_0},$$

(57)

where ${\rho }_c=9.3\times {10}^{-30}$ g cm${}^{-3}$ is the critical density, ${\Omega }_m\approx 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

$$P(<t)\sim \left[\frac{37}{29}{\left(\frac{t}{t_{\max }}\right)}^{3/37}-\frac{8}{29}{\left(\frac{t}{t_{\max }}\right)}^{3/8}\right],$$

(58)

where

$$t_{\max }\sim \frac{5}{512G^3{M}_{\mathrm{BH}}^3}\frac{{\overline{x}}^4}{f^4},$$

(59)

the average distance between the PBH is

$$\overline{x}={\left(\frac{M_{\mathrm{BH}}}{f{\rho }_{\mathrm{eq}}}\right)}^{1/3},$$

(60)

where ${\rho }_{\mathrm{eq}}$ is the average density of dark matter at the moment of transition to the dust-like stage of the universe evolution.

Because $t_{\max }\gg t_0$ is valid for LIGO/Virgo signal parameters, the first term in Equation (58) gives the dominant contribution to the merge rate. Then for a fixed f we have $R_b\propto {\overline{x}}^{-12/37}$. Because, for a non-uniform distribution of the PBH $\overline{x}\propto {(1+{\delta }_{\mathrm{PBH}})}^{-1/3}$, where ${\delta }_{\mathrm{PBH}}$ is set by Equation (54), the final expression $R_b\propto {(1+{\delta }_{\mathrm{PBH}})}^{-4/37}$ only depends very weakly on ${\delta }_{\mathrm{PBH}}$. Because it seems that the extremely large values of $\lambda \gamma$ 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.

3.3. Clusters of Primordial Black Holes

Perturbations of PBH number density Equation (54), which arose due to inhomogeneities of the scalar field, and the total perturbations Equation (56), 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. Additionally, after the onset of the cosmological stage of matter dominance, these perturbations evolve in the same way as the usual perturbations in dark matter $\propto 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 us 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 [51]

$${\sigma }_{\mathrm{cap}}\approx \frac{3}{2}\pi r_g^2{\left(\frac{c}{v}\right)}^{18/7},$$

(61)

where $r_g\equiv 2G{M}_{\mathrm{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 mass of the rest of dark matter)

$$\dot{N}\simeq (1/2)Nn{\sigma }_{\mathrm{cap}}v=9\sqrt{2}{({\Omega }_{\mathrm{PBH}}/{\Omega }_m)}^2{\left(\frac{v}{c}\right)}^{17/7}\frac{c}{R},$$

(62)

where, n is the number density of PBHs, $N=M({\Omega }_{\mathrm{PBH}}/{\Omega }_m)/M_{\mathrm{PBH}}\simeq (4\pi /3)R^{3}n$ is the number of PBHs in the cluster, $v\sim {(GM/R)}^{1/2}$, R is the radius of the cluster,

$$R\simeq {\left(\frac{3M}{4\pi {\rho }_{\mathrm{cl}}}\right)}^{1/3},$$

(63)

where, according to the spherical model of perturbation evolution, ${\rho }_{\mathrm{cl}}\simeq 18{\pi }^2\overline{\rho }\left(t_f\right)$, $\overline{\rho }\left(t\right)$ is the average cosmological density of dark matter at time t, and $t_f$ is the moment when the perturbation that is evolving according to the law $\propto t^{2/3}$ reaches the value $\simeq 2.81$.

The value ${\Omega }_{\mathrm{PBH}}/{\Omega }_m\le {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

$$R_b=\dot{N}B\frac{{\rho }_c{\Omega }_m{f}_{\mathrm{cl}}}{M},$$

(64)

where B is the gain factor due to dynamic effects, as discussed below, and the fraction of clusters in the dark matter composition

$$f_{\mathrm{cl}}=\frac{1}{\sqrt{2\pi }\Delta }\underset{{\delta }_{\mathrm{tot}}}{\overset{\infty }{\int }}d\delta exp\left(-\frac{{\delta }^2}{2{\Delta }^2}\right)=\frac{1}{2}\left[1-\mathrm{Erf}\left(\frac{{\delta }_{\mathrm{tot}}}{\sqrt{2}\Delta }\right)\right],$$

(65)

where $\Delta ={\langle {\delta }_{\mathrm{tot}}^2\rangle }^{1/2}$. The observations give $R_b\sim 10$ yr${}^{-1}$ Gpc${}^{-3}$ [15].

Numerical calculation using the above equations gives the necessary value of perturbation

$${\delta }_{\mathrm{tot}}{f}_{\mathrm{cl}}^{14/31}\simeq 0.1{\left(\frac{M}{{10}^5{M}_{\odot }}\right)}^{22/93}{\left(\frac{R_b}{10{\mathrm{yr}}^{-1}{\mathrm{Gpc}}^{-3}}\right)}^{14/31}{\left(\frac{{\Omega }_{\mathrm{PBH}}/{\Omega }_m}{0.002}\right)}^{-31/28}{\left(\frac{B}{1000}\right)}^{-14/31},$$

(66)

and $f_{\mathrm{cl}}<1$. If we ignore the dynamic effects and put $B\sim 1$, we obtain ${\delta }_{\mathrm{tot}}{f}_{\mathrm{cl}}^{14/31}\sim 2$. This means that, even if we take a sufficiently large $\Delta \sim 0.1$, it is impossible to obtain clusters of PBH that give the observed rate of gravitational bursts in this case.

However, there are two dynamical evolutionary effects that can increase the rate of PBH collisions in clusters. The first effect is the dynamical friction of PBHs in the medium of dark matter in a cluster. Because of dynamical friction, the PBH lose kinetic energy and settle closer to the center of the cluster, see Figure 2. As a result, a more compact PBH subsystem is formed inside the cluster, which mainly consists of dark matter, and the rate of merge increases. The second effect is the dynamical evolution of the PBH cluster due to two-body approaches. This evolution occurs over several dozen relaxation times and also leads to the compression of the central part of PBH cluster and a final increase in the merge rate.

The objects under consideration have a certain density profiles, which we assume to be the profiles of an isothermal sphere with density $\propto r^{-2}$. In this case, the characteristic PBH velocities are the same at any r. But under the influence of dynamical friction, PBHs lose kinetic energy and their orbit contracts. The law of compression is estimated as follows

$$r\left(t\right)=R{\left(1-\frac{t}{t_{\mathrm{df}}}\right)}^{1/2},$$

(67)

where the characteristic time

$$t_{\mathrm{df}}\sim 0.5t_{\mathrm{eq}}\frac{M}{M_{\mathrm{PBH}}}{\delta }_{\mathrm{tot}}^{-3/2}.$$

(68)

It can be seen from Equation (62) that the compression of the PBH subsystem leads to a proportional increase in the merge rate. Accordingly, B is equal to the compression ratio ${(r\left(t\right)/R)}^{-1}$. However, the numerical simulations are required in order to accurately calculate all of the processes that accompany the cluster compression. The calculation of this section should be considered as a rough estimate. The normalizing value $B\sim {10}^3$ in Equation (66) is chosen from the condition that, when the PBH cluster in the isothermal profile $M\left(r\right)\propto r$ is compressed by three orders of magnitude in radius, the PBH gravity in the central region of the object will begin to prevail over the dark matter gravity, since ${\Omega }_{\mathrm{PBH}}/{\Omega }_m\sim {10}^{-3}$. From the condition $t_{\mathrm{df}}<t_0$, we obtain

$${\delta }_{\mathrm{tot}}>0.044{\left(\frac{M}{{10}^5{M}_{\odot }}\right)}^{2/3}{\left(\frac{M_{\mathrm{PBH}}}{30M_{\odot }}\right)}^{-2/3}.$$

(69)

Two-body relaxation leads to dynamical evaporation of the outer part of the PBH cluster and further compression of its central core. In general, this will cause an additional increase in the rate of PBH merge. Thus, it can be seen that, for ${\delta }_{\mathrm{tot}}\sim \Delta \sim 0.1$, the conditions are reached under which the rate of PBH merge in clusters begins to prevail over the rate of their merge in the absence of perturbations of scalar field of the Brans-Dicke theory. The exact calculation of this effect is difficult, due to the presence of complex dynamic evolution of the PBH cluster. The possibility of reaching $\Delta \sim 0.1$, according to Equation (56), depends on the value of perturbations in the scalar field of the Brans–Dicke theory.

4. Discussion

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 that was 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 the perturbations of the scalar field can modulate the number density of PBHs in the universe, which leads 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 in 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 Section 3.3. However, a detailed study of this effect on a large-scale structure is beyond the scope of this paper. Note that only the results of the Section 3.3 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 [39] 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 only restricted by subgalactic scales. At the same time, at larger scales, as observed in cosmic microwave background, the differences with the standard picture may be below the level that is available for the current observations.