Application of peaks theory to the abundance of primordial black holes

We consider the application of peaks theory to the calculation of the number density of peaks relevant for primordial black hole (PBH) formation. For PBHs, the final mass is related to the amplitude and scale of the perturbation from which it forms, where the scale is defined as the scale at which the compaction function peaks. We therefore extend peaks theory to calculate not only the abundance of peaks of a given amplitude, but peaks of a given amplitude and scale. A simple fitting formula is given in the high-peak limit relevant for PBH formation. We also adapt the calculation to use a Gaussian smoothing function, ensuring convergence regardless of the choice of power spectrum.


I. INTRODUCTION
Primordial black holes (PBHs) are black holes which may have formed in the early universe.
Whilst no observations have been confirmed, there are several hints towards their existence [1].
They represent a viable dark matter candidate, and are a unique probe to constraint the small scale early universe. There are several different formation mechanisms, but we will focus here on PBHs formed from the collapse of large density perturbations. Shortly after the end of inflation, during the radiation dominated epoch of the universe, the cosmological horizon grows and perturbations which were super-horizon cross the horizon and can collapse to form a PBH.
In order to form a PBH, a density perturbation at horizon crossing must have an amplitude above some threshold, δ c , for gravity forces to overcome pressure forces and collapse. The density contrast δ(t, x) is the relative over-density, stated in the comoving synchronous gauge, where ρ is the energy density, and the bar denotes the background value for a flat universe,ρ(t) = 3H 2 (t)/(8π). H is the Hubble parameter, and we are using natural units c = G = 1.
The abundance of PBHs is then typically found by calculating the abundance of perturbations above this threshold value (see [2][3][4][5][6][7][8][9][10] amongst others). The simplest way of determining the abundance of such perturbations is by using a Press-Schechter-like) formalism, which has a simple condition that the density must be above the threshold. Peaks theory [11] introduces a further constraint, stating that compact objects (such as galaxies, or in our case, PBHs) form at peaks of the density. Applying this peak constraint to the condition for a perturbation to form a PBH gives us δ In the case of large-structure, it is also commonly assumed that compact objects also form at positions where the density is above a certain threshold, and the mass of the object can be determined purely by its scale -defined as the largest smoothing scale at which the perturbation is above the threshold value. However, this is not the case for PBH formation, because a larger amplitude perturbation will pull in more of the surrounding material as a PBH forms, resulting in a larger PBH mass [9,[12][13][14]. The mass of the PBH depends on both the amplitude and scale of the perturbation from which it formed (discussed further in section II).
To calculate the abundance of perturbations, we will therefore introduce another constraint requiring that perturbations have a specific scale. Previous calculations (i.e. [2,5,9,15,16]) utilising peaks theory have made the assumption that, when smoothed on a scale R, all density perturbations are exactly of this scale.
In this paper, we will extend the peaks theory approach used in previous papers to account for the fact that both scale and amplitude must be accounted for. For simplicity, we will assume Gaussian statistics, although, as has been pointed out in several recent papers related to PBH abundance [8,9,16], the density contrast δ will not be Gaussian even if the curvature perturbation ζ is. However, it is expected that large peaks in the linear, Gaussian density field can be identified with peaks in the non-linear, non-Gaussian density [16] -and so this could easily be accounted for following the methods of [9]. Primordial non-Gaussianity has also been shown to have a significant effect on PBH abundance [3,[17][18][19], which we will not consider here.
We will begin by considering how the mass of a PBH depends on the perturbation from which it forms, before deriving an expression for the number density of peaks. Finally, we will describe how this can be used to determine the abundance and mass function of PBHs.

II. PRIMORDIAL BLACK HOLE MASS
We will describe how the mass of a PBH may be determined from the initial perturbation from which it forms. Firstly, it is necessary to define the compaction function, Throughout this paper, we will consider all perturbations in the super-horizon regime [20], where the time dependence of H cancels exactly with the time dependence of δ. The compaction function C can be considered as equivalent to the volume-averaged time-independent component of the density contrast -which is typically referred to as δ R in the literature. This allows the threshold for collapse to be stated in terms of the compaction, C c . Note that this is different from the usual definition of the compaction function, equivalent to using a Gaussian smoothing window rather than a top-hat one (see [20] for more discussion of the use of this function, and why there is a factor of 2 missing from the exponential in usual definition of a Gaussian function). This is to avoid divergences in the calculation, which would later appear (in the calculation of the variances).
For a perturbation centered at x, the compaction function peaks at some scale R, and the perturbation length r m is defined as where the prime denotes a derivative with respect to R. For the same perturbation, the Gaussian window function returns slightly smaller values for r m and C than the standard top-hat. The mass of the PBH is then given by the critical-scaling relationship [9,[12][13][14], This follows the same form as the standard result utilising a top-hat window function, with different values for the constants: K ≈ 10, C c ≈ 0.25 (instead of 4 and 0.55 respectively [9]), and γ ≈ 0.36 (the details of this are given in appendix A). The horizon mass M H depends on r m as where r eq and M eq are the horizon scale and mass at the time of matter-radiation equality, respectively, where we have assumed radiation domination from PBH formation until the time of equality.  forms a vanishingly small PBH, as it has exactly the threshold density, and follows the self-similar, critical collapse described in [14]. Care, therefore, must be taken that when PBH abundance is calculated, the number density of peaks of different scales must be evaluated at the correct smoothing scale.

A. Variables
Let us now turn our attention to deriving an expression for the number density of peaks of a given height and scale. In order to derive this expression, we will introduce a large number of variables, although the final expression is much simpler. The Fourier transform of C, denoted by a hat, is where the second equality is included because this is the same expression as for the smoothed density-contrast (multiplied by R 2 H 2 which cancels out the time-dependence). The Fourier transforms of the relevant derivatives of C are We can see that C is a linear combination of C and ∇ 2 C (and C of C, ∇ 2 C and ∇ 2 ∇ 2 C). Each of these has the same form, a function of k and R, multiplied by δ R (k): The correlators of these variables can then be determined by integrating over the dimensionless and we will use the notation C A C A = σ 2 A . We will also introduce the cross-correlation coefficient Because of rotational invariance, one then gets where all cross-correlators with an odd number of gradients vanish, while where A is any of the scalars C, C , C or ∇ 2 C. These relations imply that and similar expressions hold for σ 2 RR ≡ (C ) 2 and σ 2 1R ≡ ∇C · ∇C . As we will follow the derivation used in [21], we will follow the notation used there as closely as possible for ease of reference. The variables are therefore expressed in terms of variables normalised to have unit variance: In the high-peak limit relevant for PBH formation, we will make the assumption ζ 0i 1, justified by the fact that, for large peaks, the physical location of a peak in C is not expected to move significantly under a very small change in the smoothing scale, and therefore neglect this term when it appears. Whilst ν and ζ 00 are already rotationally invariant, we can define additional rotationally invariant quantities whereζ ij ≡ ζ ij − δ ij J 1 /3. The independent, rotationally-invariant quantities will be collectively referred to as which is the same list of variables which appeared in [21], with one extra scalar variable, ζ 00 .
From equation (7), we can see that for a Gaussian smoothing window η 0 is not an independent variable but a linear combination of ν and J 1 , B. The peak constraint and number density of peaks Accounting for the fact that we wish to determine peaks of a given scale, the peak constraint becomes where D is the n-dimensional Dirac-delta function, θ H is the Heaviside step function, and || . . . || stands for the absolute value of the determinant of the Hessian matrix. This is the Jacobian determinant of the coordinate transformation from C and ∇δ (the Gaussian variables of the constraint) to the physical coordinates R and x.
This function returns the number of peaks in an infinitesimal volume d 3 x with scale between R and R + dR and height between ν and ν + dν, divided by dRdνd 3 x. The determinant in the above equation can be factorized as In the high peak limit we are interested in, the scalars ζ 00 and J 1 = −tr(ζ) are both very large, since they correlate with ν and ν =ν 1. Conversely, the traceless partζ ij and the vector ζ 0i remain of order O(1). Hence, the 3-D Hessian matrix ζ ij can be approximated by −(J 1 /3)δ ij , and the term ζ 0i ζ −1 ij ζ 0j is negligible (it is suppressed twice) compared to ζ 00 . Ignoring the φ i terms, the peak constraint then becomes Technically, we should worry about the PBHs which form from such peaks being inside the radius of larger PBHs which form later -the so-called "cloud-in-cloud" problem. However, this has been found to have a negligible effect on the PBH abundace, due to the rarity of PBHs and the very small probability of such an occurence [22,23].
The average number density can be expressed by integrating over the probability density function (PDF) of w, P (w),n pk = dwP (w)n pk .
Up until this point, the expression is completely general, and no assumptions have been made about Gaussianity. Assuming that the density contrast (and therefore the compaction) is Gaussian, the PDF can be expressed as [21] P (w)dw = N (ν, J 1 , ζ 00 )dνdJ 1 dζ 00 where N (ν, J 1 , ζ 00 ) is a 3-variate normal distribution, χ 2 n (x) is a χ 2 distribution with n degrees of freedom, Ω is vector of the angular variables, with P (Ω) its PDF. The terms inside the square brackets in the above equation are unchanged from the standard peaks theory calculation, and the integral over the 3η 2 , 5J 2 , J 3 terms and the factor 3 3/2 in equation (20) gives the function f (J 1 ) defined in reference [11], Making use of equation (16), the Dirac-delta function δ D (η 0 ) can be rewritten as and for convenience we will defineJ 1 ≡ 4σ 0ν /(R 2 σ 2 ). Finally, integrating over the remaining Dirac-delta functions and Heaviside step functions yields where N (ζ 00 |ν,J 1 ) ≡ N (ν,J 1 , ζ 00 )/N (ν,J 1 ) is the conditional PDF of ζ 00 givenν andJ 1 . We give its exact expression in Appendix C. Equation (24) can be considered as the main result of this paper. Its prefactor has dimensions of [R −4 ], consistently with the fact that we are computing a differential number density per unit of smoothing scale.
In the high-peak limit relevant for PBH formation, we can approximate f (J 1 ) ≈J 3 1 = 64σ 3 0ν 3 /(R 6 σ 3 2 ). In this same limit, the conditional PDF N (ζ 00 |ν,J 1 ) becomes highly concentrated around its mean value ζ 00 |ν,J 1 , and the integral in equation (24) tends to ανN (ν,J 1 ), where α ≡ ζ 00 |ν,J 1 /ν depends on R and on P(k) but no longer onν. We give the exact expression for α in equation (C7) (though α is close to unity except for narrow power spectra). Using the explicit expression for N (ν,J 1 ), which we derive in detail in Appendix D, we can, to high accuracy, approximate the number density of peaks as where the bars have now been dropped for simplicity, and γ 0,2 is the cross-correlation coefficient of ν and J 1 . Corrections to this result are of the same order of the terms that were already neglected from equation (18). Equation (25) can be compared to the average number density of high peaks, n hi−pk , obtained using un-modified peaks theory: which being a standard number density has dimensions of [R −3 ]. This result crucially differs from ours by one power of ν 1 and by the constant multiplying ν 2 in the exponential.
We also emphasize that our derivation relied on the assumption of Gaussian statistics only for the exact expression of f (J 1 ), equation (22). However, the limit f (J 1 ) J 3 1 holds independently of the distribution, simply because | det(ζ ij )| (J 1 /3) 3 for high peaks, and all non-scalar variables are marginalized over. Therefore, the generic number density of high peaks can be obtained from equation (25) simply by replacing α and N (ν,J 1 ) with the appropriate non-Gaussian values: the abundance of very high peaks can be easily obtained once the PDF of the scalar variables is known.

IV. DISCUSSION
Here, we have considered PBHs forming from the collapse of large-density perturbations upon horizon entry. Such PBHs form with a mass close to the horizon mass at this time -the exact mass depends on both the amplitude and scale R of the perturbation from which a PBH forms. We have therefore extended peaks theory to derive an expression for the number density of peaks of a given height ν and scale R.
Whilst we have here assumed Gaussian statistics for the density contrast δ, it has recently been discussed in a number of papers [8,9,16] that δ will not have a Gaussian distribution even if the curvature perturbation ζ does. However, the easiest way to account for this is to note that high peaks in the smoothed Gaussian field correspond to peaks in the smoothed non-Gaussian field [9,16], although the recent paper [24] may be used to describe the statistics of peak of non-Gaussian fields. In order to calculate the abundance and mass function of PBHs, the expression derived here for the number density of peaks can be substituted into the existing calculation (for example, that given in [9]) in place of the previous expression forn pk (more details for the calculation in the Gaussian case are given in appedix E). Whilst the full derivation is non-trivial and is left for future study, we can make a simple estimate of the impact of the new expression by comparing the expressions for the number density of peaks.
Whether the peak number density is increased or decreased will generally depend on the form of the power spectrum and the scale of perturbations being considered. To make a comparison, we will consider a power spectrum with a lognormal peak (the form of the power spectrum is given in equation (B5)), and consider the number density of peaks on the scale corresponding to the peak of the power spectrum, R = 1/k * (and we will take R = 1/k * = 1 for convenience). Figure 2 shows the number density of peaks changes as a function of ν for different widths ∆ of the lognormal peak. Broader power spectra typically predict a larger number density of peaks for the same ν.
For the same amplitude of the power spectrum A we also expect the variance of perturbations σ 2 0 to be larger, so perturbations of the same amplitude C will have a smaller 'relative' amplitude ν associated with it (also implying a larger number of peaks (and PBHs) would form). shown, for peaks in the power spectrum of different widths. Again, we are considering peaks with a width R corresponding to the scale of the peak of the lognormal power spectrum. For broad peaks in the power spectrum, more peaks are predicted by a factor proportional to ν. However, for narrow peaks, the new expression predicts a significantly smaller value.

NOTE ADDED
As this paper was being written, two similar papers appeared on the arXiv, although following different methodology. Reference [25] appeared first, and provides a similar argument for the inclusion of the condition that only perturbations where the compaction function peaks on a specific scale should be included when one wishes to calculate the mass function of PBHs, but does not make a comparison of the effect of including this term.
Reference [26] appeared shortly thereafter, and performed a fuller calculation. That paper concluded that previous (linear) calculations over-predict the abundance of PBHs for broad power spectra, seemingly in contrast with the results presented here. However, a direct comparison of these conclusions is complicated due to the different methods being used in each paper. We also note that, in this paper, we have utilised a Gaussian smoothing function, which ensures convergence of the variances regardless of the form of the power spectrum.
a Gaussian smoothing function is used instead.
The shape of the primordial perturbations is not known, but the "average" profile shape can be calculated from the statistics of the power spectrum. We will treat the large, rare perturbations which may form PBHs as spherically symmetric [11], and therefore describe the perturbation shape using only the radius from the centre of the perturbation r. For a wide range of power spectrum shapes, it has been shown that the central region of large perturbations (the relevant part for PBH formation) is well approximated by a sinc function [20] δ(r) = A sin(r) r , where we are ignoring the time-dependence of the perturbation, and the radial coordinate r is defined in arbitrary units (i.e. Different scale perturbations can be defined by changing the units).
Using a top-hat or Gaussian smoothing function to calculate r m gives r T H ≈ 2.74 or r G = 2 respectively. The Gaussian function gives a smaller characteristic scale for the perturbation, which then predicts a horizon mass which is smaller by a factor (r G /r T H ) 2 , Likewise, the calculated amplitude of the perturbations (stated in terms of the compaction function at its maximum) will be different depending on the smoothing function used, but are proportional in the ratio C T H ≈ 2.17C G (where the factor 2.17 is simple to find numerically).
Substituting these values into the expression for the PBH mass, equation (A1) gives corresponding to the values given in equation (4).

Appendix B: Correlation factors
In this section, we will further discuss the correlation factors which impact the calculation of peak number density. The correlation factors are normally expressed in terms of moments of the power spectrum, and we include this here for completeness. The n th moment of the power spectrum is defined as where δ R is equivalent to the compaction function, as defined in equation (6). The power spectrum P δ R can be considered equivalent to density power spectrum smoothed on a scale R, and rescaled by a factor R 4 H 4 (the rescaling can, in turn, be considered as a linear extrapolation to the amplitude at horizon entry).
The relevant correlation factors are the correlation factors of ν = C σ 0 , J 1 = C 2 σ 2 and ζ 00 = C RR σ RR . Therefore, it is desirable to express C RR in terms of spatial derivatives of C, which, in Fourier space, gives We can therefore express the variance σ 2 RR as Finally, the correlation factors are given in terms of moments of the power spectrum as, We now turn our attention to the possible range of values for the correlation factors. To achieve this, we will consider a range of shapes for the (smoothed) power spectrum parameterised as where A describes the amplitude, the (kR) 4 term describes the super-horizon growth of density perturbations, and the first exponential term is the smoothing term. The second exponential term describes a lognormal peak, with a peak at k * and width ∆. Note that, due to the k 4 growth at small k and exponential smoothing at large k, the power spectrum will inevitably be relatively narrow, regardless of the choice of ∆. Varying ∆ from 0 to ∞ describes the minimum and maximum allowed widths respectively.

Appendix C: The conditional PDF from the 3-variate normal
In order to gain a better understanding of the integral over the 3-variate normal PDF seen in Section III B, we will briefly discuss the factorization of the conditional PDF. The 3-variate normal PDF (with zero mean) is given by where X = {ν, J 1 , ζ 00 } and Σ is the covariance matrix, whose elements are given by the crosscorrelation coefficients γ AB = C A C B /(σ A σ B ) defined in equation (9). For the three variables involved in the PDF, the relevant coefficients are γ 0,2 , γ 0,RR and γ 2,RR .
To diagonalise the PDF we introduce the following normalized Gaussian variables, which depend linearly on J 1 and ζ 00 respectively, but are independent of ν (since the components that correlate with ν have been subtracted). Similarly, we define the variable z, where the correlation factor γ x,y can be expressed as .

(C4)
The variable z depends linearly on ζ 00 and is independent of both ν and x. The 3 variables ν, x and z are mutually independent and their PDF automatically factorizes. Therefore, the PDF of ν, J 1 and ζ 00 can be expressed as .

(C5)
Let us now consider the integral over the PDF which appears in equation (24). Using equation (C5), then substituting inν,J 1 = (4σ 0 /R 2 σ 2 )ν and ζ 00 into the expressions for ν, x and z gives, after some simple algebra, This is a relatively complicated formula, but in the next section we will discuss the high-peak limit relevant for PBH formation, where the benefit of expressing the PDF in this form will become clear.

Appendix D: The high-peak approximation
The expression for the average peak number-density is given by equation (24), with the function f (J 1 ) defined in equation (22) and an expression for the integral given in equation (C6). However, if ν is large, as is generally the case when considering PBH formation, this can be significantly simplified.
Firstly, the function f (J 1 ) quickly asymptotes to J 3 1 as J 1 becomes large, and so we have Secondly, we will consider the integral over v RR in (C6). In the limitν → ∞, we can take the integral to be an integral over a Dirac-delta function, centred on αν, multiplied by |v RR |.
where the factor α is given in equation (C7), and is typically of order unity, except for very narrow power spectra, where it becomes large. equality). Finally, there are several ways of defining the mass function, involving a derivative with respect to the PBH mass, which is normalised such that d(lnM PBH )f (M PBH ) = Ω PBH /Ω CDM , and which is normalised such that dM PBH ψ (M PBH ) = 1. In the case of a broad (scale-invariant) power spectrum, β is constant over M H , in which case ψ (M PBH ) ∝ M