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

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 $\mathcal{P}_\zeta$ must be a factor of $\sim2$ larger to produce the same number of PBHs as the linear model (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.

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][5][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 β, 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 (O(10 −2 ) compared to 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][15][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 ζ is typically used to study cosmological perturbations generated with the different models and to predict their observable consequences 1 . The curvature perturbation ζ appears in the metric in the comoving uniform-density gauge as where η 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 ζ 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 ζ [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 ζ was above a critical value ζ c of order unity. However, it has since been shown that the curvature perturbation ζ 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 δρ/ρ b rather than the curvature perturbation for the calculation of the abundance have since used a linear relation between the two parameters (as in the recent paper [24] for example), where ω = 1/3 is the equation of state parameter during the radiation dominated epoch of the early universe, ρ 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 1 Often R is used instead to denote the curvature perturbation, whilst ζ is used to describe the curvature perturbation only on a uniform density slicing. energy density profile. Starting for the first time from simulations of PBH formation arising from perturbations in the curvature perturbation ζ, the aims of this paper are to investigate how the fully non-linear relation between ζ and δρ/ρ 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 ζ and the abundance of PBHs at formation β.
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 ζ, and then calculating the abundance of PBHs by utilising the Gaussianity of ζ. Using the critical height of peaks in ζ 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 δ and ζ, whereas [25] compares very different methods.
The paper is organised as follows: in section II we will discuss the set up of the initial conditions in the density contrast δρ/ρ b arising from an initial curvature perturbation ζ. In section III we will discuss the criteria which should be used to determine whether an initial perturbation will eventually collapse to form a PBH. In section IV we will discuss the simulation procedure used and the numerical results obtained from the simulations. Finally, in section V we will show how the abundance of PBHs β can be obtained from the curvature perturbation power spectrum. Our findings are summarised in section VI, leaving some details of our calculations to an appendix.

II. COSMOLOGICAL PERTURBATIONS IN THE SUPER HORIZON REGIME
In this section we will first describe the general relation between the curvature perturbation ζ and the density contrast δρ/ρ b before analysing a specific parametrization of ζ that allows us to vary the profile of δρ/ρ b . This allows us, with the help of numerical simulations (see Section IV), to span almost all the possible range of values of δ 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 ζ in the literature is typically defined on a uniform density slicing, whilst the density contrast δρ/ρ 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).

A. 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 ( 1) that is conveniently identified with the ratio between the Hubble radius 1/H(t) and the length scale L characterising the perturbation where a particular value of corresponds to a particular value of the time t.
When 1 the curvature profile ζ(r) is conserved (time independent) and each spatial gradient is multiplied by , expanding the equations in a power series in up to the first non-zero order.
Putting = 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 ζ(r). The curvature profile represents the fundamental source of the perturbation, embedded into the metric (1) from the asymptotic limit, t → 0. Although this approach reproduces the time evolution of linear perturbation theory when 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 , 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 δρ/ρ b defined on a comoving uniform-cosmic time slicing and the curvature perturbation ζ is given by [23,25,29], which is exact up to the first non-zero terms in where ω ≡ p/ρ is the equation of state parameter, and is the Laplacian written assuming spherical symmetry. There are two non-linear effects contained within equation (4): the exponential term exp(−2ζ(r)) and the quadratic term of the first derivative (ζ (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 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πρ 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 ζ(r) as [23,25,29] 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 ζ (r m ) + r m ζ (r m ) = 0 .
Measuring the curvature perturbation with ζ(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) 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 and therefore according to the definition of given above, the physical length scale of the perturbation, to be called R m from here onwards, is given by 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 Rm = 4π 3 R 3 m , the amplitude of the perturbation is given by and using the explicit expression of δρ(r, t)/ρ b in terms of the curvature profile seen in (4), we get δ m ≡ δ(r m , t H ) = C(r m ), which satisfies the fundamental relation [29] for any curvature profile, ζ(r).

B. Initial conditions
We will now study the main feature of the density profile when an explicit parameterization of the curvature profile ζ(r) is specified as where A and r m respectively denote the amplitude and the scale of the perturbation. Inserting this into (4), the corresponding density profile is given by and then inserting (12) into (6) and (7), we can calculate the corresponding perturbation amplitude which gives the value of A in terms of the averaged amplitude δ m This behaviour of δ m in terms of A also shows that there is a maximum value of δ m = 2/3 corresponding to A = eα/2, that represents the transition between PBHs of type I and type II, after which formation of PBHs cannot be studied in terms of δ m but only in terms of ζ [30].
Using numerical simulations we have calculated the critical values for PBHs using initial condition in terms of ζ(r) given by (12), finding that the value of δ c is varying in terms of α (see section IV for more details). In the particular case of α = 1 we get δ c 0.55 which corresponds to a value of A 0.80. In the left frame of figure 1 we have plotted the critical profile of ζ(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 δ 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 δ c .
In figure 2 we plot the ζ-profiles (plotted in the left panel) and the corresponding density contrast (right panel) measured at horizon crossing, defined by (8), for the threshold values δ c associated to each shape. Although the ζ-profiles given by (12) are always centrally peaked, the energy density profile is centrally peaked only for α ≤ 1. In particular, only the case of α = 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 α the peak is diverging. Nevertheless the amplitude of the perturbation, measured by the averaged value δ m always remains finite.
The right panel of figure 2 also shows that for α > 1 the density contrast is off-centred, with an increasing value of the peak for larger values of α. One can see therefore that a ζ-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 (4) and the correspondence among the peaks is guaranteed only at the linear level. In general, the correspondence between the peaks in ζ and the peaks in δρ/ρ b requires the assumption that ζ (r) = 0 at r = 0, such that in the centre the only non-linear effect is given by the exponential term exp (−2ζ(r)) which reduces the amplitude.
A second issue already mentioned is that finite peaks of ζ do not always correspond to finite values of the peak of the density contrast, which happens here for α < 1. This is because the ζprofiles for α < 1 are not smooth in the centre (they are not infinitely differentiable), and the term ζ (r)/r diverges in the limit r → 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 → 0.
Because in practice a finite value of , 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 = ( √ 3H) −1 , are not completely negligible even at the initial time. For spiky shapes (which have ζ (r)/r = 0 in the centre), like those obtained from (12) with α < 1 and the effect of the transfer function might significantly change the amplitude of the peak, while the value of the averaged δ m hardly changes.
Because for every ζ-peak with a finite amplitude there is always a peak of the compaction function C evaluated at r m , with finite amplitude equal to δ m , we have used the averaged amplitude δ m to calculate the abundance of PBHs (see section V), 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.

III. 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 The left plot shows the critical profiles ζ α (r) given by (13) (13) plotted against r/r m : for α < 1 the profile has a spiky shape, with increasing steepness for decreasing value of α; for α = 1 the profiles is centrally peaked with (δρ) /ρ b → 0 when r → 0 while for α > 1 the profile has an increasing off-peak. In both plots the dashed line corresponds to the smooth centrally peaked case with α = 1.
estimate for the threshold δ c ∼ ω = 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][32][33][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 δρ/ρ b should be used rather than a metric perturbation such as the curvature perturbation ζ. 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 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 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 (equation (3)), 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 top-hat smoothed density contrast varies from 0.41 to 2/3 -a much smaller range of values. It was also discussed in section II that, for certain ζ-profiles, the peak in the density contrast may be off-centred 3 (when α > 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 • 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 II). 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 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].

A. Numerical scheme
The calculations made in this paper to calculate the threshold of PBH formation for the different initial ζ-profiles described in Section II have been made with the same code as used in [29,[31][32][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 II, specified on a space-like slice at constant initial cosmic time t i defined as a(t i )r m e ζ(rm) = 10/H, ( = 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 (δ − δ c ).

B. Threshold, scaling law and mass spectrum
In the left panel of Figure 3 we plot the threshold values δ c against the parameter α used in (12) to vary the initial curvature profile. The lowest limit δ c 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 (α 1) and a smooth tail (r 0 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 δ c = 2/3, (−r m ζ (r m ) = 1), corresponds instead to shapes with r m = r 0 , (α → ∞), 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 range of values we have been able to compute here are 0.442 δ c 0.656, (0.34 ≤ α ≤ 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 δ c 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 where M H is the mass of the cosmological horizon at horizon crossing, γ 0.36 is the critical exponent depending only on the value of ω of the equation of state and K is a parameter depending on the particular shape of the density contrast. Because this parameter will play some role in the next sectionto determine the abundance of PBHs, we have performed numerical simulations to quantify its variation, finding that it varies between 3 and 11 for α varying from 0.4 and 1.9 (corresponding to δ c varying between 0.45 and 0.65) for the profiles considered here, with a representative value of K 4 when δ c 0.55, (α = 1 in (12)).

V. 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 ζ . 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 ζ taking into account the non-linear relation between ζ and δ 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 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][46][47][48] unless the inflaton field has a non-canonical kinetic term [49][50][51][52].
The density contrast δρ/ρ b is related to the curvature perturbation ζ as in equation (4). However, the key parameter to use for PBH formation is instead the smoothed density contrast δ m . Using a top-hat window function with areal radius R = a(t) exp(ζ(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] δ m = − 2 3 r m ζ (r m ) 2 + r m ζ (r m ) .
Because ζ has a Gaussian distribution, its derivative will also have a Gaussian distribution. Therefore, equation (17) The probability density function (PDF) of δ l then follows a Gaussian distribution δ l represents the linear component of the smoothed density contrast and its variance σ 2 can be calculated by integrating the linear component of the density power spectrum as follows: whereW (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 δl (and T (k, r m )) should be evaluated 5 . For simplicity here, we assume that the relevant scale for PBH formation is given by kr m 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 bỹ 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 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 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 − 4 3 rmζ (rm) to the linear component of δr and the derivation of its variance are discussed in more detail in Appendix A. 5 We note, whilst the initial conditions of the perturbations are defined in the super-horizon regime, as in (4), 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 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 end, equation (18) is solved with δ R = δ c to find the critical amplitude of the linear component where there are 2 solutions because equation (18) is quadratic. When necessary, we will take δ c = 0.55 in this paper, the critical value of the volume averaged density contrast when ζ is taken to have a Gaussian profile shape. However, we note that only the first solution δ 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 δ c < δ < 2/3, where 2/3 is the maximum value for the density contrast given by equation (18). This corresponds to δ c,l− < δ l < 4/3.
The final mass of a PBH, M PBH , which forms during radiation domination depends on the shape and amplitude δ R of the initial perturbation, and the horizon mass at horizon reentry M H , where K depends on the profile shape and typically takes a value between 3 and 5. During radiation domination γ 0.36, which is the value we will take [58]. The horizon mass M H is proportional to the horizon scale squared r 2 m , For a random Gaussian field, the number density of sufficiently rare peaks in a comoving volume where σ is given by equation (20), and µ is given by Here, the application of equation (26) 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 (17) 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 (26) and spherical symmetry) -discussed further in appendix B.
The fraction of the energy of the universe at a peak of given height ν which collapses to form PBHs, β ν , then depends on the mass of the PBHs relative to the horizon scale and the number density of the peaks: where the time dependance of the horizon mass M H (r m ) is parameterised by the horizon scale r m .
θ(ν − ν c ) is the Heaviside step function which indicates that no PBH will form if ν 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 ν which forms PBHs, where ν c− ≡ δ c,l− /σ, see equation (23), and equations (24), (26), and (27) have been explicitly substituted into (29). 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]) where M H is the horizon mass at the time of formation and M eq is the horizon mass at the time of matter-radiation equality. We will use the value M eq = 2.8×10 17 M (using the same parameters as [60]). M min and M max are respectively the smallest and largest horizon masses at which PBHs form 6 .
The derivation of this formula assumes that the matter density Ω 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 PBH ) can then be obtained by differentiating Ω PBH with respect to the PBH mass: where the equations (20), (28) and (30) should be recast in terms of the horizon mass and PBH mass using the substitutions in equations (24) and (25). Ω 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 ζ = A s = 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 ζ = A s δ 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 (20) and (28) give scale invariant values of σ 2 ∼ 1.06A ∫ and (µ/(aH)) 2 ∼ 6.86A s respectively. Figure 4 shows the abundance of PBHs, β (equation (30)), as a function of A s , whilst figure 4 shows the mass function f (M ) (equation (32)) evaluated at M = M 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 , with corresponding horizon scale k * = k , 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 β corresponds to the total energy fraction collapsing into PBHs at all epochs, rather than integrating over lnM H as in equation (31). Whilst equation (20) will give a significant variance σ 2 for values of r m close to r (suggesting perturbations of varying scales exist), this is because δ 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 enters the horizon, equations (20) and (28) give σ 2 = (µ/(aH)) 2 ∼ 0.151A 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 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 β as when the linear relation is used. The dominant term in the calculation for the abundance of black holes, equation (30), is the exponential term exp(−ν 2 ). After integrating, this will give a dependance roughly proportional to exp(−ν 2 c− ). In order to produce (approximately) the same number of PBHs from the linear calculation as from the non-linear calculation, we therefore require ν c,L = ν c,N L -where the subscripts L and N L denote the linear and non-linear calculations respectively, where δ c,l− is given by equation (23). In both the linear and non-linear case, σ 2 for a given power spectrum shape is proportional to A s by the same constant of proportionality, σ 2 = CA s . Finally, we can write down the ratio between A L and A N L as For values 0.41 < δ c < 2/3, this factor varies approximately from 1.5 to 4. For more typical values 0.5 < δ 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 ζ, with δ 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 ζ and δ.

VI. 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 δ m , rather than the amplitude of the curvature perturbation ζ. The non-linear relation between the curvature perturbation ζ and the density contrast δρ/ρ b , given by equation (4), means that δ m will have a non-Gaussian distribution even when ζ 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 δ m should be used rather than the value δρ/ρ 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 (5)) is maximised, and note that at this scale C(r) is equal to δ m .
In this paper, we use the amplitude of δ 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 δρ/ρ 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 IV, we have considered different profiles of ζ which form density perturbations to determine a threshold value for PBH formation and how this can depend on the shape of the profile in ζ. We noted in section II that the relation between peaks in ζ and δρ/ρ b is not straightforward and may not coincide. For a given profile shape of ζ, the profile of δρ/ρ b depends also on the amplitude of the perturbation, and we have calculated a relation between the amplitude of the perturbation δ 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 δ m and ζ, given by equation (17), which can be used to fully describe the non-Gaussianity of δ m when ζ 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 β, the present-time density parameter for PBHs, Ω P BH , or the mass fraction of dark matter contained within PBHs of a given mass, f (M P BH ). 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 δ 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 ζ and δ, we find that the amplitude of the power spectrum must therefore be ∼ 2 − 3 times larger, with the amount depending on the value taken for the collapse threshold, see (34), 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 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 δ m and it does not generate a bispectrum because the derivatives of ζ 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 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 NL < 0 of ζ, see [61] for details.
The relation used between ζ and δ, equation (4), is exact in ζ, only neglecting higher-order terms in , valid because 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 ζ and δ makes the formation of PBHs less likely.
For typical values of δ c ∼ 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 ∼ 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 δρ/ρ b and ζ leads to a suppression in the PBH formation rate. Unlike these two papers, we take critical collapse into account.
with the window function W (x, r m ) given by: For our purposes, it is more convenient to express equation (A4), a convolution in real space, as a multiplication in Fourier space, where the Fourier transform of the window functionW (k, r m ) is given by equation (21) and δ(k) is given by where the linear transfer function T (k, r m ) is given by equation (22) and δ m (k) can therefore be written as Since we have assumed ζ(x) has a Gaussian distribution, it also has a Gaussian distribution in Fourier space: ζ(k) has a Gaussian distribution (being a linear combination of Gaussian variables ζ(x)). δ m (k) is then related to ζ(k) by a linear factor, meaning that δ m (k) also has a Gaussian distribution. Finally, δ m (x) is a linear combination of the Gaussian Fourier modes; hence it also has a Gaussian distribution.
Finally, we can calculate the variance σ 2 by integrating the power spectrum, 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 (12) Figure 2 shows the critical profiles of peaks in ζ and δρ/ρ b (in spherical symmetry). All of the profiles in ζ 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 δρ/ρ b . As mentioned in section III 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 → 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 ζ correspond to peaks in the smoothed density field δ 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 (ζ or δ l ) are very unlikely to be close to other large peaks, that is, that the local region surrouding 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 preservedimplying that the perturbation in the non-linear field must also be a peak (when smoothed on an appropriate scale, see above).
The correlation of peaks in ζ and the density δρ/ρ 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 (20). The smoothed density fields (linear and non-linear) therefore only inherit peaks from the ζ-field on a small range of scales. This means that, whilst a perturbation of scale r m in the ζ-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.