PBH abundance from random Gaussian curvature perturbations and a local density threshold

The production rate of primordial black holes is often calculated by considering a nearly Gaussian distribution of cosmological perturbations, and assuming that black holes will form in regions where the amplitude of such perturbations exceeds a certain threshold. A threshold $\zeta_{\rm th}$ for the curvature perturbation is somewhat inappropriate for this purpose, because it depends significantly on environmental effects, not essential to the local dynamics. By contrast, a threshold $\delta_{\rm th}$ for the density perturbation at horizon crossing seems to provide a more robust criterion. On the other hand, the density perturbation is known to be bounded above by a maximum limit $\delta_{\rm max}$, and given that $\delta_{\rm th}$ is comparable to $\delta_{\rm max}$, the density perturbation will be far from Gaussian near or above the threshold. In this paper, we provide a new plausible estimate for the primordial black hole abundance based on peak theory. In our approach, we assume that the curvature perturbation is given as a random Gaussian field with the power spectrum characterized by a single scale, while an optimized criterion for PBH formation is imposed, based on the locally averaged density perturbation. Both variables are related by the full nonlinear expression derived in the long-wavelength approximation of general relativity. We do not introduce a window function, and the scale of the inhomogeneity is introduced as a random variable in the peak theory. We find that the mass spectrum is shifted to larger mass scales by one order of magnitude or so, compared to a conventional calculation. The abundance of PBHs becomes significantly larger than the conventional one, by many orders of magnitude, mainly due to the optimized criterion for PBH formation and the removal of the suppresion associated with a window function.

The production rate of primordial black holes is often calculated by considering a nearly Gaussian distribution of cosmological perturbations, and assuming that black holes will form in regions where the amplitude of such perturbations exceeds a certain threshold.A threshold ζ th for the curvature perturbation is somewhat inappropriate for this purpose, because it depends significantly on environmental effects, not essential to the local dynamics.By contrast, a threshold δ th for the density perturbation at horizon crossing seems to provide a more robust criterion.On the other hand, the density perturbation is known to be bounded above by a maximum limit δ max at the horizon entry, and given that δ th is comparable to δ max , the density perturbation will be far from Gaussian near or above the threshold.In this paper, we provide a new plausible estimate for the primordial black hole abundance based on peak theory.In our approach, we assume that the curvature perturbation is given as a random Gaussian field with the power spectrum characterized by a single scale, while an optimized criterion for PBH formation is imposed, based on the locally averaged density perturbation around the nearly spherically symmetric high peaks.Both variables are related by the full nonlinear expression derived in the long-wavelength approximation of general relativity.We do not introduce a window function which is usually introduced to obtain the scale dependence of the spectrum.The scale of the inhomogeneity is introduced as a random variable in the peak theory, and the scale dependent PBH fraction is automatically induced.We find that the mass spectrum is shifted to larger mass scales by one order of magnitude or so, compared to a conventional calculation.The abundance of PBHs becomes significantly larger than the conventional one, by many orders of magnitude, mainly due to the optimized criterion for PBH formation and the removal of the suppression associated with a window function.

I. INTRODUCTION
Processes which may lead to the formation of primordial black holes (PBHs), along with their cosmological implications, have been extensively investigated in the literature since the pioneering work of Zel'dovich and Novikov [1] and Hawking [2].PBHs may be produced by gravitational collapse in regions with a large amplitude of density perturbations in the early Universe, and measurements or constraints on their abundance can be regarded as a probe of the primordial spectrum and the underlying inflationary model.The latest observational constraints are summarized in, e.g.Refs.[3,4].So far, PBHs have been actively studied as candidates of dark matter (e.g., see Refs.[4][5][6][7][8][9][10][11][12][13][14][15] and references therein).In addition, the recent discovery of gravitational waves emitted from binary black holes (BBHs) [16,17] has stimulated the investigation of PBH binaries and their merger rates [18][19][20][21].
In this paper, we will focus on the formation of PBHs in the radiation dominated era (see Refs. [22][23][24] for cases of the matter dominated era), due to some enhanced feature in the primordial power spectrum of density perturbations around some specific scale. 1A rough criterion for the formation of PBHs was first proposed by Carr [6], and much numerical work has been done in search of a more accurate criterion [28][29][30][31][32][33][34][35][36][37][38].The perturbation variables which are used to characterize the amplitude of the initial inhomogeneity are roughly divided into two sorts: the density perturbation and the curvature perturbation.For instance, Shibata and Sasaki [30] discussed the threshold for PBH formation by using the curvature variable which is given by the conformal factor of the spatial metric.On the other hand, in Refs.[32,33,[35][36][37][38], the threshold value is given for the averaged density perturbation at the horizon entry in the comoving slicing, and in the lowest order long-wavelength expansion.The threshold value of the density perturbation is given by δ th ≈ 0.42 − 0.66 depending on the perturbation profile.The lowest threshold value seems to correspond to the value analytically derived in Ref. [39] with significant simplification [40].
As for the curvature variable, it has been suggested in Ref. [41] that the threshold is strongly affected by environmental effects, while that of a density perturbation is not.This fact has been also numerically demonstrated in Ref. [42].One extreme example which shows the significance of the environmental effects is the estimate of the threshold of the curvature perturbation suggested in Ref. [39].There, an (irrelevant) extremely low value of ζ th ≃ 0.0862 is obtained, due to the existence of an unphysical negative density region in the environment in the specific model adopted there.However, even if we keep the positivity of the density in the environment, the threshold value of the curvature perturbation can be significantly affected [42].This can be intuitively understood if we consider the curvature perturbation as the general relativistic counterpart of the Newtonian potential, which can be shifted by an arbitrary constant.Since, at least in spherically symmetric systems, the process of PBH formation can be described in a quasi-local manner, the use of a threshold value for a quasi-local perturbation variable seems to be more appropriate (see Sec.VII and VIII in Ref. [42] for details).
A useful criterion has been proposed in Ref. [30] for spherically symmetric systems based on the so-called compaction function, which is equivalent to one half of the volume average of the density perturbation δ at the time of horizon entry [42].The criterion for PBH formation is that the maximum value of the compaction function as a function of the averaging radius lies above a specified threshold C th = δ th /2, at the time when the averaging radius enters the horizon.Such threshold value has been found to be more robust.This has been tested by considering two different families of profiles for the perturbation, and a broad range of parameters [42].In what follows, we will not further discuss the possible profile dependence of the threshold, but simply assume the existence of a typical value(see e.g.Ref. [40] for an analysis about the profile dependence).We also note that, although our framework is applicable to generic non-spherical systems, we will adopt a criterion for PBH formation by referring to the compaction function in the corresponding spherical system.This is justified because high peaks of a random Gaussian field tend to be spherical.
The main purpose of this paper is to find an estimate for the abundance of PBHs once a threshold value of the averaged density perturbation is provided.One conventional way is to apply the Press-Schechter(PS) formalism to the density perturbation by assuming that this variable is Gaussian distributed.However, due to the local dynamics of overdense regions, there is an upper limit for the value of the density perturbation at horizon crossing.This was first observed in Refs.[39,43], in the context of a simplified "3 zone model" where a spherical homogeneous overdensity is embedded in a flat Friedmann-Lemaître-Robertson-Walker(FLRW) environment.More generally, it was found [42] that for spherically symmetric perturbations with any profile, the maximum density perturbation at horizon crossing in the co-moving slicing is bounded by δ max ≈ 2/3.The argument will be reviewed in Section II. 2 Noting that δ th is in fact comparable to the maximum value δ max (above which the probability distribution should vanish) a naive application of the Gaussian distribution seems rather questionable.In addition, while the threshold is often set for the density perturbation, the statistical properties of the primordial curvature perturbation are usually better understood.Therefore, it is natural to consider the abundance of PBHs by combining the threshold of the density perturbation with the statistical properties of the curvature perturbation.Since PBH formation is a non-linear process, a non-linear relation between these perturbation variables should be taken into account.In this paper, we address the calculation of the PBH abundance by using the peak theory for the Gaussian probability distribution of the curvature perturbation, and the non-linear relation between curvature and density perturbation in the long-wavelength limit.Readers not interested in the details of the derivation can skip directly to Eq. (61), and the ensuing explanation on how to use it.
Another significant problem is the window function dependence of the mass spectrum.In Ref. [44], it is reported that the mass spectrum significantly depends on the choice of the window function in the PS formalism.For an extended curvature power spectrum, the scale dependence of the PBH fraction is introduced by a window function in the PS formalism.Therefore, without a window function, an extended mass spectrum of PBH cannot be obtained along the conventional PS formalism.In contrast, according to the peak theory, the scale of the inhomogeneity can be also introduced as a random variable.The probability distribution of the random variable is associated with the power spectrum of the curvature perturbation.Then, in our procedure, the scale dependence is naturally induced depending on the profile of the curvature power spectrum.
The plan of the paper is the following.In Section II we discuss the perturbation variables and the implementation of a threshold for PBH formation.In Section III we consider the statistical properties of the Taylor expansion coefficients of the Gaussian random field ζ around a given point.In Section IV, where we discuss the probability for PBH formation, based on the peak theory and on our implementation of the threshold on the averaged density perturbation.The results will be compared to the more conventional PS approach for the monochromatic and a simple extended curvature power spectrum in Section V. Our conclusions are summarized in Section V. Some technical aspects are discussed in the Appendices.Throughout this paper, we use the geometrized units in which both the speed of light and Newton's gravitational constant are unity, c = G = 1.

II. PERTURBATION VARIABLES AND THRESHOLD FOR PBH FORMATION
Let us consider the density perturbation in the comoving slicing, which is orthogonal to the fluid world line.In the long-wavelength approximation, the curvature perturbation ζ and the density perturbation δ with the comoving slicing are related by [42], where w is the equation of state parameter, a is the scale factor, H is the Hubble rate, △ is the Laplacian of the reference flat metric, and the spatial metric is given by with det γ being the same as the determinant of the reference flat metric.
We will be interested mainly in high peaks, which tend to be nearly spherically symmetric [45].Therefore, in this section, we introduce the criterion for PBH formation assuming spherical symmetry originally proposed in Ref. [30].Here, we basically follow and refer to the discussions and calculation in Ref. [42].
First, let us define the compaction function C as where R is the areal radius at the radius r, and δM is the excess of the Misner-Sharp mass M MS enclosed by the sphere of the radius r compared with the mass M F inside the sphere in the fiducial flat FLRW universe with the same areal radius.The Misner-Sharp mass for the comoving slicing is given by where ρ is the matter density and the prime " ′ " denotes the derivative with respect to the radial coordinate(see Sec.V-B in Ref. [42] for general slicing).Then, we find the following expressions for each variable: In the limit of the long-wavelength approximation, for comoving and constant mean curvature gauge, the compaction function is given by: where δ is the volume average of the density perturbation δ within the radius r (see Eq. (5.31) in Ref. [42]).From the definition of C, we can derive the following simple form in the comoving slicing (see also Eq. (6.33) in Ref. [42]): From this expression, it is clear that C ≤ 1/3.If we identify the time of horizon entry of the perturbation from the condition HR = 1, then, form Eq. ( 8), the corresponding averaged density perturbation is δ < δ max = 2/3, as discussed in the introduction.We can rewrite C(r) as 1 − R ′2 e 2ζ /a 2 /3.The existence of the region R ′ < 0 corresponds to the Type II PBH formation reported in Ref. [43].In what follows, we focus on the Type I cases, that is, R ′ > 0.
We will also assume that the function C is a smooth function of r for r > 0.Then, the value of C takes the maximum value C max at r m which satisfies3 We consider the following criterion for PBH formation: In the constant-mean-curvature(CMC) slice, the threshold C CMC th for PBH formation is evaluated as ≃ 0.4 ± 0.03 (see Figs. 2 and 3 or TABLE I and II in Ref. [42]).This threshold corresponds to the perturbation profiles of Refs.[30,33,35], and is found to be quite robust for a broad range of parameters.Since the relation between the density perturbation in the comoving slice (δ) and the CMC slice (δ CMC ) is given by the threshold value in the comoving slice is given by C th ≃ 0.267 which corresponds to δ th ≃ 0.533.In this paper we shall use this as a reference value.We should keep in mind, however, that the threshold value is not completely independent of profile.For instance, as mentioned in the introduction, the threshold in a 3-zone model with a homogeneous overdensity the threshold is somewhat lower.Throughout this paper, we assume the random Gaussian distribution of ζ with its power spectrum P(k) defined by the following equation: where ζ(k) is the Fourier transform of ζ and the bracket < ... > denotes the ensemble average.Each gradient moment σ n can be calculated by Focusing on a high peak and taking it as the origin of the coordinates, we introduce the amplitude µ and the curvature scale 1/k * of the peak as follows: Peaks of ζ do not necessarily correspond to peaks of δ.Nevertheless, as is shown in Appendix A, if the value of δ is comparable to the threshold value δ th at a peak, we can almost always find the associated peak of ζ well inside the horizon patch centered at the peak of δ.
In some cases, there may be a local minimum of δ where there is a maximum of ζ, although this typically happens for values of ζ above the threshold for black hole formation.This is discussed in Appendix D, where we also give the relation between k * and the inverse length scale k δ ≡ (−△δ/δ| x=0 ) 1/2 of the density perturbation at the peak.According to the peak theory [45], for a high peak, we may expect the typical form of the profile ζ can be described by using µ, k * and the two point correlation function ψ as follows: where It is worthy of note that, for k * = k c := σ 1 /σ 0 , we obtain It will be shown that regarding k * as a probability variable, we obtain k c as the mean value of k * .The profile ( 21) is introduced in the recent paper [47] as a typical profile associated with the curvature power spectrum.Here, we also take k * dependence into account, and introduce the scale dependence to the profile.Applying Eq. ( 9) to ζ, we obtain the relation between µ and C as where the smaller root is taken.Let us define the threshold value µ where rm (k * ) is the value of r m for ζ = ζ, and In Eq. ( 23), we have explicitly denoted the k * dependence of rm and g m to emphasize it.
Although we obtain the threshold of the amplitude µ as a function of k * , since our goal is to obtain the mass spectrum, we need the threshold value as a function of the PBH mass M. For this purpose, let us consider the horizon entry condition.In Eq. ( 11), we have implicitly used Eq. ( 8) with the horizon entry condition We note that this coincides with the condition 2M F (r m e − ζ(rm) ) = H −1 .Since the PBH mass is given by M = α/(2H) with α being a numerical factor, from the horizon entry condition (25), the PBH mass M can be expressed as follows: where we have used the fact H ∝ a −2 and a = a 2 eq H eq rm e −µgm with a eq and H eq being the scale factor and Hubble expansion rate at the matter radiation equality.M eq and k eq are defined by M eq = αH −1 eq /2 and k eq = a eq H eq , respectively.We have also introduced the function M (µ,k * ) (µ, k * ).The value of the numerical factor α is rather ambiguous, and we set α = 1 as a fiducial value 4 .Then, we may obtain the threshold value of µ While, from Eq. ( 26), we can describe µ as a function of M and k * as follows: As is explicitly shown for an extended power spectrum, the value of µ may be bounded below by µ min (M) for a fixed value of M.Then, for a fixed value of M, the region of µ for PBH formation can be given by

III. RANDOM GAUSSIAN DISTRIBUTION OF ζ
A key assumption is the random Gaussian distribution of ζ with its power spectrum P(k).In this section, we briefly review Ref. [45] to introduce the probability distribution for the curvature variables.Due to the random Gaussian assumption, the probability distribution of any set of linear combination of the variable ζ(x i ) is given by a multi-dimensional Gaussian probability distribution [45,54]: where the components of the matrix M are given by the correlation < V I V J > defined by Here, we consider the value of ζ and its derivatives up to the second order of the Taylor expansion of the field ζ(x i ): 4 If we take into account the critical behavior [48,49] near the PBH formation threshold, α should be given by a function of µ and k * as α = K(k * )(µ − µ th (k * )) γ with γ ≃ 0.36 [31,32,34,35,47,[50][51][52][53].Since the profile depenence of the function K(k * ) and the parameter range of the scaling behavior is not well understood, for simplicity, we just treat α as a numerical factor which takes a typical value of order 1 in this paper.
The non-zero correlations between two of ζ 0 , ζ i 1 and ζ ii 2 are given by Let us focus on the variables ζ ij 2 .There are 6 independent variables 2 ).It can be shown that, taking the principal direction of the matrix ζ ij 2 , the volume element can be rewritten as follows: where λ i are eigen values of the matrix ζ ij 2 with λ 1 ≥ λ 2 ≥ λ 3 and θ i are the Euler angles to take the principal direction.From the integration with respect to the Euler angles, the factor 2π 2 arises.
Following Ref. [45], we introduce new variables ν, η i and ξ i as follows: λ i is described in terms of ξ i as follows: Then, the probability distribution can be expressed as where and ξ 2 ≥ ξ 3 ≥ −ξ 2 and ξ 2 ≥ 0. We have abbreviated the three components variables ξ i and η i as ξ and η.In terms of µ = νσ 0 and k 2 * = ξ 1 σ 2 /µ, the probability P 1 can be expressed as where k 2 c = γσ 2 /σ 0 = σ 2 1 /σ 2 0 , and IV. PBH FRACTION TO THE TOTAL DENSITY

A. General expression
Following Ref. [45], we start by deriving an expression for the peak number density.The probability distribution can be written as Let us focus on the parameters ν and ξ to characterize each extremum.We define n ext (x, ν, ξ 1 ) as the distribution of extrema of the field ζ in the space of (x, ν, ξ 1 ), that is, n ext (x, ν, ξ 1 )∆x∆ν∆ξ 1 = number of extrema in the volume ∆x∆ν∆ξ 1 .
Then, n ext (x, ν, ξ 1 ) can be expressed as follows: where we have expressed the variables at each extremum with the subscript p.Then, x p is the position of the extremum, that is, ζ 1 = η = 0 at x = x p .Therefore, we obtain where The peak number density n pk (ν, ξ 1 ) is given by the ensemble average of n ext Θ(λ 3 ) as follows: where Θ(λ 3 ) is multiplied to pick peaks out of extrema, and the function f is given by We note that, due to the condition λ 3 > 0, we obtain ζ 2 = ξ 1 σ 2 > 0. Let us change the variables from ν = −ζ 0 /σ 0 = µ/σ 0 and ξ 1 = △ζ| r=0 /σ 2 = µk 2 * /σ 2 to variables µ and k * .Then, we obtain Since the direct observable is not k * but the PBH mass M, we further change the variable from k * to M as follows: where k * should be regarded as a function of µ and M given by solving Eq. ( 26) for k * .Here, we note that an extended power spectrum is implicitly assumed in the above expression.The monochromatic spectrum case will be independently discussed in Sec.V A. We also note that, since we relate k * to M with µ fixed, we have implicitly assumed that there is only one peak with △ζ = µk 2 * in the region corresponding to the mass M, that is, inside r = r m .If the spectrum is broad enough or has multiple peaks at far separated scales, and the typical PBH mass is relatively larger than the minimum scale given by the spectrum, we would find multiple peaks inside r = r m .Then, the PBH abundance would be overestimated because we count every peak as a candidate for a PBH formation.In order to avoid this difficulty, we simply assume that the power spectrum is characterized by a single scale k 0 and has a localized peak around the scale k 0 .Therefore, the current version of our procedure cannot be directly applied to a spectrum with a broad support or multiple scales.
The number density of PBHs is given by We also note that the scale factor a is a function of M as a = 2M 1/2 M 1/2 eq k eq /α.Then, the fraction of PBHs to the total density β 0 d ln M can be given by Here we note again that k * should be regarded as a function of µ and M. The above formula can be evaluated in principle once the form of the power spectrum is given.The PBH fraction to the total density f 0 at the equality time is given by f 0 = β 0 (M eq /M) 1/2 .Let us summarize how to use the above formula.Once the power spectrum characterized by a single scale k 0 is given, the typical profile is given by Eq. (17).Taking the radius of the maximum compaction function (9) for this profile, we obtain the function rm (k * ) and g m (k * ) = g(r m (k * ); k * ).Since k * is implicitly given as a function of µ and M as Eq. ( 26), we can express the integrand in Eq. (61) as a function of µ and M. The lower boundary µ b of the integration is given by Eq. ( 29), where µ (M ) th (M) is implicitly given by eliminating k * from Eqs. ( 23) and ( 26) as is defined in Eq. (27).
From the functional form (49) of P 1 , we may expect that the integrand of Eq. ( 61) has a non-negligible value only around k = k c and µ = µ b .Assuming µ b ≫ σ, we can obtain the following approximate form without performing the integral: Again, k * should be regarded as a function of µ and M in the above expression.Let us roughly estimate the typical width ∆ ln M in the mass spectrum from the above approximate expression.From the horizon entry condition and a rough dimensional analysis, we may estimate the relation between M and k * as M ∼ H −1 ∝ k −2 * .Then, we find ∆ ln M ∼ ∆ ln k 2 * σ 2 /k 2 c ∼ σ 0 .The probability P 1 takes the maximum value at k * = k c for a fixed value of µ, and the threshold value of µ for k * = k c is given by µ c := µ (k * ) th (k c ).Then, the value of the mass at the peak probability is given by where g c := g m (k c ).The value of β 0 can be evaluated as where ℓ := rm (k c )k c , and we have evaluated as [45] in the first line.A numerical factor of the order of unity is neglected in the second line, and the values of moments are assumed to be σ 0 ∼ σ 1 /k c ∼ σ 2 /k 2 c .Let us summarize how to use the above approximate formula.Once the power spectrum is given, the typical profile with k * = k c = σ 1 /σ 0 is given by its two point correlation function ψ(r).Calculating the compaction function and taking the radius of the maximum C for the typical profile, we can obtain the value of ℓ.The threshold value µ c can be evaluated by the formula (23) with k * being k c , where the value of δ th should be provided.Then, the simplified version of the PBH fraction (64) at the spectrum peak can be calculated.Therefore, necessary ingredients are the power spectrum and the values of δ th and α.

B. Estimation from the Press-Schechter formalism
For a comparison, we review a conventional estimate of the fraction of PBHs based on the PS formalism.In the conventional formalism, the scale dependence is introduced by a window function W (k/k M ), where Then, each gradient moment is replaced by the following expression: The conventional estimate starts from the following Gaussian distribution assumption for the density perturbation δ: where σ δ is given by the coarse-grained density contrast We note that the definition of δ and δ th used in this conventional estimate are rather vague and not necessarily identical to our definition of δ and δ th .Therefore, there is an ambiguity of which definition of the density perturbation is supposed in this formalism.Here, for simplicity, we just use the same numerical value of δ th as in our approach, in other words, we assume that the volume average of the density perturbation obeys the Gaussian probability distribution given by Eq. ( 68) with the coarse-grained density contrast (69) in the PS formalism.This Gaussian distribution and the dispersion are motivated by the linear relation between ζ and δ.The fraction β PS 0 is then evaluated as follows(see e.g.[4]): The PBH fraction to the total density f PS 0,δ at the equality time is given by f PS 0,δ = β PS 0,δ (M eq /M) 1/2 .Let us consider another way of estimation as a reference.First, referring to the linear relation δ = 4△ζ/(9k 2 * ) = 4µ/9 at the conventional horizon entry aH = k * , we change the variable for the Gaussian distribution Eq. ( 68) to µ as As in the case of β PS 0,δ , we can evaluate the PBH fraction β PS 0,µ as follows: where the threshold value µ c = µ ) is evaluated by the procedure introduced in Sec.II, where the non-linearity is taken into account.Therefore, comparing Eq. (70) and Eq. ( 72), we can extract the effect of the optimized criterion proposed in Sec.II.The PBH fraction to the total density f PS 0,µ at the equality time is given by f PS 0,µ = β PS 0,µ (M eq /M) 1/2 .As will be shown later, we obtain µ c ∼ 0.52 for the monochromatic spectrum and ∼ 0.75 for a specific model of extended spectra.Therefore, the value of µ c is typically smaller than 9δ th /4 ≈ 1.20.This simple analysis clearly indicates that the optimized criterion given in Sec.II will significantly increase the PBH fraction compared to the conventional estimate (70).

V. SPECIFIC EXAMPLES A. Monochromatic power spectrum
Let us consider the monochromatic power spectrum given by Then, the moments are calculated as It leads to k c = k 0 and γ = 1.Since k c = k 0 , from Eq. ( 21), the functional form of g(r; k 0 ) is given by Then, we can find rm (k 0 ) = ℓ/k 0 = 2.74/k 0 , µ c = 0.520 and g c = −0.141.Since the value of γ is given by 1. Taking the limit γ → 1, in the expression (46), we obtain Then, the k * integration in Eq. ( 58) can be performed, and we obtain the following expression for the peak number density: Under the condition k * = k 0 , the relation between M and µ is given by Therefore, we obtain the following PBH number density: where Θ (M − M c ) has been multiplied to extract the distribution above the threshold.Finally, we obtain At M = M c = M eq k 2 eq ℓ 2 k −2 0 e −2µcgc , the value of β 0 is given by Since the function f (x) behaves like f (x) ∼ x 3 for x ≫ 1 [45], in the limit of small σ 0 , we obtain where a numerical factor of the order of unity is neglected.The mass spectra β 0 and f 0 are depicted as functions of the PBH mass M for σ 0 = 0.08, 0.06 and 0.05 in the left panels of Figs. 1 and 2, while β 0 and f 0 = (M eq /M) 1/2 β 0 are depicted as functions of σ 0 for M = M c in the right panels of Figs. 1 and 2. As a result, we obtain an extended mass spectrum of PBHs even for the monochromatic power spectrum of the curvature perturbation.This is caused by the non-linear relation between the density perturbation and the curvature perturbation through Eq. ( 26) 6 .However, the width of the spectrum is hardly significant.
It is worthy to be mentioned that comparing Eq. (65) with Eq. (82), we find the factor 1/σ 0 of difference.This clearly shows that, the monochromatic spectrum case gets larger peak amplitude in the mass spectrum instead of the loss of the mass spectrum width.Note that the amplitude of the mass spectrum is huge compared to the conventional one β PS 0,δ for a small value of σ 0 .The main reason for this deviation comes from the optimization of the PBH formation threshold µ c .The σ −4 0 dependence of the prefactor in Eq. ( 82) also contributes to increase the fraction, but not so dramatically.We note that the strong enhancement in the abundance of PBH is a robust feature, for any chosen value of the threshold δ th .Indeed, it follows from (23) that for δ th < δ max = 2/3, we have µ c < (9/4)δ th .According to Eqs. ( 70) and (72), this implies β PS 0,µ (M) ≫ β PS 0,δ (M).In fact, for δ th < 0.623, we have µ c < (9/8)δ th (the right hand side of this inequality is actually a good linear fit to the value of µ c , which typically overestimates the actual value by less than 30%).Assuming that the probability of PBH formation is low, we can approximate the complementary error function as erfc(x) ≈ e −x 2 /( √ πx), and it follows that within this range of δ th we have where we have ignored the sub-exponential dependence.Given that the probability of PBH formation is exponentially small, this represents a very strong enhancement.Furthermore, there is a non-trivial correction to the expression for the mass in terms of the wave number at the time of horizon crossing, which increases the mass of the black holes by one order of magnitude or so relative to the expectation from linear theory.This is clear from the expression of M c given in Eq. ( 63), which contains the factor ℓ 2 e −2µcgc ≈ 8.7.Here, we have used the numerical values for µ c , ℓ and g c corresponding to the sinc profile, which is the appropriate one for the monochromatic spectrum.These values are listed in Table I of Appendix C. Hence, not only do we find more PBHs, but they are also significantly larger than naively expected.

B. An extended power spectrum
Let us consider the simple extended power spectrum given by Gradient moments are calculated as where Γ means the gamma function 7 .The functional form of ψ(r) is given by  9) for k * = 0.1k 0 , k 0 , 1.5k 0 and 2k 0 are shown in Fig. 3.The value of rm is analytically calculated as and is depicted as a function of k * in Fig. 4. The functional form of µ th defined in Eq. ( 23) value of µ min is given by Substituting the function k th * (M) into Eq.( 23), we can draw µ (M ) th as a function of M as is shown in Fig. 7.In the small mass region, namely large k * region, the second term in Eq. ( 17) dominates the profile.Then, the magnitude of k * degenerates with the overall amplitude of the inhomogeneity.Consequently, PBH can form even for very small value of µ if k * is sufficiently large.However, as will be shown below, the probability of such cases are exponentially suppressed.
Using all the functional forms shown above, we can numerically evaluate the integral (61).The results are shown in Figs. 8 and 9.
We note that the mass spectrum from the conventional PS formalism is much smaller for the same value of σ.There are two reasons for this difference.One is the smaller threshold value 0.753 < 9δ th /4 ≈ 1.2 as is also discussed for the monochromatic power spectrum in Sec.V A. In addition, for the extended spectrum case, we also find that the value of β PS 0,µ (M) ≫ β PS 0,δ (M) is also much smaller than our estimate if we use the following Gaussian window function in the PS formalism: This is because of the difference between the typical dispersions σ2 /k 2 * and σ.In Fig. 10, we plot σ and σ2 /k 2 * for the Gaussian, real-space top-hat, k-space top-hat window functions.The real-space and k-space top-hat window functions are defined as follows: Fig. 10 shows that the value of σ is larger than σ2 /k 2 * with the Gaussian window function.This difference is reflected to the many orders of magnitude difference in the PBH fraction.

C. Behavior near the spectrum peak
In order to obtain a simpler analytic formula, let us consider the typical profile with the value of k * being its mean value k c .Replacing k 2 * by k 2 c in the expression of g, we obtain Typical value g c of g m can be given by In order to simplify the k * dependence of rm , we assume the following simple form Then, from Eq. ( 26), the relation between k * and M can be approximated by For Eq. ( 62), let us perform the following replacement: The second replacement implies µ b = max {µ min (M), µ c }.In addition, since σ n ≪ 1, we can also approximate the function f as [45] Then, finally, we obtain the following analytic expression: It should be noted that, while the same peak value as Eq. ( 64) is given by the above simple formula, the curvature of the spectrum peak cannot be well approximated due to the nontrivial scale dependence of µ (k * ) th (see Fig. 5) which is not taken into account in Eq. (98).In order to take k * dependence of µ where m 0 , m 1 and m 2 are, respectively, given by with the subscript "c" denoting the value at k * = k c .Then, the threshold value µ (k * ) is given by Substituting this expression into the combination µ 2 /σ 2 in the exponential of the probability (49), we find For our specific example, given by Eq. ( 84), we find m 0 = m 2 /k 4 0 = 2/e and m 1 = 0, and the first term in the quadratic coefficient is 3/(2k 4 0 ) while the remaining correction terms are −1/k 4 0 .Therefore, the value of the effective variance increases and the spectrum peak is flatten.For our specific example, we obtain m 1 = 0.However, in general, we may expect non-vanishing value of m 1 and it may cause the shift of the spectrum peak.Defining the modified dispersion σmod as we obtain the following modified simple expression: This modified version of the simple expression gives better approximation for the shape of the spectrum around the peak as is shown in Fig. 11.

VI. SUMMARY AND DISCUSSION
Primordial black holes (PBHs) have attracted much attention not only from a theoretical point of view but also observationally.In order to make an observationally relevant prediction from a theoretical result, the estimation of the abundance of PBHs is essential.
The conventional Press-Schechter(PS) formalism assumes a Gaussian distribution of the primordial density or curvature perturbation.However, the threshold value of the curvature perturbation has an ambiguity from environmental effects [41,42].On the other hand, the existence of an upper limit for the density perturbation at the time of horizon entry [43] is in conflict with the naive assumption of a Gaussian probability distribution for such variable.In order to overcome the above issues, we have developed a formalism to estimate the abundance of PBHs by combining the Gaussian probability distribution of the curvature perturbation with a threshold value for the locally averaged density perturbation.More precisely, we consider an optimized criterion for PBH formation which is given by a threshold for the compaction function proposed in Ref. [30].The compaction function is given by the value of the mass excess above the background divided by the area radius.Our approach is based on the peak theory of Gaussian random fields [45], and takes into consideration the non-linear relation between the curvature perturbation and the density perturbation.A general expression for the fraction of PBHs to the total density of the Universe is presented.Although our current procedure is only directly applicable to a spectrum which has a localized peak around a specific scale, the expression for the PBH fraction, in principle, can be evaluated once the curvature power spectrum is given.In our program to calculate the PBH mass spectrum, there is no need to introduce a window function as far as we consider a power spectrum localized around a single scale in the Fourier space.The scale dependence of the PBH fraction, which is conventionally induced by a window function, is induced by considering a typical peak profile of the curvature perturbation characterized by the amplitude(µ) and the spatial peak curvature(k 2 * ) for a given curvature power spectrum.We note that, for each profile characterized by µ and k * , at the horizon crossing time, the PBH mass is evaluated from the mass inside the radius at which the compaction function takes the maximum value.
The case of the monochromatic power spectrum is particularly simple, and illustrative of the general case.First of all, compared to the conventional PS approach, the PBH spectrum is shifted to larger masses, by an order of magnitude or so.This is due to the mass estimation from the mass inside the maximum compaction radius, which is substantially affected by the metric perturbation.A related effect is that there is a slight spread in the mass spectrum, even when the underlying primordial spectrum is monochromatic.Such spread, however, where C is the integration constant.This constant is related to ζ 0 as e −2ζ 0 = 1/C 2 and can be absorbed into the renormalized scale factor ā as in the text or equivalently into the rescaling of the radial coordinate.Therefore, hereafter, we set C = 1 for simplicity.The value of △ζ at the origin can be evaluated as where we have identified it as µk 2 * as in the text.From the Hubble equation, the density perturbation δ CMC in the uniform Hubble slice is given by (B13) This quantity takes the maximum value 1 at µ = 6/q.For a small radius, we may approximate r * ≃ r * , and comparing the horizon entry conditions (B11) and ( 25), we can estimate the value of q as ∼ 10(see TableI in AppendixC).This relation between δ CMC H and µ explicitly shows that the value of the density perturbation at the horizon entry time has a maximum value while the value of µ is not bounded.The case µ = 6/q corresponds to the 3-hemisphere which divides the Type I and Type II PBH formation.

FIG. 1 .
FIG. 1. PBH fraction to the total density β 0 at the formation time is depicted as a function of the PBH mass M (left panel) for each value of σ 0 , and a function of σ 0 (right panel) for M = M c with k 0 = 10 5 k eq and α = 1.We also depict the conventional estimation from the PS formalism without a window function in the right panel.

FIG. 4 .
FIG. 4. rm and the radius of the inner peak of C(r).

FIG. 9 .
FIG.9.PBH fraction to the total density f 0 at the equality time is depicted as a function of the PBH mass M (left panel) for each value of σ 0 , and a function of σ 0 (right panel) for M = M c with k 0 = 10 5 k eq and α = 1.We also depict the value at the peak of the mass spectrum from the PS formalism with the Gaussian window function in the right panel.

k 2 * = k 2 c
(k * ) th into account, let us consider the Taylor expansion of µ (k * ) th around k * = k c .First, we focus on the function m(r, k 2 * ) defined by m(r, k 2 * ) := rg ′ (r; k * )., we may consider the following expansion of m(r m (k * ), k 2 * ):

From 1
the metric forms, we obtain the following relations: dr √

) 2 *
Let us define the radius r * corresponding to the scale 1/k * as r q being a numerical factor.Then, we identify the horizon entry by 1 qµ) 2 .