Abundance of Primordial Black Holes in Peak Theory for an Arbitrary Power Spectrum

We modify the procedure to estimate PBH abundance proposed in arXiv:1805.03946 so that it can be applied to a broad power spectrum such as the scale-invariant flat power spectrum. In the new procedure, we focus on peaks of the Laplacian of the curvature perturbation $\triangle \zeta$ and use the values of $\triangle \zeta$ and $\triangle \triangle \zeta $ at each peak to specify the profile of $\zeta$ as a function of the radial coordinate while the values of $\zeta$ and $\triangle \zeta$ are used in arXiv:1805.03946. The new procedure decouples the larger-scale environmental effect from the estimate of PBH abundance. Because the redundant variance due to the environmental effect is eliminated, we obtain a narrower shape of the mass spectrum compared to the previous procedure in arXiv:1805.03946. Furthermore, the new procedure allows us to estimate PBH abundance for the scale-invariant flat power spectrum by introducing a window function. Although the final result depends on the choice of the window function, we show that the $k$-space tophat window minimizes the extra reduction of the mass spectrum due to the window function. That is, the $k$-space tophat window has the minimum required property in the theoretical PBH estimation. Our procedure makes it possible to calculate the PBH mass spectrum for an arbitrary power spectrum by using a plausible PBH formation criterion with the nonlinear relation taken into account.

We modify the procedure to estimate PBH abundance proposed in Ref. [1] so that it can be applied to a broad power spectrum such as the scale-invariant flat power spectrum.In the new procedure, we focus on peaks of the Laplacian of the curvature perturbation △ζ and use the values of △ζ and △△ζ at each peak to specify the profile of ζ as a function of the radial coordinate while the values of ζ and △ζ are used in Ref. [1].The new procedure decouples the larger-scale environmental effect from the estimate of PBH abundance.Because the redundant variance due to the environmental effect is eliminated, we obtain a narrower shape of the mass spectrum compared to the previous procedure in Ref. [1].Furthermore, the new procedure allows us to estimate PBH abundance for the scale-invariant flat power spectrum by introducing a window function.Although the final result depends on the choice of the window function, we show that the k-space tophat window minimizes the extra reduction of the mass spectrum due to the window function.That is, the k-space tophat window has the minimum required property in the theoretical PBH estimation.Our procedure makes it possible to calculate the PBH mass spectrum for an arbitrary power spectrum by using a plausible PBH formation criterion with the nonlinear relation taken into account.

I. INTRODUCTION
Since Zel'dovich, Novikov and Hawking had pointed out its possibility [2,3], primordial black holes(PBHs) have continued to attract attention.They are still viable candidates for a substantial part of dark matter(see e.g.Refs.[4,5] and references therein) and a possible origin of observed binary black holes [6,7].Mass, spin or spatial distribution of PBHs provides us valuable information about relatively small scale inhomogeneity in the early universe.When we connect a PBH production scenario and observational constraints on it, theoretical estimation of the PBH distribution is inevitable.Here, we provide a plausible procedure to calculate the PBH mass spectrum for an arbitrary power spectrum based on the peak theory.
Until relatively recently, the Press-Schechter (PS) formalism is applied to the estimation of PBH abundance based on a perturbation variable such as the comoving density or the curvature perturbation.As PBHs started to draw more attention after the discovery of the black hole binary as well as gravitational waves, people have begun to seriously doubt the relevance of the PS formalism in the estimation of PBH abundance.In order to improve the estimation, one needs to resolve the following mutually related issues: PBH formation criterion, statistical treatment of non-linear variables [8], and use of a window function [9][10][11].
For the criterion of the PBH formation, there has been a long-term debate since Carr had proposed a rough criterion [12].A lot of efforts to clarify the appropriate criterion have been made through numerical and analytic treatments [13][14][15][16][17][18][19][20][21][22].One useful criterion was proposed in Ref. [15] by using the compaction function, which is equivalent to the half of the volume average of the density perturbation in the long-wavelength limit [23].In Ref. [23], through spherically symmetric numerical simulations, it was shown that the threshold of the maximum value of the compaction function gives relatively accurate criterion, which is within about 10% accuracy for a moderate shape of initial configuration.More recently, the threshold value for the volume average of the compaction function was proposed in Ref. [24], and it was shown that this variable gives the PBH formation criterion within 2% accuracy for a moderate inhomogeneity in the radiation dominated universe(see Ref. [25] for general cosmological backgrounds).These recent developments show that the use of the compaction function is crucial for an accurate estimation of PBH abundance.
Another important ingredient in the calculation of PBH abundance is the statistics of perturbation variables.Naively, we expect that the curvature perturbation would be relevant for the Gaussian distribution assumed in the PS formalism.However, the absolute value of the curvature perturbation does not have any physical meaning in a local sense because it can be absorbed into the coordinate rescaling.Therefore, setting the threshold value for the absolute value of the curvature perturbation seems irrelevant.Conversely, while setting the threshold on the compaction function would be appropriate, the compaction function cannot be a Gaussian variable even if the curvature perturbation is totally Gaussian because of their non-linear relation.Furthermore, difference of the gauge confuses the relation between perturbation variables.In Ref. [23], relations between different gauge conditions were summarized and the gauge issue has been clarified.The compaction function is expressed in terms of the curvature perturbation in the same reference [23].Then, apart from the window function, the remaining issue is how to count the number of PBHs by taking into account the nonlinear relation between the curvature perturbation and the compaction function.Although a few procedures treating the non-linear relation have been proposed [1,26], their consistency with each other has not been clear yet.
In Ref. [1], a plausible procedure to estimate PBH abundance for a narrow power spectrum of the Gaussian curvature perturbation was proposed, where the threshold for the compaction function is used and the non-linear relation is taken into account.However, this procedure cannot be directly applied to a broad spectrum(see Ref. [27] for a simple approach with linear relations).Our aim in this paper is to improve the procedure in Ref. [1] so that we can introduce a window function, and make it possible to apply to any power spectrum.
This paper is organized as follows.First, the criterion based on the compaction function is introduced in Sec.II.In Sec.III, focusing on a high peak of the Laplacian of the curvature perturbation △ζ, we characterize the typical profile of the curvature perturbation ζ around the peak by using the values of △ζ and △△ζ.This treatment allows us to decouple the environmental effect on the absolute value of the curvature perturbation, and the criterion can be expressed in a purely local manner.In Sec.IV, the procedure to estimate PBH abundance is explained, and applied to the single-scale narrow power spectrum previously presented in Ref. [1].We discuss the case of the scale-invariant flat spectrum implementing a window function in Sec.V. Sec.VI is devoted to a summary and discussion.
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. CRITERION FOR PBH FORMATION
Let us consider the spatial metric given by ds 2  3 = a 2 e −2ζ γij dx i dx j with det γ being the same as the determinant of the reference flat metric, where a and ζ are the scale factor of the background universe and the curvature perturbation, respectively.In the long-wavelength approximation, the curvature perturbation ζ and the density perturbation δ with the comoving slicing are related by [23], (2) in the radiation dominated universe, where H is the Hubble expansion rate and △ is the Laplacian of the reference flat metric.
We will be interested mainly in high peaks, which tend to be nearly spherically symmetric [28].Therefore, in this section, we introduce the criterion for PBH formation originally proposed in Ref. [15] assuming spherical symmetry.Here, we basically follow and refer to the discussions and calculation in Ref. [23].
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 enclosed by the sphere of the radius r compared with the mass inside the sphere in the fiducial flat Friedmann-Lemaître-Robertson-Walker universe with the same areal radius.
From the definition of C, we can derive the following simple form in the comoving slicing (see also Eq. ( 6.33) in Ref. [23]): We will 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 satisfies We consider the following criterion for PBH formation: In the comoving slicing, the threshold C th for PBH formation is evaluated as ≃ 0.267 (see Figs. 2 and 3 or TABLE I and II in Ref. [23]).This threshold corresponds to the perturbation profiles of Refs.[15,18,29], and is found to be quite robust for a broad range of parameters(see Ref. [24] for a more robust criterion).In this paper we shall use this value as a reference value.

III. PEAK OF △ζ AND THE SPHERICAL PROFILE
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 Hereafter we suppose that the power spectrum is given.Then the gradient moments can be calculated from the power spectrum and regarded as constants.
In this paper, we focus on high peaks of ζ 2 := △ζ, which coincide with peaks of δ with linear relation.We note that this procedure is different from the previous one proposed in Ref. [1], where peaks of −ζ is considered.
Focusing on a high peak of ζ 2 and taking it as the origin of the coordinates, we introduce the amplitude µ 2 and the curvature scale 1/k • of the peak as follows: According to the peak theory [28], for a high peak, assuming the spherical symmetry, we may expect the typical form of the profile ζ2 can be described by using µ 2 and k • as follows: ζ2 (r) with It is worthy of note that, for k It will be shown in Eq. ( 27) that regarding k • as a probability variable, we obtain σ 3 /σ 2 as the mean value of k • .Let us consider the profile ζ given by integrating Eq. (11).Integrating ζ2 , and assuming the regularity at the center, we obtain ζ(r) where ζ ∞ = ζ| r=∞ is an integration constant.Because we have We may consider either ζ 0 or ζ ∞ as a probability variable.Since the constant shift of ζ can be absorbed into the renormalization of the background scale factor, the non-zero value of ζ ∞ would be regarded as a larger-scale environmental effect.Actually, in AppendixA, we show that the mean value of ζ ∞ is 0 for a given set of µ 2 and k • .In Ref. [1], we used ζ 0 to characterize the profile of the curvature perturbation.However, the use of ζ 0 would mix the environmental effect with the local state.Therefore, in this paper, we ignore ζ ∞ by renormalizing the background scale factor as e −ζ∞ a → a to eliminate the environmental effect, and regard ζ 0 as a dependent variable on µ 2 and k • through Eq. ( 15) with ζ ∞ = 0.In order to obtain PBH abundance, we can follow the procedure proposed in Ref. [1] replacing µ and k * by µ 2 and k • .Here, we just copy the part of the procedure from Ref. [1](the flow char of our procedure can be seen in Ref. [30]).Applying Eq. ( 4) to ζ, we obtain the relation between µ 2 and C as where g(r; k • ) := ζ/µ 2 and the smaller root is taken.Let us define the threshold value µ where rm (k • ) is the value of r m for ζ = ζ, and In Eq. ( 17), we explicitly denoted the k • dependence of rm and g m to emphasize it.
In order to express the threshold value as a function of the PBH mass M, let us consider the horizon entry condition: Since the PBH mass is given by M = α/(2H) with α being a numerical factor, from the horizon entry condition (19), the PBH mass M can be expressed as follows: r2 m e −2µ2 gm =: where we have used the fact H ∝ a −2 and a = a 2 eq H eq rm e −µ 2 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.For simplicity, we set α = 1 as a fiducial value 2 .
Then we can obtain the threshold value of µ While, from Eq. ( 20), we can describe µ 2 as a function of M and k • as follows: The value of µ 2 may be bounded below by µ 2min (M) for a fixed value of M. Actually, in the specific examples in Sec.IV and V, the value of µ 2min (M) is given as follows: Then, for a fixed value of M, the region of µ for PBH formation can be given by IV. PBH ABUNDANCE From Ref. [1], we obtain the expression for the peak number density characterized by µ 2 and k • as where and In Eq. ( 25), the following replacements have been made from Eq. (58) in Ref. [1]: 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 µ 2 and M given by solving Eq. ( 20) for k • .We note that an extended power spectrum is implicitly assumed in the above expression.In the monochromatic spectrum case, the expression reduces to the same expression in Ref. [1] because It should be noted that, since we relate k • to M with µ 2 fixed, we have implicitly assumed that there is only one peak with △ζ 2 = −µ 2 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 we cannot correctly count the number of peaks in the scale of interest.In order to avoid this difficulty, we need to introduce a window function to smooth out the smaller-scale inhomogeneities.We discuss this issue in the subsequent section.In this section, we simply assume that the power spectrum is characterized by a single scale k 0 and there is no contribution from the much smaller scales k ≫ k 0 .
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 µ 2 and M. The above formula can be evaluated in principle once the form of the power spectrum is given.A crucial difference of Eq. ( 32) from Eq. (61) in Ref. [1] is that the expression does not depend on σ 0 , which has IR-log divergence for the flat scale-invariant spectrum.Thus we can consider the PBH mass spectrum for the flat spectrum without introducing IR cut-off. 3 In order to give a simpler approximate form of Eq. ( 32), we approximately perform the integral with respect to µ as follows: 3 The PBH fraction to the total density f 0 at the equality time is given by f 0 = β 0 (M eq /M ) 1/2 .We do not explicitly show the form of f 0 in this paper since the scale dependence can be more easily understood by the form of β 0 .
Since P 1 given in Eq. ( 27) has the exponential dependence, we may expect that the value of β 0 is sensitive to the exponent −µ 2 2 /2σ 2 .Therefore, assuming µ 2b = µ , we can roughly estimate the maximum value of β 0 at the top of the mass spectrum by considering the value k t of k • which minimizes the value of The value of k t cannot be given in an analytic form in general, and a numerical procedure to find the value of k t is needed.We note that the value of k t is independent of the overall constant factor of the power spectrum, and depends only on the profile of the spectrum.Substituting k t into k • in Eq. ( 33), we obtain the following rough estimate for the maximum value β 0,max : where We note that the expression (35) gives a better approximation for a smaller value of the amplitude of the power spectrum σ 2 , and may have a factor of difference from the actual maximum value for σ 0.1 due to the mass dependence of the factors other than the exponential part as can be found in the examples below.
Let us consider the extended power spectrum given by Gradient moments are calculated as where Γ means the gamma function.The result is shown in Fig. 1.Our new procedure gives a narrower and slightly higher spectrum than that obtained in Ref [1].This behavior can be understood as the environmental effect induced by the variance of ζ ∞ .Although the effect is so small that it could be practically ignored in this example, we successfully decoupled the environmental effect.

V. IMPLEMENTING A WINDOW FUNCTION
In our new procedure, a window function can be straightforwardly implemented.That is, introducing the UV cut-off scale k W , instead of Eq. ( 7), we consider the following power FIG. 1. PBH mass spectrum(left) and β approx 0,max as a function of σ(right) for the extended power spectrum (36) with k 0 = 10 5 k eq .The solid lines correspond to the spectra calculated by our new procedure with α = 1, and the dashed lines show the spectra calculated in Ref. [1].We also plot the mass spectrum β PS 0 obtained from the Press-Schechter formalism explained in Appendix for comparison.In the left panel, the dotted horizontal lines show the corresponding values of β approx 0,max .
spectrum of ζ: where Then, following the procedure given in the previous section, we can calculate PBH abundance for a given value of k W .The final PBH mass spectrum is given by the envelope curve of the mass spectra for all values of k W .We note that, for a narrow power spectrum, P W (k) → P(k) in the limit k W → ∞, and the mass spectrum results in the case without the window function irrespective of the choice of the window function.
One important issue here is the window function dependence of the final mass spectrum.In order to clarify this issue, for a sufficiently broad power spectrum, let us consider the effect of the window function for a peak number density at a fixed scale given by the wave number k • = k 0 .If k 0 ≫ k W , no peak can be found.On the other hand, if k 0 ≪ k W , we would find many smaller scale peaks inside the region of the radius ∼ 1/k 0 because the smaller scale modes with k ≫ k 0 are superposed on top of the inhomogeneity with k ∼ k 0 .Thus every peak has a sharp profile due to the superposed small-scale inhomogeneity and satisfies k • ≫ k 0 .This means that there is essentially no peak with k • = k 0 ≪ k W if k 0 ≪ k W and the original power spectrum has a sufficiently broad support in k > k 0 (Fig. 2 would be helpful for understanding).For a fixed k • = k 0 , in the both limits k 0 ≪ k W and k 0 ≫ k W , the number of peaks decreases.
In our procedure, since we take the envelope curve for all values of k W , the final estimate for the peak number density at k • = k 0 is given by the value for k W which maximizes the peak number density at k • = k 0 .For this specific value of k W , k 0 corresponds to k t introduced in Eq. (35), which basically maximizes the peak number, and generally we have k W > k t ≃ k 0 .If the window function reduces the amplitude of the power spectrum in the region of k much smaller than k W , the maximum number density of peaks with k 2. Three functions, cos(x/10) + cos(x) + cos(10x), cos(x/10) + cos(x) and cos(x/10) are plotted as functions of x.For cos(x/10) + cos(x), every peak has the scale of order 1, but the peak profiles are sharper for cos(x/10) + cos(x) + cos(10x) and broader for cos(x/10).also inevitably decreases due to the window function.Thus the final estimate for the peak number density at k • = k 0 also decreases.For this reason, we expect that a sharp cut-off of the window function would provide a larger value of the peak number density minimizing the extra reduction of the mass spectrum due to the window function.
Let us check the above discussion by considering the flat scale-invariant spectrum with a window function: We consider the following window functions: where we note that W 1 is the standard Gaussian window function 6 .For each window function, we can calculate the PBH mass spectrum with a fixed value of k W following the procedure presented in the previous section.The results are shown in Fig. 3.As is shown in Fig. 3, the result significantly depends on the window function.For the overall mass spectrum, taking the envelope curve of the mass spectra for all values of k W , we obtain the flat mass spectrum with the amplitude given by the maximum value in the plot.Therefore the k-space tophat window function gives the largest abundance as is expected.This behavior is contrary to the case of the conventional Press-Schechter (PS) formalism where the Gaussian window function gives a larger abundance than that for the k-space tophat window(see Appendix B and Ref. [9]).The σ-dependence of β approx 0,max , which gives an order-of-magnitude estimate for the maximum value of the mass spectrum, is also shown in Fig. 3.

VI. SUMMARY AND DISCUSSION
In this paper, we have improved the procedure proposed in Ref. [1] so that we can decouple the larger-scale environmental effect, which is irrelevant to the PBH formation.Thus we can eliminate the redundant variance due to the environmental effect, and obtain a narrower mass spectrum than that in Ref. [1].This new procedure also allows us to straightforwardly implement a window function and calculate the PBH abundance for an arbitrary power spectrum of the curvature perturbation.For a sufficiently narrow power spectrum, the PBH mass spectrum results in the case without the window function irrespective of the choice of the window function.That is, there is no window function dependence for a sufficiently narrow spectrum in our procedure.
It should be noted that, in Ref. [26], the authors attempted to estimate PBH abundance for a broad spectrum without a window function.In the results in Ref. [26], we can find significant enhancement of the mass spectrum in the large-mass region compared with our results.Although the reason for this discrepancy should be further investigated in the future, we make some discussion in AppendixC.
The PBH abundance for the scale-invariant flat power spectrum has been calculated in Sec.V as an example.The result largely depends on the choice of the window function.Nevertheless, we found that the k-space tophat window function has the minimum required property.Specifically, it minimizes the extra reduction of the mass spectrum due to the window function.When one estimates PBH abundance without any concrete physical process of the smoothing, the choice of the k-space tophat window function would be the best in our procedure.Finally, we emphasize that our procedure makes it possible to calculate the PBH mass spectrum for an arbitrary power spectrum by using a plausible PBH formation criterion with the nonlinear relation taken into account.where r R is a certain radius and the subscript R denotes the value at r = r R 7 .This equation is valid for super-horizon spherically symmetric perturbations.We may define δ l which is linearly related to ζ as follows: Then we obtain δ ⊜ δ l − 3 8 δ 2 l . (C3) The linear density perturbation δ l should be compared with δ R defined in Eq. ( 9) in Ref. [26] as follows: where θ is the Heviside step function, which effectively acts as the real-space tophat window function.We note that this expression is defined not only for spherically symmetric perturbations but also for general ones.Using the following linear relation: in a spherically symmetric case, we can find [40] δ at the horizon entry time defined by r R = 1/(aH).
7 r R corresponds to R in Ref. [26].
as a function of M by eliminating k • from Eqs. (20) and µ 2 = µ (k * ) 2th (k • ), and solving it for µ 2 .That is, defining k th • (M) by the inverse function of

FIG. 3 .
FIG. 3. PBH mass spectrum(left) and β approx 0,max as a function of σ(right) for the flat power spectrum with each window function.In the left panel, we set σ = 0.1 and the dotted horizontal lines show the corresponding values of β approx 0,max .

FIG. 4 .
FIG.4.The window function dependences of the PS formalism and our procedure are shown for the monochromatic spectrum (left panel) and the flat spectrum (right panel).The solid lines and dashed lines show the results for our procedure based on the peak theory and for the PS formalism, respectively.In the left panel, the values for the PS formalism are given by Eq. (B6) with k 0 = k M , and the value in our procedure is depicted as a single solid line because the fraction of PBH does not depend on the window function in our procedure.