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

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
1) Max Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, 85748 Garching bei Muenchen, Germany

2) Institut de Ciències del Cosmos, Universitat de Barcelona,
Martí i Franquès 1, 08028 Barcelona, Spain

3) Laboratoire Univers et Théories, UMR 8102 CNRS, Observatoire de Paris, Université Paris Diderot, 5 Place Jules Janssen, F-92190 Meudon, France

4) Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, United Kingdom
September 20, 2019

The formation and abundance of primordial black holes (PBHs) arising from the curvature perturbation is studied. The non-linear relation between and the density contrast means that, even when 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 must be a factor of 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.

I 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 Hawking:1971ei ; Carr:1974nx ; Carr:1975 . If a density perturbation is above a threshold , 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 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 Carr:2017jsz ; Kuhnel:2017pwq ; Bellomo:2017zsr 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 Clesse:2017bsw . For a review of PBH formation and constraints, see Green:2014faa ; Sasaki:2018dmp .

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 () 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 () 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 ( compared to )). 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 Byrnes:2018txb .

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 Drees:2011hb , axion inflation Bugaev:2013fya ; Ozsoy:2018flq , or a waterfall transition during hybrid inflation GarciaBellido:1996qt ; Lyth:2012yp ; Bugaev:2011wy , 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 consequences111Often is used instead to denote the curvature perturbation, whilst is used to describe the curvature perturbation only on a uniform density slicing.. The curvature perturbation appears in the metric in the comoving uniform-density gauge as


where is the conformal time, is the scale factor and 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 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 Shibata:1999zs , which was then used to derive a relation between the abundance of PBHs and the power spectrum Green:2004wb . At the same time Niemeyer and Jedamzik performed a numerical study of PBH formation using instead an initial perturbation of the energy density. For a long time the abundance of PBHs was calculated assuming that regions where the curvature perturbation was above a critical value 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 Young:2014ana . The effect of super-horizon modes on PBH formation is discussed in detail in Harada:2015yda . 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 Young:2014ana ) 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 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 Germani:2018jgr for example),


where is the equation of state parameter during the radiation dominated epoch of the early universe, is the background density and 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 , the aims of this paper are to investigate how the fully non-linear relation between and 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 and the abundance of PBHs at formation .

The paper is organised as follows: in section II we will discuss the set up of the initial conditions in the density contrast 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 before analysing a specific parametrization of that allows us to vary the profile of . This allows us, with the help of numerical simulations (see Section IV), to span almost all the possible range of values of . 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 Bardeen:1985tr , and the perturbation profile is therefore defined using only a radial coordinate .

Note that the curvature perturbation in the literature is typically defined on a uniform density slicing, whilst the density contrast 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 Lyth:2004gb and the references therein).

ii.1 Gradient expansion approach

In the super-horizon regime perturbations have a length scale much larger than the cosmological Hubble horizon. This approach, usually called the gradient expansion or long-wavelength approximation Shibata:1999zs ; Salopek:1990jq is based on expanding the exact solution in a power series of a small parameter () that is conveniently identified with the ratio between the Hubble radius and the length scale characterising the perturbation


where a particular value of corresponds to a particular value of the time .

When the curvature profile 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 , 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 . The curvature profile represents the fundamental source of the perturbation, embedded into the metric (1) from the asymptotic limit, . Although this approach reproduces the time evolution of linear perturbation theory when , it also allows one to consider non-linear curvature perturbations if the spacetime is sufficiently smooth on scales greater than Lyth:2004gb . 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 defined on a comoving uniform-cosmic time slicing and the curvature perturbation is given by:


where 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 and the quadratic term of the first derivative , where the prime denotes a derivative with respect to the radial coordinate . 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 Musco:2018rwt , 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 is the areal radius, is the Misner-Sharp mass within a sphere of the radius and is the background mass within the same areal radius calculated with respect a FRW Universe. In the superhorizon regime, applying the gradient expansion approximation, the compaction function is conserved and is related to as


We can then compute the length-scale of the perturbation, identified as the location where the compaction function is maximized (), which gives


Measuring the curvature perturbation with introduces an intrinsic rescaling of the comoving coordinate with respect to the background FRW 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 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 from here onwards, is given by


The perturbation amplitude can be measured as the mass excess of the energy density within the scale , measured at the horizon crossing time . 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 Musco:2018rwt for more details). Computing the mass excess as the integral of the density contrast averaged over the background volume , the amplitude of the perturbation is given by


and using the explicit expression of in terms of the curvature profile seen in (4), we get , which satisfies the fundamental relation Musco:2018rwt


for any curvature profile, .

ii.2 Initial conditions

We will now study the main feature of the density profile when an explicit parameterization of the curvature profile is specified as


where and respectively denote the amplitude and the scale of the perturbation. Inserting this into (4), the corresponding density profile is given by


and then using (12) into (6) and (7) we can calculate the corresponding perturbation amplitude as


which gives the value of in terms of the averaged amplitude


This behaviour of in terms of also shows that there is a maximum value of corresponding to which represents the transition between PBHs of type I and type II after which formation of PBHs cannot be studied in terms of but only in terms of Kopp:2010sh .

Figure 1: The left plots shows the critical curvature profiles given by (12) for ( and ), while the right plot shows the corresponding profile. The upper limit represents the maximum theoretical upper limit of .

Using numerical simulations we have calculated the critical values for PBHs using initial condition in terms of given by (12), finding that the value of is varying in terms of (see section IV for more details). In the particular case of we get which corresponds to a value of . In the left frame of figure 1 we have plotted the critical profile of as function of identifying the critical peak of the profile. In this plot we are also showing the value of corresponding to the maximum limit of . In the right hand plot of figure 1 we show the corresponding compactness function with the peak amplitude of being equal to .

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 associated to each shape. Although the -profiles given by (12) are always centrally peaked, the energy density profile is centrally peaked only for . In particular, only the case of , plotted with a dashed line, is smooth at , 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 remain always finite.

The right panel of figure 2 also shows that for 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 requires the assumption that at , such that in the centre the only non-linear effect is given by the exponential term 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 . This is because the -profiles for are not smooth in the centre (they are not infinitely differentiable), and the term diverges in the limit . 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 . Because in practice a finite value of , corresponding to a finite initial time , 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 , are not completely negligible even at the initial time. For spiky shapes (which have in the centre), like those obtained from (12) with and the effect of the transfer function might significantly change the amplitude of the peak, while the value of the averaged hardly changes.

Because for every -peak with a finite amplitude there is always a peak of the compactness function evaluated at , with finite amplitude equal to , we have used the averaged amplitude 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 , leaving a study of the effects of the non-linear transfer function to future work.

Figure 2: The left plot shows the critical profiles given by (13) plotted against for where for increasing values of the central value is decreasing. The right plot shows the corresponding critical profiles of given by equation (13) plotted against : for the profile has a spiky shape, with increasing steepness for decreasing value of ; for the profiles is centrally peaked with when while for the profile has an increasing off-peak. In both plots the dashed line corresponds to the smooth centrally peaked case with .

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 Carr:1975 , using a Jeans length argument, provided an order of magnitude estimate for the threshold at horizon crossing for a radiation fluid. Since then, there has been extensive work to determine the collapse threshold Shibata:1999zs ; Niemeyer:1999ak ; Musco:2004ak ; Polnarev:2006aa ; Musco:2012au ; Nakama:2013ica ; Harada:2015yda ; Musco:2018rwt , as well as discussions about the appropriate parameter to use to determine whether a perturbation will collapse Green:2004wb ; Young:2014ana ; Germani:2018jgr ; Yoo:2018esr . 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 Harada:2013epa .

As discussed previously, in order to determine a clear criterion to distinguish which perturbations are able to form a PBH, the density contrast 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 Bardeen:1985tr - although non-spherical symmetry has been considered Harada:2015ewt . In this section, we will discuss several ambiguities within the literature over how the criteria for collapse should be defined222A full description and investigation of these considerations is beyond the scope of this paper, although will be discussed further in a paper by the authors currently in preparation.:

  • 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 Musco:2018rwt . 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 Germani:2018jgr 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 for collapse of the top-hat smoothed density contrast varies from 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 (meaning the central value is smaller than the peak) or infinite - a problem which is avoided by using the smoothed value.

  • The choice of window function has a significant effect on the calculated abundance of PBHs (as was discussed in Ando:2018qdb ). The threshold value for collapse is typically stated in terms of the volume-averaged density contrast (as for example in Musco:2004ak ; Musco:2018rwt ), 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 but perturbations grow as - meaning the top hat window function is not typically 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 . 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 approximation is not very accurate for the large perturbations required to generate PBHs).

  • The scale of a perturbation is best stated in terms of the radius at which the compactness function is maximised (see Section II). This is different to the previously used definition Musco:2004ak , where the scale of the perturbation was defined as the radius at which the density contrast becomes negative. As was shown in Musco:2018rwt computing the density contrast at does not give a sensible parameter to compare different shapes. The averaged value of the density contrast evaluated at is characterised by the general relation given by (11): 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 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 , corresponding to a top-hat window function with radius 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).

Iv Numerical results of PBH formation

iv.1 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 Musco:2004ak ; Polnarev:2006aa ; Musco:2008hv ; Musco:2012au ; Musco:2018rwt . 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 defined as , (), while the outer edge of the grid has been placed at , 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 ().

iv.2 Threshold, scaling law and mass spectrum

In the left panel of Figure 3 we plot the threshold values against the parameter used in (12) to vary the initial curvature profile. The lowest limit corresponds to the analytic solution derived in Harada:2013epa where the gravitational effect of pressure was taken into account while pressure gradients were instead neglected. It was shown in Musco:2018rwt that this corresponds to shapes of the density contrast with a very large peak () and a smooth tail (): 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 , (), corresponds instead to shapes with , (), 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 as a function of the corresponding behaviour of : the range of values we have been able to compute here are , (). 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 Musco:2018rwt where different parameterizations of the curvature profiles has been considered, using more than one parameter, getting much closer the lower limit of .

As has been shown in previous works Niemeyer:1997mt ; Niemeyer:1999ak ; Musco:2004ak ; Musco:2008hv ; Musco:2012au the mass spectrum of PBHs follows the critical collapse, characterized by a scaling law relation


where is the mass of the cosmological horizon at horizon crossing, is the critical exponent depending only on the value of of the equation of state and 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 nuerical simulations to quantify its variation, finding that it varies between and for varying between and for the profiles considered here, with a representative value of when , ( in (12)).

Figure 3: The left plot show the value of as function of obtained with numerical simulations using an initial curvature profile given by (12), while the right panel show the behaviour of as function of the corresponding behaviour if . 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 shows to the largest possible value of , with shapes characterized by very large pressure gradients at the scale .

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 . 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 taking into account the non-linear relation between and 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 will generate an unacceptably large dark-matter isocurvature perturbation in the CMB Young:2015kda ; Tada:2015noa - although such bounds can be evaded if non-Gaussianity exists 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 Maldacena:2002vr ) does not generate isocurvature perturbations Pajer:2013ana ; Cabass:2016cgp . Even when taking ultra slow roll inflation into account, it remains uncertain whether the non-Gaussianity can have a relevant effect Bravo:2017wyw ; Bravo:2017gct ; Atal:2018neu ; Passaglia:2018ixg unless the inflaton field has a non-canonical kinetic term Shandera:2012ke ; Young:2015cyn ; Franciolini:2018vbk ; Kamenshchik:2018sig .

The density contrast 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 . Using a top-hat window function with areal radius , this is related to the curvature perturbation in radiation domination as Musco:2018rwt


Because has a Gaussian distribution, its derivative will also have a Gaussian distribution. Therefore, equation (17) can be expressed in terms of a linear Gaussian component as333It is noted that a recent paper Kawasaki:2019mbl 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 to the linear component of and the derivation of its variance are discussed in more detail in Appendix A.


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


represents the linear component of the smoothed density contrast and its variance can be calculated by integrating the linear component of the density power spectrum as follows:


where is the Fourier transform of the top-hat smoothing function, is the linear transfer function, and the smoothing scale is equal to the horizon scale. The horizon scale is used to define the time at which (and ) should be evaluated. For simplicity here, we assume that the relevant scale for PBH formation is given by , although this is not a very accurate approximation, and the exact relation between and depends on the profile of the density contrast, which depends on the shape of the power spectrum Germani:2018jgr . The Fourier transform of the top-hat smoothing function is given by


and the linear transfer function, where we consider as a time dependent measure of the horizon, is given by Josan:2009qn


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 Byrnes:2012yx ; Young:2013oia . To this end, equation (18) is solved with to find the critical amplitude of the linear component ,


where there are 2 solutions because equation (18) is quadratic. When necessary, we will take 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 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 , where is the maximum value for the density contrast given by equation (18). This corresponds to .

The final mass of a PBH, , which forms during radiation domination depends on the shape and amplitude of the initial perturbation, and the horizon mass at horizon reentry ,


where depends on the profile shape and typically takes a value between 3 and 5, and we will take where a numerical value is needed. During radiation domination , which is the value we will take Evans:1994pj . The horizon mass is proportional to the horizon scale squared .


For a random Gaussian field, the number density of sufficiently rare peaks in a comoving volume is Bardeen:1985tr


where is given by equation (20),


and is given by


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 is parameterised by the horizon scale . 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 , 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) Byrnes:2018clq


where is the horizon mass at the time of formation and is the horizon mass at the time of matter-radiation equality. We will use the value (using the same parameters as Nakama:2016gzw ). and are respectively the smallest and largest horizon masses at which PBHs form444In principle can be arbitrarily large because can be arbitrarily small. However, the largest PBH mass which can be formed when the horizon mass is is given by .. The derivation of this formula assumes that the matter density 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 can then be obtained by differentiating 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). 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 . For the narrow peak, we will take the peak to be narrow enough such that it can be treated as a Dirac-delta function, , although note that such a power spectrum is unphysical (for example, Byrnes:2018txb describes that the power spectrum cannot be steeper than ).

For the scale invariant case, equations (20) and (28) give scale invariant values of and respectively. Figure 4 shows the abundance of PBHs, (equation (30)), as a function of , whilst figure 5 shows the mass function (equation (32)) evaluated at as a function of the power spectrum amplitude.

Figure 4: The energy fraction of the universe collapsing to form PBHs at the time of formation, , is plotted against the amplitude of a scale invariant power spectrum, . The critical value has been taken as . Accounting for the non-linear relation means that a significantly smaller abundance of PBHs is calculated, by many orders of magnitude.
Figure 5: The mass function of PBHs of PBHs (defined in equation (32)) with a mass of as a function of the amplitude of the (scale invariant) power spectrum . The critical value has been taken as .

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 , with corresponding horizon scale , 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 as in equation (31). Whilst equation (20) will give a significant variance for values of close to (suggesting perturbations of varying scales exist), this is because 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 enters the horizon, equations (20) and (28) give .

Figure 6: The energy fraction of the universe collapsing to form PBHs at the time of formation, , is plotted against the amplitude of the power spectrum, where a Dirac delta power spectrum has been assumed, . The critical value has been taken as .
Figure 7: The mass function of PBHs of PBHs with a mass of as a function of the amplitude of the power spectrum, where a Dirac delta power spectrum has been assumed, . is defined in equation (32). The critical value has been taken as .

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 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 Germani:2018jgr .

However, we will here make a simple calculation to describe by what fraction the amplitude of the curvature perturbation power spectrum 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 . After integrating, this will give a dependance roughly proportional to . In order to produce (approximately) the same number of PBHs from the linear calculation as from the non-linear calculation, we therefore require - where the subscripts and denote the linear and non-linear calculations respectively,


where is given by equation (23). In both the linear and non-linear case, for a given power spectrum shape is proportional to by the same constant of proportionality, . Finally, we can write down the ratio between and as


For values , this factor varies approximately from 1.5 to 4. This means that, in order for the same number of PBHs to form, the amplitude of the power spectrum must be roughly 2–3 times greater than previously calculated with the linear model. In particular, for the commonly considered case of a Gaussian profile in , with , 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 , rather than the amplitude of the curvature perturbation . The non-linear relation between the curvature perturbation and the density contrast , given by equation (4), means that 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 should be used rather than the central peak value . The scale over which the perturbation should be averaged, , is the scale at which the compactness function (defined in equation (5)) is maximised, and note that at this scale is equal to . In this paper, we use the amplitude of 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 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 is not straightforward and may not coincide. For a given profile shape of , the profile of depends also on the amplitude of the perturbation, and we have calculated a relation between the amplitude of the perturbation 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 and , given by equation (17), which can be used to fully describe the non-Gaussianity of 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, , or the mass fraction of dark matter contained within PBHs of a given mass, . 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 is taken into account, we find that this always reduces the number of PBHs which form by many orders of magnitude, see figures 4-7. 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 the linear model, we find that the amplitude of the power spectrum must therefore be 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) 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 and it does not generate a bispectrum because the derivatives of are uncorrelated on scales much larger than the PBH scale, . 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 Tada:2015noa ; Young:2015kda ). This also implies that the “apparent” negative local cannot lead to an enhanced formation of PBHs in the case of a primordial power spectrum with a broad peak, in contrast to a the a physical value of of , see Young:2014oea for details.

Note added: While this paper was approaching completion a related work by Kawasaki and Nakatsuka Kawasaki:2019mbl 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 , we confirm their finding that the suppression of PBH formation due to the non-linear terms requires the power spectrum to have an amplitude 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. toni 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 and leads to a suppression in the PBH formation rate. Unlike these two papers, we take critical collapse into account.


We thank V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto and C. Ünal, the authors of toni , 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.

Appendix A Gaussianity and variance of

Here we will derive the variance and discuss the Gaussianity of the linear term in equation (17), . The linear relation between the density contrast and the curvature perturbation in radial coordinates is given by


where the prime denotes a derivative with respect to the radial co-ordinate , and the parameter describes the super-horizon time-evolution of the perturbation and is given in the linear approximation by , where is the perturbation scale and is the horizon scale. In this paper we consider the perturbations at the time of horizon reentry, , such that , and we will therefore drop the . The volume averaged density contrast used to determine PBH formation is given by


where in the linear approximation the volume is given by . Although this integral is normally integrated over the comoving distance , this difference represents higher-order effects neglected in the linear approximation. Combining the above equations gives


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

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


where is the location of a peak in cartesian coordinates (corresponding to ) and is the Heaviside step function. The second equality is the same integral expressed in cartesian coordinates, with the window function given by:


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


where the Fourier transform of the window function is given by equation (21) and is given by


where the linear transfer function is given by equation (22) and can therefore be written as


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

Finally, we can calculate the variance by integrating the power spectrum,



Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description