The Ineludible non-Gaussianity of the Primordial Black Hole Abundance

TL;DR

This paper demonstrates that the formation of primordial black holes inherently involves non-Gaussian statistics due to nonlinear relations, and provides methods to accurately calculate their abundance considering this non-Gaussianity.

Contribution

It introduces two techniques based on peak theory and threshold statistics to compute non-Gaussian PBH abundance, highlighting the importance of non-Gaussian effects.

Findings

Non-Gaussianity is unavoidable in PBH formation due to nonlinear density relations.

Standard perturbative approaches may underestimate PBH abundance.

Gaussian approximations can be recovered by increasing the power spectrum amplitude by a factor of 2-3.

Abstract

We study the formation of primordial black holes when they are generated by the collapse of large overdensities in the early universe. Since the density contrast is related to the comoving curvature perturbation by a nonlinear relation, the overdensity statistics is unavoidably non-Gaussian. We show that the abundance of primordial black holes at formation may not be captured by a perturbative approach which retains the first few cumulants of the non-Gaussian probability distribution. We provide two techniques to calculate the non-Gaussian abundance of primordial black holes at formation, one based on peak theory and the other on threshold statistics. Our results show that the unavoidable non-Gaussian nature of the inhomogeneities in the energy density makes it harder to generate PBHs. We provide simple (semi-)analytical expressions to calculate the non-Gaussian abundances of the primordial black holes and show that for both narrow and broad power spectra the gaussian case from threshold statistics is reproduced by increasing the amplitude of the power spectrum by a factor ${\cal O}(2\div 3)$.

Figures

Equations

$${\sigma }_{\delta }^2=\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}{W}^2(k,R_H)P_{\delta }(k),$$

$$P_{\text{G}}(\delta >{\delta }_{\mathrm{c}})=\beta (M)={\int }_{{\delta }_{\mathrm{c}}}\frac{\mathrm{d}\delta }{\sqrt{2\pi }{\sigma }_{\delta }}{e}^{-{\delta }^2/2{\sigma }_{\delta }^2}.$$

$${\nu }_{\mathrm{c}}=\frac{{\delta }_{\mathrm{c}}}{{\sigma }_{\delta }},$$

$${\beta }_{\text{G}}^{\mathrm{th}}\simeq \sqrt{\frac{1}{2\pi {\nu }_{\mathrm{c}}^2}}{e}^{-{\nu }_{\mathrm{c}}^2/2}.$$

$${\beta }_{\text{G}}^{\text{pk}}\simeq \frac{1}{3\pi }{\left(\frac{⟨k^2⟩}{3}\right)}^{3/2}{R}_H^3({\nu }_{\text{pk}}^2-1)e^{-{\nu }_{\text{pk}}^2/2}\mathrm{with}⟨k^2⟩=\frac{1}{{\sigma }_{\delta }^2}\int \frac{{\mathrm{d}}^{3}k}{(2\pi {)}^3}{k}^2{P}_{\delta }(k),$$

$${\nu }_{\text{pk}}=\frac{{\delta }_{\text{pk}}^{\mathrm{c}}}{{\sigma }_{\delta }},$$

$$\delta (\vec{x},t)=-\frac{8}{9a^2{H}^2}{e}^{-5\zeta (\vec{x})/2}{\nabla }^2{e}^{\zeta (\vec{x})/2}=-\frac{4}{9}\frac{1}{a^2{H}^2}{e}^{-2\zeta (\vec{x})}\left({\nabla }^2\zeta (\vec{x})+\frac{1}{2}{\partial }_i\zeta (\vec{x}){\partial }^i\zeta (\vec{x})\right).$$

$$\delta (\vec{x},t)=-\frac{4}{9a^2{H}^2}{\nabla }^2\zeta (\vec{x})$$

$$P_{\delta }(k,t)=\frac{16}{81}\frac{k^4}{a^4{H}^4}{P}_{\zeta }(k).$$

$$\zeta \simeq \frac{9a^2{H}^2}{4k_{\star }^2}\delta \simeq \frac{9a^2{H}^2}{4k_{\star }^2}{\delta }_{\mathrm{c}}\simeq 0.15,$$

$$\displaystyle P(\delta>\delta_{\rm c})=\Big{<}\Theta(\delta-\nu_{\rm c}\sigma_{\delta})\Big{>}=\int[D\delta(\vec{x})]P[\delta(\vec{x})]\Theta\Big{(}\delta(\vec{x})-\nu_{\rm c}\sigma_{\delta}\Big{)},$$

$$\displaystyle\Big{<}\delta(\vec{x}_{1})\cdots\delta(\vec{x}_{n})\Big{>}_{c}=\xi_{n}(\vec{x}_{1},\cdots,\vec{x}_{n}),$$

$$P(\delta >{\delta }_{\mathrm{c}})=\beta (M)=\frac{1}{\sqrt{2\pi {\nu }_{\mathrm{c}}^2}}\exp \left\{-{\nu }_{\mathrm{c}}^2/2+\sum _{n=3}^{\infty }\frac{(-1{)}^n}{n!}{\xi }_n(0)({\delta }_{\mathrm{c}}/{\sigma }_{\delta }^2{)}^n\right\},$$

$$S_n=\frac{{\xi }_n(0)}{{\left({\xi }_2(0)\right)}^{n-1}}=\frac{{\stackrel{n-\mathrm{times}}{⟨\overbrace{\delta (\vec{x})\cdots \delta (\vec{x})}⟩}}_c}{{\sigma }_{\delta }^{2(n-1)}}.$$

$${\Delta }_n=\frac{\mathrm{d}\ln \beta (M)}{\mathrm{d}\ln S_n}.$$

$$\frac{{\beta }_{\text{NG}}^{\text{th}}(M)}{{\beta }_{\text{G}}^{\text{th}}(M)}=e^{{\Delta }_n}.$$

$$|{\Delta }_n|to0.0pt\text{lower}4.0pt\sim \text{hss}\text{raise}1.0pt<1.$$

$$|{\Delta }_n|=\frac{1}{n!}{\left(\frac{{\delta }_{\mathrm{c}}}{{\sigma }_{\delta }}\right)}^2|S_n|{\delta }_{\mathrm{c}}^{n-2}.$$

$$\left|S_n\right|to0.0pt\text{lower}4.0pt\sim \text{hss}\text{raise}1.0pt<{\left(\frac{{\sigma }_{\delta }}{{\delta }_{\mathrm{c}}}\right)}^2\frac{n!}{{\delta }_{\mathrm{c}}^{n-2}}.$$

$$P_{\zeta }(k)=\frac{2{\pi }^2}{k^3}{\mathcal{P}}_{\zeta }(k)\mathrm{and}{\mathcal{P}}_{\zeta }(k)=A_s{k}_{\star }{\delta }_D(k-k_{\star }).$$

$$⟨{\delta }^2{⟩}_c$$

$$⟨{\delta }^3{⟩}_c$$

$$⟨{\delta }^4{⟩}_c$$

$$c_{\star }{k}_{\star }^2$$

$$A_{s}to0.0pt\text{lower}4.0pt\sim \text{hss}\text{raise}1.0pt>6.0(4.0)\cdot 10^{-3},$$

$$\beta to0.0pt\text{lower}4.0pt\sim \text{hss}\text{raise}1.0pt<1.3\times 10^{-9}{\left(\frac{M}{M_{\odot }}\right)}^{1/2}.$$

$${\mathcal{P}}_{\zeta }(k)=\frac{A_g}{\sqrt{2\pi }\sigma }\exp \left[-\frac{{\ln }^2\left(k/k_{\star }\right)}{2{\sigma }^2}\right].$$

$$⟨{\delta }^2{⟩}_c$$

$$⟨{\delta }^3{⟩}_c$$

$$⟨{\delta }^4{⟩}_c$$

$$A_{g}to0.0pt\text{lower}4.0pt\sim \text{hss}\text{raise}1.0pt>3.8(2.2)\cdot 10^{-3},$$

$$\delta (\vec{x},t)=-\frac{4}{9a^2{H}^2}{e}^{-2\zeta (\vec{x})}\left[{\nabla }^2\zeta (\vec{x})+\frac{1}{2}{\partial }_i\zeta (\vec{x}){\partial }^i\zeta (\vec{x})\right].$$

$$\zeta (\vec{x})=\zeta ({\vec{x}}_{\text{pk}})+{\partial }_i\zeta ({\vec{x}}_{\text{pk}})(x^i-x_{\text{pk}}^i)+\frac{1}{2}{\partial }_i{\partial }_j\zeta ({\vec{x}}_{\text{pk}})(x^i-x_{\text{pk}}^i)(x^j-x_{\text{pk}}^j).$$

Full text

1 Introduction

Since the first detection of gravitational waves originated by the merging of two $\sim 30M_{\odot}$ black holes [1], the idea that Primordial Black Holes (PBHs) might form a considerable fraction of the dark matter [2, 3, 4] has attracted again much interest [5] (see Ref. [6] for a recent review). A popular mechanism for the formation of PBHs is the scenario in which PBHs are originated from the enhancement of the curvature power spectrum at a given short length scale due to some features [6]. If the power spectrum of the curvature perturbation is enhanced during inflation to values $\sim 10^{-2}$ on small scales and subsequently transferred to radiation during the reheating process, PBHs may form from sizeable fluctuations if the latter overcome the counter effect of the radiation pressure.

Since the perturbation of fixed comoving size does not collapse till it re-enters the cosmological horizon, the size of a PBH at formation is related to the horizon length and its mass $M$ is approximately the mass contained in such a horizon volume. Fluctuations collapse immediately after horizon re-entry to form PBHs if they are sizeable enough. We indicate by $\delta$ the overdensity and by $\sigma_{\delta}^{2}$ its variance

[TABLE]

where $P_{\delta}$ is the overdensity power spectrum, $R_{H}$ being the comoving horizon length $R_{H}=1/aH$, $H$ is the Hubble rate and $a$ the scale factor. The quantity $W(k,R_{H})$ is a window function, for which we choose a top-hat in real space. Under the assumption that the density contrast is a linear quantity obeying gaussian statistics, threshold statistics (or Press-Schechter) predicts that the primordial mass fraction $\beta(M)$ of the universe stored into PBHs at the formation time is given by111In the literature sometimes this expression may be multiplied by a factor of 2 to account for the cloud-in-cloud problem [9]. There seems to be no agreement if this factor should be included for PBHs. Numerically it makes little difference.

[TABLE]

Here $\delta_{\rm c}$ is the threshold for formation of the PBHs which quantifies how large the overdensity perturbations must be and depends on the shape of the power spectrum [7, 8, 10]. By defining

[TABLE]

the Gaussian mass fraction can be well approximated by (\nu_{\rm c}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt\sim}\hss}\raise 1.0pt\hbox{>}}5)

[TABLE]

This expression for the PBH mass fraction comes about when identifying the PBHs with regions whose overdensity is above a given threshold, hence the name of threshold statistics.

Alternatively, one can identify the PBHs with the local maxima of the overdensity, and one may use peak theory [11] to compute their mass fraction. In such a case one has [12]222We differ slightly from the corresponding expression in Ref. [12]. First by a factor of 3 to account for the fact that one counts the number density of peaks at superhorizon scales, but the PBHs formed once the overdensity crosses the horizon at a slightly later time [8] (see also section 3). Secondly, by the fact that we define the mass going into PBH to be $M=(4\pi/3)\overline{\rho}R_{H}^{3}$, where $\overline{\rho}$ is the background radiation density. More importantly, we use here the definition (1.6) for the critical value $\nu_{\text{\tiny pk}}$. We will give more details in section 3. At the gaussian level, peak theory gives a PBH abundance which is systematically larger than the one provided by the threshold statistics [12].

[TABLE]

where now [8]

[TABLE]

and $\delta_{\text{\tiny pk}}^{\rm c}$ is to be identified with the critical value of the overdensity at the center of the peak above which an initial perturbation eventually collapses into a PBH [8, 10]. Notice that here we follow Refs. [8, 7] and do not introduce a window function for the peak theory. Indeed, for the examples we will discuss the window function is not strictly necessary because they are characterised by a well-defined scale in momentum space and the corresponding distribution is already smooth on length scales smaller than that characteristic scale. Also, in the case of peak theory a typical length pops out automatically, that is the scale $R_{*}$.

The gaussian expressions (1.4) and (1.5) make already manifest the essence of the problem we are going to discuss in this paper. PBHs are generated through very large, but rare fluctuations. Therefore, their mass fraction at formation is extremely sensitive to changes in the tail of the fluctuation distribution and therefore to any possible non-Gaussianity in the density contrast [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. This implies that non-Gaussianities need to be accounted for as they can alter the initial mass fraction of PBHs in a dramatic way. For instance, the presence of a primordial local non-Gaussianity in the comoving curvature perturbation can significantly alter the number density of PBHs through mode coupling [27, 28, 29, 30, 31, 32].

In this paper we will be dealing with a source of non-Gaussianity which is unavoidably generated by the non-linear relation among the overdensity $\delta(\vec{x},t)$ ($t$ is the cosmic time) and the comoving curvature perturbation $\zeta(\vec{x})$. It is important to stress that this non-linear relation makes the overdensity non-Gaussian even if the curvature perturbation is gaussian. In this sense, the non-Gaussianity we will discuss here is ineludible.

Let us briefly discuss where this non-linearity relation comes from. As we mentioned above, in the early radiation-dominated universe, the PBHs are generated when highly overdense regions gravitationally collapse directly into a black hole. Before collapse, the comoving sizes of such regions are larger than the horizon length and the separate universe approach can be applied [33]. One therefore expands at leading order in spatial gradients of the various observables, e.g. the overdensity. At this stage, the slicing and the threading of the spacetime manifold are to be fixed. For instance, the so-called comoving gauge seems appropriate as it has been adopted to perform numerical relativity simulations to describe the formation of PBHs and to calculate the threshold for PBH formation [7].

In the comoving slicing, the overdensity turns out to be [33]

[TABLE]

As the universe expands, the overdensity grows. Regions where it becomes of order unity eventually stop expanding and collapse. This happens when the comoving scale of such a region becomes of the order of the horizon scale. Even though the gradient expansion approximation breaks down, it has been used to obtain an acceptable criterion for the PBH formation (that is to compute the overdensity threshold) and this approximation has been confirmed to hold by nonlinear numerical studies [6, 34].

The standard procedure in the literature is to expand the relation (1.7) to first-order in $\zeta$

[TABLE]

and to relate the power spectrum of the overdensity to the one of the curvature perturbation by the relation

[TABLE]

The question is to what extent this is a good approximation given the fact that even tiny changes (percent level) in the square root of the overdensity variance are exponentially amplified in the PBH mass fraction.

To get the feelings of the numbers, let us roughly estimate the impact of the exponential $e^{-2\zeta(\vec{x})}$. Calling $k_{\star}$ the typical momentum of the perturbation, from Eq. (1.8) we get

[TABLE]

where we have taken the threshold $\delta_{\rm c}\simeq 0.5$ and $k_{\star}\simeq 2.7aH$ [8]. This gives $e^{-2\zeta(\vec{x})}\simeq 0.7$. This looks as a small change, but in fact it has an exponentially large effect in the mass fraction when the corresponding overdensity variance is calculated.

The goal of this paper is to deal with the intrinsically non-Gaussian nature of the overdensity onto the mass fraction of PBHs. First of all, we will provide a simple argument to convince the reader that the non-Gaussianity introduced by the non-linear relation (1.7) between the overdensity and the gaussian curvature perturbation has an impact on the PBH mass fraction which may not be accounted for by a perturbative approach. Based on this finding, we will proceed by computing the mass fraction taking into account such intrinsic non-Gaussianity. We will do so by using two methods.

Since PBHs may be thought to originate from peaks, that is, from maxima of the local overdensity, we will resort to peak theory [11] to calculate the probability of formation of the PBHs. This method is based on the fact that for high values of the overdensity at the peaks, their location can be confused with the location of the peaks in the comoving curvature perturbation as long as such peaks are sufficiently spiky, that is if their curvature (proportional to the second spatial derivatives) is large enough at the center of the peak [7].

Alternatively, we will use the non-Gaussian threshold statistics and provide an exact expression for the probability to form PBHs. Both methods indicate that the inevitable non-Gaussian nature of the overdensity makes more difficult to generate PBHs, independently from the shape of the power spectrum.

Let us also add a cautionary note. The intrinsic non-Gaussianity of the overdensity changes also the shape of the profile of the peaks which eventually give rise to PBHs upon collapse. Since the threshold depends on the shape of the overdensity, such non-Gaussianity influences as well the threshold value. This will be discussed in a separate publication [35].

The paper is organised as follows. In section 2 we offer a simple criterion to show that the intrinsic non-Gaussianity cannot be described by perturbative methods. Sections 3 and 4 will describe the two methods mentioned above. Section 5 contains our conclusions. The paper contains as well several appendices for the technical details.

2 A simple criterion to show that intrinsic non-Gaussianity matters

In order to establish if the intrinsic non-Gaussianity introduced by the non-linear relation (1.7) is relevant, we start from the non-Gaussian threshold statistics developed in Ref. [36] and refined in Ref. [24] by means of a path-integral approach. We do not report all the details here and the interested reader is refereed to those references for more details. We do not use here the window function which would introduce painful, but useless technicalities without changing the conclusions. Suffice to say that the probability of having the overdensity larger than a given threshold can be viewed as the one-point function of the threshold quantity

[TABLE]

where $\Theta(x)$ is the Heaviside function. By defining the connected correlators of the overdensity as

[TABLE]

one finds that, in the limit of large $\nu_{\rm c}$, the threshold statistics is given by [36, 24]

[TABLE]

where the label 0 means that the correlators are computed at equal points. To see under which circumstances the non-Gaussianity of the overdensity alters the predictions of the gaussian primordial abundance of PBHs in a significant way, we define dimensionless quantities, the cumulants, by the relations

[TABLE]

Following Ref. [24] we may define the fine-tuning $\Delta_{n}$ to be the response of the PBH abundance to the introduction of the $n$-th cumulant as

[TABLE]

Each cumulant allows to express the non-Gaussian PBH abundance in terms of the gaussian abundance as

[TABLE]

This implies that the PBH abundance is exponentially sensitive to the non-Gaussianity unless $\Delta_{n}$ is in absolute value smaller than unity

[TABLE]

Inspecting Eqs. (1.4) and (2.3), we see that

[TABLE]

This tells us that intrinsic non-Gaussianity in the overdensity alters exponentially the gaussian prediction for the PBH abundance unless

[TABLE]

To investigate how restrictive this condition is, we take the simplest case possible, i.e. a very narrow power spectrum for the comoving curvature perturbation which we approximate by a Dirac delta

[TABLE]

Here $A_{s}$ is the amplitude of the power spectrum and $k_{\star}$ is the characteristic scale of the power spectrum. Its relation with the cosmological horizon at formation $R_{H}$ has to be fixed running numerical simulations [8, 10]. For the case at hand, it is given by $k_{\star}\simeq 2.7/R_{H}$ (more comments on this later on). We do not report all the technical details here, which can be found in Appendix A, where we have consistently calculated the variance, the skewness $S_{3}$ and the kurtosis $S_{4}$ up to third-order in perturbation theory (in the power spectrum $P_{\zeta}$, that is up to $A_{s}^{3}$). We get

[TABLE]

One can check that the criterion (2.9) for the skewness (kurtosis) gives the lower bound

[TABLE]

where we have taken $\delta_{\rm c}=0.5$. We now impose the condition that the PBHs form at most the totality of dark matter, which provides an upper bound on their mass fraction given by

[TABLE]

For instance, for PBH masses around the interesting value of $10^{-12}M_{\odot}$ [38, 39], one would get from the gaussian mass fraction (1.4) $\beta\sim 10^{-15}$, $\nu_{\rm c}\simeq 8$ and therefore $A_{s}\simeq 3.7\cdot 10^{-4}$. This figure violates the bound required (2.12) to neglect the non-Gaussianity by one order of magnitude. More importantly, the kurtosis does not provide a bound which is much weaker than the skewness. This signals the breaking of the perturbative approach and calls for a more refined treatment.

The same conclusion can be obtained in the case where the power spectrum of the comoving curvature perturbation is parametrised by a log-normal shape of the form

[TABLE]

Using the results in Appendix A, one finds the following (for $\sigma=0.2$)

[TABLE]

The criterion (2.9) in this case results in a lower bound

[TABLE]

for the skewness and kurtosis respectively, while requiring again $\beta\sim 10^{-15}$ for $M\sim 10^{-12}M_{\odot}$ gives $A_{g}=2.5\cdot 10^{-4}$. Again we do not see signs of convergence in the perturbative approach.

3 The non-Gaussian probability from peak theory

Having shown that perturbation theory fails to provide the probability for PBH formation, we first resort to peak theory [11]. As we already mentioned in the introduction, PBHs trace the peaks of the radiation density field on superhorizon scales where the number of peaks per comoving volume is constant. Notice that we are dealing with peaks of the overdensity rather than the peaks of the curvature perturbation. This is because one cannot impose any constraint on the value of the gravitational potential (or curvature perturbation) on superhorizon scales because constant gravitational potentials cannot lead to any observable effect. Nevertheless, one can start from the following important point: large threshold peaks of the overdensity may be identified within a Hubble volume with the peaks of the curvature perturbation if the Laplacian of the curvature perturbation (that is the curvature of the peak) at the peak is large enough [7]. More in details, one can show that if the value of $\delta$ is comparable to the threshold value at a peak, one can find the associated peak of $\zeta$ well inside the horizon patch and centered at the peak of $\delta$ as long as the peaks in $\zeta$ is spiky enough. Let us elaborate about this point in the next subsection.

3.1 Spiky peaks of the curvature perturbation may be confused with peaks of the overdensity for large thresholds

The argument given in Ref. [7] is as follows. Let us consider the nonlinear expression (1.7) relating $\delta$ and $\zeta$ on superhorizon scales and in radiation domination

[TABLE]

We can expand the comoving curvature perturbation $\zeta(\vec{x})$ for points $\vec{x}$ around the peak position $\vec{x}_{\text{\tiny pk}}$ of the overdensity 333We indicate by $\partial_{i}\zeta(\vec{x}_{\text{\tiny pk}})$ the gradient $\partial_{i}\zeta(\vec{x})$ computed at $\vec{x}_{\text{\tiny pk}}$, and so on. $\delta(\vec{x},t)$

[TABLE]

Around such a peak we can also write

[TABLE]

where we neglected the second term in the square bracket since its contribution is of higher order in $\zeta$ with respect to (3.3).

Since the peak amplitude of the overdensity must be larger than some critical value $\delta^{\rm c}_{\text{\tiny pk}}$, we deduce that the curvature of the peak in $\zeta$ is bounded from above

[TABLE]

This is what we meant by saying that the peaks in $\zeta$ must be spiky enough. Now, the peak in $\zeta$ is located in $\vec{y}_{\text{\tiny pk}}$ such that $\partial_{i}\zeta(\vec{y}_{\text{\tiny pk}})=0$, or

[TABLE]

where we have used in the last passage the notation $\partial_{i}\partial_{j}\zeta(\vec{x}_{\text{\tiny pk}})=\zeta_{ij}(\vec{x}_{\text{\tiny pk}})$. Performing a rotation of the coordinate axes to be aligned with the principal axes of the constant-$\zeta$ ellipsoids gives the eigenvalues of the shear tensor $\zeta_{ij}$ to be equal to $-\sigma_{2}\lambda_{i}$, where $\sigma_{2}$ is the characteristic root-mean-square variance of the components of $\zeta_{ij}$ (that of $\partial_{i}\zeta$ is $\sigma_{1}$) and

[TABLE]

The crucial point is now that the moments $\sigma_{j}^{2}$ are typically much smaller than $(aH)^{j}$ (because of the presence of the amplitude of the power spectrum). From Eq. (3.4), we deduce that

[TABLE]

and therefore $\lambda_{i}\sim\gamma\nu\gg 1$ (the probability to have negative eigenvalues is small for large curvatures around the peak [11]) . This implies

[TABLE]

where in the last equality we have used the fact that \sigma_{1}/\sigma_{2}\simeq k_{\star}^{-1}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt\sim}\hss}\raise 1.0pt\hbox{<}}R_{H}. Therefore the high overdensity peaks in $\delta$ lie close to the peaks of the curvature perturbation (i.e. within the Hubble volume) if the latter are characterised by a large second derivatives at the origin of the peak. This statement if of course valid in the probabilistic sense.

Since some approximations have been made along the way, in Appendix B the reader can find a numerical simulation we have performed to support this result.

3.2 The calculation of the probability from peak theory

If the argument above is correct, one can associate the number of rare peaks in the overdensity with the number of peaks in the curvature perturbation which are spiky enough, see Eq. (3.7). Therefore, expanding around the peak location $\vec{x}_{\text{\tiny pk}}$ of $\zeta$ (where $\partial_{i}\zeta(\vec{x}_{\text{\tiny pk}})=0$) we can write

[TABLE]

where

[TABLE]

Since the number of peaks (if spiky enough) in $\zeta$ is approximately the number of peaks in $\delta$, we can use the expression (A.14) of Ref. [11] to find the number of peaks of the overdensity

[TABLE]

where

[TABLE]

and $f(x)$ is provided by the expression

[TABLE]

Thus the number density of non-Gaussian peaks of the overdensity above a given threshold $\delta_{{\text{\tiny pk}}}^{\rm c}$ is simply given by

$\displaystyle{\cal N}_{\text{\tiny pk}}=\int_{-\infty}^{\infty}{\rm d}\nu\int_{x_{\delta}^{\rm c}(\nu)}^{\infty}{\rm d}x\,\frac{e^{-\nu^{2}/2}}{(2\pi)^{2}R_{*}^{3}}f(x)\frac{{\rm exp}[-(x-x_{*})^{2}/2(1-\gamma^{2})]}{[2\pi(1-\gamma^{2})]^{1/2}},$

where

[TABLE]

accounts for the fact that only large enough Laplacian values at the peak of the curvature perturbation have to be accounted for, see Eq. (3.4). Notice that if we take the lower limit (3.15) at $\nu=0$, $x_{\delta}^{\rm c}(0)\simeq(9a^{2}H^{2}/4\sigma_{2})\delta^{\rm c}_{\text{\tiny pk}}$, we automatically reproduce the gaussian case. We have checked numerically that in such a case, the peak theory abundance of PBHs obtained from the number density (3.2) with $x_{\delta}^{\rm c}(0)$ reproduces the abundance (1.5) within a factor of order unity. This gives us extra confidence that identifying large threshold peaks in $\delta$ with the spiky enough peaks in $\zeta$ is a correct procedure. From the expression above one can see that the narrower is the power spectrum (that is the closer to unity is the parameter $\gamma$) the more the integrand is peaked at the value $x\simeq x_{*}\simeq\nu$.

We conclude that the non-liner relation between the curvature perturbation and the overdensity makes it harder to generate PBHs, independently from the shape of the curvature perturbation power spectrum.

From the knowledge of the number density of peaks ${\cal N}_{\text{\tiny pk}}$ we can compute the mass fraction of PBHs $\beta$ at the time of the formation $t_{\rm f}$. Since PBHs trace the peaks of the radiation density field on superhorizon scales and since the number of peaks per comoving volume is constant, the number of enough sizeable peaks on superhorizon scales provides the number of PBHs formed once the overdensity has crossed the horizon and one has properly rescaled it to the formation time [8, 10].

The next question is therefore what defines the horizon crossing. In cosmology we are used to the concept of the horizon crossing associated to a given comoving wavelength $k^{-1}$ and we say that horizon crossing takes place when $k=aH$. In the case of PBHs, the large inhomogeneities have characteristic profiles in coordinate space and therefore it is not immediate to associate to them a given wavelength or momentum. The procedure we will follow is the one adopted to define the threshold for collapse [10]. Suppose the overdensity has an average profile in real space given by [11]444As mentioned already in the introduction, we do not include here the non-Gaussianities in the average overdensity profile, whose effect we will study elsewhere [35]. As for the variance around the average profile, it is negligible for $\delta_{\text{\tiny pk}}^{\rm c}/\sigma_{\delta}\gg 1$.

[TABLE]

where $\xi_{2}(r,t)$ is the two-point correlator. One can define a scale $r_{m}$ through the relation

[TABLE]

This scale is relevant since one can show that the threshold for PBH formation is given by [8, 10]

[TABLE]

where $\delta_{\rm c}=3\overline{\delta}(r_{m},t_{m})$, since $r_{m}$ is precisely the scale at which the compaction function ${\cal C}\simeq 2\delta M/ar$ (being $\delta M$ the overmass generated by the averaged curvature perturbation) is maximised [10]. Such a maximum is located at distances larger than the cosmological horizon. It is then natural to define the “horizon crossing” as the time at which555The condition should read $e^{\zeta(r_{m})}a(t_{m})H(t_{m})r_{m}=1$, but $\zeta(r_{m})\ll\zeta(x_{\text{\tiny pk}})$ and we can safely neglect this correction. $a(t_{m})H(t_{m})r_{m}=1$. Numerical simulations must provide a relation between the scale $r_{m}$ and the characteristic momentum appearing in the power spectrum of the curvature perturbation.

The mass fraction at formation time (that is when the horizon forms) from peak theory will then be

[TABLE]

where $M(R_{H})$ is the mass of the PBH associated with the horizon size 666In case, one can take into account that the PBH mass is not precisely the expression $M(R_{H})$, but scales with the initial perturbations [37]. $R_{H}$,

[TABLE]

and $\overline{\rho}_{\rm f}$ and $\overline{\rho}_{m}$ are the background radiation energy densities at the time of formation and horizon crossing, respectively. Numerical simulations show that the ratio $a_{f}/a_{m}$ is rather independent from the shape of the power spectrum and $\sim 3$ [8]. We therefore have

[TABLE]

3.3 The log-normal power spectrum

We assume a power spectrum of the form

[TABLE]

where changing the value of $\sigma$ changes the broadness of the power spectrum. For the case at hand it turns out that [8]

[TABLE]

and one has to choose the critical value $\delta_{{\text{\tiny pk}}}^{\rm c}=1.16$ corresponding to $\delta_{\rm c}=0.51$ [8, 10].

In Fig. 1 we plot the mass fraction for various values of $\sigma$ as a function of $A_{g}$.

We see that the inclusion of the intrinsic non-Gaussian effects systematically lowers the PBH abundance (having kept fixed the amplitude of the power spectrum of the curvature perturbation). Said in other words, keeping the amplitude of the fluctuations fixed, it is more difficult to generate PBHs. This will remain true also using the threshold statistics, as we show in the next Section. Quantitatively, in the case considered, for the usual value of $\beta\sim 10^{-15}$ necessary for PBHs to be all the dark matter in the universe for masses of the order of $10^{-12}\,M_{\odot}$, we find that in the gaussian case the value of the amplitude is consistent with the one reported in Ref. [8] once the difference in the normalisation of the power spectrum is taken into account, while the non-Gaussian abundance is suppressed.

3.4 Broad power spectrum

We also consider a broad power spectrum, that is a top-hat function with amplitude $A_{t}$ as

[TABLE]

where $\Theta$ stands again for the Heaviside step function and $k_{\rm max}\gg k_{\rm min}$, such that the scale $k_{\rm min}$ in practice does not participate in the PBH formation [8]. In this case one finds $k_{\rm max}\simeq 3.5/r_{m}$, $\delta_{\rm c}$ is again 0.51, and $\delta_{\text{\tiny pk}}^{\rm c}\simeq 1.22$ [8] and the variances are obtained by putting $a_{m}H_{m}$ as the infrared cut-off since the unphysical long wavelength modes should be disregarded.

Fig. 2 shows the mass fraction as a function of $A_{t}$.777We do not introduce a window function to be able to compare with the gaussian results of Ref. [8] which are reproduced in the gaussian case. As predicted, both for narrow spectra and broad ones, the intrinsic non-Gaussianity in the overdensity makes it harder to produce PBHs.

4 The non-Gaussian probability from threshold statistics

In this section we present an alternative way to calculate the non-Gaussian probability to form PBHs which does not rely on the fact that spiky peaks of the curvature perturbation coincide with peaks of the overdensity for large thresholds. The price to pay is that we will be dealing with the threshold statistics (the threshold being identified with $\delta_{\rm c}$ [8]). This might be not a great sacrifice as regions characterised by large thresholds are likely to be regions of maxima of the overdensity [40]. The gain is that the expressions we are going to obtain are exact.

Let us consider again the curvature perturbation $\zeta(\vec{x})$ as a random field. Following the notation of the Appendix A of Ref. [11], we define

[TABLE]

The correlations of these fields are provided by the expressions

[TABLE]

These variances will be computed numerically using the Fourier transform of the top-hat window function in real space, that is

[TABLE]

The matrix $-\zeta_{ij}$ can be diagonalized with eigenvalues $\sigma_{2}\lambda_{i}$, ordered such that $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$. Thus we define

[TABLE]

Introducing again $\nu=\zeta(\vec{x})/\sigma_{0}$, the correlations become

[TABLE]

and all the others are zero. The joint gaussian probability distribution for these variables is provided by the expression (from now on we will label $\eta_{i}\equiv\zeta_{i}$)

[TABLE]

as a function of

[TABLE]

and

[TABLE]

The variables $y$ and $z$ are unconstrained and we integrate them out. With the ordering of the eigenvalues previously defined, we see that the variable $z$ lies in the range $[-y,y]$, while $y\geq 0$. The result is therefore given by888 Notice that, assuming the linear relation between $\delta$ and $\zeta$ as in Eq. (1.8), one recovers the Press-Schechter result in Eq. (1.2).

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

We can then write the $\delta$ as a function of these variables as

[TABLE]

Now we perform the change of variables:

[TABLE]

The argument of the logarithm is positive for $x_{\delta}>\vec{\eta}_{\delta}\cdot\vec{\eta}_{\delta}/2\sigma_{2}$. The Jacobian of the transformation is given by

[TABLE]

Therefore the distribution in terms of the new variables is given by

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

Finally, since the probability distribution is only a function of the modulus $\vec{\eta}_{\delta}\cdot\vec{\eta}_{\delta}=\eta_{\delta}^{2}$, one can change variable as ${\rm d}^{3}\eta_{\delta}=\eta_{\delta}^{2}\sin\theta_{\delta}{\rm d}\eta_{\delta}{\rm d}\theta_{\delta}{\rm d}\phi_{\delta}$ and perform the integration on the angles which trivially results in

[TABLE]

Finally we get

$\displaystyle\beta_{\text{\tiny NG}}^{\text{\tiny th}}$ $\displaystyle=$ $\displaystyle 4\pi\int_{\delta_{\rm c}}{\rm d}\delta\int_{0}^{\infty}{\rm d}\eta_{\delta}\,\eta_{\delta}^{2}\int_{\eta^{2}_{\delta}/2\sigma_{2}}^{\infty}{\rm d}x_{\delta}\,D(\delta,x_{\delta},{\eta}_{\delta})e^{-Q_{3}}.$

(4.24)

This is an exact result, no approximations have been made at this stage999We checked that using Eq. (4.24) gives the same numerical result obtained by computing the probability of the overdensity integrating Eq. (4.14) with the insertion of a Heaviside function of the form $\Theta\left(\delta-\delta_{\rm c}\right)$ leading to the limit of integration in the variable $x$ given by the condition $x>(9a^{2}H^{2}/4\sigma_{2}){\rm exp}(2\nu\sigma_{0})\delta_{\rm c}$..

4.1 Spiky power spectrum

In the limit of $\gamma\simeq 1$, i.e. for power spectra whose width is very narrow (typical of the PBHs), we can simplify our expressions dramatically. First of all, from Eq. (4.16) one sees that the distribution in $x_{\delta}$ becomes a Dirac delta centered in $x_{*}\simeq\nu$. We then obtain

[TABLE]

where

[TABLE]

Then, to perform the integral in ${\rm d}\eta_{\delta}$, we rewrite the Dirac delta as

[TABLE]

where

[TABLE]

and where we have chosen the positive root since $\eta_{\delta}$ is always positive. The root imposes the condition

[TABLE]

which is solved by ($W_{0}$ and $W_{-1}$ are the so-called principal and negative branches of the Lambert function)

[TABLE]

with the requirement that101010The condition (4.31) leaves a really narrow window in terms of $\delta$, $0<\delta<0.59$. This means that the Dirac delta power spectrum would not be a good choice where the threshold is larger than 0.59. Of course, such a monochromatic power spectrum is only an approximation for more physical narrow power spectra.

[TABLE]

After integrating in ${\rm d}\eta_{\delta}$, we find that the joint probability is

[TABLE]

This means that the threshold probability is

$\displaystyle\beta_{\text{\tiny NG}}^{\text{\tiny th}}$ $\displaystyle=$ $\displaystyle\int_{\delta_{\rm c}}^{\delta_{+}}{\rm d}\delta\int_{x_{-}(\delta)}^{x_{+}(\delta)}{\rm d}x_{\delta}\,\frac{54\sqrt{3}}{8\pi\sqrt{2}}\frac{a^{3}H^{3}}{\sigma_{1}^{3}}\sqrt{4x_{\delta}\frac{\sigma_{2}}{a^{2}H^{2}}-9\delta e^{2\sigma_{0}x_{\delta}}}$

(4.33)

$\displaystyle\times$ $\displaystyle\exp\left[-\frac{1}{2}x_{\delta}^{2}+2\sigma_{0}x_{\delta}-\frac{3a^{2}H^{2}}{4\sigma_{1}^{2}}\left(4x_{\delta}\frac{\sigma_{2}}{a^{2}H^{2}}-9\delta e^{2\sigma_{0}x_{\delta}}\right)\right],$

where the higher extremum of integration in $\delta$ is due to (4.31). In Fig. 3 one can find the comparison of the gaussian and non-Gaussian mass functions computed using the threshold statistics for a spiky power spectrum.

To proceed further and provide more analytical insights, we notice that the integration over $x_{\delta}$ in Eq. (4.33) is highly dominated by the lower extremum of integration $x_{-}\left(\delta\right)$. As we show in Appendix C, the integrand in this region is very well approximated by

[TABLE]

Since the integral in the second line is highly dominated by the lower extremum of integration, we can set $x_{+}(\delta)\to\infty$ and perform the integration analytically, obtaining (for $\gamma\simeq 1$)

[TABLE]

In Appendix C we show that this expression is extremely accurate for the case of a Dirac delta power spectrum of the curvature perturbation.

We can perform the final integral (4.35) by changing the variable of integration from $\delta$ to $x_{-}\left(\delta\right)$. The lower and higher extrema of integration then become, respectively, $x_{-}\left(\delta_{\rm c}\right)$ and $1/2\sigma_{0}$. The integrand is highly dominated by the region around the lower extremum, so that we can send the higher extremum to infinity. We can also evaluate all the integrand, apart from the exponential factor ${\rm exp}(-x_{-}^{2}\left(\delta\right)/2)$ at the lower extremum. The integral of this exponent can be then done analytically, and its result (the complementary error function) can be expanded in the limit of large argument. This leads to

[TABLE]

The accuracy of this result is shown in Figure 10 of Appendix C, performed for the case of a Dirac delta power spectrum of the curvature perturbation, where it is compared with a two-dimensional numerical integration of the starting expression (4.33).

4.2 Log-normal power spectrum

We assume again a power spectrum of the form

[TABLE]

Then, one can integrate Eq. (4.24) numerically to get the mass fraction. In Fig. 4 we plot the beta for various values of $\sigma$ as a function of $A_{g}$.

4.3 Broad power spectrum

We also consider a broad power spectrum, that is a top-hat with amplitude $A_{t}$ as

[TABLE]

where $\Theta$ stands for the Heaviside step function and $k_{\rm max}\gg k_{\rm min}$. Again, the parameters used are $k_{\rm max}\simeq 3.5/r_{m}$, $\delta_{\rm c}=0.51$ [8] and, to disregard unphysical long wavelength modes, variances are obtained by choosing $a_{m}H_{m}$ as the infrared cut-off. The results are presented in Fig. 5.

We conclude that threshold statistics confirms what we found in peak theory: independently from the power spectrum, non-Gaussian abundances are smaller than the gaussian ones.

We also see that the difference between the gaussian and the non-Gaussian cases in terms of the amplitude of the power spectrum is about a factor $(2\div 3)$, the same for the Dirac delta case. This is the shift one should adopt if insisting in using the gaussian expressions.

5 Conclusions

In this paper we have discussed the impact of the non-Gaussianity arising from the non-linear relation between the density contrast and the curvature perturbation when dealing with PBH abundances. We have proposed two different methods to deal with such unavoidable and intrinsic non-Gaussianity, providing simple analytical expressions for the abundance to take it into account.

The first method is based on peak theory and on the realisation that the number of peaks in the overdensity is approximately equal to the number of peaks in the curvature perturbation as long as one restricts her/himself to those peaks having large spatial second derivatives at the peak location.

The second method relies on the threshold statistics and contains no approximations. Both methods show that the intrinsic non-Gaussianity makes it harder to generate PBHs. In particular, if one insists in adopting the gaussian expression for the abundance coming from threshold statistics, one has simply to increase the amplitude of the power spectrum by a factor 111111Our findings agree with those recently obtained in Refs. [41, 42], where the PBH abundance has been derived using the averaged (over a volume of radius $r_{m}$) density fluctuation constructed out of a radial profile in the curvature perturbation $\zeta$. Adopting the volume averaged density provides a clear relation between the linear gaussian component of the peak height and the non-linear peak height. However, in order to make use consistently of the obtained critical threshold value, one needs to identify peaks in $\zeta$ with the peaks in $\delta$, which we have shown here is true for spiky peaks in $\zeta$. In Ref. [7] the authors computed the abundance using peak theory for the comoving curvature perturbation by setting a threshold on the $\zeta$, contrarily to our choice of expressing the abundance and the threshold in terms of the overdensity field. ${\cal O}(2\div 3)$.

Our findings do not alleviate the differences between peak theory and threshold statistics in the computation of the abundance, already present at the gaussian level [12].

Our results can be surely improved along some directions. It would be important to have a full non-Gaussian extension of peak theory. More importantly, the intrinsic non-Gaussianity of the overdensity is expected to change the shape of the profile of the peaks which eventually give rise to PBHs upon collapse. Since the threshold $\delta_{\text{\tiny pk}}^{\rm c}$ depends on the shape of the overdensity, such non-Gaussianity might change as well the value of $\delta_{\text{\tiny pk}}^{\rm c}$. We leave this study for a future publication [35].

Acknowledgments

We thank I. Musco for many fruitful discussions. We thank C. Byrnes, I. Musco and S. Young for sharing their draft [42] with us and for useful interactions. M.P. and C.U. thank Geneva Theoretical Physics group for their hospitality during this project. V. DL. thanks the Galileo Galilei Institute for Theoretical Physics (Florence, Italy) for the nice hospitality during the realisation of this project. A.R., V. DL., and G.F. are supported by the Swiss National Science Foundation (SNSF), project The non-Gaussian Universe and Cosmological Symmetries, project number: 200020-178787. C.U. is supported by European Structural and Investment Funds and the Czech Ministry of Education, Youth and Sports (Project CoGraDS- CZ.02.1.01/0.0/0.0/15 003/0000437) and would like to acknowledge networking support by the COST Action GWverse CA16104.

Appendix A The cumulants for a narrow power spectrum

In this Appendix we derive the relations (2.11) of the main text. We start from Eq. (1.7), that we need to expand as a power series of $\zeta$. We denote by $\delta_{n}$ the term that is of ${\cal O}\left(\zeta^{n}\right)$

[TABLE]

Using the convention

[TABLE]

for the Fourier transform of the curvature perturbation, and symmetrizing over the momenta $p_{i}$ of the Fourier modes, the above relations can be cast in the form

[TABLE]

We are interested in computing the connected $2$-, $3$- and $4$-point correlation functions of $\delta\left(0\right)=\sum_{n=1}^{\infty}\delta_{n}$, where by connected we mean terms that cannot be factorized as products of smaller-order correlation functions. Under the assumption of Gaussianity of the curvature $\zeta$, all the correlators can be broken down to the products of the two-point function of $\zeta$,

[TABLE]

The practical effect of computing a connected, rather than a full, correlator is that some of the contractions are not included. To give just one example, we have

[TABLE]

with

[TABLE]

with the omission of the $\left\langle\zeta\left(\vec{p}_{1}\right)\zeta\left(\vec{p}_{2}\right)\right\rangle\,\left\langle\zeta\left(\vec{q}_{1}\right)\zeta\left(\vec{q}_{2}\right)\right\rangle$ term.

More in general, we note that the first cumulants are related to the full correlators by

[TABLE]

It is worth noting that only the first cumulant is affected by the average of $\delta$. In fact, the expressions (A.7) show that a shift $\delta\rightarrow\delta+C$, where $C$ is a constant, only affects the first cumulant, $\left\langle\delta\right\rangle_{c}\rightarrow\left\langle\delta\right\rangle_{c}+C$, while the higher cumulants are unchanged.

Working up to cubic order in the power of $\zeta$, we compute

[TABLE]

where we have kept together terms that are of the same order in $P_{\zeta}$. We note that the last expression does not contain the contraction of $\delta_{1}^{4}\left(0\right)$ as it has no connected component.

The evaluations of the correlators in (A.8) is tedious, but straightforward. We expand the various terms according to (A.3) and we then split the correlators in sums of connected products of $\left\langle\zeta\left(\vec{p}\right)\zeta\left(\vec{q}\right)\right\rangle$. Half of the integrals over momenta are then removed with the Dirac delta functions arising from Eq. (A.4). We divide the remaining half into integrals over the magnitude of the momenta and the angles. We encounter the following nontrivial angular integrals

[TABLE]

The explicit evaluations then give

[TABLE]

In the case of a very narrow power spectrum of the curvature perturbation, that can be approximated by a Dirac delta function as in Eq. (2.10), these expressions give the results (2.11) reported in the main text.

Appendix B Spiky peaks in the curvature perturbation versus peaks in overdensity: a numerical treatment

We start from the relation between $\delta$ and $\zeta$

[TABLE]

One can simulate numerically a realisation of the gaussian random field $\zeta(\vec{x})$ in a $n$-dimensional box of dimensions $N$ which is discretised using a grid of $N^{n}$ points with a spacing $\Delta x=1$ between them in all directions. We choose to present the analysis in a 2-dimensional space $(n=2)$ since the results can be more easily depicted. We set the parameters of the perturbation assuming a narrow power spectrum described by a log-normal function as

[TABLE]

The variance of the field turns out to be $\sigma_{0}^{2}=2.5\cdot 10^{-3}$. The characteristic momentum has been chosen to be $k_{\star}=0.2/\Delta x$. The realisation of the field $\zeta(\vec{x})$ and the corresponding field $\delta_{\rm r}(\vec{x})$ can be seen in Fig. 6. There the stars indicate the location of the spiky peaks in $\zeta$ and the peaks in $\delta_{\rm r}$, showing the location correspondence. The color code is the same as in Fig. 7.

In Fig. 7 one can find an analysis of the field values obtained in the simulation. More in detail, each point of the plot represents a peak in $\zeta$ with the corresponding values of the rescaled amplitude $\nu$ and the curvature $x$. The red, cyan and yellow lines correspond to lower bounds on $x>x_{\delta}^{\rm c}(\nu)$ in terms of the absolute maximum of the density contrast $\delta_{\max}$ in the simulation as

[TABLE]

with $\delta_{\rm max}=0.4$. This bound corresponds to the condition (3.15). With red, cyan and yellow dots we highlight the points which, at the same positions, have a peak in $\delta$, with $x$ satisfying the corresponding lower limits. Green dots are peaks in $\zeta$ as well, but they do not satisfy these conditions. This shows the correspondence between peaks of $\zeta$ and peaks of $\delta$, provided the condition (3.15) is met. We expect that this correspondence will be even more satisfied when rarer events are simulated. We also checked that, by extending the simulation to three dimensions, and these findings are confirmed.

These results strongly indicates that, assuming condition (3.15), peaks in $\delta$ are located at the positions of peaks in $\zeta$.

Appendix C Analytic integration of the PBH abundance for spiky power spectra using threshold statistics

In this appendix we derive the expressions (4.34) and (4.35) of the main text. We start from Eq. (4.33). One can verify that the integration over $x_{\delta}$ of this equation is highly dominated by the lower extremum $x_{-}\left(\delta\right)$ (from now on, in this appendix, we do not write the dependence of $x_{-}$ on $\delta$ to shorten the notation). We therefore perform an expansion of the integrand for $x_{\delta}\simeq x_{-}$ that allows us to perform the integration analytically. We expand the expression in the square root and in the exponent by linearising the exponential in $x_{\delta}-x_{-}$

[TABLE]

where the second line has been obtained exploiting the fact that $x_{-}$ satisfies (exactly) $\delta{\rm e}^{2\sigma_{0}x_{-}}=4\sigma_{2}x_{-}/9a^{2}H^{2}$. We also approximate the first two terms in the exponent of Eq. (4.33) as

[TABLE]

where we have linearised the first term on the left-hand side to first order in $x_{\delta}-x_{-}$, while in the second term we simply put $x_{\delta}=x_{-}$ (since this term is highly subdominant). With these approximations, the expression (4.33) reduces to the form (4.34) written in the main text.

The integration over $x_{\delta}$ in Eq. (4.34) is highly dominated by the lower extremum of integration, and we can set $x_{+}\to\infty$. In this way the integration can be done analytically, leading to

[TABLE]

Recalling that these results are valid for $\gamma\equiv\frac{\sigma_{1}^{2}}{\sigma_{0}\,\sigma_{2}}\simeq 1$ then leads to the expression (4.35) written in the main text.

In the case of a Dirac delta power spectrum of the curvature perturbation $\zeta$, see Eq. (2.10), we have $\sigma_{i}=w\sqrt{A_{s}}\,k_{\star}^{i}$, where $w=W(k_{\star},r_{m})$. Recalling that $k_{\star}\simeq(27/10)a_{m}\,H_{m}$, the probability distribution reduces to

[TABLE]

where on the right-hand side we have defined ${\hat{x}}\equiv 2w\sqrt{A_{s}}x_{\delta}$ and ${\hat{x}}_{-}\equiv 2w\sqrt{A_{s}}x_{-}=-W_{0}\left(-50\delta/81\right)$ (which is the expression of the first root in eq. (4.30) in the present case). Figure 8 confirms the validity of this result. The probability in the figure is shown for ${\hat{x}}\simeq{\hat{x}}_{-}\simeq 0.54$ (for the value of $\delta$ chosen in the figure), while ${\hat{x}}_{+}\simeq 1.67$. We note that indeed this expression is highly dominated by the lower bound ${\hat{x}}\simeq{\hat{x}}_{-}$ (this extends also for the values of ${\hat{x}}$ not shown in the figure).

The integration over $x_{\delta}$ of this expression leads to

[TABLE]

where we stress that ${\hat{x}}_{-}$ depends on $\delta$. The higher extremum of integration is the upper bound in Eq. (4.31) written in the present context. This result is extremely accurate, as we show in Figure 9.

The expression (C.5) can be integrated, proceeding as we did in the main text to obtain the result (4.36) from (4.35). We obtain

[TABLE]

This expression also follows immediately from (4.36), in the limit of Dirac delta power spectrum of the curvature perturbation, and noting that $x_{\rm c}={\hat{x}}_{\rm c}/2\,\sigma_{0}={\hat{x}}_{\rm c}/2\,w\sqrt{A_{s}}$. The high accuracy of this result is shown in Figure 10, where we compare it with a fully numerical two-dimensional integration of the starting expression (4.33).

Appendix D Peaks versus thresholds

In the past literature PBHs have been identified either with peaks or with thresholds of the superhorizon overdensity, where by thresholds one means those regions in real space where the value of the density contrast is larger than a given threshold, in our case the critical value $\delta_{\rm c}$. Regions characterised by large thresholds of the overdensity are indeed probable to be also local extrema. We first find the average threshold statistics profile $\overline{\delta}(r)$ of the density contrast $\delta(r)$ at a given distance $r$ from the point $r=0$ (therefore without threshold) in the following way

[TABLE]

where

[TABLE]

and $\delta_{0}=\delta(0)$. If both $\delta(r)$ and $\delta_{0}$ are Gaussian variables, one can calculate the above quantity by recalling that $P(\delta(r),\delta_{0})$ is constructed in the standard way through the covariance matrix

[TABLE]

where

[TABLE]

is the two-point correlator in coordinate space. From these expressions we derive

[TABLE]

where ${\rm Erfc}(x)$ is the complementary error function. Combining the different terms we finally get

[TABLE]

Using the expansion for large values of the argument

[TABLE]

we can finally evaluate the average $\overline{\delta}(r)$ at distance $r$ from the threshold for $\nu\gg 1$

[TABLE]

Taking $\nu=\delta_{\text{\tiny pk}}/\sigma_{\delta}$ one finds

[TABLE]

which is exactly the average profile derived in peak theory [11]121212For the non-Gaussian extension of this result, see Ref. [35].. This already suggests that large thresholds overdensity should correspond to extrema. To have further evidence, we follow Ref. [40] and consider the curvature of the large threshold regions. The mean value of the second derivative of $\delta(r)$ in any random direction at $r=0$ is (by expanding the density contrast around the origin in powers of $r$ and taking the mean value of it) with $\delta(0)=\delta_{\text{\tiny pk}}$

[TABLE]

The scatter of the second derivative from its mean value is found by averaging over all ${\rm d}^{2}\delta(r)/{\rm d}r^{2}|_{r=0}$ and $\delta_{\text{\tiny pk}}$, yet keeping $\delta(0)=\delta_{\text{\tiny pk}}$,

[TABLE]

We then get

[TABLE]

where we have taken $\xi_{2}(0)\sim\sigma_{\delta}^{2}$ and assumed that the profile varies over a characteristic scale $r_{m}$. The condition to have large threshold, that is $\delta_{\text{\tiny pk}}\gg\sigma_{\delta}$, implies that large threshold regions are most likely to be local extrema, that is peaks.