Primordial Black Holes from Broad Spectra: Abundance and Clustering

TL;DR

This paper investigates the formation and distribution of primordial black holes generated by broad scalar perturbation spectra during inflation, concluding they are rare, predominantly Poissonian, and have similar mass distributions to narrow spectra.

Contribution

It provides a detailed analysis of primordial black hole formation and clustering for broad spectra, showing the rarity of their formation and their mostly Poissonian distribution.

Findings

Primordial black holes from broad spectra are rare events.

Their mass distribution is similar to that from narrow spectra.

They are not significantly clustered at formation, with a Poissonian distribution.

Abstract

A common mechanism to form primordial black holes in the early universe is by enhancing at small-scales the scalar perturbations generated during inflation. If these fluctuations have a large enough amplitude, they may collapse into primordial black holes upon horizon re-entry. Such primordial black holes may comprise the totality of the dark matter. We offer some considerations about the formation and clustering of primordial black holes when the scalar perturbations are characterised by a broad spectrum. Using the excursion set method, as well as the supreme statistics, we show that the cloud-in-cloud phenomenon, for which small mass primordial black holes may be absorbed by bigger mass ones, is basically absent. This is due to the fact that the formation of a primordial black hole is an extremely rare event. We also show that, from the point of view of mass distribution, broad and narrow spectra give similar results in the sense that the mass distribution is tilted towards a single mass. Furthermore, we argue that primordial black holes from broad spectra are not clustered at formation, their distribution is dominantly Poissonian.

Equations

$$C(r)=\frac{2\left[M(r,t)-M_b(r,t)\right]}{R(r,t)}.$$

$$\overline{\delta }(r)={\delta }_0\frac{\xi (r)}{{\sigma }^2},$$

$$\xi (r)=\int \frac{\mathrm{d}k}{k}\frac{\sin kr}{kr}{\mathcal{P}}_{\delta }(k,t),{\mathcal{P}}_{\delta }(k,t)=\frac{16}{81}{\left(\frac{k}{aH}\right)}^4{\mathcal{P}}_{\zeta }(k),$$

$${\delta }_c=3\overline{\delta }(r_m).$$

$$\beta (M(r_m))={\int }_{{\delta }_c}\frac{\mathrm{d}{\delta }_m}{\sqrt{2\pi }{\sigma }_m}{e}^{-{\delta }_m^2/2{\sigma }_m^2}\simeq 6\cdot 10^{-9}{\left(\frac{M(r_m)}{M_{\odot }}\right)}^{1/2},$$

$${\delta }_m$$

$${\sigma }_m^2$$

$$T(k,r_m)=3\frac{\sin (k{r}_m/\sqrt{3})-(k{r}_m/\sqrt{3})\cos (k{r}_m/\sqrt{3})}{(k{r}_m/\sqrt{3}{)}^3}.$$

$${\mathcal{P}}_{\zeta }={\mathcal{P}}_0{\delta }_D(\ln k-\ln k_{\star }).$$

$$\delta (\vec{x},t)={\delta }_l(\vec{x},t)+{\delta }_s(\vec{x},t)={\int }_{k_l}^{k_{\mathrm{ref}}}\frac{{\mathrm{d}}^{3}k}{{\left(2\pi \right)}^3}{\delta }_{\vec{k}}{e}^{-i\vec{k}\cdot \vec{x}}+{\int }_{k_{\mathrm{ref}}}^{k_s}\frac{{\mathrm{d}}^{3}k}{{\left(2\pi \right)}^3}{\delta }_{\vec{k}}{e}^{-i\vec{k}\cdot \vec{x}},k_l\ll k_{\mathrm{ref}}\ll k_s,$$

$${\mathcal{P}}_{\zeta }\approx {\mathcal{P}}_0(k/k_s{)}^{n_p}\Theta \left(k_s-k\right)\Theta \left(k-k_l\right),$$

$$\delta (\vec{x},R)=\int {\mathrm{d}}^3{x}^{\prime }W(|\vec{x}-{\vec{x}}^{\prime }|,R)\delta ({\vec{x}}^{\prime }),$$

$$\delta (\vec{k},R)=\tilde{W}(\vec{k},R)\delta (\vec{k})$$

$${\tilde{W}}_{\mathrm{sharp}-k}(k,k_f)=\Theta (k_f-k),$$

$$\frac{\partial {\delta }_R}{\partial R}=\zeta (R),$$

$$\zeta (R)\equiv \int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}\tilde{\delta }(\vec{k})\frac{\partial \tilde{W}(k,R)}{\partial R}.$$

$$⟨\zeta (R_1)\zeta (R_2)⟩={\int }_{-\infty }^{\infty }\mathrm{d}\ln k{\Delta }^2(k)\frac{\partial \tilde{W}(k,R_1)}{\partial R_1}\frac{\partial \tilde{W}(k,R_2)}{\partial R_2},$$

$$\frac{\partial \delta (k_f)}{\partial \ln k_f}=Q(k_f),$$

$$⟨Q({k_f}_1)Q({k_f}_2)⟩={\Delta }^2({k_f}_1){\delta }_D(\ln {k_f}_1-\ln {k_f}_2).$$

$$S(R)={\int }_{-\infty }^{\infty }\mathrm{d}\ln k{\Delta }^2(k)|\tilde{W}(k,R){|}^2,$$

$$S(k_f)={\int }_{-\infty }^{\ln k_f}\mathrm{d}\ln k{\Delta }^2(k).$$

$$\frac{\partial \delta (S)}{\partial S}=\eta (S),$$

$$⟨\eta (S_1)\eta (S_2)⟩={\delta }_D(S_1-S_2).$$

$$\frac{\partial P}{\partial S}=\frac{1}{2}\frac{{\partial }^{2}P}{\partial {\delta }^2}.$$

$${P(\delta ,S)|}_{\delta ={\delta }_c}=0,$$

$$P(\delta ,S)=\frac{1}{\sqrt{2\pi S}}\left[e^{-{\delta }^2/(2S)}-e^{-(2{\delta }_c-\delta {)}^2/(2S)}\right].$$

$$P(\delta >{\delta }_c)=1-{\int }_{-\infty }^{{\delta }_c}\mathrm{d}\delta P(\delta ,S)=\mathrm{Erfc}\left(\frac{{\delta }_c}{\sqrt{2S}}\right),$$

$$\omega (a_1)={\left(\frac{a_2}{a_1}\right)}^2\omega (a_2).$$

$$P(S_2,{\omega }_2|S_1,{\omega }_1)\mathrm{d}{S}_2$$

$$S(R_1)=S(R_2){\left(\frac{M_{H_1}}{M_{H_2}}\right)}^{-(4+n_p)/2}.$$

$${\mathcal{P}}_{\zeta }\approx {\mathcal{P}}_0(k/k_s{)}^{n_p}\Theta \left(k_s-k\right)\Theta \left(k-k_l\right),n_p>0,$$

Full text

\newdateformat

mydate\THEDAY \monthname[\THEMONTH] \THEYEAR

**Primordial Black Holes from Broad Spectra: Abundance and Clustering **

Azadeh Moradinezhad Dizgah 1,2, Gabriele Franciolini1, and Antonio Riotto1,3

1* Département de Physique Théorique and Centre for Astroparticle Physics (CAP),

Université de Genève, 24 quai E. Ansermet, CH-1211 Geneva, Switzerland*

2* Department of Physics, Harvard University, 17 Oxford Street, Cambridge, MA 02138, USA*

3* CERN, Theoretical Physics Department, Geneva, Switzerland*

Abstract

A common mechanism to form primordial black holes in the early universe is by enhancing at small-scales the scalar perturbations generated during inflation. If these fluctuations have a large enough amplitude, they may collapse into primordial black holes upon horizon re-entry. Such primordial black holes may comprise the totality of the dark matter. We offer some considerations about the formation and clustering of primordial black holes when the scalar perturbations are characterised by a broad spectrum. Using the excursion set method, as well as the supreme statistics, we show that the cloud-in-cloud phenomenon, for which small mass primordial black holes may be absorbed by bigger mass ones, is basically absent. This is due to the fact that the formation of a primordial black hole is an extremely rare event. We also show that, from the point of view of mass distribution, broad and narrow spectra give similar results in the sense that the mass distribution is tilted towards a single mass. Furthermore, we argue that primordial black holes from broad spectra are not clustered at formation, their distribution is dominantly Poissonian.

††Email: Azadeh.MoradinezhadDizgah@unige.ch, Gabriele.Franciolini@unige.ch, Antonio.Riotto@unige.ch

1 Introduction

The detection of gravitational waves produced by the merger of two $\sim 30\,M_{\odot}$ black holes [1] has started the era of gravitational wave astronomy. At the same time, it has revived the interest in the physics of Primordial Black Holes (PBHs) and the exciting hypothesis that all (or a significant part) of the dark matter of the universe might be in the form of PBHs [2] (see Ref. [3] for a review and additional references).

PBHs in the early universe might owe their origin to the small-scale enhancement of the primordial curvature perturbation $\zeta$ generated during inflation [4, 5, 6]. After inflation, during reheating, these small-scale perturbations are communicated to the radiation fluid. PBHs (of mass roughly equal to the horizon mass) are generated at horizon re-entry, if the density perturbation $\delta$ has a large enough amplitude, and hence the radiation density contrast is larger than a critical value $\delta_{c}$ [7].

The abundance of PBHs is estimated using either the peak theory [8] or the Press-Schechter (PS) approach [9]. Assuming spherical symmetry, to determine the characteristic scale and amplitude of the fluctuations that collapse to form a PBH, one defines the compaction function $C(r)$ as the ratio of mass-excess within a sphere of radius $r$ to the areal radius at $r$, [10, 7]

[TABLE]

Here $R(r,t)$ is the areal radius, $M(r,t)$ is the mass within a sphere of radius $r$ centred on the peak location, and $M_{b}(r,t)$ is the background mass within the same radius, calculated with respect to a flat FRW Universe. The criterion for collapse widely used in the recent literature states that a PBH forms if the maximum of the compaction function is above a certain threshold, found numerically [7]. It can be shown that the amplitude of the fluctuations in terms of the excess mass within a spherical volume is equivalent to the local value of compaction function [7]. Therefore the scale of fluctuations relevant for the formation of PBHs, commonly referred to as $r_{m}$, can be determined by maximizing the compaction function. Given that the PBHs form from large-amplitude fluctuations, one can apply the peak statistics [8], assuming spherical peaks (which is a good approximations since the PBHs are rare events), to calculate $r_{m}$.

Considering the average shape around a peak of height $\delta_{0}$ [11]

[TABLE]

with $\xi(r)$ being the two-point correlator in real space if working in the comoving slicing,

[TABLE]

the threshold for the formation of a PBH is defined in terms of volume-averaged mass-excess, which is three times the peak profile at $r_{m}$ [7]

[TABLE]

The corresponding critical peaks amplitude, $\delta_{0}$, has to be found numerically.

Once the critical threshold is known, the abundance of PBHs can be computed using threshold statistics. For Gaussian fluctuations222For the non-Gaussian extension, see [12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 21]. The non-Gaussian corrections are going to be neglected in this paper and the study of their impact on the present discussion is left for future work., the mass fraction of PBHs at the time of collapse, $\beta$, is given by

[TABLE]

where in the last passage we have assumed that the PBHs contribute to the totality of the dark matter, and the volume-averaged field and its variance are given by

[TABLE]

with $T(k,r_{m})$ being the radiation transfer function

[TABLE]

We have conveniently multiplied the smoothed overdensity by the the factor $(aHr_{m})^{2}$ to obtain a time-independent quantity which can be thus calculated on superhorizon scales [7, 23].

The dynamics of PBH formation, the PBH mass and mass function are rather well understood in the case of a narrow curvature perturbation power spectrum, where only one length scale $1/k_{\star}$ dominates. A benchmark power spectrum often used in the literature is the Dirac-delta type

[TABLE]

In such a case the scale $r_{m}$ is the order of $1/k_{\star}$ [24] and PBHs are formed with a mass $\sim M(1/k_{\star})$.

For broad spectra things are more blurred. Consider a power spectrum which extends over a large range of scales, say from $k_{l}$ to $k_{s}\gg k_{l}$. In a realistic universe there are many perturbations superimposed on top of each other with varying scales, amplitudes, and physical locations (i.e. different Fourier modes with varying wavelength, amplitude and phase). One therefore expects peaks with various amplitude and, above all, with different length scales, running from $\sim 1/k_{l}$ to $\sim 1/k_{s}$. Each peak profile will have its own scale $r_{m}$ associated with it.

Several questions will then arise; What is the distribution of PBH masses for broad spectra? Are the small PBHs swallowed by bigger PBHs, i.e. is there a so-called cloud-in-cloud problem, similar to that for the dark matter halo formation? Are PBHs clustered differently from a Poisson distribution, contrary to what happens for a narrow power spectrum [37]? The goal of the paper is to offer some considerations about these questions. Our basic findings are the following.

Using the excursion set method (with further support from an alternative technique based on the supreme statistics), we will show that the cloud-in-cloud phenomenon is basically absent for PBHs. This is due to the fact that the probability of forming a PBH (if we wish them to compose all or most of the dark matter) is very small, i.e. PBHs are very rare events 333If PBHs constitute only a small fraction of the dark matter, their mass fraction is smaller, and hence they are even rarer events. Therefore, our results hold even stronger.. This in practice prevents pre-existing PBHs to be swallowed by bigger ones. Also, despite the broadness of the power spectrum of the curvature perturbation, the production of PBHs is tilted towards either small-scale (small mass) or large-scale (large mass) PBHs. In this sense, the mass spectrum of PBHs does not differ from the one derived for a narrow power spectrum. This result also suggests that the clustering properties at formation will be similar for broad and narrow spectra: PBHs at their birth are likely to be distributed randomly following a Poisson distribution.

The paper is organized as follows. In section 2 we will discuss the threshold for formation of PBHs in the case of a broad spectrum. Section 3 is a brief summary of the excursion set method and its consequences when applied to the physics of PBHs formation. Section 4 supports some of the results found in the previous section through the so-called supreme statistics. Section 5 is devoted to the clustering, and section 6 contains our conclusions.

2 Broad spectra and their threshold

The first question to ask for broad spectra is what is the critical threshold. Here we follow Ref. [23] and imagine to separate the overdensity into long- (background) and short-modes (foreground)

[TABLE]

where, for simplicity, we suppose that both the long and short mode power spectra are flat, even though they can have different amplitudes. The operation of selecting the longer modes can be in principle done by a smoothing444One should not give too much emphasis on the use of the smoothing to determine the average profile shape. The physical evolution of a perturbation depends on the profile shape, but it should not depend on the smoothing function. The calculation of the abundance does depend on the choice of window function, but this just highlitghts the uncertainty in the calculation rather than suggesting different amounts of PBHs form [23]. over volumes of radius $\sim 1/k_{\rm ref}$. Since what matters is always the smallest length scale [7], computing the average profile of the short-mode perturbations delivers a value $r_{m,s}$ which is approximately given by $3.5/k_{s}$ [24]. Clearly the long-mode perturbations have a much larger characteristic scale $r_{m,l}\sim 3.5/k_{\rm ref}$.

Although the values of $r_{m}$ for short and long mode differ, it is important to note that the critical threshold value is the same for both cases and is $\delta_{c}\simeq 0.51$. In figure 1 we show the typical profile in Eq. (1.2), for a generalised broad power spectrum of the form

[TABLE]

where $\Theta$ is the Heaviside step function, $k_{s}=10^{4}\,k_{l}$ and $n_{p}=-2,-1,0,1$. The approximate sign in Eq. (2.2) reminds us that transition from 0 to $k_{l}$ and from $k_{s}$ from 0 is not sudden in real cases.

The plot shows that the profiles are nearly identical within their corresponding distance $r_{m}$ from the centre, which leads to the same critical threshold (see also [23]). We conclude that for broad spectra, peaks with different values of $r_{m}$ will eventually collapse if the volume-averaged over a sphere of radius $r_{m}$ overcomes the same threshold $\delta_{c}$.

Notice that causality prevents the large-scale perturbations longer than the Hubble horizon from affecting the formation of light PBH on short-scales. On the other side such PBHs formed at small-scales do not affect the fluctuations on large-scales after they re-enter the horizon. Indeed, and more in general, the effect of a set of self-gravitating particles in a virialized cluster does not alter the behaviour of the perturbations on much larger scales. In this sense, gravitational interactions in the cosmological setting are renormalizable.

3 The excursion set method and PBHs

Having argued that peaks with different characteristic scales $r_{m}$ have the same threshold, the next question is what are the masses which count as far as the PBH formation is concerned.

As already mentioned in the introduction, the abundance of PBHs can be determined by the Press-Schecter method, very much as done when computing the mass function of dark matter halos. An extension of it, the so-called excursion set theory [25], corrects the original argument of Press and Schechter which miscounts the number of virialized dark matter objects because of the so-called cloud-in-cloud problem. In the dark matter spherical collapse model one assumes that a region of radius $R$, with a smoothed density contrast $\delta_{R}$, collapses and virializes once $\delta_{R}$ exceeds a critical value. This procedure misses the cases in which, on a given smoothing scale $R$, $\delta_{R}$ is below the threshold, but still it happened to be above the threshold at some scale $R^{\prime}>R$. Such a configuration corresponds to a virialized object of mass $M^{\prime}>M$.

In the excursion set method the density perturbations evolve stochastically with the smoothing scale, and the problem of computing the probability of halo formation is mapped into the so-called first-passage time problem in the presence of a barrier. The cloud-in-cloud problem is avoided by counting only the trajectories which pass the threshold for the first-time.

The excursion set method is therefore well-suited to develop a better understanding of formation of PBHs form a broad spectrum, above all the cloud-in-cloud problem. For the reader not familiar with the excursion set theory, we provide a brief overview of the method. The expert reader might wish to skip this part and go directly to section 3.2.

3.1 A brief review of the excursion set method

Consider the density contrast $\delta({\vec{x}})$, smoothed on a particular scale $R$ as

[TABLE]

where $W(|{\vec{x}}-{\vec{x}}^{\prime}|,R)$ is the window function. Since the convolution in real space is a product in Fourier space, the Fourier transform of the smoothed field is

[TABLE]

with $\widetilde{W}({\vec{k}},R)$ being the Fourier transform of the real-space window function. There are various choice for the window function. As we discuss below, for a sharp filter in $k$-space, defined as

[TABLE]

with $k_{f}=1/R$, $k=|{\vec{k}}|$, and $\Theta$ being the Heaviside step function, the formulation of the excursion-set theory is the simplest.

In Eq. (3.1), the density field is independent of time. We will denote its variance as $S$ in the following (while $\sigma_{R}^{2}$ is reserved to the time-dependent variance). This is because in the framework of excursion set theory, to determine the collapsed fraction, one determines the evolution of the density field in terms of the variance at fixed time. The formation history of the collapsed objects, can be studied by retrieving the time-dependance via promoting the collapse threshold to be time-dependant, while keeping the variance fixed in time. This is due to the fact that the probability of forming virialized objects is determined in terms of the ratio $\delta_{c}/\sigma_{R}$, with $\sigma_{R}$ being the square root of the variance of the smoothed density field. We will further illustrate on this for the particular case of PBHs at the end of this section.

One proceeds by studying the “evolution” of $\delta({\vec{x}},R)$ with $R$, at a fixed value of ${\vec{x}}$, which is given by

[TABLE]

where

[TABLE]

Without loss of generality, we set $\vec{x}=0$, and to keep the notation compact, denote the smoothed field at $\vec{x}=0$ as $\delta_{R}$. Since the modes $\widetilde{\delta}({\vec{k}})$ are stochastic variables, $\zeta(R)$ is also a stochastic variable, and Eq. (3.4) has the form of a Langevin equation, with $R$ playing the role of time, and $\zeta(R)$ playing the role of noise. Statistical properties of a Gaussian distributed $\delta_{R}$, and hence $\zeta(R)$, are fully determined in terms of their connected two-point correlator, which for $\zeta(R)$ is given by

[TABLE]

where $\Delta^{2}(k)=k^{3}P_{\delta}(k)/(2\pi^{2})$. For a generic filter, the right-hand side is a function of $R_{1}$ and $R_{2}$, different from a Dirac-delta $\delta_{D}(R_{1}-R_{2})$.

For the particular case of sharp $k$-space filter, things simplify considerably, and one retrieves the Langevin equation with Dirac-delta noise. This can be clearly seen using $k_{f}=1/R$ instead of $R$, and defining $Q(k_{f})=-(1/k_{f})\zeta(k_{f})$ to get

[TABLE]

with

[TABLE]

The equations can be further simplified using the variance $S$ of the overdensity as a “pseudotime” variable,

[TABLE]

which for a sharp $k$-space filter becomes

[TABLE]

Finally, redefining $\eta(k_{F})=Q(k_{F})/\Delta^{2}(k_{F})$, Eq. (3.7) can be re-written as

[TABLE]

with

[TABLE]

Therefore, the evolution of the density field in terms of the variance (playing the role of time) obeys the Langevin equation with a Gaussian white noise. This implies that $\delta(S)$ performs a Brownian random walk, with respect to the “time” variable $S$, with no memory from one step to the other555If the window function is not a step function in momentum space, and is, for instance, a step function in real space (which gives rise to the volume-averaged field in Eq. (1)), the random walk of the smoothed overdensity is non-markovian, introducing memory effects from one step to the other (the white-noise case is markovian and no memory effects are present). This case has been studied in Refs. [26, 27, 28] using the path-integral methods, but is much more involved than the case of a white-noise. For the sake of simplicity, in this paper we restrict ourselves to the markovian case originally studied in Ref. [25].. Following Ref. [25], we refer to the evolution of $\delta$ as a function of $S$ as a “trajectory” (see a schematic picture in figure 2).

It is well-known that for a stochastic variable evolving according to Langevin equation with white noise, the distribution obeys the Fokker-Planck equation. Therefore, the distribution of $\delta_{R}$ can be determined from

[TABLE]

The Press-Schecter approach is equivalent to solving this equation, setting the boundary condition that the distribution vanishes at $\delta\rightarrow\pm\infty$, and accounting for all the trajectories that exceed the threshold of collapse $\delta_{c}$. This leads to the so-called “cloud-in-cloud” problem, which is associated with trajectories that make multiple crossings of the threshold. This problem is cured in excursion set theory approach by accounting for the first up-crossing of the trajectories. This is equivalent to calculating the lowest value of $S$ (or, equivalently, the highest value of $R$) for which the trajectory pierces the threshold. Once a trajectory passes the threshold, a virialized object forms, and this trajectory should be excluded from further consideration. Therefore, the solution of the cloud-in-cloud problem is implemented by imposing the boundary condition

[TABLE]

for which the solution of the FP equation is given by

[TABLE]

The fraction $F(S)$ of trajectories that have crossed the threshold at “time” smaller or equal to $S$ is then calculated from the fraction of trajectories that at “time” $S$ have never crossed the threshold, so

[TABLE]

which is twice the Press-Schecter prediction.

Excursion set theory can also be used to study the formation history of virialized objects. For dark matter halos, this was briefly discussed in the classic paper of Bond et al., [25] and further elaborated in Ref. [29]. As we briefly discussed earlier, the time evolution of the density field is conveniently recovered by absorbing the time-dependance into the collapse threshold, and defining a time-dependant threshold $\omega(a)\equiv\delta_{c}/D(a)$, with $D(a)$ being the linear growth factor. Therefore, the collapse is determined in terms of the linear density field evolved to the time of interest. Since for PBHs, the threshold scales like $\omega(a)\equiv\delta_{c}/a^{2}$, the barriers at two different time are related to one another as

[TABLE]

In other words, the barrier is lowered as time passes, while the variance is constant. In this way, the real time evolution is given by the time evolution of the barrier, which for PBHs is in one-to-one correspondence with the increase of the comoving Hubble length.

A great deal of insight into the formation history can be gained by considering the “two-barrier” problem [29], which allows us not only to calculate the probability to form a virialized object at a given time, but also the conditional probability that a trajectory with a first up-crossing of the barrier $\omega_{1}$ at some $S_{1}$ will have a first upcrossing of $\omega_{2}$ between $S_{2}$ and $S_{2}+{\rm d}S_{2}$ with $S_{1}\gg S_{2}$ and $\omega_{1}>\omega_{2}$ (see figure 3 for a schematic picture). This conditional probability is given by

[TABLE]

For PBHs, analogous to halos, the trajectories provide the variation of $\delta_{R}$ at fixed time. However, the condition for forming a PBH as well as its mass apply when a given scale crosses the horizon. One must generate the trajectories at a given fixed time and then evolve each value of $\delta_{R}$ backwards or forwards in time till a given scale of interest crosses the horizon [30]. Ultimately one is interested in the mass $M_{H}$ associated with a given scale when the latter enters the horizon, $R\sim 1/aH\sim a\sim M_{H}^{1/2}$ (in the radiation-dominated case), so that, still at fixed time

[TABLE]

Eq. (3.18) then gives the conditional probability that a PBH of mass $M_{1}$ present at the time $t_{1}$ has merged into a PBH of mass between $M_{2}$ and $M_{2}+{\rm d}M_{2}$ at time $t_{2}>t_{1}$. This quantity therefore seems the most suitable one to discuss the fate of the small PBHs in the broad spectrum case.

3.2 Broad spectra with a blue tilt

Having reviewed basics of the excursion set theory, let us get back to the questions of what are the relevant scales for the formation of PBHs for a broad primordial spectrum. We advocate here that for realistic power spectra with a blue tilt

[TABLE]

it is always the smallest scale (the largest momentum) that is relevant. In other words, if a broad spectrum is characterized by a range of scales running from $1/k_{l}$ to $1/k_{s}$, production of PBHs will be dominated by PBHs with mass $M(r_{m,s})$. To show this we can resort to two alternative techniques, one based on the excursion set that we discuss in this section, and the other on the so-called supreme statistics to which we devote the next section.

Let us consider the fluctuations on small and large scales which collapse to form PBHs at the time when the scale factors are $a_{s}$ and $a_{l}$ correspondingly. From what we discussed previously, the critical threshold for formation, $\delta_{c}$, is the same for all the PBHs. This means that the rescaled thresholds in Eq. (3.17) are hierarchical. Neglecting the radiation transfer function, the variances (at fixed time) as a function of the scale are also hierarchical

[TABLE]

From the excursion point of view, the separation of scales $k_{s}\gg k_{l}$ with a blue spectrum is therefore given by the hierarchies

[TABLE]

and

[TABLE]

Notice that the last condition implies that small PBHs are generated with larger abundances than the big ones, since the abundances are given by

[TABLE]

and hence,

[TABLE]

Taking the above limit in Eq. (3.18) we have

[TABLE]

Therefore, when larger modes re-enter the horizon, PBHs of higher masses form with a negligible abundance. Moreover, even if they do form, they do not engulf the previously formed lighter PBHs. In other words, the cloud-in-cloud phenomenon does not exist in practice for broad spectra with a blue tilt.

The final spectrum of masses is dominated by the small PBHs. One could argue that when longer modes enter the horizon, the PBHs could also form if $\sigma_{l}$ is sufficiently large. However, since at each horizon crossing time $\sigma_{s}\gg\sigma_{l}$, this implies that the variance of the small-scales would be dangerously in the non-perturbative regime. In this sense broad spectra with a blue tilt give predictions similar to the ones of narrow spectra.

3.3 Flat broad spectra

Next, let us consider a broad spectrum which is flat,

[TABLE]

In such a case, the hierarchies for the thresholds and the variances is maintained

[TABLE]

However (see Eq. (3.23) with $n_{p}=0$)

[TABLE]

This condition implies that small PBHs are generated with the same abundance as large ones.

The probability that a small PBHs have been incorporated into a new halo of a larger mass at $S(<S_{s})$ at later times is given by [29] (we approximate the lower limit of integration to zero, $S_{l}\simeq 0$)

[TABLE]

Assuming that $S\ll S_{s}$, $\omega(a)\ll\omega(a_{s})$, but $(S_{s}/S)=(\omega(a_{s})/\omega(a))^{2}$, it is easy to see that $X\simeq Y\simeq\omega(a)/\sqrt{2S}$ the above probability is again suppressed,

[TABLE]

where we used the fact that $\omega(a_{s})/\sqrt{S_{s}}={\cal O}(6\div 8)$, and the argument in the exponential enhancement factor is teamed by the small coefficient $\omega(a)/\omega(a_{s})\ll 1$.

This confirms the intuitive picture that, since formation of a PBH is a rather rare event, having a small PBH swallowed by a larger one is even rarer and it goes like $\beta(M(r_{m}))\beta(M(r_{m,s}))\simeq\beta^{2}(M(r_{m,s}))$, where in the last passage we have used the fact that for a broad spectrum the abundance is the same for small and big PBHs (since $\omega^{2}(a)/S=\omega^{2}(a_{s})/S_{s})$).

Having established the absence of the cloud-in-cloud problem for broad flat spectra, and since the abundance $\beta$ is the same for all masses, a natural question is then what is the final mass distribution. To answer this question we recall that

[TABLE]

Taking into account that the density of small-mass PBHs grows like non-relativistic matter, and that $\beta=\rho_{\text{\tiny PBH}}/\rho_{\rm tot}$ at the time of the PBH formation, the abundance of the small PBHs at the formation time $t_{l}$ of the large ones is

[TABLE]

Therefore, even though the probability of forming a PBH is independent of the mass, i.e. $\beta(M(r_{m,s}))=\beta(M(r_{m,l}))$, the PBHs with the smallest mass will give the largest contribution to the mass distribution. In this sense, again, the predictions are similar to the ones of a narrow spectrum.

Of course if the hierarchy between the smallest and largest momenta is not that large, the big PBHs could have a non-negligible abundance and play an interesting cosmological role.

3.4 Broad spectra with a red tilt

If the spectrum is broad, but with a red tilt, $n_{p}<0$, from the expression (3.2) one deduces that large-scale variances computed at their horizon crossing time $t_{l}$ are larger than those at small scales when they entered the horizon. The formation of PBHs will be therefore dominated by the largest scales, and PBHs will have the largest possible mass allowed by the power spectrum. In such a case, the presence of the small PBHs is highly suppressed, and the cloud-in-cloud problem is absent. In this sense, broad spectra with a red tilt also give predictions similar to the peaked power spectra.

One could argue that if the broad spectrum is red or flat, the production of the lighter PBHs could be significantly enhanced at later times by the impact of the long modes perturbations on the small ones after they have entered the horizon [31]. However, long modes cross the horizon when the short modes are well inside the horizon. By that time, the radiation pressure has already dispersed the small scales peaks (see, for example, Ref. [32] for a discussion of the dynamics of the sub-threshold peaks upon horizon re-entry). This effect is captured already at the linear level by the transfer function in Eq. (1.7) in a radiation-dominated universe, which scales like $1/(k\eta)^{2}\ll 1$ on subhorizon scales. Therefore short-scale variances promptly decay in contrast to what was suggested in Ref. [31]. The time evolution of each mode prevents different scales from sizeably impacting each other.

Another option is that PBHs are also generated at small-scales. This will require though that $\sigma_{l}\gg\sigma_{s}$ and the danger is to enter a non-linear regime before the long mode enter the horizon. We will come back to these considerations in the section devoted to the clustering.

4 The supreme statistics

In this section we return to the case of a broad spectrum as in Eq. (2.2) with a blue tilt $(n_{p}>0)$ and we wish to show with a different technique, the supreme statistics, that small PBHs are the only ones which are relevant. The reader should be aware of the fact that in this section we return to work with the physical quantities and not with the rescaled ones.

Consider a spherical peak (we consider for simplicity only spherical peaks, which is a good approximation for those rare peaks giving rise to PBHs) whose typical scale is $r_{m,l}\gg r_{m,s}$. This peak re-enters the horizon when $a_{l}H_{l}r_{m,l}\sim 1$, and its volume-averaged density contrast is

[TABLE]

Inside this large spherical peak of size $r_{m,l}$ there are about

[TABLE]

subpeaks with volume of $4\pi r_{m,s}^{3}/3$ and typical scale of $r_{m,s}\sim 1/k_{s}$. The volume-averaged density contrast characterising each subpeak is

[TABLE]

These subpeaks re-enter the horizon when $a_{s}H_{s}r_{m,s}\sim 1$, with $a_{s}H_{s}\gg a_{l}H_{l}$, therefore they re-enter the horizon much earlier than the large peak does. The critical threshold for the formation of a PBHs is the same for the large peak and its subpeaks, while the variances are different at each horizon crossing $t_{s}$ and $t_{l}$ repsectively (we do not consider here the radiation transfer function)

[TABLE]

The question we wish to answer is whether the subpeaks form PBHs more easily than the large peaks. To answer this question we can use the so-called supreme statistics (or extreme-value statistics) [33] to study the statistics of the highest subpeaks, and determine the largest value of the subpeak density contrast at the short scale. It is this highest subpeak which will determine if a collapse into a PBH takes place on the short scale. In more general terms, considering $n$ samples, each of size $m$, all drawn from the same underlying distribution (which is Gaussian in our case), the distribution of the maxima within each sample, and therefore the most probable maximum value, can be determined using the supreme statistics 666Of course in the limit $n\to\infty$ one recovers the Gaussin distribution and the largest value goes to infinity..

In our case we have $m=1$ and $n=N$. Specializing to a Gaussian parent distribution, the probability that the density in any one subpeak is less than a given value $\widehat{\delta}_{s}$ is given by the cumulative distribution function,

[TABLE]

where $P(\delta_{s}|\delta_{l})$ is the conditional probability for the density at short-scales $\delta_{s}$, given the density at large scales $\delta_{l}$. Since both $\delta_{s}$ and $\delta_{l}$ are Gaussian distributed, their conditional probability is also a Gaussian. For large values of $\delta_{s}$ and $\sigma_{s}$, this expression can be safely approximated to

[TABLE]

and hence

[TABLE]

The probability that the density of all $N$ subpeaks is less than $\widehat{\delta}_{s}$ is given by

[TABLE]

and the differential probability that a subpeak of density $\widehat{\delta}_{s}$ at a scale $r_{m,s}$ has the highest density among all the subpeaks is given by ${\rm d}P_{N}/{\rm d}\widehat{\delta}_{s}$. From this probability distribution, we can compute the mean (and variance) of the full set of subpeaks. In figure 4, we show a schematic picture of the subpeaks within a large spherical peak.

In the limit of large $\widehat{\delta}_{s}$, one can obtain a simple analytic expression for the most probable value of $\widehat{\delta}_{s}$, as was done in Ref. [34]. The maximum of ${\rm d}P_{N}/{\rm d}\widehat{\delta}_{s}$ determines the most likely value of $\widehat{\delta}_{s}$,

[TABLE]

Since by definition ${\rm d}P_{1}/{\rm d}\widehat{\delta}_{s}=P(\widehat{\delta}_{s}|\delta_{l})$, and $P_{1}$ is nearly one in the large-$\widehat{\delta}_{s}$ limit, Eq. (4.9) reduces to

[TABLE]

which using Eq. (4.7) gives

[TABLE]

Therefore, the exponential suppression is offset by the large number of subregions that can attain independent values. Since $N\simeq(r_{m,l}/r_{m,s})^{3}\gg 1$, to a good approximation we finally find [34]

[TABLE]

which, as expected, for $N\to\infty$ goes to infinity. Note that the density on small-scales does not depend on the boundary conditions, i.e. the large-scale density. Physically, this can be understood as that the very small scales fluctuations are so large that they are virtually independent of the large-scale density field.

In fact, a more refined computation for the limit of large $N$, shows that the probability distribution of the largest value $\widehat{\delta}_{s}$ is given by [36, 33, 35]

[TABLE]

A comparison with a numerical simulation of several samples of $N=10^{4}$ is shown in figure 5.

The key point is that, even in the presence of a broad peak with size $r_{m,l}$, there is on average one subpeak at a given scale $r_{m,s}$ whose value of the overdensity is given by the maximal value $\widehat{\delta}_{s}\gg\delta_{l}$. It is therefore at the smallest scale $r_{m,s}$ that the PBH forms first once the comoving horizon grows up to a value $\sim r_{m,s}$. In other words, at the re-entry of the small scale modes into the comoving Hubble horizon, if the corresponding smoothed overdensity is large enough, the PBHs will form.

One would conclude that, if $\widehat{\delta}_{s}$ is above (and close to) the threshold, one expects to form roughly one PBH in each volume $V_{s}$ (with $\widehat{\delta}_{s}$ being the maximum which is realized only in one out of the $N$ samples), that is

[TABLE]

The mass fraction composed of the PBHs with mass $M(r_{m,s})$ is therefore $(x=\widehat{\delta}_{s}/\sigma_{s}$)

[TABLE]

regardless of the value of $N$. Here the factor $f_{l}$ accounts for the number density of the larger peaks.

The result in Eq. (4) summarizes our point. Even if the initial power spectrum of the curvature perturbation is broad, if it grows with the momentum, the population of PBHs is dominated by those which form first, i.e. from the smallest scale fluctuations. Furthermore, the cloud-in-cloud probability for these small PBHs to be swallowed by larger PBHs is very tiny. In our final expression, this results in the further suppression factor $f_{l}$, which is extremely small, given that the typical value of $\delta_{l}$ is much smaller than $\delta_{c}$. Note that summing over all peaks of size $r_{m,l}$, even those not giving rise to PBHs, the factor $f_{l}$ tends to unity and one recovers the standard abundance of PBHs with mass $M(r_{m,s})$.

To summarize, our expectation is that when larger modes re-enter the horizon, the PBHs of higher masses do not start forming, neither do they engulf the previously formed lighter PBHs. This result reenforces the conclusion in the previous section, which was based on the excursion set method applied to PBHs.

5 The clustering at formation

Another relevant issue arising when PBHs are generated by broad spectra is whether or not they are clustered at the formation time. Initial clustering can affect PBH merging and accretion, leading to different spatial distributions at later times. This in turn can weaken (or strengthen) the current observational constraints on their abundances.

Being discrete objects, the two-point correlation function of PBHs have the general form

[TABLE]

where $\delta_{D}({\vec{r}})$ is the three-dimensional Dirac distribution, $\bar{n}_{\text{\tiny PBH}}$ is the average number density of PBH per comoving volume, and $\xi_{\text{\tiny PBH}}(r)$ is the reduced PBH correlation function. In Ref. [37] it was shown that for narrow spectra with amplitude ${\cal P}_{0}$, in the allowed range of PBH masses, the comoving scale $r_{\xi}$ at which $\xi_{\text{\tiny PBH}}(r)\simeq 1$ is always smaller than the average comoving separation between two PBHs, $\overline{r}=(3/4\pi\bar{n}_{\text{\tiny PBH}})^{1/3}$

[TABLE]

For instance, for the case of PBHs with $M_{\text{\tiny PBH}}\sim 30\ M_{\odot}$, the initial distances relevant for the present merger rate is \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt\sim}\hss}\raise 1.0pt\hbox{>}}4\cdot 10^{-5} Mpc [31], and it was concluded that for narrow spectra clustering is not relevant.

In the previous sections we have argued that, since PBHs are rare events, the cloud-in-cloud phenomenon is basically absent. Furthermore, the fact that the final spectrum of PBHs is tilted either towards small-scales PBHs (for flat broad spectra or blue-tilted broad spectra) or large-scale PBHs (for red-tilted broad spectra) the expectation as far as the mass function is concerned is similar to that of a narrow power spectrum.

These two considerations suggest that the clustering properties at formation will not differ from the case of narrow spectra. In other words, the PBHs will not be clustered. As already mentioned, we do not expect this conclusion to be altered by the possibility, in the case of flat or red-tilted broad spectra, that the formation of small-scale PBHs is enhanced at later times by the effect of the long modes perturbations once the latter have entered the horizon. Indeed, by that time the small-scale peaks with an amplitude below the threshold (at the time when they enter the horizon) have been already washed-out by the radiation pressure.

We use the supreme statistics again to elaborate further on the clustering for the case of a broad spectrum with a blue tilt. Having argued that for broad spectra with a blue tilt PBHs form at the smallest scale without being swallowed by the more massive PBHs, one can ask if nevertheless the correlation length of the lighter PBHs is larger than the Poisson distance. At the time of formation the average comoving separation between such PBHs is

[TABLE]

provided that only one subpeak within the spherical peak of radius $r_{m,L}$ collapses to form a PBH (this is the assumption made in the previous section). However, if more than one subpeak collapses into a PBH the correlation length will be given roughly by $r_{m,l}$. For instance, one can compute the probability distribution function of the $p$-th maximum out of a sample of size $N$, which is given by [38]

[TABLE]

where $f(x)=\exp(-x^{2}/2)/\sqrt{2\pi}$ is the Gaussian probability distribution function. It is easy to show that, in the case of the first maximum $p=1$, the previous distribution recovers the one in Eq. (4.13) in the large $N$ limit.

In order to deal with more tractable analytical formulae one can also approximate the probability $P_{p}(x)$ with the following form (we have checked that for large $N$ the results coincide) [33]

[TABLE]

where $x_{p}$ is the mode of the probability, where the latter has its maximum. For $p=1$ one recovers Eq. (4.13). As for the coefficients, we have found them numerically, see Fig. 6, and find, for instance, that $a_{1}\simeq a_{2}\simeq\log_{10}N$.

In the limit of large $N$ and small $p$, the distance between different ranked maxima, defined as

[TABLE]

is distributed as

[TABLE]

Therefore, knowing the expectation value for the first ranked maximum, we can infer the expectation value for the second and so on. In particular, the most probable $d_{m}$ is

[TABLE]

For the second maximum we therefore get

[TABLE]

Using the above distribution we can compute the probability to have two PBHs in a peak of size $r_{m,l}$. This is equivalent to computing the probability that the second maximum is above or equal the threshold $\delta_{c}$ as

[TABLE]

The result is shown in Fig. 7. We stress the dependence of $P(\widehat{\delta}_{s,2}\geq\delta_{c})$ on the variance $\sigma_{s}$ through the lower integration limit.

Suppose now that the broad spectrum extends over a range of scales such that $r_{m,l}\simeq N^{1/3}r_{m,s}\simeq N^{1/3}/a_{s}H_{s}$. In order for the distance of these two subpeaks to be smaller than the Poissonian one, the following condition must be imposed

[TABLE]

where in the last passage we have assumed that the PBHs form the dark matter. On the other hand from Fig. 7 we learn that the probability of having two subpeaks above the threshold is negligibly small unless N\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt\sim}\hss}\raise 1.0pt\hbox{>}}\beta^{-1}. This can be understood by realizing that the two maxima are very close to each other and therefore

[TABLE]

This indicates that the correlation length is never smaller than the average Poisson distance for a broad spectrum with a blue tilt, suggesting that PBHs are not clustered at formation.

6 Conclusions

In this paper we have examined the formation of PBHs that are generated by sizeable small-scale curvature perturbations originated during inflation, upon horizon re-entry. We have particularly focused on the case of broad spectra for the curvature perturbation in which one may not immediately identify a single mass-scale to associate with the PBHs. By using the excursion set method we have argued that the probability of the small-mass PBHs being absorbed by the bigger ones is negligible. This is due to the fact that the formation of PBHs is by itself a rare event. This result has been confirmed by using the supreme statistics which describes the distribution of extreme values.

Our findings indicate that the predictions from broad spectra are qualitatively the same as for narrow spectra. The mass distribution of the PBHs is dominated by a single mass (which can be associated to the maximum or the minimum of the momenta of the power spectrum depending on its tilt). In particular, we have offered indications that the PBHs at formation are not clustered. Of course the subsequent evolution of the clustering will depend upon the merging and accretion of the PBHs, and might be fully captured only by real-time numerical simulations. Another aspect worth investigating more in detail is the role of non-Gaussianities in our considerations. For instance, the presence of local-type primordial non-Gaussianity might enhance the clustering on large scales [14, 39] due to the correlation of small and large scales.

Acknowledgments

We thank V. Dejaques, I. Musco and S. Young for many useful discussions. The authors are supported by the Swiss National Science Foundation (SNSF), project The Non-Gaussian Universe and Cosmological Symmetries, project number: 200020-178787.