Non-Gaussian Formation of Primordial Black Holes: Effects on the Threshold

TL;DR

This paper studies how intrinsic non-Gaussianity affects the formation threshold of primordial black holes, finding that the threshold is relatively robust against non-linearities and small non-Gaussian effects.

Contribution

It provides a detailed analysis of the impact of non-Gaussianity on the critical threshold for primordial black hole formation, showing minimal effect from non-linearities.

Findings

The threshold $\delta_c$ changes by a few percent due to non-Gaussianity.

The threshold remains robust against small local primordial non-Gaussianity.

Non-linear relations introduce minor corrections to black hole formation criteria.

Abstract

Primordial black holes could have been formed in the early universe from sufficiently large cosmological perturbations re-entering the horizon when the Universe is still radiation dominated. These originate from the spectrum of curvature perturbations generated during inflation at small-scales. Because of the non-linear relation between the curvature perturbation $\zeta$ and the overdensity $\delta\rho$, the formation of the primordial black holes is affected by intrinsic non-Gaussianity even though the curvature perturbation is Gaussian. We investigate the impact of this non-Gaussianity on the critical threshold $\delta_c$ which measures the excess of mass of the perturbation, finding a relative change with respect to the value obtained using a linear relation between $\zeta$ and $\delta\rho$, of a few percent. This shows that the value of the critical threshold is rather robust against non-linearities. The same holds also when cosmologically interesting values of local primordial non-Gaussianity are added to the curvature perturbation.

Figures

Tables

Table 1. Table 1: Values and variations of the basic quantities.

Type	$\delta {\rho }_c$	${\delta }_c$	$x_m$	$r_0/r_m$	$\delta {\rho }_{0_G}$	$\Delta \delta {\rho }_0$	$\Delta {\delta }_c$
linear	1.218	0.516	2.744	1.145	1.218	$-$	$-$
$f_{\text{NL}}=0$	1.269	0.511	2.722	1.160	0.607	0.042	0.010
${\overline{\zeta }}_{\text{NL}}$	0.716	0.582	2.989	1.027	$-$	0.412	0.129
$f_{\text{NL}}=1$	1.317	0.507	2.620	1.178	0.507	0.081	0.017
$f_{\text{NL}}=2$	1.334	0.504	2.549	1.189	0.468	0.095	0.022
$f_{\text{NL}}=3$	1.341	0.503	2.497	1.197	0.427	0.101	0.025
$f_{\text{NL}}=4$	1.343	0.502	2.458	1.201	0.394	0.103	0.027
$f_{\text{NL}}=-1/2$	1.218	0.515	2.789	1.148	0.664	0.00	0.002
$f_{\text{NL}}=-1$	1.125	0.521	2.871	1.133	0.736	0.076	0.010
$f_{\text{NL}}=-3/2$	0.952	0.529	2.970	1.116	0.829	0.218	0.025
$f_{\text{NL}}=-2$	0.626	0.542	3.084	1.097	0.949	0.486	0.050

Equations

$$\beta \equiv {\frac{{\rho }_{\text{PBH}}}{{\rho }_{\mathrm{tot}}}|}_{\mathrm{form}}=P_{\text{G}}(\delta >{\delta }_c)={\int }_{{\delta }_c}\frac{\mathrm{d}\delta }{\sqrt{2\pi }\sigma }{e}^{-{\delta }^2/2{\sigma }^2}\simeq \sqrt{\frac{1}{2\pi }}\frac{\sigma }{{\delta }_c}{e}^{-{\delta }_c^2/2{\sigma }^2}$$

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

$$\delta \rho (\vec{x},t)=-{\left(\frac{2\sqrt{2}}{3aH}\right)}^2{e}^{-5\zeta (\vec{x})/2}{\nabla }^2{e}^{\zeta (\vec{x})/2}.$$

$$\delta \rho (\vec{x},t)\simeq -{\left(\frac{2}{3aH}\right)}^2{\nabla }^2\zeta (\vec{x}),$$

$$\mathrm{d}{s}^2=-\mathrm{d}{t}^2+a^2(t)\left[\frac{\mathrm{d}{r}^2}{1-K(r)r^2}+r^2\mathrm{d}{\Omega }^2\right]=-\mathrm{d}{t}^2+a^2(t)e^{2\zeta (\hat{r})}\mathrm{d}{\vec{x}}^2,$$

$$\{\begin{array}{rl}& r=\hat{r}{e}^{\zeta (\hat{r})},\\ & \frac{\mathrm{d}r}{\sqrt{1-K(r)r^2}}=e^{\zeta (\hat{r})}\mathrm{d}\hat{r}.\end{array}$$

$$K(r)r^2=-\hat{r}{\zeta }^{\prime }(\hat{r})\left[2+\hat{r}{\zeta }^{\prime }(\hat{r})\right],$$

$$\delta \rho \equiv \frac{\rho (r,t)-{\rho }_b(t)}{{\rho }_b(t)}=\frac{1}{a^2{H}^2}\frac{(1+\omega )}{5+3\omega }\frac{{\left[K(r)r^3\right]}^{\prime }}{r^2}=-\frac{1}{a^2{H}^2}\frac{4(1+w)}{5+3w}{e}^{-5\zeta (\hat{r})/2}{\nabla }^2{e}^{\zeta (\hat{r})/2}.$$

$$p=\omega \rho ,$$

$$\mathcal{C}\equiv 2\frac{\delta M(r,t)}{R(r,t)},$$

$$\mathcal{C}=\frac{3(1+w)}{5+3w}K(r)r^2,$$

$${\mathcal{C}}^{\prime }(r_m)=0\text{or}{\mathcal{C}}^{\prime }({\hat{r}}_m)=0$$

$$K(r_m)+\frac{r_m}{2}{K}^{\prime }(r_m)=0\mathrm{or}{\zeta }^{\prime }({\hat{r}}_m)+{\hat{r}}_m{\zeta }^{\prime \prime }({\hat{r}}_m)=0.$$

$$ϵ\equiv \frac{R_H(t)}{R_b(r_m,t)}=\frac{1}{aH{r}_m}=\frac{1}{aH{\hat{r}}_m{e}^{\zeta ({\hat{r}}_m)}},$$

$$\delta \rho ={\left(\frac{1}{aH}\right)}^2\frac{3(1+\omega )}{5+3\omega }\left(K(r)+\frac{r}{3}{K}^{\prime }(r)\right)=-{\left(\frac{1}{aH}\right)}^2\frac{2(1+\omega )}{5+3\omega }\left[{\zeta }^{\prime \prime }(\hat{r})+{\zeta }^{\prime }(\hat{r})\left(\frac{2}{\hat{r}}+\frac{1}{2}{\zeta }^{\prime }(\hat{r})\right)\right]e^{-2\zeta (\hat{r})},$$

$$U=H(t)R(1+\delta U),$$

$$R=a(t)r(1+\delta R),$$

$$\delta U=-\frac{1}{(1+\omega )}\frac{1}{r^3}\int r^2\mathrm{d}r\delta \rho ,$$

$$\delta R=-\frac{\omega }{(1+3\omega )(1+\omega )}\delta \rho +\frac{1}{1+3w}\delta U.$$

$$\delta (r_m,t_H)=\frac{4\pi }{V_{R_m}}{\int }_0^{R_m}\mathrm{d}R{R}^2\delta \rho =\frac{3}{r_m^3}{\int }_0^{r_m}\mathrm{d}r{r}^2\delta \rho ,\text{where}{V}_{R_m}=\frac{4\pi }{3}{R}_m^3.$$

$${\delta }_m\equiv \delta (r_m,t_H)=\frac{3(1+w)}{5+3w}K(r_m)r_m^2=\mathcal{C}(r_m)$$

$${\delta }_m=3\delta \rho (r_m,t_H).$$

$$r_m^3=\frac{{\int }_0^{r_m}\mathrm{d}r{r}^2\delta \rho (r,t_H)}{\delta \rho (r_m,t_H)},$$

$$\overline{\delta \rho }(r)=⟨\delta \rho (r)|\delta {\rho }_0>\nu \sigma ⟩={\int }_{-\infty }^{\infty }\mathrm{d}\delta \rho (r)\delta \rho (r)P(\delta (r)|\delta {\rho }_0>\nu \sigma ),$$

$$P(\delta \rho (r)|\delta {\rho }_0>\nu \sigma )$$

$$P(\delta \rho (r),\delta {\rho }_0)$$

$${\overrightarrow{\delta \rho }}^T$$

$$C$$

$${\xi }^{(2)}(r)=⟨\delta \rho ({\vec{x}}_1)\delta \rho (\vec{0})⟩$$

$$P(\delta \rho (r),\delta {\rho }_0>\nu \sigma )$$

$$P(\delta {\rho }_0>\nu \sigma )$$

$$\overline{\delta \rho }(r)=\frac{{\xi }^{(2)}(r)}{\sigma }\sqrt{\frac{2}{\pi }}\frac{e^{-{\nu }^2/2}}{\mathrm{Erfc}\left(\nu /\sqrt{2}\right)}.$$

$$\mathrm{Erfc}\left(x\gg 1\right)\approx \frac{e^{-x^2}}{x\sqrt{\pi }},$$

Full text

1. I Introduction

2. II Initial conditions for PBH formation

3. III The average density profile

4. III.1 The Gaussian case

5. III.2 The non-Gaussian case

6. IV The average density profile around the threshold for PBH formation

7. IV.1 The case of a peaked power spectrum

8. IV.2 Including non-Gaussianity of the power spectrum

9. V The average profile including the three-point correlation function

10. V.1 The averaged density profile from threshold statistics

11. V.2 The density profile from the averaged curvature profile $\bar{\zeta}$

12. VI Numerical results

13. VII Conclusions

I Introduction

Primordial Black Holes (PBHs) have recently received much attention starting from the discovery of the gravitational waves emitted by the merging of two $\sim 30\,M_{\odot}$ black holes ligo1 . In particular, the focus has been on the possibility that PBHs may describe the nature of dark matter we observe in the universe bird (see also revPBH and references therein).

Even though there are many ways to generate PBHs in the early universe, the mechanism which has been investigated more extensively in the recent literature is the one obtained from inflation s1 ; s2 ; s3 . During such stage of primordial acceleration, the curvature perturbation $\zeta$ may be enhanced at small-scales with respect to the large-scale perturbation $\zeta\sim 10^{-5}$ which is ultimately responsible for the CMB anisotropies. At cosmological horizon re-entry the small-scale fluctuations in the overdensity $\delta\rho$ might collapse into a PBH if they are large enough to overcome the pressure gradients: a PBH would form if the perturbation amplitude $\delta$ is larger than a given threshold $\delta_{c}$, with a mass of the order of the mass contained within the horizon volume at horizon re-entry. The mechanism of PBH formation has been investigated in details by several authors performing spherically symmetric numerical simulations Jedamzik:1999am ; Shibata:1999zs ; Hawke:2002rf ; Hawke:2002rf ; Musco:2004ak and it has been shown that the critical collapse mechanism Choptuik:1992jv arises when $\delta>\delta_{c}$, with the mass spectrum of PBHs described by a scaling law Niemeyer:1997mt ; Musco:2008hv ; Musco:2012au .

The abundance $\beta$ of PBHs at formation is exponentially sensitive to the threshold (for simplicity we give the Gaussian expression)

[TABLE]

where $\rho_{\text{\tiny PBH}}$ is the energy density collapsed into PBHs while $\rho_{\rm tot}$ is the total energy density. The expression $P_{\text{\tiny G}}(\delta>\delta_{c})$ indicates the Gaussian probability of a perturbation collapsing to a PBH if its amplitude $\delta$ is larger than a certain threshold $\delta_{c}$. Here $\sigma^{2}$ is the variance of the overdensity

[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 an appropriate window function.

Recently the investigation of the value of the threshold $\delta_{c}$ has been very intensive and it has been pointed out that the value of $\delta_{c}$ is not unique, but depends on the shape of the power spectrum of the curvature perturbation haradath ; musco ; mg . In particular the exact value, lineraly extrapolated at horizon crossing, is varying between $0.4$ and $2/3$, depending on the particular initial shape of the curvature/density profile musco , which affects the impact of the pressure gradients during the non-linear evolution of the collapse. This is closely related to the shape of the power spectrum which determines, using peak theory, the average shape of the density perturbation BBKS ; mg .

One point of particular importance is the fact that the overdensity111The notation used here for the density constrast is slightly different from what has been used in the literature. Usually papers on PBH formation are using $\delta\rho/\rho_{b}$ while other papers, coming from the cosmological community, use the simpler notation $\delta$ for the same quantity. Because with $\delta$ here we are referring to the average threshold integrated over the volume, to keep a clear distinction between the two quantities, we have decided to simplify a bit the notation calling the density contrast just as $\delta\rho$, properly defined later in (2.4). $\delta\rho$ and the curvature perturbation $\zeta$ are related to each other by a non-linear relation. In the comoving slicing, when the Universe is radiation dominated, it reads harada

[TABLE]

This implies that, even when the curvature perturbation is a Gaussian random field, the overdensity $\delta\rho$ is intrinsically and unavoidably non-Gaussian ng2 ; ng3 ; ng4 . In the presence of such ineludible non-Gaussianity the abundance of PBHs is significantly reduced compared to the Gaussian (linear) case where one approximates the relation (1.3) as

[TABLE]

and one has to reduce the amplitude of the power spectrum of $\zeta$ by a factor ${\cal O}(2\div 3)$ to have the same non-Gaussian number of PBHs starting from the Gaussian expression ng2 ; ng3 ; ng4 .

The goal of this paper is to assess the impact of the intrinsic non-Gaussianity of the overdensity onto the critical threshold $\delta_{c}$. PBHs are identified with the local maxima of the overdensity and, in order to distinguish whether a cosmological perturbation will collapse forming a PBH, it is crucial to evaluate the amplitude of the peak of the corresponding compaction function. An important input is therefore the shape of the overdensity around the peak since non-linearities would have an effect on the shape of the density, and it is reasonable to expect that the intrinsic non-Gaussianity modifies the critical threshold $\delta_{c}$. As a byproduct, our investigation will allow us to check (and in fact confirm a posteriori) the validity of the assumption made in ng3 where the abundance of the PBHs, including the effect coming from the intrinsic non-Gaussianity, has been performed adopting the critical threshold $\delta_{c}$ derived for the linear Gaussian relation between $\delta\rho$ and $\zeta$.

Our results are based on a perturbative calculation of the average profile around the peak of a perturbation and suggest that the critical threshold is rather robust against the intrinsic non-Gaussianity introduced by the non-linear relation between $\delta\rho$ and $\zeta$. This also remains true if we endow the curvature perturbation with some primordial non-Gaussianity. The relative changes of the critical threshold are of the order of few percent and they do not significantly affect the calculation of the PBH abundance. The reason for this result is based on the close relation between the shape of the density perturbation and the value of the threshold $\delta_{c}$: although the amplitude of the non-linear components is of the same order of the linear one, the effect on the shape due to the non-linear and the non-Gaussian effects are not very significant, and therefore the final shape is quite close to the one obtained using the linear approximation given by (1.4). For this reason we suggest, as other works have done (e.g. ng4 ; Young:2019osy ), that the threshold $\delta_{c}$ allows the computation of the abundance of PBHs with less uncertainties with respect of using the local critical amplitude of the peak which is more sensitive to the local features of the shape.

The paper is organised as follows. In Section II we describe how to specify initial conditions for PBH formation. Section III is devoted to the calculation of the average density profile in the presence of non-Gaussianity. In section IV we discuss average density profile around the threshold for PBH formation, assuming a particular shape of the power spectrum to derive the explicit shape of the density, which is then discussed in Section V. In Section VI we discuss the numerical results obtained with the initial conditions previously derived, and finally in Section VII we give our conclusions.

II Initial conditions for PBH formation

In order to describe the formation of PBHs, we need to consider a region of the expanding Universe with a local non-linear perturbation of the metric which, after re-entering the cosmological horizon, will collapse forming a black hole. Assuming spherical symmetry the perturbation of this region is described by the two following asymptotic forms of the metric

[TABLE]

where the equivalence between the radial and the angular parts gives

[TABLE]

In Eq. (2.1) $a(t)$ is the scale factor while $K(r)$ and $\zeta({\hat{r}})$ are the conserved comoving curvature perturbations on super-Hubble scale, converging to zero at infinity where the universe is taken unperturbed and spatially flat. Combining the two expressions of Eq. (2.2) one gets the explicit transformation between $K(r)$ and $\zeta({\hat{r}})$

[TABLE]

where $\zeta^{\prime}({\hat{r}})$ is the first derivative of $\zeta({\hat{r}})$ with respect to ${\hat{r}}$. In general $K(r)$ and $\zeta({\hat{r}})$ are identified with the average curvature profile.

The metrics given by Eq. (2.1) are asymptotic solutions of the Einstein equations, while the full solution on superhorizon scales, when the curvature profile is conserved being time independent, is obtained using the gradient expansion approximation Shibata:1999zs ; Tomita:1975kj ; Salopek:1990jq ; Polnarev:2006aa . In this regime the energy density profile can be written as a function of the curvature profile harada ; musco as

[TABLE]

Here $H(t)\equiv\dot{a}(t)/a(t)$ is the Hubble parameter while $\omega$ is the coefficient of the equation of state relating the total (isotropic) pressure $p$ to the total energy density $\rho$ as

[TABLE]

where the standard scenario for PBHs assumes a radiation dominated Universe with $\omega=1/3$. The difference between the two Lagrangian coordinates $r$ and ${\hat{r}}$ is related to the particular parameterization of the comovingcoordinate, fixed by the particular form chosen to specify the curvature perturbation into the metric, i.e. $K(r)$ or $\zeta({\hat{r}})$.Here $K^{\prime}(r)$ denotes differentiation with respect to $r$ while $\zeta^{\prime}({\hat{r}})$ and $\nabla^{2}\zeta({\hat{r}})$ denote differentiation with respect to ${\hat{r}}$.

The criterion to distinguish whether a cosmological perturbation is able to form a PBH depends on the amplitude measured at the peak of the compaction function222This was originally introduced in Shibata:1999zs without the factor 2. defined as

[TABLE]

where $R(r,t)$ is the areal radius and $\delta M(r,t)$ is the difference between the Misner-Sharp mass within a sphere of radius $R$ and background mass $M_{b}(r,t)=4\pi\rho_{b}(r,t)R^{3}(r,t)/3$ with the same areal radius but calculated with respect to a spatially flat FRW metric. In the superhorizon regime, applying the gradient expansion approximation, the compaction function is time independent, and is simply related to the curvature profile by

[TABLE]

which, using Eq. (2.3), can be written also in terms of $\zeta({\hat{r}})$. As shown in musco , the length-scale of the perturbation must be identified as the location where the compaction function is reaching its maximum

[TABLE]

which gives

[TABLE]

Given the curvature profile, the value of $r_{m}$ or ${\hat{r}}_{m}$ can be then used to define the small parameter $\epsilon$ of the gradient expansion approximation as

[TABLE]

where $R_{H}$ is the cosmological horizon and $R_{b}(r,t)=a(t)r$ is the areal radius of the background (note that in terms of ${\hat{r}}_{m}$ the curvature profile $\zeta({\hat{r}}_{m})$ is necessary to compute the background value of the areal radius, because of the difference between $r$ and ${\hat{r}}$). The explicit form of the density profile seen in Eq. (2.4), valid for small $\epsilon$, is given by

[TABLE]

where in the first equality the term $[K(r)r^{3}]^{\prime}$ of Eq. (2.4) has been written explicitly, while in the second equality $\nabla^{2}e^{\zeta({\hat{r}})/2}$ has been written in spherical symmetry. Note that to write explicitly these expressions in terms of the small parameter $\epsilon$ one needs to insert $r_{m}$ into the denominator of the term $(1/aH)$ and multiply the radial profile by $r_{m}^{2}$.

Introducing only a perturbation of the energy density field as initial condition corresponds to a combination of growing and decaying mode which would affect the evolution of the cosmological perturbation and the corresponding value of the threshold. As noticed also in Musco:2004ak ; musco , to have at initial time a perturbation behaving like a pure growing mode it is necessary to introduce also a consistent perturbation of the velocity field $U$ and the areal radius $R$, that in gradient expansion have the following form

[TABLE]

where for a pure growing mode one has

[TABLE]

We are now able to define consistently the perturbation amplitude as the mass excess of the energy density within the scale $r_{m}$ measured at horizon crossing time $t_{H}$, defined when $\epsilon=1$ ($aHr_{m}=1$). Although in this regime the gradient expansion approximation is not very accurate and the perturbation amplitude does not represent the exact value of the perturbation at the “real horizon crossing”, this provides a well defined criterion that allows one to compare consistently the amplitude of different perturbations, understanding how the threshold is varying because of the different initial curvature profiles (see musco for more details). The amplitude of the perturbation is given by the excess of mass averaged over a spherical volume of radius $R_{m}$, defined as

[TABLE]

The second equality is obtained by neglecting the higher order terms in $\epsilon$, which allows $R(r,t)$ to be approximated as $R(r,t)\simeq a(t)r$, reducing the first integral over the physical sphere of areal radius $R_{m}$ to an integral over the comoving volume of radius $r_{m}$. inserting the explicit expression of $\delta\rho$ in terms of the curvature profile into (2.16) one obtains

[TABLE]

and a simple calculation seen in musco gives the fundamental relation

[TABLE]

Inserting now Eq. (2.18) into Eq. (2.16) one can easily show that

[TABLE]

which gives an alternative way to compute the length scale $r_{m}$ of the perturbation directly from the energy density profile instead of using the curvature profile.

As shown in musco the threshold of $\delta_{m}$ for PBH formation, called $\delta_{c}$, depends crucially on the shape of the perturbation, which we parameterize in the following through the average density contrast $\overline{\delta\rho}(r)$ measured at horizon crossing $t_{H}$. This quantity inevitably receives non-Gaussian corrections, even though the comoving curvature perturbation is Gaussian. This because the relation (2.4) between the density contrast $\delta\rho$ and the comoving curvature perturbation $\zeta$ is non-linear. In the next section we will calculate the average density contrast $\overline{\delta\rho}(r)$ away from a threshold in the presence of non-Gaussianity.

III The average density profile

To the best of our knowledge, the average profile around a peak has not been calculated for the non-Gaussian case in peak theory. We will therefore resort to threshold statistics, reviewing first the calculation for the Gaussian case in Section III.1, generalizing it for the non-Gaussian case then in Section III.2. Since regions with peak amplitude $\delta\rho_{0}\gg\sigma$ correspond to local maxima to high statistically degree hoffman ; ng3 , our approach should be enough when dealing with PBHs.

III.1 The Gaussian case

Let us first recall how to compute for a Gaussian statistics the average profile $\overline{\delta\rho}(x_{1})$ of the density contrast $\delta\rho({\vec{x}}_{1})$ at a given point ${\vec{x}}_{1}$ from a threshold point located at ${\vec{x}}_{2}$ pw . We define the distance $|{\vec{x}}_{2}-{\vec{x}}_{1}|=r$. Assuming spherical symmetry we can write $\delta\rho({\vec{x}}_{1})=\delta\rho(r)$ and $\delta\rho({\vec{x}}_{2})=\delta\rho_{0}>\nu\sigma$, where $\sigma^{2}=\langle\delta\rho^{2}({\vec{x}})\rangle$ is the variance of the density contrast.

At a distance $r$ from the threshold at the origin, the average $\delta\rho$ is given by

[TABLE]

where

[TABLE]

Since $\delta\rho(r)$ and $\delta\rho_{0}$ are Gaussian variables, one can derive $P(\delta\rho(r),\delta\rho_{0})$ using the covariance matrix

[TABLE]

where

[TABLE]

denotes the two-point correlator. We then deduce that

[TABLE]

being ${\rm Erfc}(x)$ the complementary error function. Using Eq. (3.1), we then obtain

[TABLE]

Finally, exploiting the asymptotic behaviour

[TABLE]

we get that the average $\overline{\delta\rho}$ at a distance $r$ from the threshold with $\nu\gg 1$ is

[TABLE]

As expected, for large values of $\nu$, it coincides with the average profile around peaks obtained using peak theory bbks .

III.2 The non-Gaussian case

In this section we generalise the calculation of the average density profile to the case in which the density contrast is a non-Gaussian field. We start by defining more conveniently the probability

[TABLE]

where $\theta(x)$ is the standard step function, and the conditional probability is therefore

[TABLE]

To proceed, we closely follow the path-integral technique developed in blm ; pbhng . Our starting point is the density contrast $\delta\rho({\vec{x}})$ endowed with a probability distribution $P[\delta\rho({\vec{x}})]$. The corresponding partition function $Z[J]$ in the presence of an external source $J({\vec{x}})$ reads

[TABLE]

The connected $n$-point correlation functions are generated by the functional Taylor expansion of $W[J]=\ln Z[J]$ in powers of the source $J({\vec{x}})$

[TABLE]

At this stage, it is also convenient to normalise the correlators as

[TABLE]

For instance,

[TABLE]

denotes the two-point correlator evaluated at the same point.

With our formalism the average density contrast is easily found as

[TABLE]

To evaluate it, we use the following representation of the $\theta$-function

[TABLE]

and the identity

[TABLE]

This implies

[TABLE]

with

[TABLE]

Using the standard expansion for $\ln Z[J]$

[TABLE]

where

[TABLE]

we find

[TABLE]

Here the prime on the sum reminds us that the sum has to be performed by omitting the terms containing $\phi_{1}^{2}$ and $\phi_{2}^{2}$, and performing the integration over the variables $\phi_{1}$ and $\phi_{2}$

[TABLE]

we find

[TABLE]

It is relevant to point out that only connected correlators appear up to fourth-order, whereas non-connected correlators start to appear at the next fifth order. Then, the connected piece of (3.25) turns out to be

[TABLE]

where $H_{m}(x)$ are the Hermite polynomials.

The one-point non-Gaussian threshold probability for $\nu\gg 1$ is blm ; pbhng

[TABLE]

where $w^{(n)}(0)$ are the normalised $n$-point correlators calculated at the same point. The final expression of the average profile at distance $r$ from the origin and for large thresholds and up to the four-point correlator (recall that $|{\vec{x}}_{2}-{\vec{x}}_{1}|=r$) is given by

[TABLE]

Of course it reduces to the expression (3.8) once the Gaussian limit is taken. The expression (3.28) is the profile we are going to use in the following. However, we will restrict ourselves to a perturbative approach and only include the three-point correlator. Including higher-order terms is unfortunately technically quite demanding. However we will show in Section VI that the modifications of the threshold due to the three-point correlator is quite small because the final non-Gaussian shape is not very different with respect the linear Gaussian one.

IV The average density profile around the threshold for PBH formation

Having calculated the generic expression for the average profile around the threshold, we are now ready to study the problem of PBH formation. As we have already stressed, equation (2.4) is a non-linear relation between the density contrast $\delta\rho$ and the comoving curvature perturbation $\zeta$. This makes the variable $\delta\rho$ non-Gaussian even if $\zeta$ is Gaussian.

First, we will assume that the comoving curvature $\zeta$ is Gaussian so that $\zeta$ does not have an intrinsically second-order component $\zeta_{2}$, but only the linear one, which will call $\zeta_{1}$ (we will promptly extend our computation to the case in which $\zeta$ has some primordial non-Gaussianity). Furthermore, we will restrict ourselves to the case in which we keep only the three-point correlator (see more comments on this later on). Let us also notice that the description of the PBH collapse involves a non-linear relation between the areal radius $R$ and the comoving coordinate $r$, i.e. $R=r\,{\rm exp}\,\zeta$, which introduces further non-linearities.

Expanded at second-order for a linear Gaussian comoving curvature pertubation $\zeta_{1}$, the density contrast is made of a first- and a second-order piece (we assume from now on a radiation phase)

[TABLE]

In Fourier space these relations become (we use here the conventions of b )

[TABLE]

where

[TABLE]

The corresponding bispectrum turns out to be

[TABLE]

where $P_{\zeta_{1}}(k)$ and $P_{\delta\rho_{1}}(k_{2})$ are the power spectrum of the comoving curvature perturbation and of the linear density contrast, respectively. The connected two-point and three-point correlators in coordinate space are given by

[TABLE]

and

[TABLE]

so that

[TABLE]

IV.1 The case of a peaked power spectrum

In order to present analytical formulae we adopt the simplest power spectrum of the comoving curvature perturbation, corresponding to the Dirac-delta case

[TABLE]

for which we have

[TABLE]

The power spectrum (4.7) should be regarded as the limit of zero width of a more physical power spectrum Byrnes .

IV.2 Including non-Gaussianity of the power spectrum

We can generalise these findings to the case in which the comoving curvature perturbation is non-Gaussian ngreview and we standardly parametrise the non-linearities as

[TABLE]

This expression is intended to parametrise the non-linearities which arise at small scales around the scale $k_{*}$ b1 ; b2 ; b3 ; b4 . We are going to consider both positive and negative values of $f_{\text{\tiny NL}}$, keeping in mind that positive values evade observational constraints which place an upper bound on the allowed amplitude of the primordial power spectrum, and allow a cosmologically relevant population of PBHs on the relevant scales333Positive values of the non Gaussianitiy reduce the variance of the curvature perturbation and the bound from the second-order gravitational waves is relaxed, while the contrary is happening for negative values revPBH . The corresponding contribution to the bispectrum is

[TABLE]

Since $B^{\text{\tiny NL}}_{\delta}({\vec{k}}_{1},{\vec{k}}_{2},{\vec{k}}_{3})=\alpha(k_{1})\alpha(k_{2})\alpha(k_{3})B^{\text{\tiny NL}}_{\zeta}({\vec{k}}_{1},{\vec{k}}_{2},{\vec{k}}_{3})$, we have

[TABLE]

and for a peaked power spectrum of the form (4.7) we finally get

[TABLE]

V The average profile including the three-point correlation function

In the previous section we have derived the general form of the three-point correlation function related to the non-linear component of the curvature profile given by (2.4) and to the possible non-Gaussian component of the curvature power spectrum ($f_{\text{\tiny NL}}\neq 0$), considering the particular case of a peaked power spectrum, which allows to get an analytic solution. In the first part of this section we are going to analyze the energy density profile obtained when the three-point correlation function term is taken into account. Although this is just a particular example, it is nevertheless interesting, as a matter of principle, to investigate this case, computing the modification obtained on the threshold $\delta_{c}$ for PBH formation, to get a hint about the general effect of the non-linearities.

In the second part of this section we are going to analyze the energy density profile obtained from the averaged profile of the curvature perturbation $\overline{\zeta}$ of a peaked power spectrum if peak theory is applied to $\zeta$ instead of $\delta\rho$ as was done in haradath . The aim is to make a comparison of the threshold with the profile obtained with threshold statistics, showing that the energy density profile as follows from (2.4), using the averaged curvature profile $\overline{\zeta}$, is very different in general from the mean profile. In other words, the knowledge of $\overline{\zeta}$ does not give a direct way to compute the corresponding threshold. A non-Gaussian method to generalize peak theory, as the one we are using here, is necessary to compute precisely the threshold of PBH formation.

V.1 The averaged density profile from threshold statistics

Considering (3.28) up to the three-point correlation function for the power spectrum given by (4.7) in spherical symmetry and inserting Eqs. (IV.1), (4.9), (IV.2), one obtains the explicit form of the averaged density profile given by

[TABLE]

where $\hat{x}\equiv k_{*}\hat{r}$. Note that in (IV.1) $\xi^{(2)}(r)\propto\sigma^{2}$ while in (4.9) $\xi^{(3)}(r)\propto\sigma^{4}$, therefore the expansion parameter $\nu$ of (3.28) corresponds in the explicit profile given by (5.1) to an expansion around the peak amplitude of the energy density. The functions ${\cal F}_{1}({\hat{x}})$ and ${\cal F}_{2}({\hat{x}})$ are modifications of the profile coming from the three-point correlation function related respectively to the non-linear term of (4.1), and to the non-Gaussianity introduced in (4.10) These two functions read as

[TABLE]

where they have been normalized such that ${\cal F}_{1}(0)={\cal F}_{2}(0)=1$. Note that in the linear limit of a Gaussian density contrast, $\xi^{(3)}(\hat{x})=0$ and the profile is simply reduced to the sync function as it has been obtained in Ref. musco . Using now (4.2) combined with (IV.1), one gets

[TABLE]

which inserted into (5.1) gives

[TABLE]

We see that for $A=0$ the perturbation vanishes ($\delta\rho({\hat{x}})=0$). This can be renormalized with respect the central value as

[TABLE]

where

[TABLE]

Finally, the peak amplitude $\delta\rho_{0}$ of the average energy density profile is related to the amplitude of the peak in the Gaussian approximation $\delta\rho_{0_{G}}$ as

[TABLE]

Apart from the exponential correction, which we will see later at the end of Section VI that can be usually neglected, the profile given by (5.6) is a second-order expansion in terms of the Gaussian amplitude of the peak, consistently with the second order approach we are following.

We are now going to assume $x\simeq\hat{x}$, neglecting the exponential term, because there is no an analytic form of $\zeta(\hat{x})$ corresponding to (5.6), necessary to calculate precisely the value of $\hat{x}_{m}$ and the perturbation of the velocity field given by (2.14). The exponential in (5.8) is only a numerical coefficient that is not going to change the profile of (5.6), changing only the relative value of $\delta\rho_{0_{G}}$ with respect the physical value of the peak given by $\delta\rho_{0}$. The value of the threshold $\delta_{c}$ and the corresponding critical peak amplitude $\delta\rho_{0c}$ are therefore independent from the value of $A$. To simplify the treatment we are therefore going to neglect this exponential term, keeping in mind that the numerical values of $\delta\rho_{0_{G}}$ obtained in the next section should be in general associated to $\delta\rho_{0_{G}}e^{-\mathcal{S}(A)}$, where

[TABLE]

According to this we simplify (5.6) as

[TABLE]

The function ${\cal F}(x)$ is a second-order correction to the profile, measured in powers of $\delta\rho_{0_{G}}$, with respect to the linear Gaussian approximation where ${\cal F}(x)=0$. This is the final form of the profile that will be used to compute numerically in the next section the corresponding value of the threshold $\delta_{c}$ for different values of $f_{\text{\tiny NL}}$. This approximation is consistent with the second-order expansion we have used here to derive the energy density profile. However one should remember that, because the threshold of PBH formation is non-linear, in principle all the non-linear components of the curvature perturbation should be taken into account. The aim of this calculation is to check if the amplitude of the modification given by the three-point correlation function truncating (3.28) at the third order, including also a possible non-Gaussian component of the power spectrum, is small. Only in this case our approach would be consistent.

The input parameter measuring the amplitude of the perturbation is given by $\delta\rho_{0}$, with the corresponding Gaussian value computed with (5.8). The shape of the energy density profile given by (5.9) is characterized by three different functions: ${\rm sync}(x)=\sin x/x$, ${\cal F}_{1}(x)$, ${\cal F}_{2}(x)$, combined together with different coefficients to determine the final shape. In Figure 1 these functions are plotted against $x=k_{*}r$, showing that ${\cal F}_{2}(x)$ is a bit steeper than ${\cal F}_{1}(x)$ which is itself slightly steeper than ${\rm sync}(x)$. Depending on the sign of $f_{\text{\tiny NL}}$, these three functions will combine in different ways and the final non-linear shape given by (5.9) would be steeper or shallower with respect to the Gaussian shape which is described simply by the sync function. We will see later in Section VI how the threshold $\delta_{c}$ for PBH formation is changing with respect to the linear Gaussian case, varying also the value of $f_{\text{\tiny NL}}$.

An analogous calculation gives the profile of the velocity field: inserting (5.9) into (2.14) and assuming $aHr_{m}=1$, we get

[TABLE]

where the functions ${\cal G}_{0}(x),{\cal G}(x)$ are defined as

[TABLE]

The integrals of the functions ${\cal F}_{1,2}(x)$ can be computed analytically

[TABLE]

where ${\cal G}_{0}(0)={\cal G}_{1}(0)={\cal G}_{2}(0)=0$ as one would expect consistently with the boundary condition of the velocity at the centre ($U(0)=0$). We notice that the function ${\cal G}(x)$ is formally modifying the profile of the velocity field with respect to the linear Gaussian case given by ${\cal G}_{0}(x)$, as the function ${\cal F}(x)$ is doing for the energy density profile given by $\sin{x}/x$.

For the numerical implementation of this perturbation, we need to compute the value of $x_{m}$ in terms of the initial input parameters, that is the amplitude measured by the central peak $\delta\rho_{0}$, the peak of the power spectrum $A$ and the non-Gaussian component of the power spectrum measured by $f_{\text{\tiny NL}}$. The integral relation for $x_{m}$ as follows from (5.9) and (2.19) is explicitly written as

[TABLE]

and needs to be solved numerically. When ${\cal F}(x)=0$, which implies that ${\cal G}(x)=0$, one gets $x_{m_{G}}\simeq 2.74$ consistently with musco . Finally we are now able to calculate the averaged amplitude $\delta_{m}$ from the input value of the central energy density peak $\delta\rho_{0}$ inserting (5.9) into (2.18) for $x=x_{m}$.

V.2 The density profile from the averaged curvature profile $\bar{\zeta}$

In the following we are going to derive the energy density profile corresponding to the mean curvature profile $\bar{\zeta}$ obtained from the peaked power spectrum when peak theory is applied to the Gaussian variable $\zeta$ instead of the standard approach using the energy density $\delta\rho$. The two approaches in general are not equivalent because of the non-linear relation of expression (2.4): even though peaks in $\zeta$ correspond to peaks in $\delta\rho$ if they are steep enough haradath ; ng3 , the energy density profile obtained with this from $\bar{\zeta}$ does not correspond to the mean profile of the energy density. The aim here is to compare in the next section the threshold of this profile with the one obtained earlier in (5.9). The mean curvature profile $\bar{\zeta}$ corresponding to a peaked power spectrum is haradath

[TABLE]

and plugged into (2.11) gives

[TABLE]

where $\hat{x}\equiv k_{*}\hat{x}$. The overdensity at the center turns out then to be

[TABLE]

which allows to renormalize (5.13) as

[TABLE]

We can then calculate the scale ${\hat{x}}_{m}$ of the perturbation by solving equation (2.9), which is explicitly written as

[TABLE]

analogous to (5.11) when ${\cal F}(x)=0$ and its solution is ${\hat{x}}_{m}\simeq 2.74$. Because the horizon crossing is calculated in real space when $aHr_{m}=1=aH{\hat{r}}_{m}e^{\zeta({\hat{r}}_{m})}$, it is necessary to renormalize the central peak of the energy density with respect to $x_{m}={\hat{x}}_{m}e^{\zeta({\hat{x}}_{m})}$, that is

[TABLE]

Using now the expression for $\delta_{m}$ given by (2.18), we find that in terms of $\zeta({\hat{x}}_{m})$

[TABLE]

which combined with (5.12) and (5.16) leads to

[TABLE]

Replacing this into (5.17), one can calculate the peak amplitude of the energy density from the averaged perturbation amplitude $\delta_{m}$.

VI Numerical results

The averaged profiles of the density, velocity and curvature profiles analyzed in the previous section have been implemented as initial conditions, using the gradient expansion approach described in Section II to calculate the corresponding threshold of PBH formation with the same code used in Musco:2004ak ; Polnarev:2006aa ; Musco:2008hv ; Musco:2012au ; musco ; ng4 . This has been fully described previously and therefore we give only a very brief outline of it here. It is an explicit Lagrangian hydrodynamics code with the grid designed for calculations in an expanding cosmological background. The basic grid uses logarithmic spacing in a mass-type comoving coordinate, allowing it to reach out to very large radii while giving finer resolution at small radii necessary to have a good resolution of the initial perturbation. The initial data are specified on a space-like slice at constant initial cosmic time $t_{i}$ defined as $a(t_{i})r_{m}=10/H$, ($\epsilon=10^{-1}$), while the outer edge of the grid has been placed at $90R_{m}$, to ensure that there is no causal contact between it and the perturbed region during the time of the calculations. The initial data are evolved using the Misner-Sharp-Hernandez equations so as to generate a second set of initial data on an initial null slice which are then evolved using the Hernandez-Misner equations. During the evolution, the grid is modified with an adaptive mesh refinement scheme (AMR), built on top of the initial logarithmic grid, to provide sufficient resolution to follow black hole formation down to extremely small values of ($\delta-\delta_{c}$). The code has a long history and has been carefully tested in its various forms. Numerically is a second order scheme, using double precision (16 digits after the coma) keeping the numerical error less than $10^{-4}$ . Further information about tests of the code, including a convergence test, could be found in the appendix B of Musco:2008hv .

We are now going to analyze the critical average profiles given by (5.9) showing explicitly the different components that gives rise to the final profile, when $\zeta$ is a Gaussian random variable ($f_{\text{\tiny NL}}=0$) and when a non-Gaussian contribution to the field is also taken into account ($f_{\text{\tiny NL}}\neq 0$), using both positive and negative values of the non-Gaussian parameter. One can write explicitly the different components as

[TABLE]

and we notice that the non-linear components are one order or magnitude higher in terms of $\delta\rho_{0_{G}}/\rho_{b}$, consistent with our perturbative approach. Because $27/2x^{2}_{m_{G}}\simeq 1.80$, if the peak amplitude of the perturbation is small ($\delta\rho_{0_{G}}\ll 1$), then the non-linear components can be neglected and linear theory can be used with good accuracy to calculate the shape of the average density peak, while if $\delta\rho_{0_{G}}\sim 1$ or larger, as it is necessary for PBH formation musco , the non-linear components have the same amplitude of the linear one and one should take them into account.

In principle this is questioning our second-order expansion approach to compute the threshold for PBH formation $\delta_{c}$, suggesting that one should compute all the higher-order terms of (3.28) for an accurate computation, which would be extremely difficult. However, because the shapes of the all three components are similar to each other in the range of the overdensity region (see Figure 1), it will turn out that the final shape is not very different from the linear one, because is a combination of three similar shapes. For this reason, the final critical amplitude of the peak is not very different from the one calculated with the linear approximation as we can see in the right plot of Figure 2 where we are comparing the linear critical average density profile (green line) with the one obtained using the non linear approximation (red line) obtained by the combination of the linear and non-linear components, represented separately in the left plot with the magenta and cyan lines respectively. As we have argued these two components have a comparable amplitude, but the two final linear and non-linear shapes are very similar, and the threshold $\delta_{c}$ of the non linear case is about $1\%$ smaller in the linear case, while the critical amplitude of the peak is about $4\%$ larger in the non linear case with respect the linear one.

The blue line in the right plot of Figure 2 represents the density profile obtained with the alternative approach of using the averaged profile $\overline{\zeta}$ discussed in Section V.2 to compute the profile of the density contrast, which is very different from the one we have obtained with our perturbative approach. In this case the critical amplitude of the peak is significantly smaller ($\delta\rho_{0}=0.716$ in the non-linear case against $\delta\rho_{0}$ = 1.218 of the linear one) with a difference larger than $40\%$. Because the relation to compute the energy density profile from the curvature profile given by (2.4) is non-linear, the mean profile of $\zeta$ does not give the corresponding mean profile of the energy density. Until it will not be clear how this expression should be modified, it is not clear if the application of peak theory in $\zeta$ could be used to compute consistently the correct value of the threshold associated to the mean energy density profile that one needs to calculate the abundance of PBHs, as it has been done in haradath ; b4 .

In Figure 3 we analyze the critical shape of the density contrast for positive values of $f_{\text{\tiny NL}}$ between $1$ and $4$: when $f_{\text{\tiny NL}}=1,2$ (top plots), the three components (linear, non-linear and non-Gaussian) have a similar amplitude while for $f_{\text{\tiny NL}}=2,4$ (bottom plot) the amplitude of the non-Gaussian component is larger than the other two and becomes progressively more and more dominant for increasing values of $f_{\text{\tiny NL}}$. In Figure 4 the same analysis is done for negative values of $f_{\text{\tiny NL}}$ between $-1/2$ and $-2$, where now the non-Gaussian component is negative.

All these pictures show that one can interpret the final critical peak as the combination of 3 different peaks with a similar length-scale, and in the case of positive value of $f_{\text{\tiny NL}}$ the three components are all positive, with the two non linear components being slightly steeper, as one can see from Figure 1. Therefore the critical amplitude of the final peak is progressively increasing with increasing values of $f_{\text{\tiny NL}}$ with respect to the critical shape obtained using a linear approximation as one can see from the left plot of Figure 5 where all the critical shapes for $f_{\text{\tiny NL}}\geq 0$ are plotted as function of $r/r_{m}$. For negative values of $f_{\text{\tiny NL}}$ instead the non-Gaussian component is a negative peak, with opposite sign with respect to the other two components.

In the right plot of Figure 5 we are plotting all the critical shapes for $f_{\text{\tiny NL}}\leq 0$, and we see that the non-Gaussian component is basically compensating the non-linear component for $f_{\text{\tiny NL}}-1/2$, giving rise to almost the same profile of the linear case (in the overdensity region the green and the blue line are indistinguishable), while for more negative values of $f_{\text{\tiny NL}}$ the negative non-Gaussian component becomes more and more important, and an off-centered peak arises when $f_{\text{\tiny NL}}=-2$. This is in the limit of what is possible to be studied with our perturbative approach because the value of the coefficient $\delta\rho_{\text{\tiny G}}$ is obtained from equation (5.8) neglecting the exponential term (see later for a comment about this approximation) which gives a critical value of $f_{\text{\tiny NL}}$

[TABLE]

beyond which the value of $\delta\rho_{\text{\tiny G}}$ becomes an imaginary number and our perturbative approach breaks down (note that replacing the values from Table 1 one obtains indeed $f_{\text{\tiny NL}_{c}}\simeq-2.04$). This is physically saying that for large enough negative values of $f_{\text{\tiny NL}}$, the non-Gaussian component will dominate and a negative peak will arise in the center, which cannot be treated consistently with peak theory because the critical amplitude of the peak is decreasing significantly if the perturbation is not anymore centrally peaked, even if the perturbation has the same value of the threshold $\delta_{c}$. This suggests that a more proper calculation that would take into account also the contribution from the higher order correlators in equation (3.28) might significantly change the behaviour of the peak.

In Table 1 we have summarized all the numerical values of the main quantities characterizing the critical cases we have studied, divided in three parts. The first part of the table gives the value of the critical profiles we have shown in Figure 2 when $f_{\text{\tiny NL}}=0$, while the second and the third part refer instead to the cases of positive and negative values of $f_{\text{\tiny NL}}$ that we have seen separately in Figure 3 and 4, and summarized in the left and right plot of Figure 5. The first two columns of data of the table give the corresponding critical values of the peak amplitude $\delta\rho_{c}$ and of the threshold $\delta_{c}$ respectively, while the third column gives the corresponding value of $x_{m}=k_{*}r_{m}$. The fourth column gives the ratio between $r_{0}$, measuring the edge of the overdensity and the typical perturbation scale $r_{m}$. As seen in musco , this is one of the crucial parameters, together with the peak amplitude $\delta\rho$ and the mass excess $\delta$ which characterize the shape. The fifth column gives the amplitude of the corresponding Gaussian peak $\delta\rho_{\text{\tiny G}}$, which is equal to the critical amplitude of the peak for the linear case when there are no non-linear corrections to the shape, while for the other cases this value is a coefficient weighting the amplitude of the different components of the profile as we have seen in Eqs. (6.1), (6.2) and (6.3). Finally the last two columns give the percentage fractional correction of the critical value of the peak and of the average threshold with respect to the linear case.

In the left plot of Figure 6 we are plotting the values of the threshold $\delta_{c}$ obtained with the non linear corrections, as function of $f_{\text{\tiny NL}}$: $\delta_{c}$ is monotonically decreasing for increasing values of $f_{\text{\tiny NL}}$, converging to $\delta_{c}\simeq 0.5$ for large positive values of $f_{\text{\tiny NL}}$. The corresponding amplitude of the critical peak is increasing and converging to a maximum value $\delta\rho_{c}\simeq 1.35$, about $10\%$ larger than the amplitude of the critical peak obtained with the linear approximation. On the contrary, for negative values of $f_{\text{\tiny NL}}$, the threshold $\delta_{c}$ is increasing for $f_{\text{\tiny NL}}$ becoming more and more negative, while the amplitude of the critical peak is decreasing. This inverse behaviour of the critical peak amplitude increasing against the corresponding value of the threshold is due to the different amplitude of the pressure gradients modifying the shape during the collapse musco . In particular, for $f_{\text{\tiny NL}}=-2$, we find that the density contrast is not anymore centrally peaked, and an off-centered peak arises. Looking at Table 1 we see that for negative values of $f_{\text{\tiny NL}}$, the critical amplitude of the density contrast is varying more significantly with respect to the variation obtained for positive values.

Looking at the right plot of Figure 6, where we plot the relative change of $\delta_{c}$ as function of $f_{\text{\tiny NL}}$, we see that this is much more under control compared to the critical amplitude of the peak: for the positive values we have analyzed, the variation is less than $3\%$ and the convergent behaviour suggests that it will not increase significantly more than this limit, while for negative values of $f_{\text{\tiny NL}}$ the relative change becomes more and more significant, tending to diverge. For values of $f_{\text{\tiny NL}}$ down to $-3/2$ the variation is still of few percent, which is consistent with our perturbative approach, but for $f_{\text{\tiny NL}}=2$, when the peak becomes off-centered the variation is more significant (about of $5\%$) which is in the limit of what could be considered consistent with a perturbative approach. We argue therefore that it is reasonable that the inclusion of higher-order terms in the calculation of the average shape will not change significantly our results for $f_{\text{\tiny NL}}\geq 3/2$, while for more negative values a significant change could be possible.

Before concluding, we would like to comment about the approximation done in Section V.1 where we have neglected the exponential term in equation (5.8). Using $\delta\rho_{0_{G}}\sim 1$, as one can see from table 1, one obtains

[TABLE]

This shows that for values of $A$ large enough, such that the abundance of PBHs is able to explain a significant amount of dark matter, the impact of the exponential term in (5.8) is quite small and in first approximation can be neglected as we have done for simplicity. However, as we have already explained in the previous section, even for values of $A\ll 1$ such that this term should be taken into account, the correction of the value of $\delta\rho_{0_{G}}$ with respect to the peak amplitude $\delta\rho$, is affecting significantly only the relative coefficients of the profiles ${\rm Sync}(x),{\cal F}_{1}(x),{\cal F}_{2}(x)$, and not the final profile which is a linear combination of similar shapes seen in Figure 1.

VII Conclusions

In this paper we have considered the change in the critical threshold induced by the ineludible and intrinsic non-Gaussianity originated by the non-linear relation between the overdensity and the curvature perturbation in the case in which the latter is Gaussian. We have also extended our results by assuming a non-Gaussian curvature perturbation.

The impact of non-Gaussianity in the density contrast threshold, even if the curvature perturbation is Gaussian, alters the PBH abundance. Denoting the change in the threshold by $\Delta\delta_{c}$ and using the expression (1.1), we can estimate the contribution to the non-linear abundance from the shift in the threshold to be

[TABLE]

Interesting abundances of PBHs are obtained for $\delta_{c}/\sigma={\cal O}(6\div 8)$. On the other hand, our results show that the relative change $\Delta\delta_{c}/\delta_{c}$ is at the percent level (see Figure 6), leading to a change in the abundance, with respect to the result assuming the Gaussian critical threshold, which is $\lesssim{\cal O}(10^{2})$. In any case it is smaller than other uncertainties present in the estimate (e.g. the use of statistics), and therefore is not cosmologically significant.

Our results are based on a perturbative approach which is restricted to the second-order. Even though limited, we argue that our results should be rather robust against the addition of higher-order terms for f_{\text{\tiny NL}}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt\sim}\hss}\raise 1.0pt\hbox{>}}-3/2, because the critical threshold is only sensitive to the final shape of the profile of the overdensity, that determines the value of the threshold $\delta_{c}$, and not to the amplitude of the different non linear components. This will be true if the higher order of non-Gaussianity will behave as the $f_{\text{\tiny NL}}$-term, not altering significantly the final shape.

We have also seen that the relative change of the threshold $\delta_{c}$ is much more robust than the relative change of the critical amplitude of the peak $\delta\rho_{c}$, changing up to $10\%$ for positive values of $f_{\text{\tiny NL}}$, and more than $20\%$ for negative values, because the critical amplitude of the peak is much more sensitive to the local features of the shape than the threshold $\delta_{c}$ which is an averaged quantity. We expect that going beyond the perturbative approach, at least for $f_{\text{\tiny NL}}\geq-3/2$, will not alter significantly the threshold. This suggest also that the threshold $\delta_{c}$ would allow to compute the abundance of PBHs with less uncertainties than using the local critical amplitude of the peak, as was also pointed out in Young:2019osy .

It will be interesting in the future to understand to which extent the conclusions we have reached here are valid also for a more general shape of the peak of the cosmological power spectrum.

Acknowledgments

We thank V. Atal, N. Bellomo, Chris Byrnes, V. De Luca, G. Franciolini, J. Garriga, C. Germani, J. Miller, L. Verde and S. Young for useful discussions. A.R. is supported by the Swiss National Science Foundation (SNSF), project The Non-Gaussian Universe and Cosmological Symmetries, project number: 200020-178787. I.M. is supported by the Unidad de Excelencia María de Maeztu Grant No. MDM-2014-0369. A.K. is supported by the GSRT under EDEIL/67108600.