10

# The basics of primordial black hole formation and
abundance estimation

Chul-Moon Yoo<EMAIL_ADDRESS>Division of Particle and
Astrophysical Science, Graduate School of Science, Nagoya University, Nagoya
464-8602, Japan

###### Abstract

This paper is a biased review of the 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.

###### Contents

  1. 1 Introduction
  2. 2 PBH formation process: three-zone model
    1. 2.1 FLRW models
    2. 2.2 Background flat FLRW
    3. 2.3 Super-horizon inhomogeneity
    4. 2.4 Three-zone model
    5. 2.5 Jeans criterion
    6. 2.6 Refinement of the threshold
  3. 3 Cosmological long-wavelength solutions
    1. 3.1 Cosmological 3+1 decomposition
    2. 3.2 Gradient expansion
  4. 4 PBH formation criterion
  5. 5 Estimation of PBH abundance in peak theory
    1. 5.1 Step 1: Rewriting the criterion
    2. 5.2 Step 2: Derivation of the PBH mass expression
    3. 5.3 Step 3: Derivation of the mass function
    4. 5.4 Critical behavior
    5. 5.5 The monochromatic curvature power spectrum
  6. 6 Summary and further related topics
  7. A Horizon and mass
  8. B Peak number density in peak theory
  9. C Typical profile in peak theory
  10. D Equations for peaks of $\triangle\zeta$

##  1 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 recent years, the possibility of PBHs is seriously considered and
observational constraints have been actively investigated. While we have not
found evidence of the existence of PBHs, there are three interesting
candidates: dark matter [4, 5, 6], black hole binaries observed by
gravitational wave interferometers [7, 8, 9, 10] and the earth-mass objects of
microlensing events [11]. In addition, PBHs may explain a potential signal at
a pulsar timing array [12, 13, 14, 15, 16, 17, 18] and the seeds of
supermassive black holes and cosmic structures [19, 20, 21, 22, 23, 24, 25].
It is very exciting to clarify whether those candidates are PBHs or not.

There are several aspects of PBH research: early universe models which provide
a substantial number of PBHs, dynamical process of PBH formation [26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41] effects of the evaporation
of PBHs[42, 43, 44, 45, 46, 47, 48, 49, 50], associated gravitational wave
background[51, 52, 53, 54, 55, 56, 57, 58], observational consequences and
constraints [5] and so on. The purpose of this review is to give a brief
introduction to the PBH formation process and show a recently developed method
for the abundance calculation with the PBH formation criterion and the peak
statistics for the Gaussian curvature perturbation. We would like to ask
readers to refer to other recent reviews and lecture notes (e.g., Refs. [59,
60, 61, 10]) on PBHs for the scientific significance and historical reviews.

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. [62, 63]. 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 [64] 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$.

##  2 PBH formation process: three-zone model

First, in order to understand the overall picture of PBH formation, let us
review the three-zone 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.

###  2.1 FLRW models

Let us write the line element of an FLRW model as follows:

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\left[\frac{\mathrm{d}r^{2}}{1-Kr^{2}}+r^{2}\mathrm{d}\Omega^{2}\right],$
(2.1)

where $a(t)$, $K$ and $\mathrm{d}\Omega$ are the scale factor, spatial
curvature and line element of the unit two-sphere. One of the Friedmann
equations is given by

$\left(\frac{\dot{a}}{a}\right)^{2}=\frac{8\pi}{3}\rho-\frac{K}{a^{2}},$ (2.2)

where the dot ”$\dot{~}{}$” denotes the time derivative. Introducing the areal
radius $R:=a(t)r$, Eq. (2.2) can be rewritten as follows:

$\frac{\dot{R}^{2}}{2}-\frac{M}{R}=E(r):=-\frac{Kr^{2}}{2},$ (2.3)

where $M$ is the mass inside the comoving radius $r$ given by $4\pi\rho
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_{\rm max}=M/E(r)$. At the maximum expansion, the density
$\rho_{\rm max}$ and the scale factor $a_{\rm max}$ satisfy

$\frac{8\pi}{3}\rho_{\rm max}=\frac{K}{a_{\rm max}^{2}}.$ (2.4)

In contrast, the fluid continues to expand if $K\leq 0$.

###  2.2 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\rho$ with a constant $w$, where $\rho$
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 ($\Lambda$),
respectively. From the conservation law of the fluid $\mathrm{d}(\rho
a^{3})=-p\mathrm{d}a^{3}$, we obtain

$\rho\propto a^{-3(1+w)}.$ (2.5)

Then, we find the behavior of the Hubble radius as

$R_{\rm H}:=1/H\propto a^{3(1+w)/2}.$ (2.6)

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_{\rm H}\sim{\rm const.}$) is considered before the
decelerated expansion phases of our universe. The inflationary era may be
followed by an early MD phase ($R_{\rm H}\sim a^{3/2}$) supported by the
inflaton oscillation, and the universe supposed to become RD universe ($R_{\rm
H}\sim 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.

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

###  2.3 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\ll aH$, will be summarized in Sec. 3.
Here let us take a result in advance. At the leading order of the long-
wavelength perturbation ($\epsilon:=k/(aH)\ll 1$), we obtain the following
metric form:

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}\Psi(\bm{x})\delta_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j},$
(2.7)

where $\delta_{ij}$ is the Kronecker delta and $\Psi(\bm{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:

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}\frac{1}{\left(1+\frac{K}{4}|\bm{x}|^{2}\right)^{2}}\delta_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j}.$
(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.

###  2.4 Three-zone model

Motivated by the previous subsection, let us consider the three-zone model
composed of the spherical closed FLRW region $r<r_{\rm c}$ surrounded by an
under-dense region embedded in the flat FLRW background (Fig. 2).

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

This model can be an exact solution only for the dust case ($w=0$)1)1)1)The
under-dense region is necessary for the dust case to be an exact solution. .
In the case of $w\neq 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.

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

$\frac{8}{3}\pi\rho_{\rm b}=\frac{8}{3}\pi\rho-\frac{K}{a^{2}},$ (2.9)

where $\rho$ and $\rho_{\rm 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

$\delta:=\frac{\rho-\rho_{\rm b}}{\rho_{\rm b}}=\frac{K}{a^{2}H^{2}}.$ (2.10)

The value of the density perturbation $\delta_{\rm H}$ at the horizon entry
$R_{\rm c}:=ar_{\rm c}=1/H$ is given by

$\delta_{\rm H}=Kr_{\rm c}^{2}.$ (2.11)

Since the value of the comoving radius $r$ cannot be larger than $\sqrt{K}$,
we find the following geometrical upper bound [65]:

$\delta_{\rm H}<1.$ (2.12)

###  2.5 Jeans criterion

If the whole dynamics of the overdense region are 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_{\rm s}$ would be
given by

$t_{\rm s}=\frac{R_{\rm c}}{\sqrt{w}}=\frac{ar_{\rm c}}{\sqrt{w}}\propto a,$
(2.13)

where we have used the fact that the sound speed $c_{\rm s}$ is given by
$\sqrt{w}$. The free-fall time-scale $t_{\rm ff}$ is given by

$t_{\rm ff}=\frac{1}{H}=\left(\frac{8\pi}{3}\rho\right)^{-1/2}=\frac{a_{\rm
max}}{\sqrt{K}}\left(\frac{a}{a_{\rm max}}\right)^{3(1+w)/2}\propto
a^{3(1+w)/2}.$ (2.14)

Let us define $a_{\rm th}$ as the value of the scale factor satisfying $t_{\rm
s}|_{a=a_{\rm th}}=t_{\rm ff}|_{a=a_{\rm th}}$. Then if $a_{\rm max}>a_{\rm
th}$, the sound-wave propagation time scale is shorter than the free-fall time
scale at the time ${a=a_{\rm 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:

$a_{\rm max}<a_{\rm th}=a_{\rm max}\left(\frac{w}{Kr_{\rm
c}^{2}}\right)^{-1/(1+3w)}\Leftrightarrow Kr_{\rm c}^{2}=\delta_{\rm H}>w.$
(2.15)

This condition is often quoted as Carr’s condition [66]. Here it should be
noted that Carr’s condition has the ambiguity of a factor of $\mathcal{O}(1)$
coming from the order estimations of Eq. (2.13) and Eq. (2.14).

###  2.6 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_{\rm s}$ and $t_{\rm 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 $\chi$ defined by

$\sin^{2}(\sqrt{K}\chi)=Kr^{2}$ (2.16)

as follows:

$\displaystyle\mathrm{d}s^{2}$ $\displaystyle=$
$\displaystyle-\mathrm{d}t^{2}+a^{2}\left[\mathrm{d}\chi^{2}+\frac{1}{K}\sin^{2}(\sqrt{K}\chi)\mathrm{d}\Omega^{2}\right]$
(2.17) $\displaystyle=$ $\displaystyle
a^{2}\left[-\mathrm{d}\eta^{2}+\mathrm{d}\chi^{2}+\frac{1}{K}\sin^{2}(\sqrt{K}\chi)\mathrm{d}\Omega^{2}\right],$
(2.18)

where we have also introduced the conformal time defined by
$a\mathrm{d}\eta=\mathrm{d}t$. The Friedmann equation can be written as

$\left[\frac{\mathrm{d}}{\mathrm{d}\eta}\ln
a\right]^{2}=K\left[\left(\frac{a}{a_{\rm max}}\right)^{-1-3w}-1\right].$
(2.19)

Introducing $u:=(a/a_{\rm max})^{1+3w}$, we obtain

$\left(\frac{\mathrm{d}u}{\mathrm{d}\eta}\right)^{2}=K(1+3w)^{2}\left(u-u^{2}\right).$
(2.20)

The solution is given by

$u=\frac{1}{2}\left[1-\cos\left(\sqrt{K}(1+3w)\eta\right)\right],$ (2.21)

where the value of $\eta$ is restricted to

$0\leq\eta/2\leq\eta_{\rm max}:=\frac{\pi}{\sqrt{K}(1+3w)}.$ (2.22)

The times $\eta=0$ and $\eta=2\eta_{\rm max}$ correspond to the big-bang and
big-crunch, and the maximum expansion is realized at $\eta=\eta_{\rm max}$. A
spacetime diagram for the closed universe with $w=1/3$ is shown in Fig. 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.

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_{\rm s}$ is given by $c_{\rm s}=\sqrt{w}$. Then the
trajectory of the sound wave is parallel to the straight line
$\chi=\pm\sqrt{w}\eta$. 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_{\rm 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:

$\sqrt{w}\eta_{\rm max}=\chi_{\rm c},$ (2.23)

where $\chi_{\rm c}$ is the value of the coordinate $\chi$ at the edge of the
overdense region, that is, defined by

$\sin^{2}\sqrt{K}\chi_{\rm c}=\sqrt{K}r_{\rm c}.$ (2.24)

Combining Eqs. (2.11), (2.22), (2.23) and (2.24), we can derive the following
threshold value for the amplitude of the density perturbation at the horizon
entry:

$\delta_{\rm th}=\sin^{2}\left(\frac{\pi\sqrt{w}}{1+3w}\right).$ (2.25)

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,

$\delta^{\rm co}_{\rm
th}=\frac{3(1+w)}{(5+3w)}\sin^{2}\left(\frac{\pi\sqrt{w}}{1+3w}\right).$
(2.26)

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.

##  3 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 $\epsilon=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
[67] 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 [68], 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.

###  3.1 Cosmological 3+1 decomposition

First, we write the line element of spacetime as follows:

$\mathrm{d}s^{2}=-\alpha^{2}\mathrm{d}t^{2}+a^{2}\psi^{4}\tilde{\gamma}_{ij}(\mathrm{d}x^{i}+\beta^{i}\mathrm{d}t)(\mathrm{d}x^{j}+\beta^{j}\mathrm{d}t),$
(3.1)

where $a$ is the scale factor of the background universe, and the lapse
function $\alpha$, shift vector $\beta^{i}$, spatial conformal factor $\psi$
and $\tilde{\gamma}_{ij}$ are functions of $t$ and $x^{i}$. We choose
$\tilde{\gamma}_{ij}$ so that
$\tilde{\gamma}:=\det(\tilde{\gamma}_{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 $\gamma_{ij}$/$\tilde{\gamma}_{ij}$ and
$\gamma^{ij}$/$\tilde{\gamma}^{ij}$, respectively. The extrinsic curvature of
a time slice can be decomposed as

$K_{ij}=\psi^{4}a^{2}\left(\tilde{A}_{ij}+\frac{1}{3}\tilde{\gamma}_{ij}K\right),$
(3.2)

where $\tilde{A}_{ij}$ satisfies

$\tilde{\gamma}^{ij}\tilde{A}_{ij}=0.$ (3.3)

Then the Hamiltonian and momentum constraint equations can be written as
follows:

$\displaystyle\tilde{\triangle}\psi=\frac{1}{8}\tilde{\mathcal{R}}^{k}_{k}\psi-2\pi\psi^{5}a^{2}E-\frac{1}{8}\psi^{5}a^{2}\left(\tilde{A}_{ij}\tilde{A}^{ij}-\frac{2}{3}K^{2}\right),$
(3.4)
$\displaystyle\tilde{\mathcal{D}}^{j}(\psi^{6}\tilde{A}_{ij})-\frac{2}{3}\psi^{6}\tilde{\mathcal{D}}_{i}K=8\pi
J_{i}\psi^{6},$ (3.5)

where $E=n_{\mu}n_{\nu}T^{\mu\nu}$ and $J_{i}=-\gamma_{i\mu}n_{\nu}T^{\mu\nu}$
with $n_{\mu}$, $\gamma_{\mu\nu}$ and $T^{\mu\nu}$ being the unit-normal form
of the time slice, spatial metric and stress-energy tensor, respectively.
$\tilde{\triangle}$, $\tilde{\mathcal{D}}_{i}$ and $\tilde{\mathcal{R}}_{ij}$
are the Laplacian, covariant derivative and the Ricci tensor for
$\tilde{\gamma}_{ij}$, respectively. The evolution equations for geometrical
variables are given by

$\displaystyle(\partial_{t}-{\cal L}_{\beta})\tilde{\gamma}_{ij}$
$\displaystyle=$
$\displaystyle-2\alpha\tilde{A}_{ij}-\frac{2}{3}\tilde{\gamma}_{ij}{\cal
D}_{k}\beta^{k},$ (3.6) $\displaystyle(\partial_{t}-{\cal
L}_{\beta})\tilde{A}_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{a^{2}\psi^{4}}\left[\alpha\left({\cal
R}_{ij}-\frac{\gamma_{ij}}{3}{\cal
R}\right)-\left(D_{i}D_{j}\alpha-\frac{\gamma_{ij}}{3}D_{k}D^{k}\alpha\right)\right]$
(3.7)
$\displaystyle+\alpha(K\tilde{A}_{ij}-2\tilde{A}_{ik}\tilde{A}_{j}^{k})-\frac{2}{3}({\cal
D}_{k}\beta^{k})\tilde{A}_{ij}-\frac{8\pi\alpha}{a^{2}\psi^{4}}\left(S_{ij}-\frac{\gamma_{ij}}{3}S_{k}^{k}\right),$
$\displaystyle(\partial_{t}-{\cal L}_{\beta})\psi$ $\displaystyle=$
$\displaystyle-\frac{1}{2}H\psi+\frac{\psi}{6}(-\alpha K+{\cal
D}_{k}\beta^{k}),$ (3.8) $\displaystyle(\partial_{t}-{\cal L}_{\beta})K$
$\displaystyle=$
$\displaystyle\alpha\left(\tilde{A}_{ij}\tilde{A}^{ij}+\frac{1}{3}K^{2}\right)-D_{k}D^{k}\alpha+4\pi\alpha(E+S_{k}^{k}),$
(3.9)

where $D_{i}$ and ${\mathcal{D}}_{i}$ are the covariant derivatives for
$\gamma_{ij}$ and the 3-dim flat metric $f_{ij}$, respectively, and
$\mathcal{R}_{ij}$ is the Ricci tensor for $\gamma_{ij}$. $S_{ij}$ is defined
by $\gamma_{i\mu}\gamma_{j\nu}T^{\mu\nu}$.

For a perfect fluid, the stress-energy tensor can be written as

$T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu},$ (3.10)

where $u^{\mu}$ is the four-velocity of the fluid element which is normalized
as $u^{\mu}u_{\mu}=-1$. Introducing

$v^{l}:=\frac{u^{l}}{u^{t}}$ (3.11)

and

$\Gamma:=\alpha
u^{t}=\left[1-\alpha^{-2}(\beta_{k}+v_{k})(\beta^{k}+v^{k})\right]^{-1/2},$
(3.12)

we can write the hydrodynamical equations in the form

$\displaystyle\left[\psi^{6}a^{3}\left\\{(\rho+p)\Gamma^{2}-p\right\\}\right]_{,t}+\frac{1}{\sqrt{f}}\left[\sqrt{f}\psi^{6}a^{3}\left\\{(\rho+p)\Gamma^{2}-p\right\\}v^{l}\right]_{,l}$
$\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{f}}\left[\sqrt{f}\psi^{6}a^{3}p(v^{l}+\beta^{l})\right]_{,l}+\alpha\psi^{6}a^{3}pK-\alpha^{-1}\alpha_{,l}\psi^{6}a^{3}\Gamma^{2}(\rho+p)(v^{l}+\beta^{l})$
$\displaystyle+\alpha^{-1}\psi^{10}a^{5}\Gamma^{2}(\rho+p)(v^{l}+\beta^{l})(v^{m}+\beta^{m})\left(\tilde{A}_{lm}+\frac{\tilde{\gamma}_{lm}}{3}K\right),$
$\displaystyle(\Gamma\psi^{6}a^{3}(\rho+p)u_{j})_{,t}+\frac{1}{\sqrt{f}}(\sqrt{f}\Gamma\psi^{6}a^{3}(\rho+p)v^{k}u_{j})_{,k}$
$\displaystyle=$
$\displaystyle-\alpha\psi^{6}a^{3}p_{,j}+\Gamma\psi^{6}a^{3}(\rho+p)\left(-\Gamma\alpha_{,j}+u_{k}\beta^{k}_{,j}-\frac{u_{k}u_{l}}{2u^{t}}\gamma^{kl}_{,j}\right),$
(3.14)

and

$\displaystyle(\Gamma\psi^{6}a^{3}n)_{,t}+\frac{1}{\sqrt{f}}(\sqrt{f}\Gamma\psi^{6}a^{3}nv^{k})_{,k}=0.$
(3.15)

###  3.2 Gradient expansion

In practice, the gradient expansion is performed by replacing the spatial
derivative $\partial_{i}$ by $\epsilon\partial_{i}$ with a fictitious
parameter $\epsilon$. Expanding the equations of motion in a power series of
$\epsilon$, finally, we set $\epsilon=1$. Letting $k$ be the comoving
wavenumber for the scale of inhomogeneity of interest, we may consider
$|\partial_{i}|\sim k$. Since the only natural scale of the system other than
$k$ is $H$, the fictitious parameter can be regarded as $\epsilon\sim k/(aH)$.
We also have to make key assumptions for the leading order
$\epsilon\rightarrow 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 $\alpha-1=\mathcal{O}(\epsilon)$,
$\beta^{i}=\mathcal{O}(\epsilon)$ and
$\partial_{t}\tilde{\gamma}_{ij}=\mathcal{O}(\epsilon)$. In addition, we
assume that $\psi$ 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 $\partial_{t}\tilde{\gamma}_{ij}$
temporally decays and one can set
$\partial_{t}\tilde{\gamma}_{ij}=\mathcal{O}(\epsilon^{2})$ for growing mode
solutions [69].

In this review, we mainly focus on two specific time slicing conditions: the
constant mean curvature (CMC) slicing ($K=-3H$) and comoving slicing
($n^{\mu}=u^{\mu}$). Readers may refer to Ref. [34] for other gauge
conditions. The comoving slicing implies $u_{i}=u^{t}(v_{i}+\beta_{i})=0$,
then $\Gamma=1$, $E=\rho$ and $J_{i}=(\rho+p)u_{i}=0$. Since
$\tilde{A}_{ij}=\mathcal{O}(\epsilon^{2})$ from Eq. (3.6), from the momentum
(3.5) and Hamiltonian (3.4) constraints, we obtain
$\mathcal{\tilde{D}}_{i}K=\mathcal{O}(\epsilon^{3})$ and
$H^{2}=\frac{8\pi}{3}\rho+\mathcal{O}(\epsilon^{2})$. Therefore CMC, uniform
density and comoving slicings coincide to $\mathcal{O}(\epsilon)$. In
addition, from Eq. (3.1), we obtain

$2\partial_{t}\ln\psi=H-\partial_{t}\ln a+\mathcal{O}(\epsilon^{2}),$ (3.16)

where we have assumed that the pressure is homogeneous to
$\mathcal{O}(\epsilon)$. Since the right-hand side is a function of $t$,
$\partial_{t}\ln\psi$ is also homogeneous to $\mathcal{O}(\epsilon)$. Since we
have assumed that $\psi$ is identically unity at some point, we conclude

$\partial_{t}\psi=\mathcal{O}(\epsilon^{2}).$ (3.17)

Coming back to the general expressions without specifying the slicing
condition, from $\partial_{t}\psi=\mathcal{O}(\epsilon^{2})$ and the equations
of motion displayed in the previous section, one can find

$\displaystyle\psi=\mathcal{O}(\epsilon^{0}),~{}v^{i}=\mathcal{O}(\epsilon),~{}\rho=\rho_{\rm
b}(1+\mathcal{O}(\epsilon^{2})),~{}\tilde{A}_{ij}=\mathcal{O}(\epsilon^{2}),$
(3.18)
$\displaystyle\tilde{\gamma}_{ij}=f_{ij}+\mathcal{O}(\epsilon^{2}),~{}\alpha=1+\mathcal{O}(\epsilon^{2}),~{}K=3H(1+\mathcal{O}(\epsilon^{2})),$
(3.19)

where $\rho_{\rm b}$ is the background density. Then we introduce the
following perturbation variables:

$\displaystyle\psi(t,\bm{x})$ $\displaystyle=$
$\displaystyle\Psi(\bm{x})(1+\xi(t,\bm{x}))+\mathcal{O}(\epsilon^{3}),$ (3.20)
$\displaystyle\tilde{\gamma}_{ij}(t,\bm{x})$ $\displaystyle=$ $\displaystyle
f_{ij}+h_{ij}(t,\bm{x})+\mathcal{O}(\epsilon^{3}),$ (3.21) $\displaystyle
K(t,\bm{x})$ $\displaystyle=$ $\displaystyle
3H(1+\kappa(t,\bm{x}))+\mathcal{O}(\epsilon^{3}),$ (3.22)
$\displaystyle\alpha(t,\bm{x})$ $\displaystyle=$ $\displaystyle
1+\chi(t,\bm{x})+\mathcal{O}(\epsilon^{3}),$ (3.23)
$\displaystyle\beta^{i}(t,\bm{x})$ $\displaystyle=$
$\displaystyle\mathcal{O}(\epsilon),$ (3.24) $\displaystyle\rho(t,\bm{x})$
$\displaystyle=$ $\displaystyle\rho_{\rm
b}(1+\delta(t,\bm{x}))+\mathcal{O}(\epsilon^{3}),$ (3.25) $\displaystyle
v^{i}(t,\bm{x})$ $\displaystyle=$ $\displaystyle\mathcal{O}(\epsilon),$ (3.26)

where $\Psi(\bm{x})=\mathcal{O}(\epsilon^{0})$ and $\xi$, $h_{ij}$, $\kappa$,
$\chi$ and $\delta$ are the perturbation variables of
$\mathcal{O}(\epsilon^{2})$.

Hereafter, we focus on the case of the linear equation of state $p=w\rho$ and
consider the normal threading $\beta^{i}=0$ for simplicity. Readers may refer
to Ref. [34] for other general cases. The equations of motion can be reduced
to the followings:

$\displaystyle\partial_{t}\delta+3H(1+w)(\chi+\kappa)=O(\epsilon^{4}),$ (3.27)
$\displaystyle\partial_{t}[a^{3}(1+w)\rho_{\rm b}u_{j}]=-a^{3}\rho_{\rm
b}\left[w\partial_{j}\delta+(1+w)\partial_{j}\chi\right]+O(\epsilon^{5}).,$
(3.28) $\displaystyle 6\partial_{t}{\xi}-3H(\chi+\kappa)=O(\epsilon^{4}),$
(3.29) $\displaystyle\triangle\Psi=-2\pi\Psi^{5}a^{2}\rho_{\rm
b}(\delta-2\kappa)+O(\epsilon^{4}),$ (3.30) $\displaystyle
H^{-1}\partial_{t}{\kappa}=\frac{1}{2}(3w-1)\kappa-\frac{3}{2}(1+w)\chi-\frac{1}{2}\delta(1+3w)+O(\epsilon^{4}),$
(3.31) $\displaystyle\partial_{t}h_{ij}=-2\tilde{A}_{ij}+O(\epsilon^{4}),$
(3.32)
$\displaystyle\partial_{t}\tilde{A}_{ij}+3H\tilde{A}_{ij}=\frac{1}{a^{2}\Psi^{4}}\Bigl{[}-\frac{2}{\Psi}\left({\cal
D}_{i}{\cal D}_{j}\Psi-\frac{1}{3}f_{ij}{\triangle}\Psi\right)$ (3.33)
$\displaystyle\hskip 85.35826pt+\frac{6}{\Psi^{2}}\left({\cal D}_{i}\Psi{\cal
D}_{j}\Psi-\frac{1}{3}f_{ij}{\cal D}^{k}\Psi{\cal
D}_{k}\Psi\right)\Bigr{]}+O(\epsilon^{4}),$ (3.34) $\displaystyle{\cal
D}_{i}(\Psi^{6}\tilde{A}^{i}_{j})+2H\Psi^{6}{\cal
D}_{j}\kappa=8\pi\Psi^{6}(1+w)\rho_{\rm b}u_{j}+O(\epsilon^{5}),$ (3.35)

The solutions for Eqs. (3.32) and (3.34) are given by

$\displaystyle h_{ij}$ $\displaystyle=$
$\displaystyle-\frac{4}{(3w+5)(3w-1)}p_{ij}\left(\frac{1}{aH}\right)^{2}+O(\epsilon^{4}),$
(3.36) $\displaystyle\tilde{A}_{ij}$ $\displaystyle=$
$\displaystyle\frac{2}{3w+5}p_{ij}H\left(\frac{1}{aH}\right)^{2}+O(\epsilon^{4}),$
(3.37)

respectively, where

$p_{ij}(\bm{x}):=\frac{1}{\Psi^{4}}\left[-\frac{2}{\Psi}\left({\cal
D}_{i}{\cal
D}_{j}\Psi-\frac{1}{3}f_{ij}\triangle\Psi\right)+\frac{6}{\Psi^{2}}\left({\cal
D}_{i}\Psi{\cal D}_{j}\Psi-\frac{1}{3}f_{ij}{\cal D}^{k}\Psi{\cal
D}_{k}\Psi\right)\right]$ (3.38)

with $\triangle:=f^{ij}\mathcal{D}_{i}\mathcal{D}_{j}$.

To find solutions for other quantities, we need to fix the slicing condition.
Let us consider CMC slicing ($\kappa=0$). From Eq. (3.30), we obtain

$\delta=q\left(\frac{1}{aH}\right)^{2}+\mathcal{O}(\epsilon^{4}),$ (3.39)

where

$q(\bm{x})=-\frac{4}{3}\frac{\triangle\Psi}{\Psi}.$ (3.40)

Then, from Eqs. (3.31), (3.29) and (3.28), $\chi$, $\xi$ and $u_{j}$ are given
as

$\displaystyle\chi$ $\displaystyle=$
$\displaystyle-\frac{3w+1}{3(1+w)}q\left(\frac{1}{aH}\right)^{2}+\mathcal{O}(\epsilon^{4}),$
(3.41) $\displaystyle\xi$ $\displaystyle=$
$\displaystyle-\frac{1}{6(1+w)}q\left(\frac{1}{aH}\right)^{2}+\mathcal{O}(\epsilon^{4}),$
(3.42) $\displaystyle u_{j}$ $\displaystyle=$ $\displaystyle
v_{j}+\mathcal{O}(\epsilon^{5})=\frac{2}{3(1+w)(3w+5)}\partial_{j}qa\left(\frac{1}{aH}\right)^{3}+\mathcal{O}(\epsilon^{5}).$
(3.43)

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.28), we obtain

$\chi=\frac{w}{1+w}\delta+\mathcal{O}(\epsilon^{4}).$ (3.44)

Substituting this relation into Eq. (3.27) and (3.31), and combining them, we
can derive the following solutions for $\delta$ and $\kappa$:

$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\frac{3(1+w)}{3w+5}q\left(\frac{1}{aH}\right)^{2}+\mathcal{O}(\epsilon^{4}),$
(3.45) $\displaystyle\kappa$ $\displaystyle=$
$\displaystyle-\frac{1}{3w+5}q\left(\frac{1}{aH}\right)^{2}+\mathcal{O}(\epsilon^{4}).$
(3.46)

Then $\chi$ can be given as

$\chi=-\frac{3}{3w+5}q\left(\frac{1}{aH}\right)^{2}+\mathcal{O}(\epsilon^{4}).$
(3.47)

From Eq. (3.29), we obtain

$\xi=-\frac{1}{2(3w+5)}q\left(\frac{1}{aH}\right)^{2}+\mathcal{O}(\epsilon^{4}).$
(3.48)

Since the function $\Psi(\bm{x})$, and therefore $q(\bm{x})$ are shared in the
both gauge conditions, we can easily find the relation of gauge-dependent
quantities. For instance, we find

$\delta_{\rm co}=\frac{3(1+w)}{3w+5}\delta_{\rm
CMC}+\mathcal{O}(\epsilon^{4}),$ (3.49)

where $\delta_{\rm CMC}$ and $\delta_{\rm co}$ are density perturbations for
CMC and comoving gauges, respectively.

##  4 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 $\bar{\delta}$ be the amplitude of the volume average of
the density perturbation at the horizon entry, if $\bar{\delta}$ exceeds a
certain threshold value, we think that PBH formation is realized. The
threshold value has been estimated analytically and numerically by many
authors [66, 70, 71, 28, 27, 30, 68, 72, 73, 31, 32, 33].

A recent trend is to use the “compaction function” firstly 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 [36, 37, 38, 35]. 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 [74]. 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:

$\mathcal{C}_{\rm SS}(r):=\frac{\delta M_{\rm SS}}{R}~{}~{}{\rm
with}~{}~{}\delta M_{\rm SS}:=4\pi a^{3}\rho_{\rm
b}\int^{r}_{0}\mathrm{d}x\delta_{\rm
CMC}(\Psi^{2}x)^{2}\frac{\mathrm{d}}{\mathrm{d}x}(\Psi^{2}x),$ (4.1)

where $R$ is the areal radius of the sphere of the coordinate radius $r$ given
by $R=ar\Psi^{2}$. By using the solution (3.39), one can derive the following
expression [34]:

$\mathcal{C}_{\rm
SS}(r)=\frac{1}{2}\left[1-\left(1+2r\frac{\mathrm{d}}{\mathrm{d}r}\ln\Psi\right)\right].$
(4.2)

From this expression, we find that $\mathcal{C}_{\rm SS}$ is time independent
and calculated by the gauge independent quantity $\Psi$. Furthermore, $\Psi$
is related to the curvature perturbation, which is often used in the context
of cosmological perturbation in the early universe, as
$\zeta=-\frac{1}{2}\ln\Psi$ 2)2)2)We choose the sign convention so that the
spatial metric will be given by $a^{2}{\rm
e}^{-2\zeta}\tilde{\gamma}_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j}$. . Let us
introduce the radius $r_{\rm m}$ at which the compaction function takes the
maximum value $\mathcal{C}_{\rm SS}^{\rm max}=\mathcal{C}_{\rm SS}(r_{\rm
m})$. Then, we adopt the following PBH formation criterion:

$\mathcal{C}^{\rm max}_{\rm SS}\gtrsim\mathcal{C}_{\rm SS}^{\rm th}$ (4.3)

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 $\mathcal{C}_{\rm th}$ can be roughly estimated as
$0.4$ with about 10% uncertainty. For the specific profile given by
$\ln\Psi\propto\sin(r)/r$, which corresponds to the typical profile for the
monochromatic power spectrum of $\ln\Psi$, we obtain $\mathcal{C}_{\rm
th}\simeq 0.44$ from numerical simulations. Therefore we use the value
$\mathcal{C}_{\rm 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 [67].
As is shown in 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.20), we can derive the following expression:

$\displaystyle M(t,r)$ $\displaystyle=$ $\displaystyle
4\pi\int^{r}_{0}\mathrm{d}xx^{2}a^{3}\alpha\psi^{6}T^{t}_{~{}\mu}K^{\mu}$
(4.4) $\displaystyle=$ $\displaystyle 4\pi
a^{3}\int^{r}_{0}\mathrm{d}xx^{2}\psi^{4}\left\\{-[(\rho+p)u^{t}u_{t}+p]\partial_{x}(\psi^{2}x)+(\rho+p)u^{t}u_{r}\frac{x}{a}\partial_{t}(\psi^{2}a)\right\\}.$
(4.5)

In the comoving slice ($u_{r}=0$), we obtain

$\displaystyle M_{\rm co}(t,r)$ $\displaystyle=$ $\displaystyle 4\pi
a^{3}\int^{r}_{0}\mathrm{d}xx^{2}\psi^{4}\rho\left\\{\partial_{x}(\psi^{2}x)+\mathcal{O}(\epsilon^{4})\right\\}.$
(4.6)

The mass excess compared to the background FLRW with the identical areal
radius is calculated as

$\delta M_{\rm co}:=M(t,r)-M_{\rm FLRW}(t,\psi^{2}r)=4\pi a^{3}\rho_{\rm
b}\int^{r}_{0}\mathrm{d}x\left[(\Psi^{2}x)^{2}\partial_{x}(\Psi^{2}x)\delta_{\rm
co}+\mathcal{O}(\epsilon^{4})\right].$ (4.7)

Then we can define the compaction function in the comoving gauge as

$\mathcal{C}_{\rm co}:=\frac{\delta M_{\rm
co}}{R}=\frac{3(1+w)}{3w+5}\frac{1}{2}\left[1-(1+2r\partial_{r}\ln\Psi)\right]=\frac{3(1+w)}{3w+5}\mathcal{C}_{\rm
SS},$ (4.8)

where $R=a\Psi^{2}r$. Therefore $\mathcal{C}_{\rm SS}$ is the compaction
function in the comoving gauge multiplied by $\frac{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.5) cannot be neglected3)3)3)In Refs. [28, 34],
the second term in the integrand of Eq. (4.5) has been erroneously dropped,
and the authors got the wrong relation $\mathcal{C}_{\rm CMC}=\mathcal{C}_{\rm
SS}$. The original definition of $\mathcal{C}_{\rm 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.  and we find
$\mathcal{C}_{\rm CMC}\neq\mathcal{C}_{\rm SS}$. For a more concrete
description, readers may refer to Ref. [74].

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:

$\bar{\delta}:=\frac{4\pi\int^{r}_{0}\mathrm{d}x(\Psi^{2}x)^{2}\partial_{x}(\Psi^{2}x)\delta_{\rm
co}}{4\pi(\Psi^{2}r)^{3}/3}.$ (4.9)

Then we can easily find

$\bar{\delta}=\frac{2}{H^{2}R^{2}}\frac{\delta M_{\rm co}}{R}.$ (4.10)

Therefore we obtain

$\left.\bar{\delta}\right|_{HR=1}=2\mathcal{C}_{\rm co}.$ (4.11)

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.

##  5 Estimation of PBH abundance in peak theory

We start with writing the spatial metric in the following form:

$\mathrm{d}s_{3}^{2}=a^{2}{\rm
e}^{-2\zeta}\tilde{\gamma}_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j}.$ (5.1)

The curvature perturbation $\zeta$ is related to $\Psi$ as

$\zeta=-\frac{1}{2}\ln\Psi.$ (5.2)

Here we assume that $\zeta$ is the Gaussian random variable with the power
spectrum $\mathcal{P}(k)$ defined by

$\langle\tilde{\zeta}^{*}(\bm{k})\tilde{\zeta}(\bm{k}^{\prime})\rangle=\frac{2\pi^{2}}{k^{3}}\mathcal{P}(k)(2\pi)^{3}\delta(\bm{k}-\bm{k}^{\prime})$
(5.3)

with $\tilde{\zeta}(\bm{k})=\int\mathrm{d}^{3}x\zeta(\bm{x}){\rm
e}^{\mathrm{i}\bm{k}\bm{x}}$. Readers may refer to Refs. [75, 76, 77] for non-
Gaussian cases. From the random Gaussian assumption and the gradient moments
$\sigma_{n}$ given by

$\sigma_{n}^{2}=\int\mathrm{d}\ln kk^{2n}\mathcal{P}(k),$ (5.4)

as is shown in Appendices B, C and D, we can obtain the peak number density
$n_{\rm pk}^{(k_{\bullet})}(\tilde{\mu}_{2},k_{\bullet})$ and the typical peak
profile $\bar{\zeta}(r)$ characterized by the two parameters $\tilde{\mu}_{2}$
and $k_{\bullet}$ defined by

$\displaystyle\tilde{\mu}_{2}$ $\displaystyle=$
$\displaystyle\triangle\zeta(0)\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}},$ (5.5)
$\displaystyle k_{\bullet}^{2}$ $\displaystyle=$
$\displaystyle-\frac{\triangle\triangle\zeta(0)}{\tilde{\mu}_{2}}\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}},$
(5.6)

where $\triangle$ is the flat Laplacian. Here we just write down the final
expressions (D.3) and (D.5):

$\displaystyle n^{(k_{\bullet})}_{\rm
pk}(\tilde{\mu}_{2},k_{\bullet})d\tilde{\mu}_{2}dk_{\bullet}=\frac{2\cdot
3^{3/2}}{(2\pi)^{3/2}}\tilde{\mu}_{2}k_{\bullet}\frac{\sigma_{2}^{3}\sigma_{4}^{2}}{\sigma_{1}^{4}\sigma_{3}^{3}}f\left(\tilde{\mu}_{2}k_{\bullet}^{2}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}\sigma_{4}}\right)P_{1}\left(\tilde{\mu}_{2}\frac{\sigma_{2}}{\sigma_{1}^{2}},\mu_{2}k_{\bullet}^{2}\frac{\tilde{\sigma}_{2}^{2}}{\sigma_{1}^{2}\sigma_{4}};\gamma_{3}\right)\mathrm{d}\tilde{\mu}_{2}\mathrm{d}k_{\bullet},$
(5.7)
$\displaystyle\frac{\overline{\zeta}(r)}{\tilde{\mu}_{2}}=-\frac{1}{1-\gamma_{3}^{2}}\left(\psi_{1}+\frac{1}{3}R_{3}^{2}\triangle\psi_{1}\right)+k_{\bullet}^{2}\frac{1}{\gamma_{3}(1-\gamma_{3}^{2})}\frac{\sigma_{2}}{\sigma_{4}}\left(\gamma_{3}^{2}\psi_{1}+\frac{1}{3}R_{3}^{2}\triangle\psi_{1}\right),~{}$
(5.8)

where $P_{1}$, $f$, $R_{n}$ and $\psi_{n}$ are defined in Eqs. (B.19), (B.35),
(C.23) 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.

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

###  5.1 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,
$\mathcal{C}=\mathcal{C}_{\rm co}$. Focusing on the radiation-dominated
universe ($w=1/3$), the PBH formation criterion can be rewritten as

$\mathcal{C}^{\rm max}=\mathcal{C}(r_{\rm m})=\frac{2}{3}\mathcal{C}^{\rm
max}_{\rm SS}\gtrsim\mathcal{C}^{\rm th}\sim 0.293$ (5.9)

In terms of the typical profile $\bar{\zeta}$, the compaction function can be
written as

$\mathcal{C}=\frac{1}{3}\left[1-\left(1-r\partial_{r}\bar{\zeta}\right)^{2}\right].$
(5.10)

Assuming type I fluctuations defined by
$\partial_{r}R>0\Rightarrow\partial_{r}({\rm e}^{-\zeta}r)>0$, we obtain

$\left(\partial_{r}\bar{\zeta}+r\partial_{r}^{2}\zeta\right)|_{r=r_{\rm
m}}=0.$ (5.11)

Defining the function $g(r;k_{\bullet})$ by

$g(r;k_{\bullet}):=\frac{\bar{\zeta}(r)}{\tilde{\mu}_{2}},$ (5.12)

the condition (5.11) for $r_{\rm m}$ can be rewritten as

$(\partial_{r}g+r\partial_{r}^{2}g)|_{r=r_{\rm m}}=0.$ (5.13)

This equation implicitly gives the value of $r_{\rm m}$ as a function of
$k_{\bullet}$: $r_{\rm m}=r_{\rm m}(k_{\bullet})$. The value $\mathcal{C}^{\rm
max}$ is then given as a function of $\tilde{\mu}_{2}$ and $k_{\bullet}$
through $r_{\rm m}(k_{\bullet})$ as

$\tilde{\mu}_{2}=\frac{1-\sqrt{1-3\mathcal{C}_{\rm max}}}{r_{\rm
m}\partial_{r}g(r_{\rm m})}.$ (5.14)

Through this equation, we can obtain the threshold value of $\tilde{\mu}_{2}$
as a function of $k_{\bullet}$ as

$\tilde{\mu}_{\rm
2th}^{(k_{\bullet})}(k_{\bullet}):=\frac{1-\sqrt{\mathcal{C}_{\rm th}}}{r_{\rm
m}(k_{\bullet})\partial_{r}g_{\rm m}(k_{\bullet})},$ (5.15)

where $\partial_{r}g_{\rm m}=\partial_{r}g(r_{\rm
m}(k_{\bullet});k_{\bullet})$.

###  5.2 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 $\tilde{\mu}_{2}$, $k_{\bullet}$ 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

$aH=\frac{a}{R}=\frac{{\rm e}^{\tilde{\mu}_{2}g_{\rm m}}}{r_{\rm m}}.$ (5.16)

In the radiation-dominated universe, adopting the matter-radiation equality as
a reference time, we find

$a=a_{\rm eq}\left(\frac{H_{\rm eq}}{H}\right)^{1/2}.$ (5.17)

From Eqs. (5.16) and (5.17), at the horizon entry, we obtain

$H^{-1}=a_{\rm eq}^{2}H_{\rm eq}r_{\rm m}^{2}{\rm e}^{-2\tilde{\mu}_{2}g_{\rm
m}}=H_{\rm eq}^{-1}k_{\rm eq}^{2}r_{\rm m}^{2}{\rm e}^{-2\tilde{\mu}_{2}g_{\rm
m}},$ (5.18)

where $k_{\rm eq}=a_{\rm eq}H_{\rm eq}$. Then the PBH mass can be estimated as

$M=\frac{1}{2H}=\frac{1}{2}H_{\rm eq}^{-1}k_{\rm eq}^{2}r_{\rm m}^{2}{\rm
e}^{-2\tilde{\mu}_{2}g_{\rm m}}=M_{\rm eq}k_{\rm eq}^{2}r_{\rm m}^{2}{\rm
e}^{-2\tilde{\mu}_{2}g_{\rm m}},$ (5.19)

where we have introduced $M_{\rm eq}=1/(2H_{\rm eq})$. We have not taken into
account the mass scaling associated with the critical behavior [78, 79] 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 [65].

###  5.3 Step 3: Derivation of the mass function

Eq. (5.24) gives us the relation between $M$, $k_{\bullet}$ and
$\tilde{\mu}_{2}$. Therefore we can change the independent variables in the
peak number density (5.7) from $\tilde{\mu}_{2}$ and $k_{\bullet}$ to
$\tilde{\mu}_{2}$ and $M$ as follows:

$\displaystyle n^{(M)}_{\rm pk}(\tilde{\mu}_{2},M)d\tilde{\mu}_{2}d\ln M$
$\displaystyle=$ $\displaystyle\frac{2\cdot
3^{3/2}}{(2\pi)^{3/2}}\tilde{\mu}_{2}k_{\bullet}\frac{\sigma_{2}^{3}\sigma_{4}^{2}}{\sigma_{1}^{4}\sigma_{3}^{3}}f\left(\tilde{\mu}_{2}k_{\bullet}^{2}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}\sigma_{4}}\right)$
(5.21) $\displaystyle
P_{1}\left(\tilde{\mu}_{2}\frac{\sigma_{2}}{\sigma_{1}^{2}},\tilde{\mu}_{2}k_{\bullet}^{2}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}\sigma_{4}};\gamma_{3}\right)\left|\frac{\mathrm{d}\ln
M}{\mathrm{d}k_{\bullet}}\right|^{-1}\mathrm{d}\tilde{\mu}_{2}\mathrm{d}\ln
M,$

where $k_{\bullet}$ should be regarded as a function of $\tilde{\mu}_{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.21) for $\tilde{\mu}_{2}$.

As is explicitly shown in Ref. [63] for some specific cases, the threshold
value of $\tilde{\mu}_{2}$ given by Eq. (5.15) can be converted to the
threshold value $\tilde{\mu}_{\rm 2th}^{(M)}(M)$ for a given value of $M$ by
using Eqs. (5.15) and (5.19). In addition, the value of $\tilde{\mu}_{2}$ may
be bounded below for a fixed value of $M$. Then, letting the minimum possible
value of $\tilde{\mu}_{2}$ be $\tilde{\mu}_{\rm 2min}(M)$, the lower limit of
the integral for $\tilde{\mu}_{2}$ is given by $\tilde{\mu}_{\rm
2b}:=\max\left\\{\tilde{\mu}_{\rm 2min}(M),\tilde{\mu}_{\rm
2th}^{(M)}(M)\right\\}$. That is, the PBH number density $n_{\rm BH}(M)$ and
the mass spectrum $f_{0}(M)$ at equality time are given by

$\displaystyle n_{\rm BH}\mathrm{d}\ln M$ $\displaystyle=$
$\displaystyle\left[\int^{\infty}_{\tilde{\mu}_{\rm
2b}}\mathrm{d}\tilde{\mu}_{2}n^{(M)}_{\rm
pk}(\tilde{\mu}_{2},M)\right]\mathrm{d}\ln M,$ (5.22) $\displaystyle
f_{0}\mathrm{d}\ln M$ $\displaystyle=$ $\displaystyle\frac{Mn_{\rm BH}}{\rho
a^{3}}\frac{a_{\rm eq}}{a}\mathrm{d}\ln M=\frac{4\pi}{3}\frac{M}{M_{\rm
eq}}k_{\rm eq}^{-3}n_{\rm BH}\mathrm{d}\ln M.$ (5.23)

###  5.4 Critical behavior

It is well known that, at least for idealized spherically-symmetric
situations, the mass of the black hole scales as $\propto(b-b_{\rm
th})^{\gamma}$ with $b$ and $b_{\rm 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 [27, 30, 72, 73, 26, 80, 81, 82, 83, 84, 76]. If we take the
critical behavior into account, the PBH mass would be given by

$M=K(k_{\bullet})(\tilde{\mu}_{2}-\tilde{\mu}_{\rm
2th}^{(k_{\bullet})}(k_{\bullet}))^{\gamma}M_{\rm eq}k_{\rm eq}^{2}r_{\rm
m}^{2}{\rm e}^{-2\tilde{\mu}_{2}g_{\rm m}},$ (5.24)

where $\gamma$ is the universal exponent given by $\gamma\simeq 0.36$ for
radiation fluid, and $K(k_{\bullet})$ is a profile-dependent 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 $\gamma$ to
be 0 in Eq. (5.24).

###  5.5 The monochromatic curvature power spectrum

In this paper, we show only the case of the monochromatic curvature power
spectrum. Readers may refer to Refs. [62, 63] for extended spectra. Let us
write the monochromatic power spectrum as

$\mathcal{P}(k)=\sigma_{0}^{2}k_{0}\delta(k-k0)=\sigma_{0}\delta(\ln k-\ln
k_{0}).$ (5.25)

Then the gradient moments are calculated as

$\sigma_{n}^{2}=\sigma_{0}^{2}k_{0}^{2n}\Rightarrow\gamma_{n}=1.$ (5.26)

The typical profile is given as

$g(r;k_{\bullet})=\frac{\bar{\zeta}}{\tilde{\mu}_{2}}=-\psi_{1}=-\frac{\sin(k_{0}r)}{k_{0}r}.$
(5.27)

From the condition (5.11), we can find

$r_{\rm m}\simeq 2.74/k_{0},~{}g_{\rm m}=g(r_{\rm m})\simeq-0.141.$ (5.28)

From the condition (5.9), we can find

$\tilde{\mu}_{\rm 2th}\simeq 0.615.$ (5.29)

Then the PBH mass can be given by a function of $\tilde{\mu}_{2}$ as

$M=(\tilde{\mu}_{2}-\tilde{\mu}_{\rm 2th})^{\gamma}M_{\rm eq}k_{\rm
eq}^{2}r_{\rm m}^{2}{\rm e}^{-2\tilde{\mu}_{2}g_{\rm m}},$ (5.30)

where $\gamma=0.36$ and $0$ with and without taking into account the critical
behavior, respectively. The minimum value of the PBH mass $M_{\rm min}$ is
given by $0$ and $M_{\rm eq}k_{\rm eq}^{2}r_{\rm m}^{2}{\rm
e}^{-2\tilde{\mu}_{\rm 2th}g_{\rm m}}$ for $\gamma=0.36$ and $0$,
respectively.

For the number density, we reconsider the expression (5.7). In the limit
$\gamma_{3}\rightarrow 0$, the probability distribution
$P_{1}(\nu,\xi_{1};\gamma_{3})$ can be written as

$\lim_{\gamma_{3}\rightarrow
1}P_{1}(\nu,\xi_{1};\gamma_{3})=\lim_{\gamma_{3}\rightarrow
1}\frac{2}{2\pi}\frac{1}{\sqrt{1-\gamma_{3}}}\left[-\frac{1}{2}\left(\nu^{2}+\frac{(\xi_{1}-\gamma_{3}\nu)^{2}}{1-\gamma_{3}^{2}}\right)\right]=\delta(\xi_{1}-\nu)\frac{1}{\sqrt{2\pi}}\exp\left[-\frac{1}{2}\nu^{2}\right].$
(5.31)

Substituting $\nu=\mu_{2}/\sigma_{2}$ and
$\xi_{1}=\mu_{2}k_{\bullet}^{2}/\sigma_{4}$ and using
$\sigma_{n}=\sigma_{0}k_{0}$, we obtain

$P_{1}(\frac{\tilde{\mu}_{2}}{\sigma_{0}},\frac{\tilde{\mu}_{2}k_{\bullet}^{2}}{\sigma_{0}k_{0}^{2}};1)=\frac{\sigma_{0}k_{0}}{2\tilde{\mu}_{2}}\delta(k_{\bullet}-k_{0})\frac{1}{\sqrt{2\pi}}\exp\left[-\frac{1}{2}\frac{\tilde{\mu}_{2}^{2}}{\sigma_{0}^{2}}\right].$
(5.32)

Integrating the number density for $k_{\bullet}$, we obtain

$n_{\rm
pk}^{(\tilde{\mu}_{2})}\mathrm{d}\tilde{\mu}_{2}:=\int\mathrm{d}k_{\bullet}n_{\rm
pk}^{(k_{\bullet})}(\tilde{\mu}_{2},k_{\bullet})\mathrm{d}\tilde{\mu}_{2}=\frac{3^{3/2}}{(2\pi)^{2}}\frac{1}{\sigma_{0}}k_{0}^{3}f\left(\frac{\tilde{\mu}_{2}}{\sigma_{0}}\right)\exp\left(-\frac{\tilde{\mu}_{2}^{2}}{2\sigma_{0}^{2}}\right).$
(5.33)

The variable $\tilde{\mu}_{2}$ can be converted into $M$ through Eq. (5.30).
Then we obtain

$\displaystyle n_{\rm BH}(M)d\ln M$ $\displaystyle=$ $\displaystyle
n^{(\tilde{\mu}_{2})}_{\rm pk}\left|\frac{\mathrm{d}\ln
M}{\mathrm{d}\tilde{\mu}_{2}}\right|^{-1}\Theta(M-M_{\rm min})\mathrm{d}\ln M$
(5.34) $\displaystyle=$
$\displaystyle\frac{3^{3/2}}{(2\pi)^{2}}\frac{1}{\sigma_{0}}k_{0}^{3}f\left(\frac{\tilde{\mu}_{2}}{\sigma_{0}}\right)\exp\left(-\frac{\tilde{\mu}_{2}^{2}}{2\sigma_{0}^{2}}\right)\left|\frac{\mathrm{d}\ln
M}{\mathrm{d}\tilde{\mu}_{2}}\right|^{-1}\Theta(M-M_{\rm min})\mathrm{d}\ln
M,~{}~{}$ (5.35)

where $\tilde{\mu}_{2}$ is regarded as a function of $M$ through (5.30), and

$\left|\frac{\mathrm{d}\ln
M}{\mathrm{d}\tilde{\mu}_{2}}\right|=\left|\frac{\gamma}{\tilde{\mu}_{2}-\tilde{\mu}_{\rm
2th}}-2g_{\rm
m}\right|\simeq\left|\frac{\gamma}{\tilde{\mu}_{2}-0.615}-0.282\right|.$
(5.36)

The PBH fraction $f_{0}$ is given by

$f_{0}\simeq\frac{\sqrt{3}}{\pi}\frac{1}{\sigma_{0}}\frac{M}{M_{\rm
eq}}\frac{k_{0}^{3}}{k_{\rm
eq}^{3}}f\left(\frac{\tilde{\mu}_{2}}{\sigma_{0}}\right)\exp\left(-\frac{\tilde{\mu}_{2}^{2}}{\sigma_{0}^{2}}\right)\left|\frac{\gamma}{\tilde{\mu}_{2}-0.615}-0.282\right|^{-1}\Theta\left(M-M_{\rm
min}\right).$ (5.37)

Since small values of the PBH mass can be realized for
$\tilde{\mu}_{2}\sim\tilde{\mu}_{\rm 2th}$, the behavior of $f_{0}$ in the
small PBH mass limit is given by $f_{0}\propto M^{1+\gamma^{-1}}$. We show the
PBH fraction $f_{0}(M)$ for $k_{0}/k_{\rm eq}=10^{5}$ as a function of
$M/M_{0}$ in Fig. 5 with $M_{0}=M_{\rm eq}k_{\rm eq}^{2}/k_{0}^{2}$.

Figure 5: The PBH fractions $f_{0}$ are plotted as functions of $M/M_{0}$ for
$\sigma_{0}=0.08$, $0.10$ and $0.12$ with $k_{0}=10^{5}k_{\rm eq}$. The right
spiky plots show the PBH fractions without the critical behavior, and the left
broader plots show those with the critical behavior.

The total fraction $f_{\rm tot}=\int f_{0}\mathrm{d}\ln M$ is shown as
functions of $\sigma_{0}$ for the cases with and without critical behavior in
Fig. 6

Figure 6: The total fraction $f_{\rm tot}$ as a function of $\sigma_{0}$.

##  6 Summary and further related topics

In this review, we have run through the basic picture of PBH formation,
cosmological long-wavelength 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 three-zone 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 which 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. [68]) 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\rho$, the cosmological
long-wavelength solution should be reconsidered. For example, in Ref. [85],
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 [74]. More
concrete discussions will be reported elsewhere [74]. 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. [62, 63] mainly
focusing on the case of the monochromatic curvature power spectrum. However,
different procedures based on peak theory have been also proposed [83, 86, 87,
88]. Here we note that, while the variation of the peak curvature scale is
treated as an independent statistical variable in Refs. [62, 86, 63], the
typical profile for the mean peak curvature is only considered in Refs. [83,
87]. 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.

There are still many interesting open issues we could not mention here such as
effects of non-Gaussianity[89, 90, 80, 91, 92, 93, 94, 95, 96, 97, 98, 96, 99,
100, 75, 101, 102, 103, 76, 84, 77] PBH in matter dominated era[104, 105, 106,
107, 108], non-spherical PBH formation [109], distribution of PBH spin[110,
111, 107, 112, 113, 114, 115], PBH formation in the QCD phase transition epoch
[116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 77, 41, 126] and PBH from
iso-curvature perturbations[127, 85], and so on. PBH formation in modified
gravity theories would require more careful treatments in all aspects. We hope
this review will be helpful for researchers to study those issues in the
future.

## Acknowledgment

This review is based on the lecture given at Asia-Pacific School and Workshop
on Gravitation and Cosmology 2022. This work was supported in part by JSPS
KAKENHI Grant Numbers JP19H01895 , JP20H05850 and JP20H05853.

##  Appendix A Horizon and mass

Throughout this section, we assume spherical symmetry for simplicity and
review the references [67, 128, 129, 130]. In the case of the standard PBH
formation scenario in RD, this assumption can be justified as will be
discussed in Appendix B.

First, let us write the line element of a spherically symmetric spacetime as
follows:

$\mathrm{d}s^{2}=-{\rm e}^{2\phi(t,r)}\mathrm{d}t^{2}+{\rm
e}^{\lambda}(t,r)\mathrm{d}r^{2}+R^{2}(t,r)\mathrm{d}\Omega^{2}.$ (A.1)

Introducing the null coordinates satisfying

$\sqrt{2}{\rm e}^{-f(t,r)}\mathrm{d}\xi_{\pm}={\rm e}^{\phi}\mathrm{d}t\pm{\rm
e}^{\lambda/2}\mathrm{d}r,$ (A.2)

we can rewrite the line element as

$\mathrm{d}s^{2}=-2{\rm
e}^{-f(t,r)}\mathrm{d}\xi_{+}\mathrm{d}\xi_{-}+R^{2}(t,r)\mathrm{d}\Omega^{2}.$
(A.3)

In this general spherically symmetric spacetime, existence of causal horizons,
such as cosmological and black hole horizons, can be characterized by
marginally trapped surfaces. Defining the null expansion $\theta_{\pm}$ as

$\theta_{\pm}:=\frac{2}{R}\frac{\partial R}{\partial\xi_{\pm}},$ (A.4)

we can define a marginally trapped surface as a surface satisfying

$\theta_{+}\theta_{-}=0.$ (A.5)

For instance, assuming $\partial_{t}$ and $\partial_{r}$ as future-directed
and outward vector fields, respectively, we find that $\theta_{+}=0$ for the
black hole horizon in the Schwarzschild black hole and $\theta_{-}=0$ for the
cosmological horizon in the spatially flat FLRW universe models.

The definition of the Misner-Sharp mass is originally given by [67]

$M(t,r)=\frac{1}{2}R\left[1+{\rm e}^{-2\phi}(\partial_{t}R)^{2}-{\rm
e}^{-\lambda}(\partial_{r}R)^{2}\right].$ (A.6)

As can be easily checked, we can also write the Misner-Sharp mass in the
following different forms:

$\displaystyle M(t,r)$ $\displaystyle=$
$\displaystyle\frac{1}{2}R\left(1-g^{\mu\nu}\partial_{\mu}R\partial_{\nu}R\right)$
(A.7) $\displaystyle=$ $\displaystyle\frac{1}{2}R+\frac{1}{4}{\rm
e}^{f}R^{3}\theta_{+}\theta_{-}.$ (A.8)

From this expression, we find

$R(t,r)=2M(t,r)~{}~{}~{}\text{if}~{}~{}~{}\theta_{+}\theta_{-}=0.$ (A.9)

That is, half of the horizon radius is equal to the Misner-Sharp mass enclosed
by the horizon in general spherically symmetric spacetimes.

It would not be easy to understand the expression (A.6) describes the mass of
the system. In order to see the relation between the stress-energy tensor and
the Misner-Sharp mass, let us write down the Einstein equations in terms of
the Misner-Sharp mass and the stress-energy tensor. We obtain two independent
equations as follows:

$\displaystyle\partial_{r}M$ $\displaystyle=$ $\displaystyle 4\pi
R^{2}\left(T^{t}_{~{}t}\partial_{r}R-T^{t}_{~{}r}\partial_{t}R\right),$ (A.10)
$\displaystyle\partial_{t}M$ $\displaystyle=$ $\displaystyle 4\pi R^{2}{\rm
e}^{-\lambda}\left(T^{r}_{~{}t}\partial_{r}R-T^{r}_{~{}r}\partial_{t}R\right).$
(A.11)

For instance, in the case of a perfect fluid, we can take the comoving
coordinate for $t$ and $r$, that is, $T^{t}_{~{}r}=0$. Then, from Eq. (A.10),
we obtain

$M=4\pi\int T^{t}_{~{}t}R^{2}\mathrm{d}R.$ (A.12)

This expression would be much more familiar and understandable than Eq. (A.6).
Nevertheless, the expression (A.12) contains the integral and looks non-local
quantity in contrast to the expression (A.6). Eq. (A.6) is often more useful
than the expression (A.12) because it is composed of only locally defined
geometrical quantities and independent of the matter components.

Another way to define the quasi-local mass in spherical symmetry associated
with a conserved current has been proposed by Kodama [128]. The mass defined
by Kodama is equivalent to the Misner-Sharp mass as will be shown soon. Let us
define the following vector, often called Kodama vector, as follows:

$K^{A}:=-\epsilon^{AB}\partial_{B}R,$ (A.13)

where $A$ and $B$ run over $t$ and $r$, and $\epsilon^{AB}$ is the
2-dimensional totally anti-symmetric tensor given by

$\epsilon_{AB}={\rm e}^{\phi+\lambda/2}\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right),~{}~{}\epsilon^{AB}={\rm
e}^{-\phi-\lambda/2}\left(\begin{array}[]{cc}0&-1\\\ 1&0\end{array}\right).$
(A.14)

The sign on the left-hand side is a convention that makes the Kodama vector
future-directed for $\partial_{r}R>0$. This vector in the 2-dimensional
spacetime spanned by $t$ and $r$ can be trivially promoted to the vector field
in the 4-dimensional spacetime as

$K^{\mu}=K^{t}\left(\frac{\partial}{\partial
t}\right)^{\mu}+K^{r}\left(\frac{\partial}{\partial r}\right)^{\mu}.$ (A.15)

Then the Kodama current $S^{\mu}$ is defined by

$S^{\mu}:=-T^{\mu}_{~{}\nu}K^{\nu},$ (A.16)

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:

$\displaystyle S^{A}=-T^{A}_{~{}B}K^{B}$ $\displaystyle=$ $\displaystyle{\rm
e}^{-\phi-\lambda/2}\frac{1}{4\pi
R^{2}}\left[\partial_{r}M(\partial_{r})^{A}-\partial_{t}M(\partial_{r})^{A}\right]$
(A.17) $\displaystyle=$ $\displaystyle-\frac{1}{4\pi
R^{2}}\epsilon^{AB}\partial_{B}M.$ (A.18)

Then we obtain

$\nabla_{\mu}S^{\mu}=\frac{1}{R^{2}{\rm
e}^{\phi+\lambda/2}}\partial_{A}(R^{2}{\rm
e}^{\phi+\lambda/2}S^{A})=-\frac{1}{4\pi R^{2}}\frac{\epsilon^{AB}}{{\rm
e}^{\phi+\lambda/2}}\partial_{A}\partial_{B}M=0,$ (A.19)

where we have used $\partial_{\mu}(\epsilon^{AB}/{\rm e}^{\phi+\lambda/2})=0$.
Therefore we can define the following conserved charge (Kodama mass):

$\displaystyle M_{\rm K}$ $\displaystyle:=$ $\displaystyle-\int
S^{\mu}n_{\mu}\mathrm{d}\Sigma$ (A.20) $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\int\frac{1}{R^{2}}\epsilon^{AB}\partial_{B}Mn_{A}\mathrm{d}\Sigma=\int\partial_{r}M\mathrm{d}r=M,$
(A.21)

where $\mathrm{d}\Sigma$ is the volume element of the $t=\text{const.}$
hypersurface whose unit normal form is given by $n_{\mu}=-{\rm
e}^{\phi}(\mathrm{d}t)_{\mu}$. This equation explicitly shows the equivalence
between the Kodama mass and the Misner-Sharp mass.

##  Appendix B Peak number density in peak theory

In this and the next section, we briefly review the necessary part of the peak
theory [64] for the calculation of PBH abundance. Let us consider the Gaussian
variable with the power spectrum $\mathcal{P}(k)$ defined by

$\langle\tilde{\zeta}^{*}(\bm{k})\tilde{\zeta}(\bm{k}^{\prime})\rangle=\frac{2\pi^{2}}{k^{3}}\mathcal{P}(k)(2\pi)^{3}\delta(\bm{k}-\bm{k}^{\prime})$
(B.1)

with $\tilde{\zeta}(\bm{k})=\int\mathrm{d}^{3}x\zeta(\bm{x}){\rm
e}^{\mathrm{i}\bm{k}\bm{x}}$. The gradient moments $\sigma_{n}$ are given by

$\sigma_{n}^{2}=\int\mathrm{d}\ln kk^{2n}\mathcal{P}(k).$ (B.2)

The probability distribution of linear combinations of $\zeta(\bm{x})$ denoted
by $V_{I}~{}(I=1,2,\cdots,n)$ is given by the multivariate Gaussian
distribution:

$P(V_{I})\mathrm{d}^{n}V=(2\pi)^{-n/2}\left|\det\mathcal{M}\right|^{-1/2}\exp\left[-\frac{1}{2}V_{I}(\mathcal{M}^{-1})^{IJ}V_{J}\right]\mathrm{d}^{n}V,$
(B.3)

where

$\mathcal{M}_{IJ}=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k^{\prime}}{(2\pi)^{3}}\langle\tilde{V}^{*}_{I}(\bm{k})\tilde{V}_{J}(\bm{k}^{\prime})\rangle.$
(B.4)

For example, let us consider the probability distribution of $-\zeta$ and
$\triangle\zeta$ at a spatial point, where $\triangle$ is the flat Laplacian.
The correlation matrix can be calculated as

$\mathcal{M}_{IJ}=\left(\begin{array}[]{cc}\sigma_{0}^{2}&\sigma_{1}^{2}\\\
\sigma_{1}^{2}&\sigma_{2}^{2}\end{array}\right).$ (B.5)

Then the probability distribution of $-\zeta$ and $\triangle\zeta$ is given by

$\displaystyle
P(-\zeta,\triangle\zeta)\mathrm{d}(-\zeta)\mathrm{d}(\triangle\zeta)=\frac{1}{2\pi\sqrt{\sigma_{0}^{2}\sigma_{2}^{2}-\sigma_{1}^{4}}}$
(B.6) $\displaystyle\hskip
28.45274pt\times\exp\left[-\frac{1}{2(\sigma_{0}^{2}\sigma_{2}^{2}-\sigma_{1}^{4})}\left\\{\sigma_{0}^{2}(-\zeta)^{2}+\sigma_{2}^{2}(\triangle\zeta)^{2}-2\sigma_{1}^{2}(-\zeta)(\triangle\zeta)\right\\}\right]\mathrm{d}(-\zeta)\mathrm{d}(\triangle\zeta).$
(B.7)

Similarly to the example of $-\zeta$ and $\triangle\zeta$, let us consider the
probability distribution function of the coefficients in Taylor expansion of
$\zeta$ up through the second order:

$\zeta=\zeta_{0}+\sum_{i}\zeta^{i}_{1}x_{i}+\sum_{i,j}\frac{1}{2}\zeta_{2}^{ij}x_{i}x_{j}+\mathcal{O}(x^{3}).$
(B.8)

There are 10 independent variables $\zeta_{0}$, $\zeta_{1}^{i}~{}(i=1,2,3)$
and $\zeta_{2}^{ij}~{}(i,j=1,2,3)$. Non-zero correlations between two of these
coefficients are summarized as follows:

$\displaystyle\sigma_{0}^{2}$ $\displaystyle=$ $\displaystyle\int\mathrm{d}\ln
k\mathcal{P}(k)=\langle\zeta_{0}\zeta_{0}\rangle,$ (B.9)
$\displaystyle\sigma_{1}^{2}$ $\displaystyle=$ $\displaystyle\int\mathrm{d}\ln
kk^{2}\mathcal{P}(k)=-3\langle\zeta_{0}\zeta_{2}^{ii}\rangle=3\langle\zeta_{1}^{i}\zeta_{1}^{i}\rangle,$
(B.10) $\displaystyle\sigma_{2}^{2}$ $\displaystyle=$
$\displaystyle\int\mathrm{d}\ln
kk^{4}\mathcal{P}(k)=-5\langle\zeta_{2}^{ii}\zeta_{2}^{ii}\rangle=15\langle\zeta_{2}^{ii}\zeta_{2}^{jj}\rangle=15\langle\zeta_{2}^{ij}\zeta_{2}^{ij}\rangle~{}{\rm
with}~{}i\neq j.$ (B.11)

The 2nd order variables can be transformed into the six independent variables
$\lambda_{i}$ and $\theta_{i}$ with $i=1,2,3$, where $\lambda_{1}$,
$\lambda_{2}$ and $\lambda_{3}$ are the three eigenvalues of the matrix
$\zeta_{2}^{ij}$ with $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$ and
$\theta_{1}$, $\theta_{2}$ and $\theta_{3}$ are the Euler angles to take the
principal direction. Then it can be shown that the volume element
$\mathrm{d}^{6}\zeta_{2}$ can be given by

$\mathrm{d}^{6}\zeta_{2}=(\lambda_{1}-\lambda_{2})(\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda_{3})\mathrm{d}^{3}\lambda\sin\theta_{1}\mathrm{d}^{3}\theta.$
(B.12)

The integration with respect to the Euler angles gives the factor $2\pi$.

We introduce the following 7 independent variables by using the remaining 7
variables:

$\displaystyle\nu$ $\displaystyle=$
$\displaystyle-\frac{\zeta_{0}}{\sigma_{0}},$ (B.13) $\displaystyle\eta_{i}$
$\displaystyle=$ $\displaystyle\frac{\zeta^{i}_{1}}{\sigma_{1}},$ (B.14)
$\displaystyle\xi_{1}$ $\displaystyle=$
$\displaystyle\frac{\lambda_{1}+\lambda_{2}+\lambda_{3}}{\sigma_{2}},$ (B.15)
$\displaystyle\xi_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{\lambda_{1}-\lambda_{3}}{\sigma_{2}},$ (B.16)
$\displaystyle\xi_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{\lambda_{1}-2\lambda_{2}+\lambda_{2}}{\sigma_{2}}.$
(B.17)

Then the probability distribution is transformed into the following form:

$P(\nu,\bm{\xi},\bm{\eta})\mathrm{d}\nu\mathrm{d}^{3}\xi\mathrm{d}^{3}\eta=P_{1}(\nu,\xi_{1};\gamma_{1})P_{2}(\xi_{2},\xi_{3})P_{3}(\bm{\eta})\mathrm{d}\nu\mathrm{d}^{3}\xi\mathrm{d}^{3}\eta,$
(B.18)

where

$\displaystyle P_{1}(\nu,\xi_{1};\gamma_{1})\mathrm{d}\nu\mathrm{d}\xi_{1}$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\frac{1}{\sqrt{1-\gamma_{1}^{2}}}\exp\left[-\frac{1}{2}\left(\nu^{2}+\frac{(\xi_{1}-\gamma_{1}\nu)^{2}}{1-\gamma_{1}^{2}}\right)\right]\mathrm{d}\nu\mathrm{d}\xi_{1},$
(B.19) $\displaystyle
P_{2}(\xi_{2},\xi_{3})\mathrm{d}\xi_{2}\mathrm{d}\xi_{3}$ $\displaystyle=$
$\displaystyle\frac{5^{5/2}3^{2}}{(2\pi)^{1/2}}\xi_{2}(\xi_{2}^{2}-\xi_{3}^{3})\exp\left[-\frac{5}{2}(3\xi_{2}^{2}+\xi_{3}^{3})\right]\mathrm{d}\xi_{2}\mathrm{d}\xi_{3},$
(B.20) $\displaystyle P_{3}(\bm{\eta})\mathrm{d}^{3}\eta$ $\displaystyle=$
$\displaystyle\frac{3^{3/2}}{(2\pi)^{3/2}}\exp\left[-\frac{3}{2}\left(\eta_{1}^{2}+\eta_{2}^{2}+\eta_{3}^{2}\right)\right]\mathrm{d}^{3}\eta$
(B.21)

with

$\gamma_{n}=\frac{\sigma_{n}^{2}}{\sigma_{n-1}\sigma_{n+1}}.$ (B.22)

It would be worthy to note that there are no correlations between
$(\nu,\xi_{1})$ and $(\xi_{2},\xi_{3})$. Therefore even if we consider a
rarely high peak with $\nu\sim\xi_{1}\gg 1$, the values of $\xi_{2}$ and
$\xi_{3}$ are expected to be $\sim\mathcal{O}(1)\ll\xi_{1}$. Since $\xi_{2}$
and $\xi_{3}$ are responsible for the asphericity of the peak, we can conclude
that rarely high peaks tend to be more spherically symmetric.

Let us remind 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 $\nu$ and
$\xi_{1}$. Therefore our goal in this section is to obtain the mean number
density characterized by the parameters $\nu$ and $\xi_{1}$:
$\bar{n}_{peak}(\nu,\xi_{1})$. As a first step, let us consider the number
density distribution of extrema $n_{\rm ext}(\bm{x},\nu,\xi_{1})$ in the space
spanned by $(\bm{x},\nu,\xi_{1})$. By definition, we have

$\displaystyle n_{\rm
ext}(\bm{x},\nu,\xi_{1})\Delta\bm{x}\Delta\nu\Delta\xi_{1}$ $\displaystyle=$
$\displaystyle\text{number of extrema in }\Delta\bm{x}\Delta\nu\Delta\xi_{1}$
(B.23) $\displaystyle=$
$\displaystyle\sum_{p}\delta(\bm{x}-\bm{x}_{p})\delta(\nu-\nu_{p})\delta(\xi_{1}-\xi_{1p})\Delta\bm{x}\Delta\nu\Delta\xi_{1},$
(B.24)

where $p$ labels the extrema so that the $p$-th extremum will be characterized
by $(\bm{x}_{p},\nu_{p},\xi_{1p})$. The integration over the space spanned by
$(\bm{x},\nu,\xi_{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 $(\bm{x}_{p},\nu_{p},\xi_{1p})$. Since we are
interested in the average distribution without the explicit dependence on
$(\bm{x}_{p},\nu_{p},\xi_{1p})$, let us take the volume average and the
ensemble average for $\nu_{p}$ and $\xi_{1p}$ of the distribution as follows:

$\displaystyle\bar{n}_{\rm ext}(\nu,\xi_{1})\Delta\nu\Delta\xi_{1}$
$\displaystyle:=$
$\displaystyle\frac{1}{V}\int\mathrm{d}^{3}x\prod_{p}\int\mathrm{d}\nu_{p}\mathrm{d}\xi_{1p}P_{1}(\nu_{p},\xi_{1p};\gamma_{1})n_{\rm
ext}(\bm{x},\nu,\xi_{1})\Delta\nu\Delta\xi_{1}$ (B.25) $\displaystyle=$
$\displaystyle
P_{1}(\nu,\xi_{1};\gamma_{1})\Delta\nu\Delta\xi_{1}\frac{1}{V}\int\mathrm{d}^{3}x\sum_{p}\delta(\bm{x}-\bm{x}_{p}).$
(B.26)

To dispose the $\delta(\bm{x}-\bm{x}_{p})$ and the volume integral, let us
consider rewriting the delta function. Since we are focusing on extrema, at
$\bm{x}=\bm{x}_{p}$, $\bm{\eta}$ vanishes. Therefore we find

$\displaystyle\delta(\bm{\eta}(\bm{x}))=\sum_{p}\left|\det\left(\frac{\partial\bm{\eta}}{\partial\bm{x}}\right)\right|_{p}^{-1}\delta(\bm{x}-\bm{x}_{p})=\left|\det\left(\frac{\partial\bm{\eta}}{\partial\bm{x}}\right)\right|^{-1}\sum_{p}\delta(\bm{x}-\bm{x}_{p})$
(B.27) $\displaystyle\Leftrightarrow$
$\displaystyle\sum_{p}\delta(\bm{x}-\bm{x}_{p})=\sigma_{1}^{-3}\left|\det\zeta_{2}^{ij}(\bm{x})\right|\delta(\bm{\eta})=\sigma_{1}^{-3}\left|\lambda_{1}\lambda_{2}\lambda_{3}\right|\delta(\bm{\eta}),$
(B.28)

where $\left|\lambda_{1}\lambda_{2}\lambda_{3}\right|$ can be written as

$\left|\lambda_{1}\lambda_{2}\lambda_{3}\right|=\frac{\sigma_{2}^{3}}{27}\left|\left((\xi_{1}+\xi_{2})^{2}-9\xi_{2}^{2}\right)\left(\xi_{1}-2\xi_{3}\right)\right|$
(B.29)

in terms of $\bm{\xi}$. Then we obtain

$\bar{n}_{\rm
ext}(\nu,\xi_{1})=P_{1}(\nu,\xi_{1};\gamma_{1})\frac{\sigma_{2}^{3}}{3^{3}\sigma_{1}^{3}}\frac{1}{V}\int\mathrm{d}\bm{x}\left|\left((\xi_{1}+\xi_{2}(\bm{x}))^{2}-9\xi_{2}^{2}(\bm{x})\right)\left(\xi_{1}-2\xi_{3}(\bm{x})\right)\right|\delta(\bm{\eta}(\bm{x})),$
(B.30)

where we have fixed the value of $\xi_{1}$ and take $\xi_{1}$ as an
independent variable of $\bm{x}$. Here we replace the volume average with the
ensemble average for $\xi_{2}$, $\xi_{3}$ and $\bm{\eta}$ as follows:

$\displaystyle\bar{n}_{\rm
ext}(\nu,\xi_{1})=P_{1}(\nu,\xi_{1};\gamma_{1})\frac{\sigma_{2}^{3}}{3^{3}\sigma_{1}^{3}}\int\mathrm{d}\xi_{2}\mathrm{d}\xi_{3}\mathrm{d}^{3}\eta
P_{2}(\xi_{2},\xi_{3})P_{3}(\bm{\eta})\left|\left((\xi_{1}+\xi_{2})^{2}-9\xi_{2}^{2}\right)\left(\xi_{1}-2\xi_{3}\right)\right|\delta(\bm{\eta})$
(B.31)
$\displaystyle~{}~{}=P_{1}(\nu,\xi_{1};\gamma_{1})\frac{\sigma_{2}^{3}}{3^{3}\sigma_{1}^{3}}P_{3}(\bm{\eta}=0)\int\mathrm{d}\xi_{2}\mathrm{d}\xi_{3}P_{2}(\xi_{2},\xi_{3})\left|\left((\xi_{1}+\xi_{2})^{2}-9\xi_{2}^{2}\right)\left(\xi_{1}-2\xi_{3}\right)\right|.$
(B.32)

The peak number density can be obtained by multiplying the step function
$\Theta(\lambda_{3})$. Then the integration for $\xi_{2}$ and $\xi_{3}$ can be
explicitly done [64] and we obtain

$n_{\rm
peak}(\nu,\xi_{1})=\frac{3^{3/2}}{(2\pi)^{3/2}}\frac{\sigma_{2}^{3}}{\sigma_{1}^{3}}f(\xi_{1})P_{1}(\nu,\xi_{1};\gamma_{1}),$
(B.33)

where

$\displaystyle f(\xi_{1})$ $\displaystyle:=$
$\displaystyle\frac{1}{2}\xi_{1}(\xi_{1}^{2}-3)\left({\rm
erf}\left[\frac{1}{2}\sqrt{\frac{5}{2}}\xi_{1}\right]+{\rm
erf}\left[\sqrt{\frac{5}{2}}\xi_{1}\right]\right)$ (B.35)
$\displaystyle+\sqrt{\frac{2}{5\pi}}\left\\{\left(\frac{8}{5}+\frac{31}{4}\xi_{1}^{2}\right)\exp\left[-\frac{5}{8}\xi_{1}^{2}\right]+\left(-\frac{8}{5}+\frac{1}{2}\xi_{1}^{2}\right)\exp\left[-\frac{5}{2}\xi_{1}^{2}\right]\right\\}.$

For convenience in the calculation of PBH abundance, let us change the
variables from $(\xi_{1},\nu)$ to $(\mu_{0},k_{*})$ defined by

$\displaystyle\mu_{0}$ $\displaystyle=$
$\displaystyle-\zeta_{0}=\sigma_{0}\nu,$ (B.36)
$\displaystyle\mu_{0}k_{*}^{2}$ $\displaystyle=$
$\displaystyle\triangle\zeta|_{\rm peak}=\sigma_{2}\xi_{1}.$ (B.37)

Then the peak number density can be transformed into the following form:

$n^{(k_{*})}_{\rm pk}(\mu_{0},k_{*})d\mu_{0}dk_{*}:=n_{\rm
peak}(\nu,\xi_{1})d\nu d\xi_{1}=\frac{2\cdot
3^{3/2}}{(2\pi)^{3/2}}\mu_{0}k_{*}\frac{\sigma_{2}^{2}}{\sigma_{0}\sigma_{1}^{3}}f\left(\frac{\mu_{0}k_{*}^{2}}{\sigma_{2}}\right)P_{1}\left(\frac{\mu_{0}}{\sigma_{0}},\frac{\mu_{0}k_{*}^{2}}{\sigma_{2}};\gamma_{1}\right)d\mu_{0}dk_{*}.$
(B.38)

##  Appendix C Typical profile in peak theory

The peak theory also provides a typical profile characterized by the peak
parameter $\nu$ and $\xi_{1}$ [64]. To derive the expression of the typical
profile, let us consider the conditional probability of taking the value
$\zeta(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

$P_{Y|X}=\frac{P_{X\cap Y}}{P_{X}},$ (C.1)

where $P_{X\cap 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

$q_{\mu}=(q_{1},q_{2},\cdots,q_{N})=(\alpha_{1},\cdots,\alpha_{n},\beta_{1},\cdots,\beta_{N-n}),$
(C.2)

where $\alpha_{a}~{}(a=1,\cdots,n)$ and $\beta_{i}~{}(i=1,\cdots,N-n)$
correspond to the variables for $X$ and $Y$, respectively. Letting $q^{\mu}$
denotes the components of the transposition of the row vector $q$ as
$q^{\mu}=(q^{\rm t})^{\mu}$, the probability distribution of $X\cup Y$ can be
written as

$P_{X\cup Y}(q_{\mu})=\frac{1}{(2\pi)^{N/2}(\det
M)^{1/2}}\exp\left[-\frac{1}{2}\sum_{\mu,\nu}q_{\mu}\left(M^{-1}\right)^{\mu}_{~{}\nu}q^{\nu}\right],$
(C.3)

where $M$ is the $N\times N$ correlation matrix. We decompose the correlation
matrix in the following form:

$M=\left(\begin{array}[]{cc}A^{a}_{~{}b}&B^{a}_{~{}j}\\\
C^{i}_{~{}b}&D^{i}_{~{}j}\end{array}\right).$ (C.4)

Then the inverse matrix of $M$ is given by

$M^{-1}=\left(\begin{array}[]{cc}A^{-1}+A^{-1}B(M/A)^{-1}CA^{-1}&-A^{-1}B(M/A)^{-1}\\\
-(M/A)^{-1}CA^{-1}&(M/A)^{-1}\end{array}\right),$ (C.5)

where $M/A$ is the Schure complement of the block $A$ given by

$M/A=D-CA^{-1}B.$ (C.6)

By using the expression of $M^{-1}$, we find

$q_{\mu}\left(M^{-1}\right)^{\mu}_{~{}\nu}q^{\nu}=\alpha A^{-1}\alpha^{\rm
t}+(\beta-\alpha A^{-1}B)(M/A)^{-1}(\beta-\alpha A^{-1}B)^{\rm t},$ (C.7)

where we have used the fact $B=C^{\rm t}$. Since the conditional probability
is given by Eq. (C.1), we find

$P_{Y|X}=\frac{1}{(2\pi)^{(N-n)/2}(\det
M/A)^{1/2}}\exp\left[-\frac{1}{2}(\beta-\alpha A^{-1}B)(M/A)^{-1}(\beta-\alpha
A^{-1}B)^{\rm t}\right].$ (C.8)

Therefore the mean values and the correlation matrix for $\beta$ under the
condition $X$ are given by $\alpha A^{-1}B$ and $M/A$.

We are interested in the probability of $\beta=\nu(r)=\zeta(r)/\sigma_{0}$
under the condition $\alpha_{a}=(\nu(0),\xi_{1}(0),\bm{\eta}(0)=0)$. To
simplify the expressions we write $\nu(0)$, $\xi_{1}(0)$ and $\bm{\eta}(0)$
without the argument as $\nu$, $\xi_{1}$ and $\bm{\eta}$, but the argument of
$\nu(r)$ will be kept. Then the 4 blocks of the correlation matrix can be
calculated as

$\displaystyle A^{a}_{~{}b}$ $\displaystyle=$
$\displaystyle<\alpha^{a}\alpha^{b}>=\left(\begin{array}[]{ccc}1&\gamma_{1}&0\\\
\gamma_{1}&1&0\\\ 0&0&\frac{1}{3}I_{3\times 3}\end{array}\right),$ (C.12)
$\displaystyle B^{a}_{~{}1}$ $\displaystyle=$ $\displaystyle
C^{1}_{~{}a}=\langle\alpha^{a}\nu(r)\rangle,$ (C.13) $\displaystyle D$
$\displaystyle=$ $\displaystyle\langle\nu(r)\nu(r)\rangle=1.$ (C.14)

Therefore we obtain

$\displaystyle\left.\left(\alpha
A^{-1}B\right)\right|_{\nu,\xi_{1},\bm{\eta}=0}$ $\displaystyle=$
$\displaystyle\frac{1}{1-\gamma_{1}^{2}}(\nu-\xi_{1}\gamma_{1})\langle\nu\nu(r)\rangle+\frac{1}{\gamma_{1}^{2}}(-\nu\gamma_{1}+\xi_{1})\langle\xi_{1}\nu(r)\rangle,$
(C.15) $\displaystyle M/A=D-CA^{-1}B$ $\displaystyle=$ $\displaystyle
1-\frac{1}{1-\gamma_{1}^{2}}\langle\nu\nu(r)\rangle^{2}+\frac{2\gamma_{1}}{1-\gamma_{1}^{2}}\langle\xi_{1}\nu(r)\rangle\langle\nu\nu(r)\rangle$
(C.17)
$\displaystyle-\frac{1}{1-\gamma_{1}^{2}}\langle\xi_{1}\nu(r)\rangle^{2}-3\sum_{i=1}^{3}\langle\eta_{i}\nu(r)\rangle^{2}.$

The correlations are calculated as

$\displaystyle\langle\nu\nu(r)\rangle$ $\displaystyle=$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k^{\prime}}{(2\pi)^{3}}\langle\tilde{\nu}^{*}(\bm{k})e^{-i\bm{k}\cdot\bm{r}}~{}\tilde{\nu}(\bm{k}^{\prime})\rangle=\sigma_{0}^{-2}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\frac{\mathrm{d}^{3}k^{\prime}}{(2\pi)^{3}}e^{-i\bm{k}\cdot\bm{r}}\langle\tilde{\zeta}(\bm{k})\tilde{\zeta}(\bm{k}^{\prime})\rangle$
(C.18) $\displaystyle=$
$\displaystyle\sigma_{0}^{-2}\int\frac{\mathrm{d}^{3}k}{4\pi}e^{-i\bm{k}\cdot\bm{r}}\frac{1}{k^{3}}\mathcal{P}(k)=\sigma_{0}^{-2}\int\mathrm{d}\ln
k\frac{\sin(kr)}{kr}\mathcal{P}(k)=:\psi_{0}(r),$ (C.19)
$\displaystyle\langle\xi_{1}\nu(r)\rangle$ $\displaystyle=$
$\displaystyle-\frac{\sigma_{0}}{\sigma_{2}}\triangle\psi_{0}(r),$ (C.20)
$\displaystyle\sum_{i=1}^{3}\langle\eta_{i}\nu(r)\rangle^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma_{1}^{2}}\left(\partial_{r}\psi_{0}(r)\right)^{2}.$
(C.21)

Then the mean profile is given by

$\langle\nu(r)|\nu,\xi_{1},\bm{\eta}=0\rangle=\frac{\nu}{1-\gamma_{1}^{2}}\left(\psi+\frac{1}{3}R_{1}^{2}\triangle\psi_{0}\right)-\frac{\xi_{1}}{\gamma_{1}(1-\gamma_{1}^{2})}\left(\gamma_{1}^{2}\psi+\frac{1}{3}R_{1}^{2}\triangle\psi_{0}\right),$
(C.22)

where

$R_{n}^{2}=\frac{3\sigma_{n}^{2}}{\sigma_{n+1}^{2}}.$ (C.23)

The variance is given by

$\displaystyle\langle(\Delta\nu(r))^{2}|\nu,\xi_{1},\bm{\eta}=0\rangle$
$\displaystyle=$ $\displaystyle
1-\frac{1}{1-\gamma_{1}^{2}}\psi_{0}^{2}-\frac{1}{\gamma_{1}^{2}(1-\gamma_{1}^{2})}\left(2\gamma_{1}^{2}\psi_{0}+\frac{R_{1}^{2}}{3}\triangle\psi_{0}\right)\frac{R_{1}^{2}}{3}\triangle\psi_{0}$
(C.25)
$\displaystyle-\frac{5}{\gamma_{1}^{2}}\left(\frac{\partial_{r}\psi_{0}}{r}-\frac{\triangle\psi_{0}}{3}\right)^{2}R_{1}^{4}-\frac{R_{1}^{2}}{\gamma_{1}^{2}}(\partial_{r}\psi_{0})^{2}.$

The variance is roughly estimated as
$\Delta\nu(r)\sim\Delta\zeta/\sigma_{0}=\mathcal{O}(1)$, that is,
$\Delta\zeta\sim\sigma_{0}$. Therefore the typical profile is a good
approximation as long as $\zeta(r)\gg\sigma_{0}$ for a rarely high peak.
Changing the independent variables from $(\nu,\xi_{1})$ to $(\mu_{0},k_{*})$,
we obtain the following form of the typical profile of $\bar{\zeta}(r)$:

$\frac{\langle\zeta(r)|\mu_{0},k_{*},\bm{\eta}=0\rangle}{\mu_{0}}=-\frac{1}{1-\gamma_{1}^{2}}\left(\psi_{0}+\frac{1}{3}R_{1}^{2}\triangle\psi_{0}\right)+k_{*}^{2}\frac{1}{\gamma_{1}(1-\gamma_{1}^{2})}\frac{\sigma_{0}}{\sigma_{2}}\left(\gamma_{1}^{2}\psi_{0}+\frac{1}{3}R_{1}^{2}\triangle\psi_{0}\right).$
(C.26)

##  Appendix D Equations for peaks of $\triangle\zeta$

In practice, peaks of $\triangle\zeta$ are more important for the calculation
of PBH abundance. Thus we consider the parameters $\mu_{2}=\triangle\zeta(0)$
and $k_{\bullet}^{2}:=-\frac{\triangle\triangle\zeta(0)}{\mu_{2}}$ instead of
$\mu_{0}$ and $k_{*}^{2}$. With this modification associated with
$-\zeta\rightarrow\triangle\zeta$, to obtain the number density and typical
profile of peaks of $\triangle\zeta$, we need the following replacements:

$\psi_{0}\rightarrow\psi_{2},~{}\mu_{0}\rightarrow\mu_{2},~{}k_{*}\rightarrow
k_{\bullet},~{}\sigma_{n}\rightarrow\sigma_{n+2}~{}(\gamma_{1}\rightarrow\gamma_{3},~{}R_{1}\rightarrow
R_{3}),$ (D.1)

where 4)4)4)The definition of $\psi_{2}$ is different from that in Refs. [63,
62], but the same as that in Ref. [76].

$\psi_{n}:=\frac{1}{\sigma_{n}}\int\mathrm{d}\ln
k~{}k^{2n}\frac{\sin(kr)}{kr}\mathcal{P}(k).$ (D.2)

Then, introducing the non-dimensional quantity
$\tilde{\mu}_{2}=\mu_{2}\sigma_{1}^{2}/\sigma_{2}^{2}$, we obtain the number
density and the typical profile as follows:

$\displaystyle n^{(k_{\bullet})}_{\rm
pk}(\tilde{\mu}_{2},k_{\bullet})d\tilde{\mu}_{2}dk_{\bullet}=\frac{2\cdot
3^{3/2}}{(2\pi)^{3/2}}\tilde{\mu}_{2}k_{\bullet}\frac{\sigma_{2}^{3}\sigma_{4}^{2}}{\sigma_{1}^{4}\sigma_{3}^{3}}f\left(\tilde{\mu}_{2}k_{\bullet}^{2}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}\sigma_{4}}\right)P_{1}\left(\tilde{\mu}_{2}\frac{\sigma_{2}}{\sigma_{1}^{2}},\tilde{\mu}_{2}k_{\bullet}^{2}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}\sigma_{4}};\gamma_{3}\right)d\tilde{\mu}_{2}dk_{\bullet},$
(D.3)
$\displaystyle\frac{\overline{\triangle\zeta}(r)}{\tilde{\mu}_{2}}=\frac{1}{1-\gamma_{3}^{2}}\left(\psi_{2}+\frac{1}{3}R_{3}^{2}\triangle\psi_{2}\right)-k_{\bullet}^{2}\frac{1}{\gamma_{3}(1-\gamma_{3}^{2})}\frac{\sigma_{2}}{\sigma_{4}}\left(\gamma_{3}^{2}\psi_{2}+\frac{1}{3}R_{3}^{2}\triangle\psi_{2}\right).~{}$
(D.4)

The profile of $\bar{\zeta}$ associated with $\overline{\triangle\zeta}$ can
be obtained by integrating (D.4) as follows:

$\frac{\overline{\zeta}(r)}{\tilde{\mu}_{2}}=-\frac{1}{1-\gamma_{3}^{2}}\left(\psi_{1}+\frac{1}{3}R_{3}^{2}\triangle\psi_{1}\right)+k_{\bullet}^{2}\frac{1}{\gamma_{3}(1-\gamma_{3}^{2})}\frac{\sigma_{2}}{\sigma_{4}}\left(\gamma_{3}^{2}\psi_{1}+\frac{1}{3}R_{3}^{2}\triangle\psi_{1}\right)+\frac{\sigma_{2}}{\sigma_{1}^{2}}\zeta_{\infty}.$
(D.5)

It can be shown that $\zeta_{\infty}$ can be also regarded as a Gaussian
probability variable (see Ref. [63]).

## References

  * [1] Y. B. Zel’dovich and I. D. Novikov,  Soviet Ast. 10, 602 (1967), The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model.
  * [2] S. Hawking,  Mon. Not. Roy. Astron. Soc. 152, 75 (1971), Gravitationally collapsed objects of very low mass.
  * [3] B. J. Carr and S. Hawking,  Mon. Not. Roy. Astron. Soc. 168, 399 (1974), Black holes in the early Universe.
  * [4] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama,  Phys. Rev. D81, 104019 (2010), arXiv:0912.5297, New cosmological constraints on primordial black holes.
  * [5] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama,  (2020), arXiv:2002.12778, Constraints on Primordial Black Holes.
  * [6] B. Carr and F. Kuhnel,  (2020), arXiv:2006.02838, Primordial Black Holes as Dark Matter: Recent Developments.
  * [7] S. Bird et al.,  Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464, Did LIGO detect dark matter?
  * [8] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama,  Phys. Rev. Lett. 117, 061101 (2016), arXiv:1603.08338, Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914.
  * [9] S. Clesse and J. García-Bellido,  (2017), arXiv:1711.10458, Seven Hints for Primordial Black Hole Dark Matter.
  * [10] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama,  Class. Quant. Grav. 35, 063001 (2018), arXiv:1801.05235, Primordial black holes-perspectives in gravitational wave astronomy.
  * [11] H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki,  Phys. Rev. D 99, 083503 (2019), arXiv:1901.07120, Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events.
  * [12] NANOGrav, Z. Arzoumanian et al.,  Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496, The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background.
  * [13] V. Vaskonen and H. Veermäe,  Phys. Rev. Lett. 126, 051303 (2021), arXiv:2009.07832, Did NANOGrav see a signal from primordial black hole formation?
  * [14] V. De Luca, G. Franciolini, and A. Riotto,  Phys. Rev. Lett. 126, 041303 (2021), arXiv:2009.08268, NANOGrav Data Hints at Primordial Black Holes as Dark Matter.
  * [15] K. Kohri and T. Terada,  Phys. Lett. B 813, 136040 (2021), arXiv:2009.11853, Solar-Mass Primordial Black Holes Explain NANOGrav Hint of Gravitational Waves.
  * [16] S. Sugiyama et al.,  Phys. Lett. B 814, 136097 (2021), arXiv:2010.02189, Testing Stochastic Gravitational Wave Signals from Primordial Black Holes with Optical Telescopes.
  * [17] G. Domènech and S. Pi,  Sci. China Phys. Mech. Astron. 65, 230411 (2022), arXiv:2010.03976, NANOGrav hints on planet-mass primordial black holes.
  * [18] K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida,  Phys. Rev. Lett. 126, 131301 (2021), arXiv:2011.01270, NANOGrav Results and LIGO-Virgo Primordial Black Holes in Axionlike Curvaton Models.
  * [19] M. Kawasaki, A. Kusenko, and T. T. Yanagida,  Phys. Lett. B 711, 1 (2012), arXiv:1202.3848, Primordial seeds of supermassive black holes.
  * [20] K. Kohri, T. Nakama, and T. Suyama,  Phys. Rev. D 90, 083514 (2014), arXiv:1405.5999, Testing scenarios of primordial black holes being the seeds of supermassive black holes by ultracompact minihalos and CMB $\mu$-distortions.
  * [21] T. Nakama, T. Suyama, and J. Yokoyama,  Phys. Rev. D 94, 103522 (2016), arXiv:1609.02245, Supermassive black holes formed by direct collapse of inflationary perturbations.
  * [22] B. Carr and J. Silk,  Mon. Not. Roy. Astron. Soc. 478, 3756 (2018), arXiv:1801.00672, Primordial Black Holes as Generators of Cosmic Structures.
  * [23] P. D. Serpico, V. Poulin, D. Inman, and K. Kohri,  Phys. Rev. Res. 2, 023204 (2020), arXiv:2002.10771, Cosmic microwave background bounds on primordial black holes including dark matter halo accretion.
  * [24] C. Ünal, E. D. Kovetz, and S. P. Patil,  Phys. Rev. D 103, 063519 (2021), arXiv:2008.11184, Multimessenger probes of inflationary fluctuations and primordial black holes.
  * [25] K. Kohri, T. Sekiguchi, and S. Wang,  Phys. Rev. D 106, 043539 (2022), arXiv:2201.05300, Cosmological 21-cm line observations to test scenarios of super-Eddington accretion on to black holes being seeds of high-redshifted supermassive black holes.
  * [26] J. C. Niemeyer and K. Jedamzik,  Phys. Rev. Lett. 80, 5481 (1998), arXiv:astro-ph/9709072, Near-critical gravitational collapse and the initial mass function of primordial black holes.
  * [27] J. C. Niemeyer and K. Jedamzik,  Phys. Rev. D59, 124013 (1999), arXiv:astro-ph/9901292, Dynamics of primordial black hole formation.
  * [28] M. Shibata and M. Sasaki,  Phys. Rev. D60, 084002 (1999), arXiv:gr-qc/9905064, Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity.
  * [29] I. Hawke and J. M. Stewart,  Classical and Quantum Gravity 19, 3687 (2002), The dynamics of primordial black-hole formation.
  * [30] I. Musco, J. C. Miller, and L. Rezzolla,  Class. Quant. Grav. 22, 1405 (2005), arXiv:gr-qc/0412063, Computations of primordial black hole formation.
  * [31] 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)].
  * [32] T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama,  JCAP 1401, 037 (2014), arXiv:1310.3007, Identifying the most crucial parameters of the initial curvature profile for primordial black hole formation.
  * [33] T. Nakama,  JCAP 1410, 040 (2014), arXiv:1408.0955, The double formation of primordial black holes.
  * [34] 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 formation.
  * [35] I. Musco,  Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127, Threshold for primordial black holes: Dependence on the shape of the cosmological perturbations.
  * [36] A. Escrivà,  Phys. Dark Univ. 27, 100466 (2020), arXiv:1907.13065, Simulation of primordial black hole formation using pseudo-spectral methods.
  * [37] A. Escrivà, C. Germani, and R. K. Sheth,  Phys. Rev. D101, 044022 (2020), arXiv:1907.13311, Universal threshold for primordial black hole formation.
  * [38] A. Escrivà, C. Germani, and R. K. Sheth,  (2020), arXiv:2007.05564, Analytical thresholds for black hole formation in general cosmological backgrounds.
  * [39] A. Escrivà,  Universe 8, 66 (2022), arXiv:2111.12693, PBH Formation from Spherically Symmetric Hydrodynamical Perturbations: A Review.
  * [40] A. Escrivà, E. Bagui, and S. Clesse,  (2022), arXiv:2209.06196, Simulations of PBH formation at the QCD epoch and comparison with the GWTC-3 catalog.
  * [41] G. Franciolini, I. Musco, P. Pani, and A. Urbano,  (2022), arXiv:2209.05959, From inflation to black hole mergers and back again: Gravitational-wave data-driven constraints on inflationary scenarios with a first-principle model of primordial black holes across the QCD epoch.
  * [42] D. Toussaint, S. B. Treiman, F. Wilczek, and A. Zee,  Phys. Rev. D 19, 1036 (1979), Matter-antimatter accounting, thermodynamics, and black-hole radiation.
  * [43] J. D. Barrow, E. J. Copeland, E. W. Kolb, and A. R. Liddle,  Phys. Rev. D 43, 984 (1991), Baryogenesis in extended inflation. 2. Baryogenesis via primordial black holes.
  * [44] D. Baumann, P. J. Steinhardt, and N. Turok,  (2007), arXiv:hep-th/0703250, Primordial Black Hole Baryogenesis.
  * [45] K. Inomata, M. Kawasaki, K. Mukaida, T. Terada, and T. T. Yanagida,  Phys. Rev. D 101, 123533 (2020), arXiv:2003.10455, Gravitational Wave Production right after a Primordial Black Hole Evaporation.
  * [46] T. Papanikolaou, V. Vennin, and D. Langlois,  JCAP 03, 053 (2021), arXiv:2010.11573, Gravitational waves from a universe filled with primordial black holes.
  * [47] G. Domènech, C. Lin, and M. Sasaki,  JCAP 04, 062 (2021), arXiv:2012.08151, Gravitational wave constraints on the primordial black hole dominated early universe,  [Erratum: JCAP 11, E01 (2021)].
  * [48] N. Bhaumik and R. K. Jain,  Phys. Rev. D 104, 023531 (2021), arXiv:2009.10424, Small scale induced gravitational waves from primordial black holes, a stringent lower mass bound, and the imprints of an early matter to radiation transition.
  * [49] D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott, and R. Petrossian-Byrne,  (2020), arXiv:2004.00618, Hot Gravitons and Gravitational Waves From Kerr Black Holes in the Early Universe.
  * [50] G. Domènech, V. Takhistov, and M. Sasaki,  Phys. Lett. B 823, 136722 (2021), arXiv:2105.06816, Exploring evaporating primordial black holes with gravitational waves.
  * [51] K. N. Ananda, C. Clarkson, and D. Wands,  Phys. Rev. D 75, 123518 (2007), arXiv:gr-qc/0612013, The Cosmological gravitational wave background from primordial density perturbations.
  * [52] D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki,  Phys. Rev. D 76, 084019 (2007), arXiv:hep-th/0703290, Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations.
  * [53] R. Saito and J. Yokoyama,  Phys. Rev. Lett. 102, 161101 (2009), arXiv:0812.4339, Gravitational wave background as a probe of the primordial black hole abundance,  [Erratum: Phys.Rev.Lett. 107, 069901 (2011)].
  * [54] R. Saito and J. Yokoyama,  Prog. Theor. Phys. 123, 867 (2010), arXiv:0912.5317, Gravitational-Wave Constraints on the Abundance of Primordial Black Holes,  [Erratum: Prog.Theor.Phys. 126, 351–352 (2011)].
  * [55] H. Assadullahi and D. Wands,  Phys. Rev. D 79, 083511 (2009), arXiv:0901.0989, Gravitational waves from an early matter era.
  * [56] J. R. Espinosa, D. Racco, and A. Riotto,  JCAP 09, 012 (2018), arXiv:1804.07732, A Cosmological Signature of the SM Higgs Instability: Gravitational Waves.
  * [57] K. Kohri and T. Terada,  Phys. Rev. D 97, 123532 (2018), arXiv:1804.08577, Semianalytic calculation of gravitational wave spectrum nonlinearly induced from primordial curvature perturbations.
  * [58] G. Domènech,  Universe 7, 398 (2021), arXiv:2109.01398, Scalar Induced Gravitational Waves Review.
  * [59] A. Escrivà, F. Kuhnel, and Y. Tada,  (2022), arXiv:2211.05767, Primordial Black Holes.
  * [60] P. Villanueva-Domingo, O. Mena, and S. Palomares-Ruiz,  Front. Astron. Space Sci. 8, 87 (2021), arXiv:2103.12087, A brief review on primordial black holes as dark matter.
  * [61] C. T. Byrnes and P. S. Cole,  (2021), arXiv:2112.05716, Lecture notes on inflation and primordial black holes.
  * [62] C.-M. Yoo, T. Harada, J. Garriga, and K. Kohri,  PTEP 2018, 123E01 (2018), arXiv:1805.03946, Primordial black hole abundance from random Gaussian curvature perturbations and a local density threshold.
  * [63] C.-M. Yoo, T. Harada, S. Hirano, and K. Kohri,  (2020), arXiv:2008.02425, Abundance of Primordial Black Holes in Peak Theory for an Arbitrary Power Spectrum.
  * [64] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay,  ApJ304, 15 (1986), The statistics of peaks of Gaussian random fields.
  * [65] M. Kopp, S. Hofmann, and J. Weller,  Phys. Rev. D83, 124025 (2011), arXiv:1012.4369, Separate Universes Do Not Constrain Primordial Black Hole Formation.
  * [66] B. J. Carr,  Astrophys. J. 201, 1 (1975), The Primordial black hole mass spectrum.
  * [67] C. W. Misner and D. H. Sharp,  Phys. Rev. 136, B571 (1964), Relativistic equations for adiabatic, spherically symmetric gravitational collapse.
  * [68] A. G. Polnarev and I. Musco,  Class. Quant. Grav. 24, 1405 (2007), arXiv:gr-qc/0605122, Curvature profiles as initial conditions for primordial black hole formation.
  * [69] D. H. Lyth, K. A. Malik, and M. Sasaki,  JCAP 05, 004 (2005), arXiv:astro-ph/0411220, A General proof of the conservation of the curvature perturbation.
  * [70] D. K. Nadezhin, I. D. Novikov, and A. G. Polnarev,  Soviet Ast. 22, 129 (1978), The hydrodynamics of primordial black hole formation.
  * [71] I. D. Novikov and A. G. Polnarev,  Soviet Ast. 24, 147 (1980), The Hydrodynamics of Primordial Black Hole Formation - Dependence on the Equation of State.
  * [72] I. Musco, J. C. Miller, and A. G. Polnarev,  Class. Quant. Grav. 26, 235001 (2009), arXiv:0811.1452, Primordial black hole formation in the radiative era: Investigation of the critical nature of the collapse.
  * [73] I. Musco and J. C. Miller,  Class. Quant. Grav. 30, 145009 (2013), arXiv:1201.2379, Primordial black hole formation in the early universe: critical behaviour and self-similarity.
  * [74] T. Harada, C.-M. Yoo, and Y. Koga,  In preparation.
  * [75] C.-M. Yoo, J.-O. Gong, and S. Yokoyama,  JCAP 09, 033 (2019), arXiv:1906.06790, Abundance of primordial black holes with local non-Gaussianity in peak theory.
  * [76] N. Kitajima, Y. Tada, S. Yokoyama, and C.-M. Yoo,  JCAP 10, 053 (2021), arXiv:2109.00791, Primordial black holes in peak theory with a non-Gaussian tail.
  * [77] A. Escrivà, Y. Tada, S. Yokoyama, and C.-M. Yoo,  JCAP 05, 012 (2022), arXiv:2202.01028, Simulation of primordial black holes with large negative non-Gaussianity.
  * [78] M. W. Choptuik,  Phys. Rev. Lett. 70, 9 (1993), Universality and scaling in gravitational collapse of a massless scalar field.
  * [79] T. Koike, T. Hara, and S. Adachi,  Phys. Rev. Lett. 74, 5170 (1995), arXiv:gr-qc/9503007, Critical behavior in gravitational collapse of radiation fluid: A Renormalization group (linear perturbation) analysis.
  * [80] J. Yokoyama,  Phys. Rev. D58, 107502 (1998), arXiv:gr-qc/9804041, Cosmological constraints on primordial black holes produced in the near critical gravitational collapse.
  * [81] 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.
  * [82] 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.
  * [83] C. Germani and I. Musco,  Phys. Rev. Lett. 122, 141302 (2019), arXiv:1805.04087, Abundance of Primordial Black Holes Depends on the Shape of the Inflationary Power Spectrum.
  * [84] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto,  Phys. Lett. B 820, 136602 (2021), arXiv:2105.07810, The formation probability of primordial black holes.
  * [85] C.-M. Yoo, T. Harada, S. Hirano, H. Okawa, and M. Sasaki,  Phys. Rev. D 105, 103538 (2022), arXiv:2112.12335, Primordial black hole formation from massless scalar isocurvature.
  * [86] C. Germani and R. K. Sheth,  Phys. Rev. D 101, 063520 (2020), arXiv:1912.07072, Nonlinear statistics of primordial black holes from Gaussian curvature perturbations.
  * [87] I. Musco, V. De Luca, G. Franciolini, and A. Riotto,  Phys. Rev. D 103, 063538 (2021), arXiv:2011.03014, Threshold for primordial black holes. II. A simple analytic prescription.
  * [88] T. Suyama and S. Yokoyama,  PTEP 2020, 023E03 (2020), arXiv:1912.04687, A novel formulation of the PBH mass function.
  * [89] J. S. Bullock and J. R. Primack,  Phys. Rev. D 55, 7423 (1997), arXiv:astro-ph/9611106, NonGaussian fluctuations and primordial black holes from inflation.
  * [90] P. Ivanov,  Phys. Rev. D 57, 7145 (1998), arXiv:astro-ph/9708224, Nonlinear metric perturbations and production of primordial black holes.
  * [91] J. C. Hidalgo,  (2007), arXiv:0708.3875, The effect of non-Gaussian curvature perturbations on the formation of primordial black holes.
  * [92] C. T. Byrnes, E. J. Copeland, and A. M. Green,  Phys. Rev. D 86, 043512 (2012), arXiv:1206.4188, Primordial black holes as a tool for constraining non-Gaussianity.
  * [93] E. V. Bugaev and P. A. Klimai,  Int. J. Mod. Phys. D 22, 1350034 (2013), arXiv:1303.3146, Primordial black hole constraints for curvaton models with predicted large non-Gaussianity.
  * [94] S. Young, D. Regan, and C. T. Byrnes,  JCAP 02, 029 (2016), arXiv:1512.07224, Influence of large local and non-local bispectra on primordial black hole abundance.
  * [95] T. Nakama, J. Silk, and M. Kamionkowski,  Phys. Rev. D 95, 043511 (2017), arXiv:1612.06264, Stochastic gravitational waves associated with the formation of primordial black holes.
  * [96] K. Ando, K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida,  Phys. Rev. D 97, 123512 (2018), arXiv:1711.08956, Primordial black holes for the LIGO events in the axionlike curvaton model.
  * [97] G. Franciolini, A. Kehagias, S. Matarrese, and A. Riotto,  JCAP 1803, 016 (2018), arXiv:1801.09415, Primordial Black Holes from Inflation and non-Gaussianity.
  * [98] R.-g. Cai, S. Pi, and M. Sasaki,  Phys. Rev. Lett. 122, 201101 (2019), arXiv:1810.11000, Gravitational Waves Induced by non-Gaussian Scalar Perturbations.
  * [99] V. Atal and C. Germani,  Phys. Dark Univ. 24, 100275 (2019), arXiv:1811.07857, The role of non-gaussianities in Primordial Black Hole formation.
  * [100] S. Passaglia, W. Hu, and H. Motohashi,  Phys. Rev. D 99, 043536 (2019), arXiv:1812.08243, Primordial black holes and local non-Gaussianity in canonical inflation.
  * [101] M. Taoso and A. Urbano,  JCAP 08, 016 (2021), arXiv:2102.03610, Non-gaussianities for primordial black hole formation.
  * [102] V. Atal, J. Garriga, and A. Marcos-Caballero,  JCAP 09, 073 (2019), arXiv:1905.13202, Primordial black hole formation with non-Gaussian curvature perturbations.
  * [103] V. Atal, J. Cid, A. Escrivà, and J. Garriga,  JCAP 05, 022 (2020), arXiv:1908.11357, PBH in single field inflation: the effect of shape dispersion and non-Gaussianities.
  * [104] M. Yu. Khlopov and A. G. Polnarev,  Phys. Lett. 97B, 383 (1980), PRIMORDIAL BLACK HOLES AS A COSMOLOGICAL TEST OF GRAND UNIFICATION.
  * [105] A. G. Polnarev and M. Y. Khlopov,  Sov. Phys. Usp. 28, 213 (1985), COSMOLOGY, PRIMORDIAL BLACK HOLES, AND SUPERMASSIVE PARTICLES.
  * [106] T. Harada, C.-M. Yoo, K. Kohri, K.-i. Nakao, and S. Jhingan,  Astrophys. J. 833, 61 (2016), arXiv:1609.01588, Primordial black hole formation in the matter-dominated phase of the Universe.
  * [107] T. Harada, C.-M. Yoo, K. Kohri, and K.-I. Nakao,  (2017), arXiv:1707.03595, Spins of primordial black holes formed in the matter-dominated phase of the Universe.
  * [108] T. Kokubu, K. Kyutoku, K. Kohri, and T. Harada,  Phys. Rev. D 98, 123024 (2018), arXiv:1810.03490, Effect of Inhomogeneity on Primordial Black Hole Formation in the Matter Dominated Era.
  * [109] C.-M. Yoo, T. Harada, and H. Okawa,  (2020), arXiv:2004.01042, Threshold of Primordial Black Hole Formation in Nonspherical Collapse.
  * [110] T. Chiba and S. Yokoyama,  PTEP 2017, 083E01 (2017), arXiv:1704.06573, Spin Distribution of Primordial Black Holes.
  * [111] V. De Luca, G. Franciolini, P. Pani, and A. Riotto,  JCAP 04, 052 (2020), arXiv:2003.02778, The evolution of primordial black holes and their final observable spins.
  * [112] M. Mirbabayi, A. Gruzinov, and J. Noreña,  JCAP 03, 017 (2020), arXiv:1901.05963, Spin of Primordial Black Holes.
  * [113] M. He and T. Suyama,  Phys. Rev. D 100, 063520 (2019), arXiv:1906.10987, Formation threshold of rotating primordial black holes.
  * [114] M. M. Flores and A. Kusenko,  Phys. Rev. D 104, 063008 (2021), arXiv:2106.03237, Spins of primordial black holes formed in different cosmological scenarios.
  * [115] Y. Koga, T. Harada, Y. Tada, S. Yokoyama, and C.-M. Yoo,  (2022), arXiv:2208.00696, Effective inspiral spin distribution of primordial black hole binaries.
  * [116] K. Jedamzik,  Phys. Rev. D 55, 5871 (1997), arXiv:astro-ph/9605152, Primordial black hole formation during the QCD epoch.
  * [117] C. Schmid, D. J. Schwarz, and P. Widerin,  Phys. Rev. D 59, 043517 (1999), arXiv:astro-ph/9807257, Amplification of cosmological inhomogeneities from the QCD transition.
  * [118] P. Widerin and C. Schmid,  (1998), arXiv:astro-ph/9808142, Primordial black holes from the QCD transition?
  * [119] T. Boeckel, S. Schettler, and J. Schaffner-Bielich,  Prog. Part. Nucl. Phys. 66, 266 (2011), arXiv:1012.3342, The Cosmological QCD Phase Transition Revisited.
  * [120] J. L. G. Sobrinho, P. Augusto, and A. L. Gonçalves,  Mon. Not. Roy. Astron. Soc. 463, 2348 (2016), arXiv:1609.01205, New thresholds for Primordial Black Hole formation during the QCD phase transition.
  * [121] C. T. Byrnes, M. Hindmarsh, S. Young, and M. R. S. Hawkins,  JCAP 08, 041 (2018), arXiv:1801.06138, Primordial black holes with an accurate QCD equation of state.
  * [122] B. Carr, S. Clesse, J. García-Bellido, and F. Kühnel,  Phys. Dark Univ. 31, 100755 (2021), arXiv:1906.08217, Cosmic conundra explained by thermal history and primordial black holes.
  * [123] B. Carr, S. Clesse, and J. García-Bellido,  Mon. Not. Roy. Astron. Soc. 501, 1426 (2021), arXiv:1904.02129, Primordial black holes from the QCD epoch: Linking dark matter, baryogenesis and anthropic selection.
  * [124] S. Clesse and J. Garcia-Bellido,  (2020), arXiv:2007.06481, GW190425, GW190521 and GW190814: Three candidate mergers of primordial black holes from the QCD epoch.
  * [125] F. Gao and I. M. Oldengott,  Phys. Rev. Lett. 128, 131301 (2022), arXiv:2106.11991, Cosmology Meets Functional QCD: First-Order Cosmic QCD Transition Induced by Large Lepton Asymmetries.
  * [126] J. I. Juan, P. D. Serpico, and G. Franco Abellán,  JCAP 07, 009 (2022), arXiv:2204.07027, The QCD phase transition behind a PBH origin of LIGO/Virgo events?
  * [127] S. Passaglia and M. Sasaki,  (2021), arXiv:2109.12824, Primordial Black Holes from CDM Isocurvature.
  * [128] H. Kodama,  Prog. Theor. Phys. 63, 1217 (1980), Conserved Energy Flux for the Spherically Symmetric System and the Back Reaction Problem in the Black Hole Evaporation.
  * [129] S. A. Hayward,  Phys. Rev. D 49, 6467 (1994), General laws of black hole dynamics.
  * [130] S. A. Hayward,  Phys. Rev. D 53, 1938 (1996), arXiv:gr-qc/9408002, Gravitational energy in spherical symmetry.