Abundance of Primordial Black Holes in Peak Theory for an Arbitrary Power Spectrum
Chul-Moon Yoo, Tomohiro Harada, Shin'ichi Hirano, Kazunori Kohri
aa r X i v : . [ a s t r o - ph . C O ] O c t RUP-20-25KEK-Cosmo-261KEK-TH-2245
Abundance of Primordial Black Holes in Peak Theoryfor an Arbitrary Power Spectrum
Chul-Moon Yoo, ∗ Tomohiro Harada, † Shin’ichi Hirano, ‡ and Kazunori Kohri
3, 4, 5, § Division of Particle and Astrophysical Science, Graduate School of Science,Nagoya University, Nagoya 464-8602, Japan Department of Physics, Rikkyo University, Toshima, Tokyo 171-8501, Japan Institute of Particle and Nuclear Studies, KEK,1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan The Graduate University for Advanced Studies (SOKENDAI),1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan Kavli Institute for the Physics and Mathematics of the Universe (WPI),University of Tokyo, Kashiwa 277-8583, Japan
We modify the procedure to estimate PBH abundance proposed in Ref. [1] sothat it can be applied to a broad power spectrum such as the scale-invariant flatpower spectrum. In the new procedure, we focus on peaks of the Laplacian of thecurvature perturbation △ ζ and use the values of △ ζ and △△ ζ at each peak to specifythe profile of ζ as a function of the radial coordinate while the values of ζ and △ ζ areused in Ref. [1]. The new procedure decouples the larger-scale environmental effectfrom the estimate of PBH abundance. Because the redundant variance due to theenvironmental effect is eliminated, we obtain a narrower shape of the mass spectrumcompared to the previous procedure in Ref. [1]. Furthermore, the new procedureallows us to estimate PBH abundance for the scale-invariant flat power spectrumby introducing a window function. Although the final result depends on the choiceof the window function, we show that the k -space tophat window minimizes theextra 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 PBHestimation. Our procedure makes it possible to calculate the PBH mass spectrumfor an arbitrary power spectrum by using a plausible PBH formation criterion withthe nonlinear relation taken into account. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] I. INTRODUCTION
Since Zel’dovich, Novikov and Hawking had pointed out its possibility[2, 3], primordialblack holes(PBHs) have continued to attract attention. They are still viable candidates fora substantial part of dark matter(see e.g. Refs. [4, 5] and references therein) and a possibleorigin of observed binary black holes[6, 7]. Mass, spin or spatial distribution of PBHsprovides us valuable information about relatively small scale inhomogeneity in the earlyuniverse. When we connect a PBH production scenario and observational constraints on it,theoretical estimation of the PBH distribution is inevitable. Here, we provide a plausibleprocedure to calculate the PBH mass spectrum for an arbitrary power spectrum based onthe peak theory.Until relatively recently, the Press-Schechter (PS) formalism is applied to the estimationof PBH abundance based on a perturbation variable such as the comoving density or thecurvature perturbation. As PBHs started to draw more attention after the discovery of theblack hole binary as well as gravitational waves, people have begun to seriously doubt therelevance of the PS formalism in the estimation of PBH abundance. In order to improvethe estimation, one needs to resolve the following mutually related issues: PBH formationcriterion, statistical treatment of non-linear variables [8], and use of a window function [9–11].For the criterion of the PBH formation, there has been a long-term debate since Carrhad proposed a rough criterion[12]. A lot of efforts to clarify the appropriate criterion havebeen made through numerical and analytic treatments[13–22]. One useful criterion wasproposed in Ref. [15] by using the compaction function, which is equivalent to the half of thevolume 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 ofthe maximum value of the compaction function gives relatively accurate criterion, which iswithin about 10% accuracy for a moderate shape of initial configuration. More recently, thethreshold 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% accuracyfor a moderate inhomogeneity in the radiation dominated universe(see Ref. [25] for generalcosmological backgrounds). These recent developments show that the use of the compactionfunction is crucial for an accurate estimation of PBH abundance.Another important ingredient in the calculation of PBH abundance is the statistics ofperturbation variables. Naively, we expect that the curvature perturbation would be relevantfor the Gaussian distribution assumed in the PS formalism. However, the absolute value ofthe curvature perturbation does not have any physical meaning in a local sense because itcan be absorbed into the coordinate rescaling. Therefore, setting the threshold value for theabsolute value of the curvature perturbation seems irrelevant. Conversely, while setting thethreshold on the compaction function would be appropriate, the compaction function cannotbe a Gaussian variable even if the curvature perturbation is totally Gaussian because oftheir non-linear relation. Furthermore, difference of the gauge confuses the relation betweenperturbation variables.In Ref. [23], relations between different gauge conditions were summarized and the gaugeissue has been clarified. The compaction function is expressed in terms of the curvatureperturbation in the same reference [23]. Then, apart from the window function, the re-maining issue is how to count the number of PBHs by taking into account the nonlinearrelation between the curvature perturbation and the compaction function. Although a fewprocedures treating the non-linear relation have been proposed[1, 26], their consistency witheach other has not been clear yet.In Ref. [1], a plausible procedure to estimate PBH abundance for a narrow power spec-trum of the Gaussian curvature perturbation was proposed, where the threshold for thecompaction function is used and the non-linear relation is taken into account. However, thisprocedure cannot be directly applied to a broad spectrum(see Ref. [27] for a simple approachwith linear relations). Our aim in this paper is to improve the procedure in Ref. [1] so thatwe 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 functionis introduced in Sec. II. In Sec. III, focusing on a high peak of the Laplacian of the curvatureperturbation △ ζ , we characterize the typical profile of the curvature perturbation ζ aroundthe peak by using the values of △ ζ and △△ ζ . This treatment allows us to decouple theenvironmental effect on the absolute value of the curvature perturbation, and the criterioncan be expressed in a purely local manner. In Sec. IV, the procedure to estimate PBHabundance is explained, and applied to the single-scale narrow power spectrum previouslypresented in Ref. [1]. We discuss the case of the scale-invariant flat spectrum implementinga 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 lightand Newton’s gravitational constant are unity, c = G = 1. II. CRITERION FOR PBH FORMATION
Let us consider the spatial metric given byd s = a e − ζ ˜ γ ij d x i d x j (1)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, respec-tively. In the long-wavelength approximation, the curvature perturbation ζ and the densityperturbation δ with the comoving slicing are related by [23], δ = −
89 1 a H e ζ/ △ (cid:0) e − ζ/ (cid:1) (2)in the radiation dominated universe, where H is the Hubble expansion rate and △ is theLaplacian of the reference flat metric.We will be interested mainly in high peaks, which tend to be nearly spherically symmet-ric [28]. Therefore, in this section, we introduce the criterion for PBH formation originallyproposed in Ref. [15] assuming spherical symmetry. Here, we basically follow and refer tothe discussions and calculation in Ref. [23].First, let us define the compaction function C as C := δMR , (3)where R is the areal radius at the radius r , and δM is the excess of the Misner-Sharpmass enclosed by the sphere of the radius r compared with the mass inside the sphere inthe 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 (seealso Eq. (6.33) in Ref. [23]): C ( r ) = 13 h − (1 − rζ ′ ) i . (4)We will assume that the function C is a smooth function of r for r >
0. Then, the valueof C takes the maximum value C max at r m which satisfies C ′ ( r m ) = 0 ⇔ ( ζ ′ + rζ ′′ ) | r = r m = 0 . (5)We consider the following criterion for PBH formation: C max > C th ≡ δ th . (6)In the comoving slicing, the threshold C th for PBH formation is evaluated as ≃ .
267 (seeFigs. 2 and 3 or TABLE I and II in Ref. [23]). This threshold corresponds to the pertur-bation profiles of Refs. [15, 18, 29], and is found to be quite robust for a broad range ofparameters(see Ref. [24] for a more robust criterion). In this paper we shall use this valueas a reference value.
III. PEAK OF △ ζ AND THE SPHERICAL PROFILE
Throughout this paper, we assume the random Gaussian distribution of ζ with its powerspectrum P ( k ) defined by the following equation: < ˜ ζ ∗ ( k ) ˜ ζ ( k ′ ) > = 2 π k P ( k )(2 π ) δ ( k − k ′ ) , (7)where ˜ ζ ( k ) is the Fourier transform of ζ and the bracket < ... > denotes the ensembleaverage. Each gradient moment σ n can be calculated by σ n := Z d kk k n P ( k ) . (8)Hereafter we suppose that the power spectrum is given. Then the gradient moments can becalculated from the power spectrum and regarded as constants.In this paper, we focus on high peaks of ζ := △ ζ , which coincide with peaks of δ withlinear relation. We note that this procedure is different from the previous one proposed inRef. [1], where peaks of − ζ is considered.Focusing on a high peak of ζ and taking it as the origin of the coordinates, we introducethe amplitude µ and the curvature scale 1 /k • of the peak as follows: µ = ζ | r =0 , (9) k • = − △△ ζ | r =0 µ . (10)According to the peak theory[28], for a high peak, assuming the spherical symmetry, wemay expect the typical form of the profile ¯ ζ can be described by using µ and k • as follows:¯ ζ ( r ) σ = µ /σ − γ (cid:18) ψ + 13 R • △ ψ (cid:19) − µ k • /σ γ (1 − γ ) (cid:18) γ ψ + 13 R • △ ψ (cid:19) , (11)with γ = σ / ( σ σ ), R • = √ σ /σ and ψ n ( r ) = 1 σ Z d kk k n sin( kr ) kr P ( k ) . (12)It is worthy of note that, for k • = σ /σ , we obtain¯ ζ ( r ; σ /σ ) = µ ψ ( r ) . (13)It will be shown in Eq. (27) that regarding k • as a probability variable, we obtain σ /σ asthe mean value of k • .Let us consider the profile ¯ ζ given by integrating Eq. (11). Integrating ¯ ζ , and assumingthe regularity at the center, we obtain¯ ζ ( r ) σ = − µ /σ (1 − γ ) (cid:18) ψ + 13 R • △ ψ (cid:19) + µ k • /σ γ (1 − γ ) (cid:18) γ ψ + 13 R • △ ψ (cid:19) + ζ ∞ σ , (14)where ζ ∞ = ¯ ζ | r = ∞ is an integration constant. Because we have ψ | r =0 = σ /σ and △ ψ | r =0 = −
1, we obtain ζ := ¯ ζ | r =0 = − µ σ σ − σ σ + ( σ − σ σ ) k • σ σ − σ + ζ ∞ . (15)We may consider either ζ or ζ ∞ as a probability variable. Since the constant shift of ζ can be absorbed into the renormalization of the background scale factor, the non-zero valueof ζ ∞ would be regarded as an larger-scale environmental effect. Actually, in AppendixA,we show that the mean value of ζ ∞ is 0 for a given set of µ and k • . In Ref. [1], we used ζ to characterize the profile of the curvature perturbation. However, the use of ζ wouldmix the environmental effect with the local state. Therefore, in this paper, we ignore ζ ∞ Notations of the amplitude µ and the curvature scale k • are chosen so that they will be distinguishedfrom µ and k ∗ in Ref. [1]. by renormalizing the background scale factor as e − ζ ∞ a → a to eliminate the environmentaleffect, and regard ζ as a dependent variable on µ 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 µ and k • . Here, we just copy the part of the procedure from Ref. [1](theflow char of our procedure can be seen in Ref. [30]). Applying Eq. (4) to ¯ ζ , we obtain therelation between µ and C as µ = 1 − √ − C rg ′ , (16)where g ( r ; k • ) := ¯ ζ/µ and the smaller root is taken. Let us define the threshold value µ ( k • )2th as µ ( k • )2th ( k • ) = 1 − √ − C th ¯ r m ( k • ) g ′ m ( k • ) , (17)where ¯ r m ( k • ) is the value of r m for ζ = ¯ ζ , and g m ( k • ) := g (¯ r m ( k • ); k • ) . (18)In Eq. (17), we explicitly denoted the k • dependence of ¯ r m and g m to emphasize it.In order to express the threshold value as a function of the PBH mass M , let us considerthe horizon entry condition: aH = aR (¯ r m ) = 1¯ r m e µ g m . (19)Since the PBH mass is given by M = α/ (2 H ) with α being a numerical factor, from thehorizon entry condition (19), the PBH mass M can be expressed as follows: M = 12 αH − = 12 αa ¯ r m e − µ g m = M eq k ¯ r e − µ g m =: M ( µ ,k • ) ( µ , k • ) , (20)where we have used the fact H ∝ a − and a = a H eq ¯ r m e − µ g m with a eq and H eq being thescale factor and Hubble expansion rate at the matter-radiation equality. M eq and k eq aredefined by M eq = αH − / k eq = a eq H eq , respectively. For simplicity, we set α = 1 as afiducial value .Then we can obtain the threshold value of µ ( M )2th ( M ) as a function of M by eliminating k • from Eqs. (20) and µ = µ ( k ∗ )2th ( k • ), and solving it for µ . That is, defining k th • ( M ) by theinverse function of M = M ( µ ,k • ) ( µ ( k • )2th ( k • ) , k • ), we obtain the threshold value of µ ( M )2th for afixed value of M as µ ( M )2th ( M ) := µ ( k • )2th ( k th • ( M )) . (21)While, from Eq. (20), we can describe µ as a function of M and k • as follows: µ = µ ( M,k • ) ( M, k • ) := − g m ln (cid:18) k ¯ r MM eq (cid:19) . (22) In order to take into account the critical behavior[31, 32], α should be given by a function of µ and k • as α = K ( k • )( µ − µ ( k • )) γ with γ ≃ .
36 [16, 17, 19, 29, 33–37] and K ( k • ) being some function of k • ,which would be profile dependent. The value of µ may be bounded below by µ ( M ) for a fixed value of M . Actually, in thespecific examples in Sec. IV and V, the value of µ ( M ) is given as follows: µ = µ ( M,k • ) ( M, . (23)Then, for a fixed value of M , the region of µ for PBH formation can be given by µ > µ := max n µ ( M ) , µ ( M )2th ( M ) o . (24) IV. PBH ABUNDANCE
From Ref. [1], we obtain the expression for the peak number density characterized by µ and k • as n ( k • )pk ( µ , k • )d µ d k • = 2 · / (2 π ) / µ k • σ σ σ f (cid:18) µ k • σ (cid:19) P (cid:18) µ σ , µ k • σ (cid:19) d µ d k • , (25)where f ( x ) = 12 x ( x − erf " r x + erf "r x + r π (cid:26)(cid:18)
85 + 314 x (cid:19) exp (cid:20) − x (cid:21) + (cid:18) −
85 + 12 x (cid:19) exp (cid:20) − x (cid:21)(cid:27) , (26)and P (cid:18) µ σ , µ k • σ (cid:19) = µ k • πσ σ p − γ exp (cid:20) − µ σ ( k • ) (cid:21) (27)with 1˜ σ ( k • ) := 1 σ + 1 σ (1 − γ ) (cid:18) k • − σ σ (cid:19) . (28)In Eq. (25), the following replacements have been made from Eq. (58) in Ref. [1]: µ → µ , k ∗ → k • , σ n → σ n +2 , γ → γ . (29)Since the direct observable is not k • but the PBH mass M , we further change the variablefrom k • to M as follows: n ( M )pk ( µ , M )d µ d M := n ( k • )pk ( µ , k • )d µ d k • = 3 / (2 π ) / σ σ σ µ k • f (cid:18) µ k • σ (cid:19) P (cid:18) µ σ , µ k • σ (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) dd k • ln ¯ r m − µ dd k • g m (cid:12)(cid:12)(cid:12)(cid:12) − d µ d ln M, (30)where k • should be regarded as a function of µ and M given by solving Eq. (20) for k • .We note that an extended power spectrum is implicitly assumed in the above expression. Inthe monochromatic spectrum case, the expression reduces to the same expression in Ref. [1]because σ n = σ for P ( k ) = σ k δ ( k − k ).It should be noted that, since we relate k • to M with µ fixed, we have implicitly assumedthat there is only one peak with △ ζ = − µ k • 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 separatedscales, and the typical PBH mass is relatively larger than the minimum scale given by thespectrum, we would find multiple peaks inside r = r m . Then we cannot correctly count thenumber of peaks in the scale of interest. In order to avoid this difficulty, we need to introducea window function to smooth out the smaller-scale inhomogeneities. We discuss this issuein the subsequent section. In this section, we simply assume that the power spectrum ischaracterized by a single scale k and there is no contribution from the much smaller scales k ≫ k .The number density of PBHs is given by n BH d ln M = (cid:20)Z ∞ µ d µ n ( M )pk ( µ , M ) (cid:21) M d ln M. (31)We also note that the scale factor a is a function of M as a = 2 M / M / k eq /α . Then thefraction of PBHs to the total density β d ln M can be given by β d ln M = M n BH ρa d ln M = 4 π αn BH k − (cid:18) MM eq (cid:19) / d ln M = 2 · / αk − (2 π ) / σ σ σ (cid:18) MM eq (cid:19) / "Z ∞ µ d µµ k • f (cid:18) µ k • σ (cid:19) P (cid:18) µ σ , µ k • σ (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) dd k • ln ¯ r m − µ dd k • g m (cid:12)(cid:12)(cid:12)(cid:12) − d ln M. (32)Here we note again that k • should be regarded as a function of µ and M . The aboveformula can be evaluated in principle once the form of the power spectrum is given. Acrucial difference of Eq. (32) from Eq. (61) in Ref. [1] is that the expression does not dependon σ , which has IR-log divergence for the flat scale-invariant spectrum. Thus we canconsider the PBH mass spectrum for the flat spectrum without introducing IR cut-off. In order to give a simpler approximate form of Eq. (32), we approximately perform theintegral with respect to µ as follows: β d ln M ≃ · / αk − (2 π ) / σ σ σ (cid:18) MM eq (cid:19) / " ˜ σ ( k • ) k • f (cid:18) µ k • σ (cid:19) P (cid:18) µ σ , µ k • σ (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) dd k • ln ¯ r m − µ dd k • g m (cid:12)(cid:12)(cid:12)(cid:12) − µ = µ d ln M. (33)Since P given in Eq. (27) has the exponential dependence, we may expect that the value of β is sensitive to the exponent − µ / σ . Therefore, assuming µ = µ ( M )2th = µ ( k • )2th ( k th • ( M )), The PBH fraction to the total density f at the equality time is given by f = β ( M eq /M ) / . We do notexplicitly show the form of f in this paper since the scale dependence can be more easily understood bythe form of β . we can roughly estimate the maximum value of β at the top of the mass spectrum byconsidering the value k t of k • which minimizes the value of µ k • ( k • ) / ˜ σ , namely, k t := argmin k • h µ ( k • )2th ( k • ) / ˜ σ ( k • ) i . (34)The value of k t cannot be given in an analytic form in general, and a numerical procedureto find the value of k t is needed. We note that the value of k t is independent of the overallconstant 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 maximumvalue β , max : β , max ≃ β approx0 , max := 2 · / αk − (2 π ) / σ σ σ (cid:18) M t M eq (cid:19) / " ˜ σ ( k • ) k • f (cid:18) µ k • σ (cid:19) P (cid:18) µ σ , µ k • σ (cid:19) (cid:12)(cid:12)(cid:12)(cid:12) dd k • ln ¯ r m − µ dd k • g m (cid:12)(cid:12)(cid:12)(cid:12) − k • = k t , µ = µ ( k • )2th ( k t ) , (35)where M t := M ( µ ,k • ) ( µ ( k • )2th ( k t ) , k t ). We note that the expression (35) gives a better approx-imation for a smaller value of the amplitude of the power spectrum σ , and may have afactor of difference from the actual maximum value for σ & . P ( k ) = 3 r π σ (cid:18) kk (cid:19) exp (cid:18) − k k (cid:19) . (36)Gradient moments are calculated as σ n = 2 n +1 n √ π Γ (cid:18)
32 + n (cid:19) σ k n , (37)where Γ means the gamma function. The result is shown in Fig. 1. Our new procedure givesa narrower and slightly higher spectrum than that obtained in Ref [1]. This behavior can beunderstood as the environmental effect induced by the variance of ζ ∞ . Although the effectis so small that it could be practically ignored in this example, we successfully decoupledthe environmental effect. V. IMPLEMENTING A WINDOW FUNCTION
In our new procedure, a window function can be straightforwardly implemented. Thatis, introducing the UV cut-off scale k W , instead of Eq. (7), we consider the following powerspectrum of ζ : P W ( k ) = P ( k ) W ( k ; k W ) , (38) The expansion around k • = σ /σ like in Ref. [1] does not work well because of the large k • dependence of µ ( k • )2th . That is, the peak of the exponent − µ / ˜ σ significantly deviates from k • = σ /σ and the Taylorexpansion is not as effective as in Ref. [1]. argmin x f ( x ) = { x | ∀ y ( f ( y ) ≥ f ( x )) } . FIG. 1. PBH mass spectrum(left) and β approx0 , max as a function of σ (right) for the extended powerspectrum (36) with k = 10 k eq . The solid lines correspond to the spectra calculated by our newprocedure with α = 1, and the dashed lines show the spectra calculated in Ref. [1]. We also plotthe mass spectrum β PS0 obtained from the Press-Schechter formalism explained in Appendix forcomparison. In the left panel, the dotted horizontal lines show the corresponding values of β approx0 , max . where W ( k ; k W ) is a window function satisfying W ( k ; k W ) ≤ W ( k ; k W ) = 0 for k ≫ k W . Then, following the procedure given in the previous section, we can calculate PBHabundance for a given value of k W . The final PBH mass spectrum is given by the envelopecurve 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 withoutthe 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 theeffect of the window function for a peak number density at a fixed scale given by the wavenumber k • = k . If k ≫ k W , no peak can be found. On the other hand, if k ≪ k W ,we would find many smaller scale peaks inside the region of the radius ∼ /k because thesmaller scale modes with k ≫ k are superposed on top of the inhomogeneity with k ∼ k .Thus every peak has a sharp profile due to the superposed small scale inhomogeneity andsatisfies k • ≫ k . This means that there is essentially no peak with k • = k ≪ k W if k ≪ k W and the original power spectrum has a sufficiently broad support in k > k (Fig. 2would be helpful for understanding). For a fixed k • = k , in the both limits k ≪ k W and k ≫ k W , the number of peaks decreases.In our procedure, since we take the envelope curve for all values of k W , the final estimatefor the peak number density at k • = k is given by the value for k W which maximizes the peaknumber density at k • = k . For this specific value of k W , k corresponds to k t introduced inEq. (35), which basically maximizes the peak number, and generally we have k W > k t ≃ k .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 • = k ≃ k t < k W also inevitably decreases due to the window function. Thus the final estimate for the peaknumber density at k • = k also decreases. For this reason, we expect that a sharp cut-off ofthe window function would provide a larger value of the peak number density minimizing1 FIG. 2. Three functions, cos( x/
10) + cos( x ) + cos(10 x ), cos( x/
10) + cos( x ) and cos( x/
10) areplotted as functions of x . For cos( x/
10) + cos( x ), every peak has the scale of order 1, but the peakprofiles are sharper for cos( x/
10) + cos( x ) + cos(10 x ) and broader for cos( x/ 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 awindow function: P W ( k ; k W ) = σ W ( k ; k W ) . (39)We consider the following window functions: W n ( k/k W ) = exp (cid:18) − k n k nW (cid:19) , (40) W k TH ( k/k W ) = Θ( k W − k ) , (41)where we note that W is the standard Gaussian window function . For each windowfunction, we can calculate the PBH mass spectrum with a fixed value of k W followingthe procedure presented in the previous section. The results are shown in Fig. 3. As isshown in Fig. 3, the result significantly depends on the window function. For the overallmass spectrum, taking the envelope curve of the mass spectra for all values of k W , weobtain 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) formalismwhere the Gaussian window function gives a larger abundance than that for the k -spacetophat window(see Appendix B and Ref. [9]). The σ -dependence of β approx0 , max , which gives anorder-of-magnitude estimate for the maximum value of the mass spectrum, is also shown inFig. 3. The real-space tophat window function leads divergent gradient moments for the scale-invariant flat spec-trum, so that we practically cannot use it. FIG. 3. PBH mass spectrum(left) and β approx0 , max as a function of σ (right) for the flat power spectrumwith each window function. In the left panel, we set σ = 0 . β approx0 , max . VI. SUMMARY AND DISCUSSION
In this paper, we have improved the procedure proposed in Ref. [1] so that we can decouplethe larger-scale environmental effect, which is irrelevant to the PBH formation. Thus wecan eliminate the redundant variance due to the environmental effect, and obtain narrowermass spectrum than that in Ref. [1]. This new procedure also allows us to straightforwardlyimplement a window function and calculate the PBH abundance for an arbitrary powerspectrum of the curvature perturbation. For a sufficiently narrow power spectrum, the PBHmass spectrum results in the case without the window function irrespective of the choiceof the window function. That is, there is no window function dependence for a sufficientlynarrow spectrum in our procedure.It should be noted that, in Ref. [26], the authors attempted to estimate PBH abundancefor a broad spectrum without a window function. In the results in Ref. [26], we can findsignificant enhancement of the mass spectrum in the large-mass region compared with ourresults. 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 inSec. 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 requiredproperty. Specifically, it minimizes the extra reduction of the mass spectrum due to the win-dow function. When one estimates PBH abundance without any concrete physical processof the smoothing, the choice of the k -space tophat window function would be the best inour procedure. Finally, we emphasize that our procedure makes it possible to calculate thePBH mass spectrum for an arbitrary power spectrum by using a plausible PBH formationcriterion with the nonlinear relation taken into account.3 ACKNOWLEDGEMENTS
We thank Jaume Garriga for his contribution to a part of this work provided during discus-sion in the previous work[1]. We also thank Shuichiro Yokoyama and Cristiano Germani forhelpful comments and stimulating discussion. This work was supported by JSPS KAKENHIGrant Numbers JP19H01895 (C.Y., T.H. and S.H.), JP19K03876 (T.H.) and JP17H01131(K.K.), and MEXT KAKENHI Grant Nos. JP19H05114 (K.K.) and JP20H04750 (K.K.).
Appendix A: Random Gaussian distribution of ζ Due to the random Gaussian assumption, the probability distribution of any set of lin-ear combination of the variable ζ ( x i ) is given by a multi-dimensional Gaussian probabilitydistribution[28, 38]: P ( V I )d n V I = (2 π ) − n/ | det M| − / exp (cid:20) − V I (cid:0) M − (cid:1) IJ V J (cid:21) d n V, (A1)where the components of the matrix M are given by the correlation < V I V J > defined by < V I V J > := Z d k (2 π ) d k ′ (2 π ) < ˜ V ∗ I ( k ) ˜ V J ( k ′ ) > (A2)with ˜ V I ( k ) = R d xV I ( x )e i kx .The non-zero correlations between two of ν = − ζ /σ , ξ = △ ζ /σ and ω = −△△ ζ /σ aregiven by < νν > = < ξξ > = < ωω > = 1 , (A3) < νξ > = γ := σ / ( σ σ ) , (A4) < νω > = γ := σ / ( σ σ ) , (A5) < ξω > = γ := σ / ( σ σ ) . (A6)Then, the probability distribution function for these variables are given as follows: P ( ν, ξ, ω ) = (2 π ) − / | D | − / exp " − D n (1 − γ ) ν + (1 − γ ) ξ + (1 − γ ) ω − γ − γ γ ) νξ − γ − γ γ ) ων − γ − γ γ ) ξω o , (A7)where D = det M = 1 − γ − γ − γ + 2 γ γ γ (A8)with M = γ γ γ γ γ γ . (A9)4We re-express the probability P as a probability distribution function ˜ P of ζ , µ and k • ,that is˜ P ( ζ , µ , k • )d ζ d µ d k • = P ( ν, ξ, ω )d ν d ξ d ω = 2 µ k • σ σ σ P (cid:18) ζ σ , µ σ , µ k • σ (cid:19) d ζ d µ d k • . (A10)Then, the conditional probability p ( ζ ) with fixed µ and k • is given by p ( ζ ) = (cid:18) − γ πDσ (cid:19) / exp (cid:20) − − γ Dσ (cid:0) ζ − ¯ ζ (cid:1) (cid:21) = (cid:18) − γ πDσ (cid:19) / exp (cid:20) − − γ Dσ ζ ∞ (cid:21) , (A11)where ¯ ζ = − µ ( σ σ − σ σ ) + ( σ − σ σ ) k • σ σ − σ . (A12) Appendix B: Estimation and window function dependencein the Press-Schechter formalism
For a comparison, we review a conventional estimate of the fraction of PBHs based onthe PS formalism. In the conventional formalism, the scale dependence is introduced by awindow function W ( k/k M ), where k M = k eq ( M eq /M ) / . (B1)Then, each gradient moment is replaced by the following expression:ˆ σ n ( k M ) = Z d kk k n P ( k ) W ( k/k M ) . (B2)The conventional estimate starts from the following Gaussian distribution assumption forthe density perturbation ¯ δ : P δ (¯ δ )d¯ δ = 1 √ πσ δ exp (cid:18) −
12 ¯ δ σ δ (cid:19) d¯ δ, (B3)where σ δ is given by the coarse-grained density contrast σ δ ( k M ) = 49 ˆ σ ( k M ) k M . (B4)Here, for simplicity, we use the same numerical value of δ th as in our approach, in otherwords, we assume that the volume average of the density perturbation obeys the Gaussianprobability distribution given by Eq. (B3) with the coarse-grained density contrast (B4) inthe PS formalism. This Gaussian distribution and the dispersion are motivated by the linearrelation between ζ and δ . The fraction β PS0 is then evaluated as follows(see e.g. [39]): β PS0 ( M ) = 2 α Z ∞ δ th d¯ δP δ (¯ δ ) = α erfc (cid:18) δ th √ σ δ ( k M ) (cid:19) = α erfc (cid:18) δ th k M √ σ ( k M ) (cid:19) . (B5)5Let us check the window function dependence in the PS formalism for the cases of themonochromatic spectrum P ( k ) = σ k δ ( k − k ) and the flat spectrum P ( k ) = σ . For themonochromatic spectrum, the faction β PS0 , mono is given by β PS0 , mono = efrc (cid:18) δ th k M √ σk W ( k /k M ) (cid:19) . (B6)From this expression, the existence and significance of the window function dependence isclear (see the left panel in Fig. 4). On the other hand, as is stated in the first paragraph inSec. V, for a narrow power spectrum, there is essentially no window function dependence inour procedure. In particular, for the monochromatic spectrum case, the fraction reduces tothe result given in Ref. [1] for an arbitrary window function satisfying the properties listedin the same paragraph.For the flat spectrum, the value of the second gradient moment is given by σk M / √ σk M / k -space tophat window functions, respectively. Therefore, inthe PS formalism, the abundance is larger for the Gaussian window function differently fromour case shown in Fig. 3. In the right panel of Fig. 4, the window function dependencesin the PS formalism and our procedure are shown. In both cases, the window functiondependence is significant at a similar extent.It should be noted that in the PS formalism, the reason for the larger abundance for theGaussian window is the contribution from the high- k modes through the tail of the Gaussianfunction. Therefore, it is clear that the longer the tail of the window function is, the largerabundance becomes. Of course, we cannot take the long-tail limit because the windowfunction becomes irrelevant in this limit. Contrary to the PS formalism, in our procedure,the sharpest cutoff in the k -space gives the largest abundance and the extra reduction dueto the window function can be minimized in this well-defined limit. FIG. 4. The window function dependences of the PS formalism and our procedure are shown forthe monochromatic spectrum (left panel) and the flat spectrum (right panel). The solid lines anddashed 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 = k M ,and the value in our procedure is depicted as a single solid line because the fraction of PBH doesnot depend on the window function in our procedure. Appendix C: Discrepancy between our results and those in Ref. [26]
As is noted in Sec. VI, there can be seen a qualitative difference between our PBH massspectrum and that in Ref. [26] in the large-mass region. First, we briefly review the basicidea used in Ref. [26](see also Ref. [40]).In order to clearly distinguish the equations which are valid only for spherically symmetriccases from generally valid equations, we use the notation ⊜ for the equality with sphericalsymmetry. Let us start with the relation between the non-linear volume-averaged densityperturbation ¯ δ , the compaction function C and the curvature perturbation[23]:¯ δ ⊜ C ⊜ δMR ⊜ − r R ζ ′ R [2 + r R ζ ′ R ] , (C1)where r R is a certain radius and the subscript R denotes the value at r = r R . This equationis valid for super-horizon spherically symmetric perturbations. We may define δ l which islinearly related to ζ as follows: δ l ⊜ − r R ζ ′ R . (C2)Then we obtain ¯ δ ⊜ δ l − δ . (C3)The linear density perturbation δ l should be compared with δ R defined in Eq. (9) inRef. [26] as follows: δ R ( r R ) = 34 πr R Z d x δρρ θ ( r R − | ~x − ~x | ) , (C4)where θ is the Heviside step function, which effectively acts as the real-space tophat win-dow function. We note that this expression is defined not only for spherically symmetricperturbations but for general ones. Using the following linear relation: δρρ = −
49 1 a H △ ζ ⊜ −
49 1 a H r ∂ r ( r ∂ r ζ ) (C5)in a spherically symmetric case, we can find [40] δ R ( r R ) ⊜ r R Z r R drr ( ζ ′′ + 2 r ζ ′ ) ⊜ − r R ζ ′ R ⊜ δ l (C6)at the horizon entry time defined by r R = 1 / ( aH ).In Ref. [26], the relation with spherical symmetry 2 C = δ l − δ is extended to the generalrelation 2 C = δ R − δ R , and the PBH formation criterion for C is expressed in terms of δ R ,and used to estimate PBH abundance. δ R is equivalent to the linear density perturbationwith the real-space tophat window function. However, it should be noted that the real-spacetophat window function is naturally introduced so that the relation Eq. (C6) can be satisfied,and not introduced by hand as a window function. r R corresponds to R in Ref. [26]. /k . First, we note that, in Ref. [26], the value of r R is chosen such thatd δ R / d r R = 0 and C ( ~x, r R ) takes a maximal value at ~x = ~x , where C is regarded as a functionof ~x and r R . The scale of the relevant region to PBH formation is given by r R , which can besignificantly different from 1 /k . However, the relevance of the present criterion given interms of the compaction function C is not clear for the outer maxima. In Fig. 2 of Ref. [41],the result of a numerical simulation for a spherically symmetric and oscillatory initial profileis shown. The initial profile in the simulation corresponds to the most probable profile forthe delta-functional power spectrum peaked at 1 /k . The most probable profile is given bya peak at the center surrounded by repeated concentric overdense and underdense regions.More precisely, the most probable profile of the curvature perturbation is given by a sincfunction, where the compaction function is oscillatory with respect to the distance fromthe center. We can find that the PBH formation criterion is satisfied for the multiple radiicorresponding to the maxima in Fig. 2 of Ref. [41]. The resultant PBH, however, has themass corresponding to the typical scale 1 /k , whereas no PBH of larger mass scales is formed.This result suggests that the present criterion is relevant only for the innermost maximumof the compaction function but not for the outer maxima. If the present criterion is simplyapplied to not only the innermost but also outer maxima, the abundance of primordial blackholes of large-mass scales could be overestimated. [1] C.-M. Yoo, T. Harada, J. Garriga, and K. Kohri, PTEP , 123E01 (2018),arXiv:1805.03946, Primordial black hole abundance from random Gaussian curvature pertur-bations and a local density threshold .[2] Y. B. Zel’dovich and I. D. Novikov, Soviet Ast. , 602 (1967), The Hypothesis of CoresRetarded during Expansion and the Hot Cosmological Model .[3] S. Hawking, Mon. Not. Roy. Astron. Soc. , 75 (1971),
Gravitationally collapsed objects ofvery low mass .[4] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, (2020), arXiv:2002.12778,
Constraints onPrimordial Black Holes .[5] B. Carr and F. Kuhnel, (2020), arXiv:2006.02838,
Primordial Black Holes as Dark Matter:Recent Developments .[6] Virgo, LIGO Scientific, B. P. Abbott et al. , Phys. Rev. Lett. , 061102 (2016),arXiv:1602.03837,
Observation of Gravitational Waves from a Binary Black Hole Merger .[7] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. , 061101(2016), arXiv:1603.08338,
Primordial Black Hole Scenario for the Gravitational-Wave EventGW150914 .[8] V. De Luca et al. , JCAP , 048 (2019), arXiv:1904.00970, The Ineludible non-Gaussianityof the Primordial Black Hole Abundance . In our procedure, the relevant scale is r m ∼ /k . [9] K. Ando, K. Inomata, and M. Kawasaki, Phys. Rev. D , 103528 (2018), arXiv:1802.06393, Primordial black holes and uncertainties in the choice of the window function .[10] S. Young, Int. J. Mod. Phys. D , 2030002 (2019), arXiv:1905.01230, The primordial blackhole formation criterion re-examined: Parametrisation, timing and the choice of window func-tion .[11] K. Tokeshi, K. Inomata, and J. Yokoyama, (2020), arXiv:2005.07153,
Window function de-pendence of the novel mass function of primordial black holes .[12] B. J. Carr, Astrophys. J. , 1 (1975),
The Primordial black hole mass spectrum .[13] D. K. Nadezhin, I. D. Novikov, and A. G. Polnarev, Soviet Ast. , 129 (1978), The hydro-dynamics of primordial black hole formation .[14] I. D. Novikov and A. G. Polnarev, Soviet Ast. , 147 (1980), The Hydrodynamics of Pri-mordial Black Hole Formation - Dependence on the Equation of State .[15] M. Shibata and M. Sasaki, Phys. Rev.
D60 , 084002 (1999), arXiv:gr-qc/9905064,
Black holeformation in the Friedmann universe: Formulation and computation in numerical relativity .[16] J. C. Niemeyer and K. Jedamzik, Phys. Rev.
D59 , 124013 (1999), arXiv:astro-ph/9901292,
Dynamics of primordial black hole formation .[17] I. Musco, J. C. Miller, and L. Rezzolla, Class. Quant. Grav. , 1405 (2005),arXiv:gr-qc/0412063, Computations of primordial black hole formation .[18] A. G. Polnarev and I. Musco, Class. Quant. Grav. , 1405 (2007), arXiv:gr-qc/0605122, Curvature profiles as initial conditions for primordial black hole formation .[19] I. Musco, J. C. Miller, and A. G. Polnarev, Class. Quant. Grav. , 235001 (2009),arXiv:0811.1452, Primordial black hole formation in the radiative era: Investigation of thecritical nature of the collapse .[20] A. G. Polnarev, T. Nakama, and J. Yokoyama, JCAP , 027 (2012), arXiv:1204.6601,
Self-consistent initial conditions for primordial black hole formation .[21] T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama, JCAP , 037 (2014),arXiv:1310.3007,
Identifying the most crucial parameters of the initial curvature profile forprimordial black hole formation .[22] T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev.
D88 , 084051 (2013), arXiv:1309.4201,
Threshold of primordial black hole formation , [Erratum: Phys. Rev.D89,no.2,029903(2014)].[23] T. Harada, C.-M. Yoo, T. Nakama, and Y. Koga, Phys. Rev.
D91 , 084057 (2015),arXiv:1503.03934,
Cosmological long-wavelength solutions and primordial black hole forma-tion .[24] A. Escriv`a, C. Germani, and R. K. Sheth, Phys. Rev.
D101 , 044022 (2020), arXiv:1907.13311,
Universal threshold for primordial black hole formation .[25] A. Escriv`a, C. Germani, and R. K. Sheth, (2020), arXiv:2007.05564,
Analytical thresholds forblack hole formation in general cosmological backgrounds .[26] C. Germani and R. K. Sheth, Phys. Rev. D , 063520 (2020), arXiv:1912.07072,
Nonlinearstatistics of primordial black holes from Gaussian curvature perturbations .[27] T. Suyama and S. Yokoyama, PTEP , 023E03 (2020), arXiv:1912.04687,
A novel formu-lation of the PBH mass function .[28] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. , 15 (1986),
The statistics of peaks of Gaussian random fields .[29] I. Musco and J. C. Miller, Class. Quant. Grav. , 145009 (2013), arXiv:1201.2379, Primordialblack hole formation in the early universe: critical behaviour and self-similarity .[30] C.-M. Yoo, J.-O. Gong, and S. Yokoyama, JCAP , 033 (2019), arXiv:1906.06790, Abundanceof primordial black holes with local non-Gaussianity in peak theory .[31] M. W. Choptuik, Phys. Rev. Lett. , 9 (1993), Universality and scaling in gravitationalcollapse of a massless scalar field .[32] T. Koike, T. Hara, and S. Adachi, Phys. Rev. Lett. , 5170 (1995), arXiv:gr-qc/9503007, Critical behavior in gravitational collapse of radiation fluid: A Renormalization group (linearperturbation) analysis .[33] J. C. Niemeyer and K. Jedamzik, Phys. Rev. Lett. , 5481 (1998), arXiv:astro-ph/9709072, Near-critical gravitational collapse and the initial mass function of primordial black holes .[34] J. Yokoyama, Phys. Rev.
D58 , 107502 (1998), arXiv:gr-qc/9804041,
Cosmological constraintson primordial black holes produced in the near critical gravitational collapse .[35] A. M. Green and A. R. Liddle, Phys. Rev.
D60 , 063509 (1999), arXiv:astro-ph/9901268,
Critical collapse and the primordial black hole initial mass function .[36] F. Kuhnel, C. Rampf, and M. Sandstad, Eur. Phys. J.
C76 , 93 (2016), arXiv:1512.00488,
Effects of Critical Collapse on Primordial Black-Hole Mass Spectra .[37] C. Germani and I. Musco, Phys. Rev. Lett. , 141302 (2019), arXiv:1805.04087,
Abundanceof Primordial Black Holes Depends on the Shape of the Inflationary Power Spectrum .[38] A. G. Doroshkevich, Astrofizika , 581 (1970), The space structure of perturbations and theorigin of rotation of galaxies in the theory of fluctuation. [39] B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev.
D94 , 083504 (2016), arXiv:1607.06077,
Primordial Black Holes as Dark Matter .[40] S. Young, I. Musco, and C. T. Byrnes, (2019), arXiv:1904.00984,