# Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation

Sam Young, Ilia Musco, Christian T. Byrnes

► To cite this version:

Sam Young, Ilia Musco, Christian T. Byrnes. Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation. Journal of Cosmology and Astroparticle Physics, 2019, 11, pp.012. 10.1088/1475-7516/2019/11/012. hal-02116545

HAL Id: hal-02116545

https://hal.science/hal-02116545v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN ACCESS

# Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation

To cite this article: Sam Young *et al* JCAP11(2019)012

View the [article online](#) for updates and enhancements.

You may also like

- [Search for Gamma-Ray Emission from Local Primordial Black Holes with the Fermi Large Area Telescope](#)
M. Ackermann, W. B. Atwood, L. Baldini et al.

- [The Effect of Primordial Black Holes and Streaming Motions on Structure Formation](#)
F. Atrio-Barandela

- [How Do Primordial Black Holes Change the Halo Mass Function and Structure?](#)
Saiyang Zhang, Volker Bromm et al.
# Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation

Sam Young,ª Ilia Muscoᵇ.c and Christian T. Byrnesᵈ

E-mail: syoung@mpa-garching.mpg.de, iliamusco@icc.ub.edu, C.Byrnes@sussex.ac.uk

Received May 10, 2019

Revised August 1, 2019

Accepted October 14, 2019

Published November 13, 2019

**Abstract.** The formation and abundance of primordial black holes (PBHs) arising from the curvature perturbation $\zeta$ is studied. The non-linear relation between $\zeta$ and the density contrast $\delta$ means that, even when $\zeta$ has an exactly Gaussian distribution, significant non-Gaussianities affecting PBH formation must be considered. Numerical simulations are used to investigate the critical value and the mass of PBHs which form, and peaks theory is used to calculate the mass fraction of the universe collapsing to form PBHs at the time of formation. A formalism to calculate the total present day PBH abundance and mass function is also derived. It is found that the abundance of PBHs is very sensitive to the non-linear effects, and that the power spectrum $P_{\zeta}$ must be a factor of $O(2)$ larger to produce the same number of PBHs as if using the linear relation between $\zeta$ and $\delta$ (where the exact value depends on the critical value for a region to collapse and form a PBH). This also means that the derived constraints on the small-scale power spectrum from constraints on the abundance of PBHs are weaker by the same factor.

**Keywords:** cosmological perturbation theory, primordial black holes

**ArXiv ePrint:** 1904.00984

<figure><img src="image_3.png" /></figure>

<footnote>ªMax Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, 85748 Garching bei Muenchen, Germany</footnote>

<footnote>ᵇInstitut de Ciències del Cosmos, Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain</footnote>

<footnote>ᶜLaboratoire Univers et Théories, UMR 8102 CNRS, Observatoire de Paris, Université Paris Diderot, 5 Place Jules Janssen, F-92190 Meudon, France</footnote>

<footnote>ᵈDepartment of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, U.K.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Cosmological perturbations in the super horizon regime</td><td>3</td></tr><tr><td></td><td>2.1 Gradient expansion approach</td><td>4</td></tr><tr><td></td><td>2.2 Initial conditions</td><td>5</td></tr><tr><td>3</td><td>Criterion for collapse</td><td>7</td></tr><tr><td>4</td><td>Numerical results of PBH formation</td><td>9</td></tr><tr><td></td><td>4.1 Numerical scheme</td><td>9</td></tr><tr><td></td><td>4.2 Threshold, scaling law and mass spectrum</td><td>10</td></tr><tr><td>5</td><td>Calculation of PBH abundance</td><td>11</td></tr><tr><td>6</td><td>Summary</td><td>16</td></tr><tr><td>A</td><td>Gaussianity and variance of $-\frac{4}{3}r_m\zeta'(r_m)$</td><td>18</td></tr><tr><td>B</td><td>Correspondence of large peaks</td><td>19</td></tr></table>

## 1 Introduction

Primordial black holes (PBHs) could be formed from the gravitational collapse of large curvature perturbations created during cosmological inflation shortly after re-entering the cosmological horizon at early times [1–3]. If a density perturbation is above a threshold $\delta_c$, then an apparent horizon will form during the collapse, otherwise it will quickly disperse into the surrounding local environment. The mass of the resulting PBH is strongly related to the scale and amplitude of the perturbation from which it formed, with more massive black holes forming from larger-scale perturbations which enter the horizon at a later time. PBHs can theoretically form with any mass, and can provide a natural explanation for any observed black holes with masses which are not easily explained by the standard picture of black hole formation from collapsing stars. PBHs which formed with a mass below $10^{15}g$ would have evaporated by today (ignoring the possible accretion after formation), but more massive PBHs would persist into the present epoch.

PBHs still represent a viable dark matter candidate, although there are numerous constraints on the abundance of PBHs of varying masses (see [4–6] for recent discussions of the constraints for a broad mass spectrum), and clusters of PBHs could explain the early formation of super-massive black holes found in the centres of galaxies. In recent years, there have been several interesting observations which may hint towards the existence of PBHs [7]. For a review of PBH formation and constraints, see [8, 9].

There have been many attempts to detect them by their indirect effects on the universe. Ignoring the possible observations described above, no positive detection has been made. However, the non-detection of PBHs constrains their abundance. The abundance of PBHs is typically stated in terms of $\beta$, the energy fraction of the universe which goes into PBHs at the time of formation. The abundance of small PBHs ($< 10^{15}g$) in the early universe
which would have evaporated by today can be constrained by looking for the effects of the radiation from their evaporation, whilst the abundance of more massive PBHs (> $10^{15}g$) can be constrained by their gravitational effects. Because PBHs of different masses form from different scale perturbations, the constraints on different mass PBHs can be used to place a constraint on different scales of the primordial power spectrum in the early universe — though these constraints are sensitive to primordial non-Gaussianity in the early universe. Because PBHs form on scales much smaller than those observable in the cosmic microwave background (CMB) or large-scale structure of the universe (LSS), they therefore place the only available constraints on the small-scale power spectrum — and can be used to constrain models of inflation. In order for an interesting number of PBHs to form however, the power spectrum must become orders of magnitude larger on small scales than is observed in the CMB or LSS ($\mathcal{O}(10^{-2})$ compared to $\mathcal{O}(10^{-9})$). Therefore the derived constraints on the power spectrum are much weaker, but they span a much larger range of scales, including scales which cannot be probed by any other method [10].

There are many different models for cosmological inflation which would predict a large number of PBHs to form in the early universe. For example, the running mass model [11], axion inflation [12, 13], or a waterfall transition during hybrid inflation [14–16], a quartic action during inflation or a variable sound speed [17], amongst many others. A metric perturbation in the form of the curvature perturbation $\zeta$ is typically used to study cosmological perturbations generated with the different models and to predict their observable consequences.<fnref n="1" /> The curvature perturbation $\zeta$ appears in the metric in the comoving uniform-density gauge as

$$ds^2 = -d\eta^2 + a^2(t)e^{2\zeta}dX^2, \qquad (1.1)$$

where $\eta$ is the conformal time, $a(t)$ is the scale factor and $X$ represents the three comoving spatial coordinates. In order to translate the constraints on PBH abundance into constraints on models of inflation (or alternatively to predict PBH abundances from a given model) it is desirable to relate the primordial curvature perturbation power spectrum $P_\zeta$ to the abundance of PBHs of different masses.

The formation of PBHs from a non linear metric perturbation was initially studied by Shibata and Sasaki [18], which was then used to derive a relation between the abundance of PBHs and the power spectrum $P_\zeta$ [19]. Around the same time Niemeyer and Jedamzik performed a numerical study of PBH formation using instead an initial perturbation of the energy density [20, 21]. For a long time the abundance of PBHs was calculated assuming that regions where the curvature perturbation $\zeta$ was above a critical value $\zeta_c$ of order unity. However, it has since been shown that the curvature perturbation $\zeta$ is not a suitable parameter to use to determine whether a region will form a PBH or not, due to the effect of super-horizon modes on the calculation, and that the density contrast should be used instead [22]. The effect of super-horizon modes on PBH formation is discussed in detail in [23]. Nonetheless, in the following years it has been typical to use the curvature perturbation directly to calculate the abundance — which is valid in the case that an approximation is being used (as described in [22]) or in the case that a narrow peak in the power spectrum is being considered (so that large perturbations only exist at one scale). Papers which have used the density contrast $\delta\rho/\rho_b$ rather than the curvature perturbation for the calculation of the abundance have since

<footnote n="1">¹Often $\mathcal{R}$ is used instead to denote the curvature perturbation, whilst $\zeta$ is used to describe the curvature perturbation only on a uniform density slicing.</footnote>
used a linear relation between the two parameters (as in the recent paper [24] for example),

$$ \frac{\delta\rho}{\rho_b} = -\frac{2(1+\omega)}{5+3\omega} \left(\frac{1}{aH}\right)^2 \nabla^2\zeta, \qquad (1.2) $$

where $\omega = 1/3$ is the equation of state parameter during the radiation dominated epoch of the early universe, $\rho_b$ is the background density and $(aH)^{-1}$ is the comoving cosmological horizon scale. However, this expression ignores the non-linear relation between the curvature and the energy density profile. Starting for the first time from simulations of PBH formation arising from perturbations in the curvature perturbation $\zeta$, the aims of this paper are to investigate how the fully non-linear relation between $\zeta$ and $\delta\rho/\rho_b$ affects the calculation of the abundance of primordial black holes, and to derive the most accurate relation to date between the primordial curvature perturbation power spectrum $P_\zeta$ and the abundance of PBHs at formation $\beta$.

We note that reference [25] introduced a new calculation which included the non-linear effect. The abundance of PBHs was calculated using a method based on estimating the critical heights of peaks in the curvature perturbations $\zeta$, and then calculating the abundance of PBHs by utilising the Gaussianity of $\zeta$. Using the critical height of peaks in $\zeta$ is only valid for a narrow power spectrum, such that perturbations exist only on one scale (as they note in section 4), whereas the method presented here can be applied to power spectra of any width. Reference [25] stated that the abundance of PBHs they calculated was greatly increased compared to some previous calculations, whereas it is shown here that considering the non-linear relation decreases the abundance compared to the linear calculation. The apparent contradiction in conclusions can be attributed to differences in what is considered as the “standard calculation” between that paper and ours. We compare our results using the same method from using the linear or non-linear relations between $\delta$ and $\zeta$, whereas [25] compares very different methods.

The paper is organised as follows: in section 2 we will discuss the set up of the initial conditions in the density contrast $\delta\rho/\rho_b$ arising from an initial curvature perturbation $\zeta$. In section 3 we will discuss the criteria which should be used to determine whether an initial perturbation will eventually collapse to form a PBH. In section 4 we will discuss the simulation procedure used and the numerical results obtained from the simulations. Finally, in section 5 we will show how the abundance of PBHs $\beta$ can be obtained from the curvature perturbation power spectrum. Our findings are summarised in section 6, leaving some details of our calculations to an appendix.

## 2 Cosmological perturbations in the super horizon regime

In this section we will first describe the general relation between the curvature perturbation $\zeta$ and the density contrast $\delta\rho/\rho_b$ before analysing a specific parametrization of $\zeta$ that allows us to vary the profile of $\delta\rho/\rho_b$. This allows us, with the help of numerical simulations (see section 4), to span almost all the possible range of values of $\delta_c$. Throughout this paper, we assume that perturbations large enough to form PBHs are spherically symmetric. This is justified because such peaks must be extremely rare [26], and the perturbation profile is therefore defined using only a radial coordinate $r$.

Note that the curvature perturbation $\zeta$ in the literature is typically defined on a uniform density slicing, whilst the density contrast $\delta\rho/\rho_b$ is defined on a constant cosmic time slicing. However, in the super-horizon regime described in the following section the difference between these two gauges is a higher-order correction which can be neglected (see [27] and the references therein).
## 2.1 Gradient expansion approach

In the super-horizon regime perturbations have a length scale much larger than the cosmological Hubble horizon. In this regime it is possible to have an analytic treatment, usually called the gradient expansion or long-wavelength approximation [18, 28], based on expanding the exact solution in a power series of a small parameter ($\epsilon \ll 1$) that is conveniently identified with the ratio between the Hubble radius $1/H(t)$ and the length scale $L$ characterising the perturbation

$$ \epsilon(t) \equiv \frac{1}{H(t)L}, \qquad (2.1) $$

where a particular value of $\epsilon$ corresponds to a particular value of the time $t$.

When $\epsilon \ll 1$ the curvature profile $\zeta(r)$ is conserved (time independent) and each spatial gradient is multiplied by $\epsilon$, expanding the equations in a power series in $\epsilon$ up to the first non-zero order. Putting $\epsilon = 1$, which defines the horizon crossing time, one obtains the spatial dependence of the different matter/geometrical variables, related to the right/left hand side of the Einstein equations, in terms of the conserved curvature profile $\zeta(r)$. The curvature profile represents the fundamental source of the perturbation, embedded into the metric (1.1) from the asymptotic limit, $t \to 0$. Although this approach reproduces the time evolution of linear perturbation theory when $\epsilon \ll 1$, it also allows one to consider non-linear curvature perturbations if the spacetime is sufficiently smooth on scales greater than $L$ [27]. This is equivalent to pressure gradients being a higher-order correction in $\epsilon$, which corresponds to a self-similar growth of the perturbation, conserving the spatial profile.

In the gradient expansion the non-linear relation between the density contrast $\delta\rho/\rho_b$ defined on a comoving uniform-cosmic time slicing and the curvature perturbation $\zeta$ is given by [23, 25, 29], which is exact up to the first non-zero terms in $\epsilon$

$$ \frac{\delta\rho}{\rho_b}(r,t) = -\frac{4(1+\omega)}{5+3\omega} \left(\frac{1}{aH}\right)^2 e^{-5\zeta(r)/2} \nabla^2 e^{\zeta(r)/2}, \qquad (2.2) $$

where $\omega \equiv p/\rho$ is the equation of state parameter, and

$$ \nabla^2 = \frac{d^2}{dr^2} + \frac{2}{r} \frac{d}{dr} $$

is the Laplacian written assuming spherical symmetry. There are two non-linear effects contained within equation (2.2): the exponential term $\exp(-2\zeta(r))$ and the quadratic term of the first derivative $(\zeta'(r))^2$, where the prime denotes a derivative with respect to the radial coordinate $r$. Because PBHs form from large perturbations, the effect of the non-linear components is comparable with the linear term and should not be neglected.

As explained in [29], whether a cosmological perturbation is able to form a PBH depends on the amplitude measured at the peak of the compaction function defined as

$$ C(r,t) \equiv 2 \frac{M(r,t) - M_b(r,t)}{R(r,t)}, \qquad (2.3) $$

where $R(r,t)$ is the areal radius, $M(r,t)$ is the Misner-Sharp mass within a sphere of the radius $R$ and $M_b(r,t) = 4\pi\rho_b(r,t)R^3(r,t)/3$ is the background mass within the same areal radius calculated with respect to a Friedmann-Lemaître-Robertson-Walker (FLRW) universe.
In the superhorizon regime, applying the gradient expansion approximation, the compaction function is conserved and is related to $\zeta(r)$ as [23, 25, 29]

$$C(r) = -f(w)r\zeta'(r)[2 + r\zeta'(r)], \quad f(\omega) = \frac{3(1+\omega)}{5+3\omega}. \qquad (2.4)$$

We can then compute the length-scale of the perturbation, identified as the location $r_m$ where the compaction function is maximized ($C'(r_m) = 0$), which gives

$$\zeta'(r_m) + r_m \zeta''(r_m) = 0. \qquad (2.5)$$

Measuring the curvature perturbation with $\zeta(r)$ introduces an intrinsic rescaling of the comoving coordinate with respect to the background FLRW solution, because the exponential term appearing in the metric (1.1) introduces a local perturbation of the scale factor which depends on the local value of the curvature. This implies that the horizon crossing time $t_H$ is defined in real space when

$$a(t_H)H(t_H)r_m e^{\zeta(r_m)} = 1, \qquad (2.6)$$

and therefore according to the definition of $\epsilon$ given above, the physical length scale of the perturbation, to be called $R_m$ from here onwards, is given by

$$R_m(r, t) = a(t)r_m e^{\zeta(r_m)}. \qquad (2.7)$$

The perturbation amplitude can be measured as the mass excess of the energy density within the scale $R_m$, measured at the horizon crossing time $t_H$. Although in this regime the gradient expansion approximation is not very accurate, this represent a well defined criterion that allows a consistent comparison between the amplitude of different perturbations (see [29] for more details). Computing the mass excess as the integral of the density contrast averaged over the background volume $V_{R_m} = \frac{4\pi}{3}R_m^3$, the amplitude of the perturbation is given by

$$\delta(r_m, t_H) \equiv \frac{3}{R_m^3} \int_0^{R_m} \frac{\delta\rho}{\rho_b} R^2 dR = \frac{3}{(r_m e^{\zeta(r_m)})^3} \int_0^{r_m} \frac{\delta\rho}{\rho_b}(r, t_H)(re^{\zeta(r)})^2(re^{\zeta(r)})' dr, \qquad (2.8)$$

and using the explicit expression of $\delta\rho(r, t)/\rho_b$ in terms of the curvature profile seen in (2.2), we get $\delta_m \equiv \delta(r_m, t_H) = C(r_m)$, which satisfies the fundamental relation [29]

$$\delta_m = 3 \frac{\delta \rho}{\rho_b}(r_m, t_H) \qquad (2.9)$$

for any curvature profile, $\zeta(r)$.

## 2.2 Initial conditions

We will now study the main feature of the density profile when an explicit parameterization of the curvature profile $\zeta(r)$ is specified as

$$\zeta_{\alpha}(r) = \mathcal{A} \exp \left[ - \left( \frac{r}{r_m} \right)^{2\alpha} \right], \qquad (2.10)$$

where $\mathcal{A}$ and $r_m$ respectively denote the amplitude and the scale of the perturbation. Inserting this into (2.2), the corresponding density profile is given by

$$\frac{\delta \rho_{\alpha}}{\rho_{b}} = \left(\frac{1}{aH}\right)^{2} \frac{4}{3} f(\omega) \alpha \left(\frac{r}{r_{m}}\right)^{2\alpha} \left[ (2\alpha + 1) - \alpha \left(\frac{r}{r_{m}}\right)^{2\alpha} (2 + \zeta_{\alpha}(r)) \right] \frac{\zeta_{\alpha}(r)}{r^{2} e^{2\zeta_{\alpha}(r)}}, \qquad (2.11)$$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The left plots shows the critical curvature profiles given by (2.10) for $\alpha = 1$ ($A_c \simeq 0.8$ and $\delta_c \simeq 0.55$), while the right plot shows the corresponding $C(r)$ profile. The upper limit $\delta_{c,max} = 2/3$ represents the maximum theoretical upper limit of $\delta_c$.</figcaption></figure>

and then inserting (2.10) into (2.4) and (2.5), we can calculate the corresponding perturbation amplitude

$$ -r_m \zeta'_{\alpha}(r_m) = \frac{2A\alpha}{e} \Rightarrow \delta_m = 4f(\omega) \frac{A\alpha}{e} \left(1 - \frac{A\alpha}{e}\right), \quad (2.12) $$

which gives the value of $A$ in terms of the averaged amplitude $\delta_m$

$$ A = \frac{e}{2\alpha} \left( 1 - \sqrt{1 - \frac{\delta_m}{f(\omega)}} \right). \quad (2.13) $$

This behaviour of $\delta_m$ in terms of $A$ also shows that there is a maximum value of $\delta_m = 2/3$ corresponding to $A = e\alpha/2$, that represents the transition between PBHs of type I and type II, after which formation of PBHs cannot be studied in terms of $\delta_m$ but only in terms of $\zeta$ [30].

Using numerical simulations we have calculated the critical values for PBHs using initial condition in terms of $\zeta(r)$ given by (2.10), finding that the value of $\delta_c$ is varying in terms of $\alpha$ (see section 4 for more details). In the particular case of $\alpha = 1$ we get $\delta_c \simeq 0.55$ which corresponds to a value of $A \simeq 0.80$. In the left frame of figure 1 we have plotted the critical profile of $\zeta(r)$ as function of $r/r_m$ identifying the critical peak $A_c$ of the profile. In this plot we are also showing the value of $A$ corresponding to the maximum limit of $\delta_{c,max} = 2/3$. In the right hand plot of figure 1 we show the corresponding compaction function $C(r)$ with the peak amplitude of $C(r_m)$ being equal to $\delta_c$.

In figure 2 we plot the $\zeta$-profiles (plotted in the left panel) and the corresponding density contrast (right panel) measured at horizon crossing, defined by (2.6), for the threshold values $\delta_c$ associated to each shape. Although the $\zeta$-profiles given by (2.10) are always centrally peaked, the energy density profile is centrally peaked only for $\alpha \le 1$. In particular, only the case of $\alpha = 1$, plotted with a dashed line, is smooth at $r = 0$, giving a finite value of the peak of the density contrast, while for smaller values of $\alpha$ the peak is diverging. Nevertheless the amplitude of the perturbation, measured by the averaged value $\delta_m$ always remains finite.
The right panel of figure 2 also shows that for $\alpha > 1$ the density contrast is off-centred, with an increasing value of the peak for larger values of $\alpha$. One can see therefore that a $\zeta$-profile with a centrally finite peak does not always corresponds to the same type of peaks in the density contrast, because of the non-linear expression given by (2.2) and the correspondence among the peaks is guaranteed only at the linear level. In general, the correspondence between the peaks in $\zeta$ and the peaks in $\delta\rho/\rho_b$ requires the assumption that $\zeta'(r) = 0$ at $r = 0$, such that in the centre the only non-linear effect is given by the exponential term $\exp(-2\zeta(r))$ which reduces the amplitude.

A second issue already mentioned is that finite peaks of $\zeta$ do not always correspond to finite values of the peak of the density contrast, which happens here for $\alpha < 1$. This is because the $\zeta$-profiles for $\alpha < 1$ are not smooth in the centre (they are not infinitely differentiable), and the term $\zeta'(r)/r$ diverges in the limit $r \to 0$. Such peaks are of course unphysical and this divergence can be removed with a transfer function or a smoothing function which removes the small scale power, but there is lack of knowledge in the literature about which is the correct form of the non-linear transfer function to be used for a radiation fluid. Note that the transfer function should always be taken into account, because strictly speaking the curvature is exactly conserved only for $\epsilon \to 0$. Because in practice a finite value of $\epsilon$, corresponding to a finite initial time $t_i$, needs to be chosen, the effects of the pressure gradients within a region of the size of the initial sound horizon, which in radiation is $R_s = (\sqrt{3}H)^{-1}$, are not completely negligible even at the initial time. For spiky shapes (which have $\zeta'(r)/r \neq 0$ in the centre), like those obtained from (2.10) with $\alpha < 1$ and the effect of the transfer function might significantly change the amplitude of the peak, while the value of the averaged $\delta_m$ hardly changes.

Because for every $\zeta$-peak with a finite amplitude there is always a peak of the compaction function $C$ evaluated at $r_m$, with finite amplitude equal to $\delta_m$, we have used the averaged amplitude $\delta_m$ to calculate the abundance of PBHs (see section 5), using the linear transfer function to correct the value of the peak of the density contrast at $r = 0$, leaving a study of the effects of the non-linear transfer function to future work.

### 3 Criterion for collapse

In order for a perturbation to collapse into a PBH, the density must exceed some critical threshold. The original work by Carr [3], using a Jeans length argument, provided an order of magnitude estimate for the threshold $\delta_c \sim \omega = 1/3$ at horizon crossing for a radiation fluid. Since then, there has been extensive work to determine the collapse threshold [18, 21, 23, 29, 31–34], as well as discussions about the appropriate parameter to use to determine whether a perturbation will collapse [19, 22, 24, 25]. The collapse threshold is typically obtained from simulating the evolution of a perturbation as it reenters the cosmological horizon, although analytic attempts have been made, neglecting the effects of pressure gradients [35].

As discussed previously, in order to determine a clear criterion to distinguish which perturbations are able to form a PBH, the density contrast $\delta\rho/\rho_b$ should be used rather than a metric perturbation such as the curvature perturbation $\zeta$. There has been much ambiguity in the literature about how this critical amplitude is calculated and used (especially between the different communities of relativists modelling PBH formation and cosmologists calculating the abundance of PBHs). It is the aim of this section of the paper to discuss how this should be defined. Spherical symmetry is typically assumed when modelling PBH formation, again justified by the fact that such peaks are large and rare [26] — although
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The left plot shows the critical profiles $\zeta_{\alpha}(r)$ given by (2.11) plotted against $r/r_m$ for $\alpha = 0.5, 0.75, 1.0, 1.25, 1.50, 2.0$ where for increasing values of $\alpha$ the central value $A_c$ is decreasing. The right plot shows the corresponding critical profiles of $\delta\rho_{\alpha}(r)/\rho_b$ given by equation (2.11) plotted against $r/r_m$: for $\alpha < 1$ the profile has a spiky shape, with increasing steepness for decreasing value of $\alpha$; for $\alpha = 1$ the profiles is centrally peaked with $(\delta\rho')'/\rho_b \to 0$ when $r \to 0$ while for $\alpha > 1$ the profile has an increasing off-peak. In both plots the dashed line corresponds to the smooth centrally peaked case with $\alpha = 1$.</figcaption></figure>

non-spherical symmetry have been considered [36]. In this section, we will discuss several ambiguities within the literature over how the criteria for collapse should be defined:<fnref n="2" />

* *The time at which PBH abundance should be calculated.* The threshold for collapse is normally stated in terms of the time-independent component of the density contrast, during the linear regime whilst a perturbation is super-horizon [29]. In the linear regime, the density contrast grows proportionally to the parameter $\epsilon$ (equation (2.1)), which is the ratio of the perturbation scale to the horizon scale at a given scale. Taking the time-independent component is therefore equivalent to setting this parameter to unity, which has resulted in many papers treating this as the value of the density contrast at horizon crossing. Ideally, the abundance of PBHs should be calculated by considering the perturbations on super-horizon scales, long before they reenter the horizon.

* *Should the peak value of the density contrast be used, or the smoothed density contrast?* The peak value at the centre of a density perturbation was used in a recent paper [24] to determine the abundance of PBHs. This is valid when the distribution is already smooth on scales smaller than the scale being considered (as was considered in that paper), which is generally not the case unless a power spectrum with a very narrow peak is being considered. In order to investigate PBH formation over a wider range of scales, it is necessary to use a smoothing function. We also consider the fact that (for perturbations of arbitrary profile shapes), the threshold value for collapse of the central peak varies from 2/3 to infinity, whilst the critical value for collapse of the

<footnote n="2">²A full description and investigation of these considerations is beyond the scope of this paper, although is discussed further in [37].</footnote>
top-hat smoothed density contrast varies from 0.41 to 2/3 — a much smaller range of values. It was also discussed in section 2 that, for certain $\zeta$-profiles, the peak in the density contrast may be off-centred<fnref n="3" /> (when $\alpha > 1$, meaning the central value is smaller than the peak) or infinite — a problem which is avoided by using the smoothed value (see appendix B).

* The choice of window function has a significant effect on the calculated abundance of PBHs (as was discussed in [38]). The threshold value for collapse is typically stated in terms of the volume-averaged density contrast (as for example in [29, 31]), which corresponds to a top-hat window function — suggesting a top-hat function should be used. However, in the super-horizon regime, the smoothing function decreases as $k^2$ but perturbations grow as $k^2$ — meaning the top hat window function is typically not efficient enough at removing small-scale perturbations. For this reason a Gaussian window function is often used in the literature. This however has the drawback of changing the perturbation shape, introducing an unphysical change in the value of the threshold $\delta_c$. For this reason, we will follow the standard approach of using a top-hat smoothing function in this paper, whilst treating perturbations as if they are still linear at horizon entry, and employ the linear transfer function for sub-horizon perturbations to reduce the effects of small-scale perturbations (although note that the linear transfer function might not very accurate for the large amplitude perturbations required to generate PBHs).

* The scale of a perturbation $r_m$ is best stated in terms of the radius at which the compactness function $C(r)$ is maximised (see section 2). This is different to the previously used definition $r_0$ [31], where the scale of the perturbation was defined as the radius at which the density contrast becomes negative. As was shown in [29] computing the density contrast at $r_0$ does not give a sensible parameter to compare different shapes. The averaged value of the density contrast evaluated at $r_m$ is characterised by the general relation given by (2.9): thus it relates the local value of the density contrast with the smoothed averaged value, independently of any possible choice of the curvature profile. For this reason, measuring the amplitude of the perturbation at $r_m$ is a consistent way to quantify the effect of the curvature profile on the threshold.

In this paper, the criteria for a perturbation to collapse to form a PBH will be stated in terms of the volume-averaged density perturbation (averaged over a sphere of radius $r_m$, corresponding to a top-hat window function with radius $r_m$ centred on the peak of the perturbation) at the time of horizon reentry where the perturbation is taken to behave linearly up to this point (although this is not assumed in the simulations). The formation criterion, and its effect on the calculated abundance of PBHs obtained is discussed in more detail in [37].

# 4 Numerical results of PBH formation

## 4.1 Numerical scheme

The calculations made in this paper to calculate the threshold of PBH formation for the different initial $\zeta$-profiles described in section 2 have been made with the same code as used

<footnote n="3">³It was shown in [29] that off-centered profiles evolve to have a central peak, with almost the same value of the averaged amplitude $\delta$. The analysis, therefore, can be reduced to the centrally peaked values without loss of generality.</footnote>
in [29, 31–33, 39]. 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. The initial data follow from the initial condition seen in section 2, specified on a space-like slice at constant initial cosmic time $t_i$ defined as $a(t_i)r_m e^{\zeta(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 is then evolved using the Misner-Sharp-Hernandez equations so as to generate a second set of initial data on an initial null slice which is 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$).

## 4.2 Threshold, scaling law and mass spectrum

In the left panel of figure 3 we plot the threshold values $\delta_c$ against the parameter $\alpha$ used in (2.10) to vary the initial curvature profile. The lowest limit $\delta_c \simeq 0.41$ corresponds to the analytic solution derived in [35] where the gravitational effect of pressure was taken into account while pressure gradients were instead neglected. It was shown in [29] that this corresponds to shapes of the density contrast with a very large peak ($\alpha \ll 1$) and a smooth tail ($r_0 \gg r_m$): in this configuration most of the matter is already within in the initial cosmological horizon, and only a negligible amount of matter is spread away by the pressure gradients before the black hole is formed, without modifying the shape significantly during the collapse. The maximum value of $\delta_c = 2/3$, ($-r_m\zeta'(r_m) = 1$), corresponds instead to shapes with $r_m = r_0$, ($\alpha \to \infty$), and the density contrast is very steep. In this case the pressure gradients are very large and a substantial modification of the matter configuration is produced during the collapse.

The right panel show the same results of $\delta_c$ as a function of the corresponding behaviour of $r_0/r_m$: the range of values we have been able to compute here are $0.442 \lesssim \delta_c \lesssim 0.656$, ($0.34 \le \alpha \le 2$). Beyond this range the initial profile of the density contrast becomes too extreme, making the numerical simulations unstable. A more detailed analysis of the relationship between the density contrast and the morphology of the initial curvature profile can be found in [29] where different parameterizations of the curvature profiles has been considered, using more than one parameter, getting much closer to the lower limit of $\delta_c \simeq 0.41$.

As has been shown in previous works [20, 21, 31, 33, 39] the mass spectrum of PBHs follows the critical collapse, characterized by a scaling law relation

$$ \frac{M_{\text{PBH}}}{M_H} = \mathcal{K}(\delta - \delta_c)^\gamma \qquad (4.1) $$

where $M_H$ is the mass of the cosmological horizon at horizon crossing, $\gamma \simeq 0.36$ is the critical exponent depending only on the value of $\omega$ of the equation of state and $\mathcal{K}$ is a parameter depending on the particular shape of the density contrast. Because this parameter will play some role in the next section to determine the abundance of PBHs, we have performed numerical simulations to quantify its variation, finding that it varies between 3 and 11 for $\alpha$ varying from 0.4 and 1.9 (corresponding to $\delta_c$ varying between 0.45 and 0.65) for the profiles considered here, with a representative value of $\mathcal{K} \simeq 4$ when $\delta_c \simeq 0.55$, ($\alpha = 1$ in (2.10)).
<figure><img src="image_6.png" /><figcaption>Figure 3. The left plot shows the value of $\delta_c$ as function of $\alpha$ obtained with numerical simulations using an initial curvature profile given by (2.10), while the right panel shows the behaviour of $\delta_c$ as function of the corresponding behaviour if $r_0/r_m$. The bottom dashed horizontal line indicates the lowest limit of the threshold, obtained analytically assuming that the pressure gradients during the collapse are negligible. The upper dashed horizontal line shows to the largest possible value of $\delta_c$, with shapes characterized by very large pressure gradients at the scale $r_m$.</figcaption></figure>

## 5 Calculation of PBH abundance

In this section we will discuss how the abundance of PBHs can be calculated from the primordial curvature perturbation power spectrum $P_\zeta$. The abundance of PBHs will be stated in terms of the energy fraction of the universe (which will be) contained in PBHs at the time of formation, taken for simplicity to be the time of horizon reentry. In principle the time taken for a PBH to form depends slightly on the amplitude of the perturbation collapsing. A formalism for deriving the mass function from a given power spectrum $P_\zeta$ taking into account the non-linear relation between $\zeta$ and $\delta_m$ will also be derived.

We will assume throughout that the curvature perturbation has a Gaussian distribution, partly for simplicity and also motivated by the fact that any local-type non-Gaussianity with $f_{\text{NL}} \gtrsim \mathcal{O}(10^{-3})$ will generate an unacceptably large dark-matter isocurvature perturbation in the CMB [40, 41] — although such bounds can be evaded if non-Gaussianity only couples scales smaller than those observable in the CMB or LSS. In addition, we note that the non-Gaussianity present in single-field inflation (e.g. the Maldacena consistency relation [42]) does not generate isocurvature perturbations [43, 44]. Even when taking ultra slow roll inflation into account, it remains uncertain whether the non-Gaussianity can have a relevant effect [45–48] unless the inflaton field has a non-canonical kinetic term [49–52].

The density contrast $\delta\rho/\rho_b$ is related to the curvature perturbation $\zeta$ as in equation (2.2). However, the key parameter to use for PBH formation is instead the smoothed density contrast $\delta_m$. Using a top-hat window function with areal radius $R = a(t) \exp(\zeta(r_m)) r_m$, the amplitude of (spherically symmetric) peaks in the smoothed density contrast is related to the curvature perturbation in radiation domination as [29]

$$ \delta_m = -\frac{2}{3} r_m \zeta'(r_m) [2 + r_m \zeta'(r_m)] . \qquad (5.1) $$
Because $\zeta$ has a Gaussian distribution, its derivative will also have a Gaussian distribution. Therefore, equation (5.1) can be expressed in terms of a linear Gaussian component $\delta_l = -\frac{4}{3}r_m\zeta'(r_m)$ as<fnref n="4" />

$$ \delta_m = \delta_l - \frac{3}{8}\delta_l^2. \qquad (5.2) $$

The probability density function (PDF) of $\delta_l$ then follows a Gaussian distribution

$$ P(\delta_l) = \frac{1}{\sqrt{2\pi}\sigma^2} \exp\left(-\frac{\delta_l^2}{2\sigma^2}\right). \qquad (5.3) $$

$\delta_l$ represents the linear component of the smoothed density contrast and its variance $\sigma^2$ can be calculated by integrating the linear component of the density power spectrum as follows:

$$ \sigma^2 = \langle \delta_l^2 \rangle = \int_0^\infty \frac{dk}{k} P_{\delta_l}(k, r_m) = \frac{16}{81} \int_0^\infty \frac{dk}{k} (kr_m)^4 \tilde{W}^2(k, r_m) T^2(k, r_m) P_\zeta(k), \qquad (5.4) $$

where $\tilde{W}(k, r_m)$ is the Fourier transform of the top-hat smoothing function, $T(k, r_m)$ is the linear transfer function, and the smoothing scale $r_m$ is equal to the horizon scale. The horizon scale $r_m$ is used to define the time at which $P_{\delta_l}$ (and $T(k, r_m)$) should be evaluated.<fnref n="5" /> For simplicity here, we assume that the relevant scale for PBH formation is given by $kr_m \simeq 1$, although this is not a very accurate approximation, and the exact relation between $k$ and $r_m$ depends on the profile of the density contrast, which depends on the shape of the power spectrum [24].

The Fourier transform of the top-hat smoothing function is given by

$$ \tilde{W}(k, r_m) = 3 \frac{\sin(kr_m) - kr_m \cos(kr_m)}{(kr_m)^3}, \qquad (5.5) $$

and the linear transfer function, where we consider $r_m$ as a time dependent measure of the horizon, is given by [55]

$$ T(k, r_m) = 3 \frac{\sin\left(\frac{kr_m}{\sqrt{3}}\right) - \frac{kr_m}{\sqrt{3}}\cos\left(\frac{kr_m}{\sqrt{3}}\right)}{\left(\frac{kr_m}{\sqrt{3}}\right)^3}. \qquad (5.6) $$

The most straight-forward method to calculate the abundance of PBHs from a non-Gaussian distribution is to work instead with the Gaussian component of the perturbations [56, 57]. To this end, equation (5.2) is solved with $\delta_R = \delta_c$ to find the critical amplitude of the linear component $\delta_{c,l\pm}$,

$$ \delta_{c,l\pm} = \frac{4}{3} \left( 1 \pm \sqrt{\frac{2 - 3\delta_c}{2}} \right), \qquad (5.7) $$

<footnote n="4">⁴It is noted that a recent paper [53] performed a more detailed analysis derived from the skewness of the distribution to achieve the same result (as can be seen from combining equation (30) and (31) in that paper). The correspondence of $-\frac{4}{3}r_m\zeta'(r_m)$ to the linear component of $\delta_r$ and the derivation of its variance are discussed in more detail in appendix A.</footnote>

<footnote n="5">⁵We note, whilst the initial conditions of the perturbations are defined in the super-horizon regime, as in (2.2), we will follow the standard approach used in the literature and evaluate PBH abundance at the time a perturbation enters the horizon (as is done in [53] for example), although note this is somewhat *ad hoc*. A recent paper [54] addressed this point more directly, and estimated the effect of a non-linear transfer function on the density contrast at horizon crossing. Further investigation of the effect of the loss of accuracy from using the linear transfer function is left for future work.</footnote>
where there are 2 solutions because equation (5.2) is quadratic. When necessary, we will take $\delta_c = 0.55$ in this paper, the critical value of the volume averaged density contrast when $\zeta$ is taken to have a Gaussian profile shape. However, we note that only the first solution $\delta_{c,l-}$ corresponds to a physical solution — the second solution corresponds to type 2 perturbations. We will therefore take that a PBH will form in regions where $\delta_c < \delta < 2/3$, where $2/3$ is the maximum value for the density contrast given by equation (5.2). This corresponds to $\delta_{c,l-} < \delta_l < 4/3$.

The final mass of a PBH, $M_{\text{PBH}}$, which forms during radiation domination depends on the shape and amplitude $\delta_R$ of the initial perturbation, and the horizon mass at horizon reentry $M_{\text{H}}$,

$$M_{\text{PBH}} = \mathcal{K} M_{\text{H}} (\delta_m - \delta_c)^\gamma = \mathcal{K} M_{\text{H}} \left( \left[ \delta_l - \frac{3}{8} \delta_l^2 \right] - \delta_c \right)^\gamma, \quad (5.8)$$

where $\mathcal{K}$ depends on the profile shape and typically takes a value between 3 and 5. During radiation domination $\gamma \simeq 0.36$, which is the value we will take [58]. The horizon mass $M_{\text{H}}$ is proportional to the horizon scale squared $r_m^2$,

$$\frac{M_{\text{H}}}{M_{\odot}} \approx \left(\frac{r_m}{r_{\odot}}\right)^2 = \left(\frac{k}{k_{\odot}}\right)^{-2}. \quad (5.9)$$

For a random Gaussian field, the number density of sufficiently rare peaks in a comoving volume is [26]

$$\mathcal{N} = \frac{\mu^3}{4\pi^2\sigma^3}\nu^3 \exp\left(-\frac{\nu^2}{2}\right), \quad (5.10)$$

where $\sigma$ is given by equation (5.4),

$$\nu \equiv \delta_l / \sigma \quad (5.11)$$

and $\mu$ is given by

$$\mu^2 = \int_0^\infty \frac{dk}{k} \mathcal{P}_{\delta_l}(k, r_m) k^2 = \frac{16}{81} \int_0^\infty \frac{dk}{k} (kr_m)^4 \tilde{W}^2(k, r_m) T^2(k, r_m) \mathcal{P}_\zeta(k) k^2. \quad (5.12)$$

Here, the application of equation (5.10) relies on the assumption that peaks in the smoothed linear density field correspond to peaks in the smoothed non-linear density field, such that equation (5.1) can be applied to calculate the height of peaks in the non-linear field. In general, this will not be true, but will be valid in the case that only sufficiently rare and large peaks are considered (which is the same assumption required for the validity of (5.10) and spherical symmetry) — discussed further in appendix B.

The fraction of the energy of the universe at a peak of given height $\nu$ which collapses to form PBHs, $\beta_\nu$, then depends on the mass of the PBHs relative to the horizon scale and the number density of the peaks:

$$\beta_{\nu} = \frac{M_{\text{PBH}}(\nu)}{M_{\text{H}}(r_m)} \mathcal{N}(\nu) \theta(\nu - \nu_c), \quad (5.13)$$

where the time dependence of the horizon mass $M_{\text{H}}(r_m)$ is parameterised by the horizon scale $r_m$. $\theta(\nu - \nu_c)$ is the Heaviside step function which indicates that no PBH will form if $\nu$ is below the critical threshold. The total energy fraction of the universe contained within PBHs at a single time of formation is given by integrating over the range of $\nu$ which forms PBHs,

$$\beta = \int_{\nu_{c-}}^{\frac{4}{3\sigma}} d\nu \frac{\mathcal{K}}{3\pi} \left( \nu\sigma - \frac{3}{8}(\nu\sigma)^2 - \delta_c \right)^{\gamma} \left( \frac{\mu}{aH\sigma} \right)^3 \nu^3 \exp\left(-\frac{\nu^2}{2}\right), \quad (5.14)$$
where $\nu_{c-} \equiv \delta_{c,l-}/\sigma$, see equation (5.7), and equations (5.8), (5.10), and (5.11) have been explicitly substituted into (5.13). The total energy fraction of the universe contained within PBHs today can be approximated by integrating over all times at which PBHs form, parameterised here by the horizon mass (more details of this integral can be found in [59])

$$ \Omega_{\text{PBH}} = \int_{M_{\text{min}}}^{M_{\text{max}}} \mathrm{d}(\ln M_{\text{H}}) \left( \frac{M_{\text{eq}}}{M_{\text{H}}} \right)^{1/2} \beta(M_{\text{H}}), \qquad (5.15) $$

where $M_{\text{H}}$ is the horizon mass at the time of formation and $M_{\text{eq}}$ is the horizon mass at the time of matter-radiation equality. We will use the value $M_{\text{eq}} = 2.8 \times 10^{17} M_{\odot}$ (using the same parameters as [60]). $M_{\text{min}}$ and $M_{\text{max}}$ are respectively the smallest and largest horizon masses at which PBHs form.<fnref n="6" /> The derivation of this formula assumes that the matter density $\Omega_m$ during radiation domination grows proportionally to the scale factor until the time of matter-radiation equality, whereupon the universe immediately becomes matter dominated.

The mass function $f(M_{\text{PBH}})$ can then be obtained by differentiating $\Omega_{\text{PBH}}$ with respect to the PBH mass:

$$ f(M_{\text{PBH}}) = \frac{1}{\Omega_{\text{CDM}}} \frac{\mathrm{d}\Omega_{\text{PBH}}}{\mathrm{d}(\ln M_{\text{PBH}})}, \qquad (5.16) $$

where the equations (5.4), (5.12) and (5.14) should be recast in terms of the horizon mass and PBH mass using the substitutions in equations (5.8) and (5.9). $\Omega_{\text{CDM}}$ will be taken as 0.27 where necessary in this paper.

Broad and narrow peaks in the power spectrum are often considered when calculating the PBH abundance. To give a concrete example, we will consider the two extreme cases — a scale invariant power spectrum, and an extremely narrow peak. For the scale invariant case, we will take $P_{\zeta} = A_s = \text{constant}$. For the narrow peak, we will take the peak to be narrow enough such that it can be treated as a Dirac-delta function, $P_{\zeta} = A_s \delta_D(\ln(k/k_*))$, although note that such a power spectrum is unphysical (for example, [10] describes that the power spectrum cannot be steeper than $k^4$, at least in the context of single-field inflation).

For the scale invariant case, equations (5.4) and (5.12) give scale invariant values of $\sigma^2 \sim 1.06 A_f$ and $(\mu/(aH))^2 \sim 6.86 A_s$ respectively. Figure 4 shows the abundance of PBHs, $\beta$ (equation (5.14)), as a function of $A_s$, whilst figure 4 shows the mass function $f(M)$ (equation (5.16)) evaluated at $M = M_{\text{PBH}}$ as a function of the power spectrum amplitude.

In the case of the narrowly peaked spectrum, we assume without loss of generality that the power spectrum peaks at a scale corresponding to a horizon mass of $1M_{\odot}$, with corresponding horizon scale $k_* = k_{\odot}$, in order to allow a direct comparison to the broad power spectrum. We will also consider that PBH formation occurs only at the horizon scale corresponding exactly to the peak in the power spectrum, and that $\beta$ corresponds to the total energy fraction collapsing into PBHs at all epochs, rather than integrating over $\ln M_{\text{H}}$ as in equation (5.15). Whilst equation (5.4) will give a significant variance $\sigma^2$ for values of $r_m$ close to $r_{\odot}$ (suggesting perturbations of varying scales exist), this is because $\delta_m$ does not immediately become zero when the smoothing scale is not set exactly equal to the scale of the perturbation. However, the scale of all perturbations here is fixed, and so will all enter the horizon at the same scale. In this case, evaluated as $k_{\odot}$ enters the horizon, equations (5.4) and (5.12) give $\sigma^2 = (\mu/(aH))^2 \sim 0.151 A_s$.

What can be learned from these figures is that accounting for the non-linear relationship between the density contrast and the curvature perturbation will always serve to reduce the

<footnote n="6">⁶In principle $M_{\text{max}}$ can be arbitrarily large because $\delta_m - \delta_c$ can be arbitrarily small. However, the largest PBH mass which can be formed when the horizon mass is $M_{\text{H}}$ is given by $M_{\text{PBH}} = K M_{\text{H}} (2/3 - \delta_c)^{\gamma}$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left: the energy fraction of the universe collapsing to form PBHs at the time of formation, $\beta$, is plotted against the amplitude of a scale invariant power spectrum, $P_{\zeta} = A_s$. Right: the mass function of PBHs $f(M_{\odot})$ of PBHs (defined in equation (5.16)) with a mass of $1M_{\odot}$ as a function of the amplitude of $A_s$. Accounting for the non-linear relation means that a significantly smaller abundance of PBHs is calculated, by many orders of magnitude. For $\delta_c = 0.5$ then $K = 5.5$ and for $\delta_c = 0.6$ then $K = 3.5$. The abundance of PBHs is dominated by the uncertainty in $\delta_c$ (which has an exponential effect), rather than $K$ (which has a linear effect).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: the energy fraction of the universe collapsing to form PBHs at the time of formation, $\beta$, is plotted against the amplitude of the power spectrum. Right: the mass function of PBHs $f(M_{\odot})$ of PBHs with a mass of $1M_{\odot}$ as a function of the amplitude of the power spectrum. In both plots, a Dirac delta power spectrum has been assumed, $P_{\zeta} = A_s\delta_D(\ln(k/k_*))$. Using the profile shape predicted from a Dirac-delta power spectrum, $K = 4$, and $\delta_c = 0.51$.</figcaption></figure>

calculated abundance of PBHs, compared to using the linear model. It can be seen from comparing the figures from the scale invariant power spectrum to the narrowly peaked power spectrum, that the abundance of PBHs is strongly dependent on the shape of the power spectrum rather than simply the amplitude — and so constraints on the power spectrum from constraints on PBH abundance should be calculated on a case-by-case basis for different inflationary models which predict different shapes for the power spectrum. This fact has been well known for some time and was investigated in more detail recently by [24].

However, we will here make a simple calculation to describe by what fraction the amplitude of the curvature perturbation power spectrum $A_s$ needs to increase in order to obtain the same value for $\beta$ as when the linear relation is used. The dominant term in the calculation for the abundance of black holes, equation (5.14), is the exponential term $\exp(-\nu^2)$. After integrating, this will give a dependence roughly proportional to $\exp(-\nu_{c-}^2)$. In order to produce (approximately) the same number of PBHs from the linear calculation as from the non-linear calculation, we therefore require $\nu_{c,L} = \nu_{c,NL}$ — where the subscripts L and
NL denote the linear and non-linear calculations respectively,

$$ \frac{\delta_{c,L}}{\sigma_L} = \frac{\delta_{c,l-}}{\sigma_{NL}}, \qquad (5.17) $$

where $\delta_{c,l-}$ is given by equation (5.7). In both the linear and non-linear case, $\sigma^2$ for a given power spectrum shape is proportional to $A_s$ by the same constant of proportionality, $\sigma^2 = CA_s$. Finally, we can write down the ratio between $A_L$ and $A_{NL}$ as

$$ \frac{A_{NL}}{A_L} = \left( \frac{\delta_{c,l-}}{\delta_c} \right)^2 = \left( \frac{4 \left( 1 - \sqrt{\frac{2-3\delta_c}{2}} \right)}{3\delta_c} \right)^2. \qquad (5.18) $$

For values $0.41 < \delta_c < 2/3$, this factor varies approximately from 1.5 to 4. For more typical values $0.5 < \delta_c < 0.6$, in order for the same number of PBHs to form, the amplitude of the power spectrum must be 1.78-2.31 times greater than previously calculated with the linear model. In particular, for the commonly considered case of a Gaussian profile in $\zeta$, with $\delta_c = 0.55$, the power spectrum needs to be enhanced by a factor of 2.0 in order to get the same number of PBHs forming if one neglected the non-linear relationship between $\zeta$ and $\delta$.

# 6 Summary

In the radiation dominated epoch following the end of inflation, perturbations can collapse to form PBHs if the perturbation amplitude is large enough. Whether or not a PBH will form depends on the amplitude of the density contrast $\delta_m$, rather than the amplitude of the curvature perturbation $\zeta$. The non-linear relation between the curvature perturbation $\zeta$ and the density contrast $\delta\rho/\rho_b$, given by equation (2.2), means that $\delta_m$ will have a non-Gaussian distribution even when $\zeta$ is perfectly Gaussian, which has a significant effect on the number of PBHs which form in the early universe.

We have discussed briefly the formation criterion which should be used to determine whether a perturbation will collapse or not, and in this paper, we argue that the volume-averaged (smoothed) density contrast of a peak $\delta_m$ should be used rather than the value $\delta\rho/\rho_b$ evaluated at $r=0$. The scale over which the perturbation should be averaged, $r_m$, is the scale at which the compaction function $C(r)$ (defined in equation (2.3)) is maximised, and note that at this scale $C(r)$ is equal to $\delta_m$. In this paper, we use the amplitude of $\delta_m$ measured at the cosmological horizon entry to determine whether a PBH will form — although point out that this is somewhat inconsistent as the expression for $\delta\rho/\rho_b$ is valid in the super-horizon regime and the perturbation will not evolve linearly until horizon entry when they have a large amplitude. Further consideration of these factors is left for future study.

Making use of numerical simulations described in section 4, we have considered different profiles of $\zeta$ which form density perturbations to determine a threshold value for PBH formation and how this can depend on the shape of the profile in $\zeta$. We noted in section 2 that the relation between peaks in $\zeta$ and $\delta\rho/\rho_b$ is not straightforward and may not coincide. For a given profile shape of $\zeta$, the profile of $\delta\rho/\rho_b$ depends also on the amplitude of the perturbation, and we have calculated a relation between the amplitude of the perturbation $\delta_m$ and the mass of the PBH formed accounting for this (often referred to as the extended mass function, rather than assuming monochromatic formation of PBHs at a given epoch).
There is a simple relation between $\delta_m$ and $\zeta$, given by equation (5.1), which can be used to fully describe the non-Gaussianity of $\delta_m$ when $\zeta$ is taken as Gaussian. Making use of this relation we have derived a formalism to derive the abundance of PBHs. The abundance of PBHs may be expressed either in terms of the energy fraction collapsing into black holes at the time of formation $\beta$, the present-time density parameter for PBHs, $\Omega_{PBH}$, or the mass fraction of dark matter contained within PBHs of a given mass, $f(M_{PBH})$. These expressions are calculated utilising the theory of peaks and accounting for the extended mass function of PBHs.

When the non-Gaussianity of the density $\delta_m$ is taken into account, we find that this always reduces the number of PBHs which form by many orders of magnitude, see figures 4–5. We reproduce the known result that the abundance of PBHs depends upon both the amplitude and shape of the curvature perturbation power spectrum. In order for a comparable number of PBHs to form compared to using the linear relation between $\zeta$ and $\delta$, we find that the amplitude of the power spectrum must therefore be $\sim 2-3$ times larger, with the amount depending on the value taken for the collapse threshold, see (5.18), but otherwise being almost independent of the shape of the power spectrum.

Finally we note that the non-linear relation for the smoothed density contrast in terms of the spatial derivative of the curvature perturbation makes a clear analogy to local non-Gaussianity, with a negative value of (local) $f_{\text{NL}}$ that suppresses PBH formation (see also the “Note added” below). However, the analogy is potentially misleading since this non-Gaussian term only affects the one-point function of $\delta_m$ and it does not generate a bispectrum because the derivatives of $\zeta$ are uncorrelated on scales much larger than the PBH scale, $r_m$. Therefore it does not correspond to a coupling between long and short-wavelength modes and hence does not generate a large-scale dark matter isocurvature perturbation (which would be observationally ruled out [40, 41]). This also implies that the “apparent” negative local $f_{\text{NL}}$ cannot lead to an enhanced formation of PBHs in the case of a primordial power spectrum with a broad peak, in contrast to a physical value of $f_{\text{NL}} < 0$ of $\zeta$, see [61] for details.

The relation used between $\zeta$ and $\delta$, equation (2.2), is exact in $\zeta$, only neglecting higher-order terms in $\epsilon$, valid because $\epsilon \ll 1$ in the super-horizon limit. Therefore, we have quantified the effect of the complete non-linear relationship upon the abundance of PBHs.

**Note added:** while this paper was approaching completion a related work by Kawasaki and Nakatsuka [53] appeared on the arXiv. Our broad conclusions are in agreement with their paper, principally that the non-linear relationship between $\zeta$ and $\delta$ makes the formation of PBHs less likely. For typical values of $\delta_c \sim 0.5$, we confirm their finding that the suppression of PBH formation due to the non-linear terms requires the power spectrum to have an amplitude $\sim 1.4^2$ times greater in order to form the same number of PBHs.

Produced and appearing in parallel with our paper, the paper by De Luca et al. [62] studies the same topic. Although there are some significant differences in our methodology, our results are in broad agreement, both showing that the non-linear relation between $\delta\rho/\rho_b$ and $\zeta$ leads to a suppression in the PBH formation rate. Unlike these two papers, we take critical collapse into account.

## Acknowledgments

We thank V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto and C. Ünal, the authors of [62], for sharing their draft of a paper which was written in parallel to this one and
helpful discussions. We thank Nicola Bellomo and Eiichiro Komatsu for helpful comments on a draft of this paper and we thank Jaume Garriga, Cristiano Germani, Marcello Musso and Licia Verde for helpful discussions.

SY is funded by a Humboldt Research Fellowship for Postdoctoral Researchers. IM is supported by the Unidad de Excelencia María de Maeztu Grant No. MDM-2014-0369. CB is funded by a Royal Society University Research Fellowship.

A Gaussianity and variance of $-\frac{4}{3}r_m\zeta'(r_m)$

Here we will derive the variance and discuss the Gaussianity of the linear term in equation (5.1), $\delta_l = -\frac{4}{3}r_m\zeta'(r_m)$. The linear relation between the density contrast $\delta\rho/\rho_b$ and the curvature perturbation $\zeta$ in radial coordinates $r$ is given by

$$ \frac{\delta\rho}{\rho_b}(r,t) = -\frac{4}{9}\epsilon^2(t)r_m^2 \left( \zeta''(r) + \frac{2}{r}\zeta'(r) \right), \qquad (A.1) $$

where the prime denotes a derivative with respect to the radial co-ordinate $r$, and the parameter $\epsilon$ describes the super-horizon time-evolution of the perturbation and is given in the linear approximation by $\epsilon = (r_m a H)^{-1}$, where $r_m$ is the perturbation scale and $aH$ is the horizon scale. In this paper we consider the perturbations at the time of horizon reentry, $t_H$, such that $r_m = (aH)^{-1} \Rightarrow \epsilon(t_H) = 1$, and we will therefore drop the $\epsilon(t)$. The volume averaged density contrast used to determine PBH formation is given by

$$ \delta_m = \frac{1}{V} \int_0^{r_m} 4\pi r^2 \frac{\delta\rho}{\rho_b}(r, t_H) dr, \qquad (A.2) $$

where in the linear approximation the volume is given by $V = \frac{4}{3}\pi r_m^3$. Although this integral is normally integrated over the comoving distance $r_{com} \equiv re^{\zeta(r)}$, this difference represents higher-order effects neglected in the linear approximation. Combining the above equations gives

$$ \delta_m = -\frac{4}{3r_m} \int_0^{r_m} dr r^2 \left( \zeta''(r) + \frac{2}{r} \zeta'(r) \right) = -\frac{4}{3} r_m \zeta'(r_m) = \delta_l, \qquad (A.3) $$

which is the linear term seen in equation (5.1).

We now consider the variance of $\delta_m$, noticing that equation (A.2) can be considered as a top-hat window function centred on the peak convolved with the density contrast

$$ \delta_m(X) = \frac{1}{V} \int_0^\infty dr 4\pi r^2 \frac{\delta\rho}{\rho_b}(r, t_H) \theta(r_m - r) = \int dx^3 \frac{\delta\rho}{\rho_b}(x, t_H) W(X - x, r_m), \qquad (A.4) $$

where $X$ is the location of a peak in cartesian coordinates (corresponding to $r = 0$) and $\theta(r)$ is the Heaviside step function. The second equality is the same integral expressed in cartesian coordinates, with the window function $W(x, r_m)$ given by:

$$ W(x, r_m) = \frac{\theta(r_m - |x|)}{\frac{4}{3}\pi r_m^3}. \qquad (A.5) $$

For our purposes, it is more convenient to express equation (A.4), a convolution in real space, as a multiplication in Fourier space,

$$ \delta_m(k) = \tilde{W}(k, r_m)\delta(k), \qquad (A.6) $$
where the Fourier transform of the window function $\tilde{W}(k, r_m)$ is given by equation (5.5) and $\delta(k)$ is given by

$$ \delta(k) = -\frac{4}{9}T(k, r_m)(kr_m)^2\zeta(k). \qquad (A.7) $$

where the linear transfer function $T(k, r_m)$ is given by equation (5.6) and $\delta_m(k)$ can therefore be written as

$$ \delta_m(k) = -\frac{4}{9}\tilde{W}(k, r_m)T(k, r_m)(kr_m)^2\zeta(k). \qquad (A.8) $$

Since we have assumed $\zeta(x)$ has a Gaussian distribution, it also has a Gaussian distribution in Fourier space: $\zeta(k)$ has a Gaussian distribution (being a linear combination of Gaussian variables $\zeta(x)$). $\delta_m(k)$ is then related to $\zeta(k)$ by a linear factor, meaning that $\delta_m(k)$ also has a Gaussian distribution. Finally, $\delta_m(x)$ is a linear combination of the Gaussian Fourier modes; hence it also has a Gaussian distribution.

Finally, we can calculate the variance $\sigma^2$ by integrating the power spectrum,

$$ \sigma^2 = \langle \delta_l^2 \rangle = \int_0^\infty \frac{dk}{k} \mathcal{P}_{\delta_l}(k, r_m) = \frac{16}{81} \int_0^\infty \frac{dk}{k} (kr_m)^4 \tilde{W}^2(k, r_m) T^2(k, r_m) \mathcal{P}_\zeta(k). \qquad (A.9) $$

## B Correspondence of large peaks

We will here consider the correspondence between peaks in the various fields considered within the context of this paper — the curvature perturbation field $\zeta$, the non-linear density contrast field $\delta\rho/\rho_b$, the non-linear smoothed density field $\delta_m$ and the linear smoothed density field $\delta_l$. For the purposes of this discussion, type 2 perturbations will not be considered (corresponding to $\delta_m > 2/3$, as the abundance of such perturbations is exponentially suppressed (even compared to the exponentially small number of perturbations which form PBHs).

Figure 2 shows the critical profiles of peaks in $\zeta$ and $\delta\rho/\rho_b$ (in spherical symmetry). All of the profiles in $\zeta$ have a central peak at $r = 0$, whilst the density field can have an off-centred peak or a divergence at the centre — and it cannot necessarily therefore be stated that a PBH forms at peaks in $\delta\rho/\rho_b$. As mentioned in section 3 this problem is overcome by using the smoothed density contrast — figure 6 shows the same profiles smoothed on a scale $r_m$. The off-centred peaks are no longer seen because the smoothing scale $r_m$ will by definition larger than the radius at which the density peaks, and, being a feature smaller than the smoothing radius, is therefore removed in the process of smoothing. The divergences at the centre (while unphysical) are nonetheless removed because the divergence, whilst they represent infinite density as $r \to 0$, they do not represent infinite mass — and the integral during the smoothing therefore converges. It can therefore be stated that, for isolated spherically-symmetric type 1 perturbations, the peaks in $\zeta$ correspond to peaks in the smoothed density field $\delta_m$.

In order to calculate the abundance, it has been assumed that peaks in the (Gaussian) linear smoothed density field correspond to peaks in the (non-Gaussian) non-linear field. In general, this is not expected to be true — however, for the very large and rare peaks from which PBHs form, it is expected that this should be a valid approximation. This is primarily due to the fact that the large peaks in the Gaussian fields ($\zeta$ or $\delta_l$) are very unlikely to be close to other large peaks, that is, that the local region surrounding large peaks contains only significantly smaller perturbations. This is expected to result in the large peak being (approximately) spherically symmetric peak [26]. When the non-linear smoothed density field is calculated, this spherical symmetry is preserved — implying that the perturbation in the non-linear field must also be a peak (when smoothed on an appropriate scale, see above).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Profile shapes in the density contrast after a top-hat smoothing function has been applied. The profiles are calculated from a curvature perturbation profile given by equation (2.10), and each profile is smoothed on a scale $r_m$. The values of $\alpha$ are $0.5, 0.75, 1.0, 1.25, 1.50, 2.0$, where increasing $\alpha$ corresponds to more negative troughs seen in the figure. The amplitude in each case is $A = 0.5$.</figcaption></figure>

The correlation of peaks in $\zeta$ and the density $\delta\rho/\rho_b$ was investigated in [62], concluding that “one can associate the number of rare peaks in the overdensity with the number of peaks in the curvature perturbation which are spiky enough”, validating the assumptions applied here.

Furthermore, only modes which are similar to the smoothing scale (which is considered equal to the horizon scale) have a non-negligible effect on the smoothed density field — larger-scale and smaller-scale modes are suppressed by the $k^4$ term and the smoothing and transfer functions respectively in equation (5.4). The smoothed density fields (linear and non-linear) therefore only inherit peaks from the $\zeta$-field on a small range of scales. This means that, whilst a perturbation of scale $r_m$ in the $\zeta$-field is likely to have many scaller-scale peaks on top of it (and may also sit on top of a much larger scale peak), this is not the case for the smoothed density fields.

## References

[1] S. Hawking, *Gravitationally collapsed objects of very low mass*, Mon. Not. Roy. Astron. Soc. **152** (1971) 75 [INSPIRE].

[2] B.J. Carr and S.W. Hawking, *Black holes in the early Universe*, Mon. Not. Roy. Astron. Soc. **168** (1974) 399 [INSPIRE].

[3] B.J. Carr, *The primordial black hole mass spectrum*, Astrophysical J. **201** (1975) 1.

[4] B. Carr et al., *Primordial black hole constraints for extended mass functions*, Phys. Rev. D **96** (2017) 023514 [arXiv:1705.05567] [INSPIRE].

[5] F. Kühnel and K. Freese, *Constraints on primordial black holes with extended mass functions*, Phys. Rev. D **95** (2017) 083508 [arXiv:1701.07223] [INSPIRE].

[6] N. Bellomo, J.L. Bernal, A. Raccanelli and L. Verde, *Primordial black holes as dark matter: converting constraints from monochromatic to extended mass distributions*, JCAP **01** (2018) **004** [arXiv:1709.07467] [INSPIRE].
[7] S. Clesse and J. García-Bellido, *Seven hints for primordial black hole dark matter*, Phys. Dark Univ. 22 (2018) 137 [arXiv:1711.10458] [INSPIRE].

[8] A.M. Green, *Primordial black holes: sirens of the early Universe*, Fundam. Theor. Phys. 178 (2015) 129 [arXiv:1403.1198].

[9] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, *Primordial black holes — Perspectives in gravitational wave astronomy*, Class. Quant. Grav. 35 (2018) 063001 [arXiv:1801.05235] [INSPIRE].

[10] C.T. Byrnes, P.S. Cole and S.P. Patil, *Steepest growth of the power spectrum and primordial black holes*, JCAP 06 (2019) 028 [arXiv:1811.11158] [INSPIRE].

[11] M. Drees and E. Erfani, *Running-mass inflation model and primordial black holes*, JCAP 04 (2011) 005 [arXiv:1102.2340] [INSPIRE].

[12] E. Bugaev and P. Klimai, *Axion inflation with gauge field production and primordial black holes*, Phys. Rev. D 90 (2014) 103501 [arXiv:1312.7435] [INSPIRE].

[13] O. Özsoy, S. Parameswaran, G. Tasinato and I. Zavala, *Mechanisms for primordial black hole production in string theory*, JCAP 07 (2018) 005 [arXiv:1803.07626] [INSPIRE].

[14] J. García-Bellido, A.D. Linde and D. Wands, *Density perturbations and black hole formation in hybrid inflation*, Phys. Rev. D 54 (1996) 6040 [astro-ph/9605094] [INSPIRE].

[15] D.H. Lyth, *The hybrid inflation waterfall and the primordial curvature perturbation*, JCAP 05 (2012) 022 [arXiv:1201.4312] [INSPIRE].

[16] E. Bugaev and P. Klimai, *Formation of primordial black holes from non-Gaussian perturbations produced in a waterfall transition*, Phys. Rev. D 85 (2012) 103504 [arXiv:1112.5601] [INSPIRE].

[17] G. Ballesteros, J. Beltran Jimenez and M. Pieroni, *Black hole formation from a general quadratic action for inflationary primordial fluctuations*, JCAP 06 (2019) 016 [arXiv:1811.03065] [INSPIRE].

[18] M. Shibata and M. Sasaki, *Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity*, Phys. Rev. D 60 (1999) 084002 [gr-qc/9905064] [INSPIRE].

[19] A.M. Green, A.R. Liddle, K.A. Malik and M. Sasaki, *A new calculation of the mass fraction of primordial black holes*, Phys. Rev. D 70 (2004) 041502 [astro-ph/0403181] [INSPIRE].

[20] J.C. Niemeyer and K. Jedamzik, *Near-critical gravitational collapse and the initial mass function of primordial black holes*, Phys. Rev. Lett. 80 (1998) 5481 [astro-ph/9709072] [INSPIRE].

[21] J.C. Niemeyer and K. Jedamzik, *Dynamics of primordial black hole formation*, Phys. Rev. D 59 (1999) 124013 [astro-ph/9901292] [INSPIRE].

[22] S. Young, C.T. Byrnes and M. Sasaki, *Calculating the mass fraction of primordial black holes*, JCAP 07 (2014) 045 [arXiv:1405.7023] [INSPIRE].

[23] T. Harada, C.-M. Yoo, T. Nakama and Y. Koga, *Cosmological long-wavelength solutions and primordial black hole formation*, Phys. Rev. D 91 (2015) 084057 [arXiv:1503.03934] [INSPIRE].

[24] C. Germani and I. Musco, *Abundance of primordial black holes depends on the shape of the inflationary power spectrum*, Phys. Rev. Lett. 122 (2019) 141302 [arXiv:1805.04087] [INSPIRE].

[25] C.-M. Yoo, T. Harada, J. Garriga and K. Kohri, *Primordial black hole abundance from random Gaussian curvature perturbations and a local density threshold*, PTEP 2018 (2018) 123E01 [arXiv:1805.03946] [INSPIRE].
[26] J.M. Bardeen, J.R. Bond, N. Kaiser and A.S. Szalay, *The statistics of peaks of gaussian random fields*, *Astrophys. J.* **304** (1986) 15 [INSPIRE].

[27] D.H. Lyth, K.A. Malik and M. Sasaki, *A general proof of the conservation of the curvature perturbation*, *JCAP* **05** (2005) 004 [astro-ph/0411220] [INSPIRE].

[28] D.S. Salopek and J.R. Bond, *Nonlinear evolution of long wavelength metric fluctuations in inflationary models*, *Phys. Rev. D* **42** (1990) 3936 [INSPIRE].

[29] I. Musco, *The threshold for primordial black holes: dependence on the shape of the cosmological perturbations*, arXiv:1809.02127 [INSPIRE].

[30] M. Kopp, S. Hofmann and J. Weller, *Separate universes do not constrain primordial black hole formation*, *Phys. Rev. D* **83** (2011) 124025 [arXiv:1012.4369] [INSPIRE].

[31] I. Musco, J.C. Miller and L. Rezzolla, *Computations of primordial black hole formation*, *Class. Quant. Grav.* **22** (2005) 1405 [gr-qc/0412063] [INSPIRE].

[32] A.G. Polnarev and I. Musco, *Curvature profiles as initial conditions for primordial black hole formation*, *Class. Quant. Grav.* **24** (2007) 1405 [gr-qc/0605122] [INSPIRE].

[33] I. Musco and J.C. Miller, *Primordial black hole formation in the early universe: critical behaviour and self-similarity*, *Class. Quant. Grav.* **30** (2013) 145009 [arXiv:1201.2379] [INSPIRE].

[34] T. Nakama, T. Harada, A.G. Polnarev and J. Yokoyama, *Identifying the most crucial parameters of the initial curvature profile for primordial black hole formation*, *JCAP* **01** (2014) 037 [arXiv:1310.3007] [INSPIRE].

[35] T. Harada, C.-M. Yoo and K. Kohri, *Threshold of primordial black hole formation*, *Phys. Rev. D* **88** (2013) 084051 [Erratum ibid. D **89** (2014) 029903] [arXiv:1309.4201] [INSPIRE].

[36] T. Harada and S. Jhingan, *Spherical and nonspherical models of primordial black hole formation: exact solutions*, *PTEP* **2016** (2016) 093E04 [arXiv:1512.08639] [INSPIRE].

[37] S. Young, *The primordial black hole formation criterion re-examined: parameterisation, timing and the choice of window function*, arXiv:1905.01230 [INSPIRE].

[38] K. Ando, K. Inomata and M. Kawasaki, *Primordial black holes and uncertainties in the choice of the window function*, *Phys. Rev. D* **97** (2018) 103528 [arXiv:1802.06393] [INSPIRE].

[39] I. Musco, J.C. Miller and A.G. Polnarev, *Primordial black hole formation in the radiative era: investigation of the critical nature of the collapse*, *Class. Quant. Grav.* **26** (2009) 235001 [arXiv:0811.1452] [INSPIRE].

[40] S. Young and C.T. Byrnes, *Signatures of non-Gaussianity in the isocurvature modes of primordial black hole dark matter*, *JCAP* **04** (2015) 034 [arXiv:1503.01505] [INSPIRE].

[41] Y. Tada and S. Yokoyama, *Primordial black holes as biased tracers*, *Phys. Rev. D* **91** (2015) 123534 [arXiv:1502.01124] [INSPIRE].

[42] J.M. Maldacena, *Non-Gaussian features of primordial fluctuations in single field inflationary models*, *JHEP* **05** (2003) 013 [astro-ph/0210603] [INSPIRE].

[43] E. Pajer, F. Schmidt and M. Zaldarriaga, *The observed squeezed limit of cosmological three-point functions*, *Phys. Rev. D* **88** (2013) 083502 [arXiv:1305.0824] [INSPIRE].

[44] G. Cabass, E. Pajer and F. Schmidt, *How gaussian can our Universe be?*, *JCAP* **01** (2017) 003 [arXiv:1612.00033] [INSPIRE].

[45] R. Bravo, S. Mooij, G.A. Palma and B. Pradenas, *A generalized non-Gaussian consistency relation for single field inflation*, *JCAP* **05** (2018) 024 [arXiv:1711.02680] [INSPIRE].

[46] R. Bravo, S. Mooij, G.A. Palma and B. Pradenas, *Vanishing of local non-Gaussianity in canonical single field inflation*, *JCAP* **05** (2018) 025 [arXiv:1711.05290] [INSPIRE].
[47] V. Atal and C. Germani, *The role of non-gaussianities in primordial black hole formation*, Phys. Dark Univ. (2018) 100275 [arXiv:1811.07857].

[48] S. Passaglia, W. Hu and H. Motohashi, *Primordial black holes and local non-Gaussianity in canonical inflation*, Phys. Rev. D 99 (2019) 043536 [arXiv:1812.08243] [INSPIRE].

[49] S. Shandera, A.L. Erickcek, P. Scott and J.Y. Galarza, *Number counts and non-Gaussianity*, Phys. Rev. D 88 (2013) 103506 [arXiv:1211.7361] [INSPIRE].

[50] S. Young, D. Regan and C.T. Byrnes, *Influence of large local and non-local bispectra on primordial black hole abundance*, JCAP 02 (2016) 029 [arXiv:1512.07224] [INSPIRE].

[51] G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, *Primordial black holes from inflation and non-Gaussianity*, JCAP 03 (2018) 016 [arXiv:1801.09415] [INSPIRE].

[52] A.Y. Kamenshchik, A. Tronconi, T. Vardanyan and G. Venturi, *Non-canonical inflation and primordial black holes production*, Phys. Lett. B 791 (2019) 201 [arXiv:1812.02547] [INSPIRE].

[53] M. Kawasaki and H. Nakatsuka, *Effect of nonlinearity between density and curvature perturbations on the primordial black hole formation*, Phys. Rev. D 99 (2019) 123501 [arXiv:1903.02994] [INSPIRE].

[54] A. Kalaja et al., *From primordial black holes abundance to primordial curvature power spectrum (and back)*, JCAP 10 (2019) 010 [arXiv:1908.03596] [INSPIRE].

[55] A.S. Josan, A.M. Green and K.A. Malik, *Generalised constraints on the curvature perturbation from primordial black holes*, Phys. Rev. D 79 (2009) 103520 [arXiv:0903.3184] [INSPIRE].

[56] C.T. Byrnes, E.J. Copeland and A.M. Green, *Primordial black holes as a tool for constraining non-Gaussianity*, Phys. Rev. D 86 (2012) 043512 [arXiv:1206.4188] [INSPIRE].

[57] S. Young and C.T. Byrnes, *Primordial black holes in non-Gaussian regimes*, JCAP 08 (2013) 052 [arXiv:1307.4995] [INSPIRE].

[58] C.R. Evans and J.S. Coleman, *Observation of critical phenomena and selfsimilarity in the gravitational collapse of radiation fluid*, Phys. Rev. Lett. 72 (1994) 1782 [gr-qc/9402041] [INSPIRE].

[59] C.T. Byrnes, M. Hindmarsh, S. Young and M.R.S. Hawkins, *Primordial black holes with an accurate QCD equation of state*, JCAP 08 (2018) 041 [arXiv:1801.06138] [INSPIRE].

[60] T. Nakama, J. Silk and M. Kamionkowski, *Stochastic gravitational waves associated with the formation of primordial black holes*, Phys. Rev. D 95 (2017) 043511 [arXiv:1612.06264] [INSPIRE].

[61] S. Young and C.T. Byrnes, *Long-short wavelength mode coupling tightens primordial black hole constraints*, Phys. Rev. D 91 (2015) 083521 [arXiv:1411.4620] [INSPIRE].

[62] V. De Luca et al., *The ineludible non-Gaussianity of the primordial black hole abundance*, JCAP 07 (2019) 048 [arXiv:1904.00970] [INSPIRE].