The basics of primordial black hole formation and abundance estimation

This paper is a biased review of primordial black hole (PBH) formation and abundance estimation. We first review the three-zone model for PBH formation to help an intuitive understanding of the PBH formation process. Then, for more accurate analyses, we introduce necessary tools such as cosmological long-wavelength solutions, the definition of the mass and compaction function in a spherically symmetric spacetime and peak theory. Combining all these tools, we calculate the PBH mass spectrum for the case of the monochromatic curvature power spectrum as a demonstration.


Introduction
Primordial black holes (PBHs) keep attracting the attention of researchers in cosmology, high energy and gravitational physics since they have been proposed [1][2][3].The specific black hole formation process is fascinating for people who are interested in general relativity and gravitational physics in the first place.The cosmological consequence can be sufficiently significant to make cosmologists pay attention to PBHs.Generation of the primordial fluctuations which may cause an abundant population of PBHs is closely related to the underlying inflationary scenario and the associated high-energy physics model.
In Sec. 2, in order to provide an intuitive understanding of the PBH formation process, we introduce the three-zone model based on Ref. [31].The cosmological long-wavelength approximation, which is a necessary ingredient for PBH study, is reviewed in Sec. 3

following
Refs. [28,34].Assuming spherical symmetry, the PBH formation criterion is discussed in Sec. 4 based on the compaction function introduced in Ref. [28] for the first time.Our method to estimate the PBH abundance in peak theory is reviewed in Sec. 5 following Refs.[66,67].We calculate the PBH mass spectrum for the monochromatic curvature power spectrum as a demonstration.Brief reviews for a conserved mass in spherically symmetric spacetimes and peak theory [68] are attached as Appendices.
Throughout this paper, general relativity is assumed, and we use the geometrized units in which both the speed of light and Newton's gravitational constant are set to unity, c = G = 1.

PBH formation process: three-zone model
First, in order to understand the overall picture of PBH formation, let us review the threezone model discussed in Ref. [31].The three-zone model describes the collapsing system in a cosmological background with a patchwork of Friedmann-Lemaître-Robertson-Walker (FLRW) models.

FLRW models
Let us write the line element of an FLRW model as follows: where a(t), K and dΩ are the scale factor, spatial curvature and line element of the unit two-sphere.One of the Friedmann equations is given by where the dot "˙" denotes the time derivative.Introducing the areal radius R := a(t)r, Eq. (2.2) can be rewritten as follows: Ṙ2 where M is the mass inside the comoving radius r given by 4πρR 3 /3.This equation can be regarded as the energy conservation law for a fluid element on the comoving radius r.The first and second terms correspond to kinetic energy and gravitational potential energy, respectively.The time-independent total energy for the motion of the fluid element is given by the right-hand side.Since the second term is inversely proportional to the areal radius R, if K > 0, initially expanding fluid will turn around at the maximum radius R max = M/E(r).At the maximum expansion, the density ρ max and the scale factor a max satisfy 8π 3 ρ max = K a 2 max . (2.4) In contrast, the fluid continues to expand if K ≤ 0.

Background flat FLRW
We assume that the background universe is given by a spatially flat (K = 0) FLRW universe.
For the matter content, we suppose the perfect fluid with a linear equation of state given by p = wρ with a constant w, where ρ and p are the energy density and the pressure.The constant w is given by 0, 1/3, and −1 for the matter-dominated (MD) universe, radiation-dominated (RD) universe, and a cosmological constant (Λ), respectively.From the conservation law of the fluid d(ρa 3 ) = −pda 3 , we obtain ρ ∝ a −3(1+w) . (2.5) Then, we find the behavior of the Hubble radius as Let us review the rough sketch of the standard evolution of the background universe and the perturbations on it (see Fig. 1).In the standard scenario, an inflationary era (R H ∼ const.) is considered before the decelerated expansion phases of our universe.The inflationary era may be followed by an early MD phase (R H ∼ a 3/2 ) supported by the inflaton oscillation, and the universe supposed to become RD universe (R H ∼ a 2 ) after the completion of the reheating.During this early period, all relevant perturbations have super-horizon scales.As will be reviewed in subsequent sections, the amplitude of the curvature perturbation, which will be defined as the inhomogeneity of the spatial volume, is conserved on super-horizon scales.Thus the fluctuations generated in the inflationary era are "frozen" after the horizon exit, and start to grow after the horizon entry.

Super-horizon inhomogeneity
What we mainly discuss in this paper is the black hole formation after the horizon entry during the RD era.Therefore perturbations on super-horizon scales are considered for initial fluctuations of the subsequent evolution of an over-dense region.The evolution of super-horizon perturbations, whose comoving wavenumber k satisfies k ≪ aH, will be summarized in Sec. 3.
Here let us take a result in advance.At the leading order of the long-wavelength perturbation (ϵ := k/(aH) ≪ 1), we obtain the following metric form: where δ ij is the Kronecker delta and Ψ(x) is an arbitrary function of the spatial coordinates x i .It would be fruitful to compare Eq. (2.7) with the metric of the FLRW universe in the isotropic coordinates: (2.8) At the leading order of the long-wavelength approximation, the conformal factor of the spatial metric can be an arbitrary function, that is, the long-wavelength inhomogeneity may be intuitively regarded as the position-dependent spatial curvature.Since a sufficiently overdense region is expected to have a moment of maximum expansion and start to contract, such a region may be considered to have locally positive curvature.

Three-zone model
Motivated by the previous subsection, let us consider the three-zone model composed of the spherical closed FLRW region r < r c surrounded by an under-dense region embedded in the flat FLRW background (Fig. 2).This model can be an exact solution only for the dust case (w = 0) 1) .In the case of w ̸ = 0, the pressure gradient force diverges on the discontinuous spheres, and a singular shell must be introduced to realize the configuration as a solution of the Einstein equations.Nevertheless, the three-zone model is useful to intuitively understand the PBH formation process and roughly estimate the threshold for the amplitude of the density perturbation at the horizon entry.
1) The under-dense region is necessary for the dust case to be an exact solution.Let us take uniform Hubble time slices and compare the Friedmann equations (2.2) in the overdense region and the flat background.On a uniform Hubble time slice, since the value of the Hubble expansion rate is identical, we obtain where ρ and ρ b denote the density in the overdense and background regions, respectively.The overdense region is described by the closed FLRW with the spatial curvature K. From Eq. (2.9), the density perturbation can be calculated as (2.10) The value of the density perturbation δ H at the horizon entry R c := ar c = 1/H is given by Since the value of the comoving radius r cannot be larger than √ K, we find the following geometrical upper bound [69]: (2.12)

Jeans criterion
If the whole dynamics of the overdense region is described by the closed FLRW universe, nothing prevents the contraction of the closed universe and arbitrarily small amplitude of the positive density perturbation results in black hole formation.However, in reality, the pressure gradient would work as a preventing factor against the contraction.In order to estimate the effect of the pressure gradient, let us apply the Jeans criterion to the three-zone model.The Jeans criterion states that, if the free-fall time scale of the system is shorter than the sound-wave propagation time scale, the gravitational collapse cannot be prevented by the pressure gradient.For the overdense region, the sound-wave propagation time-scale t s would be given by where we have used the fact that the sound speed c s is given by √ w.The free-fall time-scale t ff is given by (2.14) Let us define a th as the value of the scale factor satisfying t s | a=a th = t ff | a=a th .Then if a max > a th , the sound-wave propagation time scale is shorter than the free-fall time scale at the time a = a max , and the gravitational contraction would be prevented by the pressure gradient effect.
For the gravitational collapse and black hole formation to be realized, the following condition must be satisfied: This condition is often quoted as Carr's condition [70].Here it should be noted that Carr's condition has the ambiguity of a factor of O(1) coming from the order estimations of Eq. (2.13) and Eq.(2.14).

Refinement of the threshold
Let us fix the ambiguous factor contained in Carr's condition (2.15) based on a well-motivated way [31] and make the condition more reliable.The basic idea is straightforward.We strictly define t s and t ff and explicitly calculate these time-scales based on the three-zone model.To do that, let us rewrite the metric of the overdense region in terms of the new radial coordinate χ defined by sin 2 ( as follows: where we have also introduced the conformal time defined by adη = dt.The Friedmann equation can be written as (2.19) Introducing u := (a/a max ) 1+3w , we obtain du dη The solution is given by u where the value of η is restricted to The times η = 0 and η = 2η max correspond to the big-bang and big-crunch, and the maximum expansion is realized at η = η max .A spacetime diagram for the closed universe with w = 1/3 is shown in Fig. 3. Let us reinterpret the Jeans criterion.The propagation time scale of the sound wave can be defined as the time required for the sound wave to propagate from the edge of the overdense region to the center.The free-fall time scale can be simply replaced by the collapsing time of the closed universe describing the overdense region.Therefore if the sound wave emanated from the edge of the overdense region at the time of the maximum expansion reaches the center before the big crunch, we conclude that the gravitational collapse will be prevented by the pressure gradient effect.
The sound speed c s is given by c s = √ w.Then the trajectory of the sound wave is parallel to the straight line χ = ± √ wη.From Fig. 3, we find that the Jeans criterion gives the threshold of the radius of the closed universe describing the overdense region for PBH formation.This is quite reasonable because the edge radius r c is related to the amplitude of the density perturbation at the horizon entry by Eq. (2.11).The threshold value is calculated by the following equation: where χ c is the value of the coordinate χ at the edge of the overdense region, that is, defined by sin It should be noted that this value is for the density perturbation in the uniform Hubble gauge.
In the comoving gauge, as will be shown later, the factor 3(1 + w)/(5 + 3w) must be multiplied, namely, In this section, we have refined Carr's condition based on the three-zone model.However, our estimation does not go much beyond order estimation.For more rigorous quantitative estimation, we need to rely on numerical simulations, and the threshold value somewhat depends on the profile of the inhomogeneity.Nevertheless, it has been reported that the threshold (2.26) approximately reproduces the lower bound of the threshold value [35].Therefore, in terms of results, the three-zone model roughly provides the least possible effect of the pressure gradient against PBH formation.

Cosmological long-wavelength solutions
To investigate PBH formation more accurately, we need to investigate the inhomogeneities which are initially super-horizon scales.Such inhomogeneities are described by growing mode solutions in the long-wavelength approximation.The expansion of the long-wavelength approximation with the expansion parameter ϵ = k/(aH) is often called gradient expansion.In the context of PBH formation, there are two main formulations of gradient expansion.One is the gradient expansion of the equations in the Misner-Sharp formulation [71] in which the Einstein equations for a fluid system in spherical symmetry are written down with the comoving coordinates.The other is based on the cosmological 3+1 decomposition [28].The relation between these two formulations is investigated in Ref. [34].Although the Misner-Sharp formulation has been adopted for PBH studies by many authors including Polnarev and Musco [72], in this paper, we review the formulation adopted in Ref. [28] because the latter allows more general gauge conditions.We follow and quote the discussions and calculations in Refs.[28,34] in this section.Let us start with writing down the Einstein equations in the form of the cosmological 3+1 decomposition for a perfect fluid component.We just display the necessary equations for later discussions and ask readers to refer to Ref. [34] for more details.

Cosmological 3+1 decomposition
First, we write the line element of spacetime as follows: where a is the scale factor of the background universe, and the lapse function α, shift vector β i , spatial conformal factor ψ and γij are functions of t and x i .We choose γij so that γ := det(γ ij ) = f := det(f ij ), with f ij being a time-independent metric of the flat three-space.The Latin indices run over 1 to 3 and we drop and raise the Latin indices of quantities without/with tilde by γ ij /γ ij and γ ij /γ ij , respectively.The extrinsic curvature of a time slice can be decomposed as where Ãij satisfies γij Ãij = 0. (3.3) Then the Hamiltonian and momentum constraint equations can be written as follows: where E = n µ n ν T µν and J i = −γ iµ n ν T µν with n µ , γ µν and T µν being the unit-normal form of the time slice, spatial metric and stress-energy tensor, respectively.△, Di and Rij are the Laplacian, covariant derivative and the Ricci tensor for γij , respectively.The evolution equations for geometrical variables are given by where D i and D i are the covariant derivatives for γ ij and the 3-dim flat metric f ij , respectively, and R ij is the Ricci tensor for γ ij .S ij is defined by γ iµ γ jν T µν .For a perfect fluid, the stress-energy tensor can be written as where u µ is the four-velocity of the fluid element which is normalized as u µ u µ = −1.Introducing we can write the hydrodynamical equations in the form and

Gradient expansion
In practice, the gradient expansion is performed by replacing the spatial derivative ∂ i by ϵ∂ i with a fictitious parameter ϵ.Expanding the equations of motion in a power series of ϵ, finally, we set ϵ = 1.Letting k be the comoving wavenumber for the scale of inhomogeneity of interest, we may consider |∂ i | ∼ k.Since the only natural scale of the system other than k is H, the fictitious parameter can be regarded as ϵ ∼ k/(aH).We also have to make key assumptions for the leading order ϵ → 0 to perform the gradient expansion.At the leading order, we assume that the universe is given by a spatially flat FLRW universe.More specifically, we require In addition, we assume that ψ is identically unity somewhere in the universe, so that a(t) is the scale factor for that part of the universe.Under these assumptions, it has been shown that ∂ t γij temporally decays and one can set for growing mode solutions [73].
In this review, we mainly focus on two specific time slicing conditions: the constant mean curvature (CMC) slicing (K = −3H) and comoving slicing (n µ = u µ ).Readers may refer to Ref. [34] for other gauge conditions.The comoving slicing implies u i = u t (v i + β i ) = 0, then Γ = 1, E = ρ and J i = (ρ + p)u i = 0. Since Ãij = O(ϵ 2 ) from Eq. (3.6), from the momentum (3.5) and Hamiltonian (3.4) constraints, we obtain Di K = O(ϵ 3 ) and H 2 = 8π 3 ρ + O(ϵ 2 ).Therefore CMC, uniform density and comoving slicings coincide to O(ϵ).In addition, from Eq. (3.13), we obtain where we have assumed that the pressure is homogeneous to O(ϵ).Since the right-hand side is a function of t, ∂ t ln ψ is also homogeneous to O(ϵ).Since we have assumed that ψ is identically unity at some point, we conclude Coming back to the general expressions without specifying the slicing condition, from ∂ t ψ = O(ϵ 2 ) and the equations of motion displayed in the previous section, one can find where ρ b is the background density.Then we introduce the following perturbation variables: where Ψ(x) = O(ϵ 0 ) and ξ, h ij , κ, χ and δ are the perturbation variables of O(ϵ 2 ).
Hereafter, we focus on the case of the linear equation of state p = wρ and consider the normal threading β i = 0 for simplicity.Readers may refer to Ref. [34] for other general cases.The equations of motion can be reduced to the following: The solutions for Eqs.(3.31) and (3.32) are given by respectively, where To find solutions for other quantities, we need to fix the slicing condition.Let us consider CMC slicing (κ = 0).From Eq. (3.29), we obtain where q(x) = − 4 3 Then, from Eqs. (3.30), (3.28) and (3.27), χ, ξ and u j are given as One can check that the solutions solve the other equations not used to derive the solutions.
For the comoving slicing (u j = 0), from Eq. (3.27), we obtain Substituting this relation into Eq.(3.26) and (3.30), and combining them, we can derive the following solutions for δ and κ: Then χ can be given as From Eq. (3.28), we obtain Since the function Ψ(x), and therefore q(x) are shared in the both gauge conditions, we can easily find the relation of gauge-dependent quantities.For instance, we find where δ CMC and δ co are density perturbations for CMC and comoving gauges, respectively.

PBH formation criterion
In this section, we consider spherically symmetric systems.In earlier days, the amplitude of the density perturbation averaged within a certain radius is used to give an analytic PBH formation criterion for convenience.More specifically, letting δ be the amplitude of the volume average of the density perturbation at the horizon entry, if δ exceeds a certain threshold value, we think that PBH formation is realized.The threshold value has been estimated analytically and numerically by many authors [27, 28, 30-33, 70, 72, 74-77].
A recent trend is to use the "compaction function" proposed in Ref. [28].There are several efforts to refine the simple treatment of the compaction function for a more accurate and convenient criterion of PBH formation [35][36][37][38].In this review, however, to avoid complicated discussions, we do not go into details about the recent refinements, but just review the simple use of the compaction function.That is, if the maximum value of the compaction function exceeds a certain threshold value, we think that PBH formation is realized.In this procedure, the threshold value depends on the profile of the initial perturbation.Nevertheless, we expect that, for moderate profiles, the deviation due to the profile dependence would be about 10%.Regarding the compaction function, recently, we found a misunderstanding that has not been realized for a long time [78].This misunderstanding will be also resolved.
Let us start with defining Shibata-Sasaki's original compaction function in terms of the cosmological long-wavelength perturbations as follows: where R is the areal radius of the sphere of the coordinate radius r given by R = arΨ2 .By using the solution (3.37), one can derive the following expression [34]: From this expression, we find that C SS is time independent and calculated by the gauge independent quantity Ψ.Furthermore, Ψ is related to the curvature perturbation, which is often used in the context of cosmological perturbation in the early universe, as ζ = − 1 2 ln Ψ 2) .Let us introduce the radius r m at which the compaction function takes the maximum value C max SS = C SS (r m ).Then, we adopt the following PBH formation criterion: with the uncertainty coming from the profile dependence.As is shown in Figs. 2 and 3 in Ref. [34] for two specific one-parameter families of initial profiles, the value of C th can be roughly estimated as 0.4 with about 10% uncertainty.For the specific profile given by ln Ψ ∝ sin(r)/r, which corresponds to the typical profile for the monochromatic power spectrum of ln Ψ, we obtain C th ≃ 0.44 from numerical simulations.Therefore we use the value C th = 0.44 as a reference value in this paper.
In order to understand the compaction function more deeply, let us review the relation between the compaction function and the conserved mass.In spherical symmetry, we have a well-behaved conserved mass called Misner-Sharp mass [71].As is shown in the Appendix, the Misner-Sharp mass is equivalent to the Kodama mass defined by the charge associated with a conserved current.In addition, the value of the Misner-Sharp mass on a marginally trapped surface is given by half of the areal radius of the marginally trapped surface.Therefore the Misner-Sharp mass is closely related to horizon formation, and relevant to be used for horizon formation criterion.
From the expression of the Kodama current and mass given in Eqs.(A. 16) and (A.19), we can derive the following expression: In the comoving slice (u r = 0), we obtain The mass excess compared to the background FLRW with the identical areal radius is calculated as Then we can define the compaction function in the comoving gauge as where R = aΨ 2 r.Therefore C SS is the compaction function in the comoving gauge multiplied by 3(1+w)  3w+5 .However, it should be noted that, for the CMC slice, the contribution from the second term in the integrand of Eq. (4.4) cannot be neglected 3) and we find C CMC ̸ = C SS .For a more concrete description, readers may refer to Ref. [78].
To see the relation between the compaction function and the volume average of the density perturbation, let us define the volume average of the density perturbation in the comoving slice as follows: Then we can easily find δ = 2 Therefore we obtain δ HR=1 = 2C co .(4.10) That is, the compaction function can be also regarded as the amplitude of the volume average of the density perturbation at the horizon entry.This relation also supports the relevance of the compaction function to be used for the criterion of PBH formation.

Estimation of PBH abundance in peak theory
We start with writing the spatial metric in the following form: The curvature perturbation ζ is related to Ψ as Here we assume that ζ is the Gaussian random variable with the power spectrum P(k) defined by 3) In Refs.[28,34], the second term in the integrand of Eq.
where △ is the flat Laplacian.Here we just write down the final expressions (D.3) and (D.5): ) where P 1 , f , R n and ψ n are defined in Eqs.(B.18), (B.30), (C.18) and (D.2), respectively, and we omitted the last term in (D.4), which can be absorbed into the renormalization of the scale factor.
A flow chart to reach the PBH mass spectrum by using the peak number density (5.7), typical profile (5.8) and the PBH formation criterion (4.3) is shown in Fig. 4. Let us explain each step of the flow chart.

Step 1: Rewriting the criterion
First, let us consider the compaction function for the typical profile.Hereafter we will work in the comoving gauge and express the compaction function in the comoving gauge without the subscript, that is, C = C co .Focusing on the radiation-dominated universe (w = 1/3), the PBH formation criterion can be rewritten as In terms of the typical profile ζ, the compaction function can be written as Step 1 Step 2 Step 3 Through this equation, we can obtain the threshold value of μ2 as a function of k • as where

Step 2: Derivation of the PBH mass expression
In order to obtain the number density for a fixed PBH mass, we need to know the relation between the parameters μ2 , k • and the PBH mass M .The typical mass of a PBH can be roughly estimated by the horizon mass of 1/(2H) at the horizon entry.The horizon entry condition is given by (5.16) In the radiation-dominated universe, adopting the matter-radiation equality as a reference time, we find a = a eq H eq H 1/2 . (5.17) From Eqs. (5.16) and (5.17), at the horizon entry, we obtain eq H eq r 2 m e −2μ 2 gm = H −1 eq k 2 eq r 2 m e −2μ 2 gm , (5.18) where k eq = a eq H eq .Then the PBH mass can be estimated as where we have introduced M eq = 1/(2H eq ).We have not taken into account the mass scaling associated with the critical behavior [82,83] in the expression (5.19).It should be also noted that the mass estimation given in this paper cannot be applied to Type II PBH formation [69].

Step 3: Derivation of the mass function
Eq. (5.19) gives us the relation between M , k • and μ2 .Therefore we can change the independent variables in the peak number density (5.7) from μ2 and k • to μ2 and M as follows: where k • should be regarded as a function of μ2 and M through Eq. (5.19).To obtain the PBH number density for a given value of the mass M , we need to integrate Eq. (5.20) for μ2 .
As is explicitly shown in Ref. [67] for some specific cases, the threshold value of μ2 given by Eq. (5.15) can be converted to the threshold value μ(M) 2th (M ) for a given value of M by using Eqs.(5.15) and (5.19).In addition, the value of μ2 may be bounded below for a fixed value of M .Then, letting the minimum possible value of μ2 be μ2min (M ), the lower limit of the integral for μ2 is given by μ2b := max μ2min (M ), μ(M) 2th (M ) .That is, the PBH number density n BH (M ) and the mass spectrum f 0 (M ) at equality time are given by (5.21) eq n BH d ln M. (5.22)

Critical behavior
It is well known that, at least for idealized spherically-symmetric situations, the mass of the black hole scales as ∝ (b−b th ) γ with b and b th being one parameter characterizing the initial data and the threshold value of it, respectively.PBH formation is no exception, and its realization and the effects on the PBH abundance have been investigated [26,27,30,76,77,80,[84][85][86][87][88].If we take the critical behavior into account, the PBH mass would be given by where γ is the universal exponent given by γ ≃ 0.36 for radiation fluid, and K(k • ) is a profiledependent numerical factor, which will be set to be 1 hereafter for simplicity.In practice, we can consider the case without critical behavior setting the value of γ to be 0 in Eq. (5.23).

The monochromatic curvature power spectrum
In this paper, we show only the case of the monochromatic curvature power spectrum.Readers may refer to Refs.[66,67] for extended spectra.Let us write the monochromatic power spectrum as Then the gradient moments are calculated as (5.25) The typical profile is given as From the condition (5.9), we can find μ2th ≃ 0.615. (5.28) Then the PBH mass can be given by a function of μ2 as where γ = 0.36 and 0 with and without taking into account the critical behavior, respectively.The minimum value of the PBH mass M min is given by 0 and M eq k 2 eq r 2 m e −2μ 2th gm for γ = 0.36 and 0, respectively.
For the number density, we reconsider the expression (5.7).In the limit γ 3 → 0, the probability distribution P 1 (ν, ξ 1 ; γ 3 ) can be written as lim (5.30) Substituting ν = µ 2 /σ 2 and ξ 1 = µ 2 k 2 • /σ 4 and using σ n = σ 0 k 0 , we obtain Integrating the number density for k • , we obtain (5.32) The variable μ2 can be converted into M through Eq. (5.29).Then we obtain where μ2 is regarded as a function of M through (5.29), and The PBH fraction f 0 is given by Since small values of the PBH mass can be realized for μ2 ∼ μ2th , the behavior of f 0 in the small PBH mass limit is given by f 0 ∝ M 1+γ −1 .We show the PBH fraction f 0 (M ) for k 0 /k eq = 10 5 as a function of M/M 0 in Fig. 5 with M 0 = M eq k 2 eq /k 2 0 .The total fraction f tot = f 0 d ln M is Figure 5: The PBH fractions f 0 are plotted as functions of M/M 0 for σ 0 = 0.08, 0.10 and 0.12 with k 0 = 10 5 k eq .The right spiky plots show the PBH fractions without the critical behavior, and the left broader plots show those with the critical behavior.
shown as functions of σ 0 for the cases with and without critical behavior in Fig. 6 6 Summary and further related topics In this review, we have run through the basic picture of PBH formation, cosmological longwavelength solutions, PBH formation criterion and estimation of PBH abundance in peak theory.In this last section, let us summarize and make notes on each topic.First, in Sec. 2, we have considered the three-zone model.It should be noted that the threezone model cannot be an exact solution for the Einstein field equations with non-vanishing pressure of the fluid, and just gives us intuition about the PBH formation process.Nevertheless, it is interesting that the three-zone model provides us a quantitative threshold value that is roughly equal to the lower bound of the threshold value of the density perturbation [35].Similar ideas to the three-zone model would be very useful as a foothold when we extend our focus to PBH formations in other settings.For more accurate analyses, the cosmological long-wavelength solutions are needed to be considered.In Sec. 3, we just quoted the equations and solutions from Refs.[28,34].Readers can refer to Refs.[28,34] for more general and concrete derivations and discussions.We note that in Sec. 3, we have not assumed spherical symmetry.In this sense, the formulation given in Sec. 3 is more general than the gradient expansion in the Misner-Sharp formulation (e.g., Ref. [72]) which is widely used to construct the initial data for numerical simulations.When one considers other matter components than the perfect fluid with the linear equation of state p = wρ, the cosmological long-wavelength solution should be reconsidered.For example, in Ref. [89], cosmological long-wavelength solutions for massless scalar isocurvature modes have been investigated.
In Sec. 4, assuming spherical symmetry, we have considered the criterion of PBH formation.Although we need numerical simulations for an accurate criterion, some of the approximate analytic criteria are still practically useful.As a useful and relatively simple criterion, we have introduced the criterion based on the compaction function.The definition of the compaction function is carefully introduced in Sec. 4 because we recently found a misunderstanding that has not been realized for a long time [78].More concrete discussions will be reported elsewhere [78].Since the compaction function is closely related to the conserved mass in a spherically symmetric spacetime, we have provided a short review of the conserved mass and horizons in spherically symmetric spacetimes in Appendix A.
One of the goals of the theoretical PBH study is the estimation of PBH abundance.Since PBHs would be associated with rarely high peaks of a perturbation variable, the peak statistics are relevant to be used.We have attached a short minimal review of the peak theory in Appendix B and C. We have explained the abundance estimation following Refs.[66,67] mainly focusing on the case of the monochromatic curvature power spectrum.However, different procedures based on peak theory have been also proposed [87,[90][91][92].Here we note that, while the variation of the peak curvature scale is treated as an independent statistical variable in Refs.[66,67,90], the typical profile for the mean peak curvature is only considered in Refs.[87,91].We may expect that those different procedures predict similar PBH mass spectra for a narrow curvature power spectrum.However, the predictions may be different from each other for a broad curvature power spectrum.This issue is still open and a direct comparison between those different procedures is needed to make a deeper understanding and a consensus.
where we have fixed the sign on the right-hand side such that S t > 0 for K t > 0.
Let us check that the Kodama current is conserved.By using the Einstein equations, we can rewrite the Kodama current (A.16) as follows: (A.17) Then we obtain where we have used ∂ µ (ϵ AB /e ϕ+λ/2 ) = 0. Therefore we can define the following conserved charge (Kodama mass): where dΣ is the volume element of the t = const.hypersurface whose unit normal form is given by n µ = −e ϕ (dt) µ .This equation explicitly shows the equivalence between the Kodama mass and the Misner-Sharp mass.

B Peak number density in peak theory
In this and the next section, we briefly review the necessary part of the peak theory [68] for the calculation of PBH abundance.Let us consider the Gaussian variable with the power spectrum P(k) defined by The probability distribution of linear combinations of ζ(x) denoted by given by the multivariate Gaussian distribution: where For example, let us consider the probability distribution of −ζ and △ζ at a spatial point, where △ is the flat Laplacian.The correlation matrix can be calculated as Then the probability distribution of −ζ and △ζ is given by Similarly to the example of −ζ and △ζ, let us consider the probability distribution function of the coefficients in Taylor expansion of ζ up through the second order: There are 10 independent variables ζ 0 , ζ i 1 (i = 1, 2, 3) and ζ ij 2 (i, j = 1, 2, 3).Non-zero correlations between two of these coefficients are summarized as follows: The 2nd order variables can be transformed into the six independent variables λ i and θ i with i = 1, 2, 3, where λ 1 , λ 2 and λ 3 are the three eigenvalues of the matrix ζ ij 2 with λ 1 ≥ λ 2 ≥ λ 3 and θ 1 , θ 2 and θ 3 are the Euler angles to take the principal direction.Then it can be shown that the volume element d 6 ζ 2 can be given by The integration with respect to the Euler angles gives the factor 2π.We introduce the following 7 independent variables by using the remaining 7 variables: Then the probability distribution is transformed into the following form: P (ν, ξ, η)dνd 3 ξd 3 η = P 1 (ν, ξ 1 ; γ 1 )P 2 (ξ 2 , ξ 3 )P 3 (η)dνd 3 ξd 3 η, (B.17 where It would be worthy to note that there are no correlations between (ν, ξ 1 ) and (ξ 2 , ξ 3 ).Therefore even if we consider a rarely high peak with ν ∼ ξ 1 ≫ 1, the values of ξ 2 and ξ 3 are expected to be ∼ O(1) ≪ ξ 1 .Since ξ 2 and ξ 3 are responsible for the asphericity of the peak, we can conclude that rarely high peaks tend to be more spherically symmetric.
Let us remember that what we are seeking is the mean number density of extrema characterized by the amplitude and curvature of the peak.The corresponding variables to the amplitude and the peak curvature would be ν and ξ 1 .Therefore our goal in this section is to obtain the mean number density characterized by the parameters ν and ξ 1 : npeak (ν, ξ 1 ).As a first step, let us consider the number density distribution of extrema n ext (x, ν, ξ 1 ) in the space spanned by (x, ν, ξ 1 ).By definition, we have where p labels the extrema so that the p-th extremum will be characterized by (x p , ν p , ξ 1p ).
The integration over the space spanned by (x, ν, ξ 1 ) trivially gives the total number of extrema.This expression of the number density distribution is the specific distribution for a given set of the parameters (x p , ν p , ξ 1p ).Since we are interested in the average distribution without the explicit dependence on (x p , ν p , ξ 1p ), let us take the volume average and the ensemble average for ν p and ξ 1p of the distribution as follows: next (ν, ξ 1 )∆ν∆ξ To dispose the δ(x − x p ) and the volume integral, let us consider rewriting the delta function.
Since we are focusing on extrema, at x = x p , η vanishes.Therefore we find where |λ 1 λ 2 λ 3 | can be written as in terms of ξ.Then we obtain next (ν, ξ 1 ) = P 1 (ν, ξ 1 ; γ 1 ) σ 3 2 where we have fixed the value of ξ 1 and take ξ 1 as an independent variable of x.Here we replace the volume average with the ensemble average for ξ 2 , ξ 3 and η as follows: The peak number density can be obtained by multiplying the step function Θ(λ 3 ).Then the integration for ξ 2 and ξ 3 can be explicitly done [68] and we obtain where For convenience in the calculation of PBH abundance, let us change the variables from (ξ 1 , ν) to (µ 0 , k * ) defined by Then the peak number density can be transformed into the following form:

C Typical profile in peak theory
The peak theory also provides a typical profile characterized by the peak parameter ν and ξ 1 [68].To derive the expression of the typical profile, let us consider the conditional probability of taking the value ζ(r) at the radius r from a peak.Before going into details, we define some notations about conditional probabilities.We represent the probability for the realization of the event X as P X .The conditional probability for the realization of Y under the condition X is denoted by P Y |X and given by where P X∩Y is the probability for the realization of X and Y .Let us consider the multi-variate Gaussian probability for the probability variables given by where α a (a = 1, • • • , n) and β i (i = 1, • • • , N − n) correspond to the variables for X and Y , respectively.Letting q µ denotes the components of the transposition of the row vector q as q µ = (q t ) µ , the probability distribution of X ∪ Y can be written as where M is the N ×N correlation matrix.We decompose the correlation matrix in the following form: By using the expression of M −1 , we find where we have used the fact B = C t .Since the conditional probability is given by Eq. (C.1), we find Therefore the mean values and the correlation matrix for β under the condition X are given by αA −1 B and M/A.We are interested in the probability of β = ν(r) = ζ(r)/σ 0 under the condition α a = (ν(0), ξ 1 (0), η(0) = 0).To simplify the expressions we write ν(0), ξ 1 (0) and η(0) without the argument as ν, ξ 1 and η, but the argument of ν(r) will be kept.Then the 4 blocks of the correlation matrix can be calculated as The variance is given by The variance is roughly estimated as ∆ν(r) ∼ ∆ζ/σ 0 = O(1), that is, ∆ζ ∼ σ 0 .Therefore the typical profile is a good approximation as long as ζ(r) ≫ σ 0 for a rarely high peak.Changing the independent variables from (ν, ξ 1 ) to (µ 0 , k * ), we obtain the following form of the typical profile of ζ(r):

D Equations for peaks of △ζ
In practice, peaks of △ζ are more important for the calculation of PBH abundance.Thus we consider the parameters µ 2 = △ζ(0) and k 2 • := − △△ζ(0) µ 2 instead of µ 0 and k 2 * .With this modification associated with −ζ → △ζ, to obtain the number density and typical profile of peaks of △ζ, we need the following replacements: where 4)   ψ n := 1 σ n d ln k k 2n sin(kr) kr P(k).(D.2) Then, introducing the dimensionless quantity μ2 = µ 2 σ 2 1 /σ 2 2 , we obtain the number density and the typical profile as follows: The profile of ζ associated with △ζ can be obtained by integrating (D.4) as follows: It can be shown that ζ ∞ can be also regarded as a Gaussian probability variable (see Ref. [67]).
4) The definition of ψ 2 is different from that in Refs.[66,67], but the same as that in Ref. [80].

Figure 1 :
Figure 1: The rough sketch of the standard evolution of the background universe and the perturbations on it.

Figure 2 :
Figure 2: Schematic figure of the three-zone model.

Figure 3 :
Figure 3: Spacetime diagram for the closed universe with w = 1/3 describing the overdense region and the trajectory of the sound wave propagation from the edge to the center emanated at the time of the maximum expansion.
(4.4)  has been erroneously dropped, and the authors got the wrong relation C CMC = C SS .The original definition of C SS given in Eq. (4.1) would have been provided based on this wrong calculation neglecting the difference between the CMC and comoving slice conditions.
with ζ(k) = d 3 xζ(x)e ikx .Readers may refer to Refs.[79][80][81] for non-Gaussian cases.From the random Gaussian assumption and the gradient moments σ n given by σ 2 n = d ln kk 2n P(k), (5.4) as is shown in Appendices B, C and D, we can obtain the peak number density n (k•) pk (μ 2 , k • ) and the typical peak profile ζ(r) characterized by the two parameters μ2 and k • defined by μ2 = △ζ(0)

Figure 4 :
Figure 4: A flow chart to reach the PBH mass spectrum.

Figure 6 :
Figure 6: The total fraction f tot as a function of σ 0 .