Abundance of primordial black holes
with local non-Gaussianity in peak theory

Chul-Moon Yooa, Jinn-Ouk Gongb and Shuichiro Yokoyamac,d

aDivision of Particle and Astrophysical Science, Graduate School of Science
Nagoya University, Nagoya 464-8602, Japan
bKorea Astronomy and Space Science Institute, Daejeon 34055, Korea
cKobayashi-Maskawa Institute, Nagoya University, Nagoya 464-8602, Japan
dKavli Institute for the Physics and Mathematics of the Universe, The
University of Tokyo
Kashiwa, Chiba 277-8583, Japan

We discuss the effect of local type non-Gaussianity on the abundance of
primordial black holes (PBH) based on the peak theory. We provide the PBH
formation criterion based on the so-called compaction function and use the
peak theory statistics associated with the curvature perturbation with the
local type non-Gaussianity. Providing a method to estimate the PBH abundance,
we demonstrate the effects of non-Gaussianity. It is explicitly shown that the
value of non-linear parameter $|f_{\rm NL}|\sim 1$ induces a similar effect to
a few factors of difference in the amplitude of the power spectrum.

##  1 Introduction

Primordial black holes (PBHs) repeatedly come up as a hot topic in cosmology
and astrophysics since the pioneering works of Zel’dovich and Novikov [1] and
Hawking [2]. For instance, PBHs are still fascinating candidate for dark
matter (e.g., see [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] and references
therein), and observational constraints are frequently updated (see e.g. [15,
10]). Number of interesting scenarios which produce substantial number of PBHs
are proposed (see e.g. a recent review [16] and references therein).

One of inevitable issues is the estimation of PBH abundance. Simplest
conventional way to estimate PBH abundance is the Press-Schechter (PS)
formalism [17], where the Gaussian distribution of a perturbation variable and
its threshold for PBH formation are assumed. The validity of these assumptions
has been under discussion for a long time (see e.g. [4, 18, 19, 21, 20, 22,
23, 24, 25, 26, 27, 28, 29, 30, 31, 32]). Recently, a more plausible procedure
for PBH formation in a radiation dominated era is proposed in [33], where an
enhanced feature is assumed in the primordial power spectrum around some
specific scale which enters the horizon in a radiation dominated era. Another
important assumption is the random Gaussian distribution of the curvature
perturbation. Then, a reliable procedure for the estimation of PBH abundance
is derived based on the peak theory [34, 35] of the curvature perturbation, by
taking the non-linear effect between the curvature perturbation and the
density perturbation into account (see e.g. [36, 37, 38] for the significance
of the non-linearity).

The purpose of this article is to show the impact of the local type non-
Gaussianity of the curvature perturbation extending the procedure given in
[33]. We note that in [33], the intrinsic non-Gaussianity that originates from
the non-linear relation between the curvature perturbation and the density
perturbation is taken into account. However, throughout [33] the curvature
perturbation is assumed to be a random Gaussian variable. In this article, we
provide a method to introduce the local type non-Gaussianity of the curvature
perturbation into the procedure provided in [33]. The effect of non-
Gaussianity on the PBH abundance and the spatial distribution of PBHs have
been discussed in [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
54, 55, 56, 57]. However, since a new plausible procedure is proposed, it is
important to check how the previously reported effect of the non-Gaussianity
works in the new procedure. Providing a method to estimate PBH abundance
including the local type non-Gaussianity, we demonstrate the effect of the
non-Gaussianity of the curvature perturbation.

This article is organized as follows. In Section 2, the typical profile of the
curvature perturbation is derived for a given scale and amplitude of the peak
with the non-linear parameter $f_{\rm NL}$. The general implicit expression
for the PBH fraction is derived in Section 3. Two specific primordial power
spectra are shown as examples to present the impact of non-Gaussianity in
Section 4. Section 5 is devoted to a short summary. Some technical detail is
given in the appendices. Throughout this article, we use the geometrized units
in which both the speed of light and the Newton’s gravitational constant are
unity, $c=G=1$.

##  2 Typical peak profile with local non-Gaussianity

In this article, we consider the spatial metric given by

$ds_{3}^{2}=a^{2}e^{-2\zeta}\tilde{\gamma}_{ij}dx^{i}dx^{j}\,,$ (1)

and now we assume that the reference spatial metric $\tilde{\gamma}_{ij}$ is
flat, and that the curvature perturbation $\zeta$ is given by the following
form with the local type non-Gaussianity:

$\zeta=\zeta_{\rm G}-\frac{3}{5}f_{\rm NL}\left(\zeta_{\rm
G}^{2}-\left\langle\zeta_{\rm G}^{2}\right\rangle\right)\,,$ (2)

where the brackets denote the ensemble average, and $\zeta_{\rm G}$ is the
random Gaussian distribution with the power spectrum ${\cal P}(k)$ defined by

$\left\langle\tilde{\zeta}_{\rm G}^{*}(\bm{k})\tilde{\zeta}_{\rm
G}(\bm{k}^{\prime})\right\rangle=(2\pi)^{3}\delta(\bm{k}-\bm{k}^{\prime})\frac{2\pi^{2}}{k^{3}}{\cal
P}(k)\,,$ (3)

where $\tilde{\zeta}_{\rm G}(\bm{k})$ is the Fourier transform of $\zeta_{\rm
G}$. Note that we follow the same sign convention of $\zeta$ as [58, 33], so
the negative sign of the non-linear term in (2) keeps $f_{\rm NL}$ consistent
with the conventional definition, e.g. adopted in the Planck 2018 report on
non-Gaussianity [59].

Let us focus on a high peak of $\zeta_{\rm G}$. Here we assume that all peaks
of $\zeta_{\rm G}$ are peaks of $\zeta$. The validity of this assumption is
discussed in Appendix A. According to the peak theory [35], for a high peak at
$\bm{r}=\bm{0}$, we obtain the probability distribution for the spherical
profile of $\zeta_{\rm G}(r)$ as follows:

$P(\zeta_{\rm
G}(r)|\mu,k_{*})=\frac{1}{\sqrt{2\pi}\sigma_{\zeta}}\exp\left\\{-\frac{1}{2\sigma_{\zeta}^{2}}\left[\zeta_{\rm
G}(r)-\bar{\zeta}_{\rm G}(r)\right]^{2}\right\\}\,,$ (4)

where $\mu$ and $k_{*}$ are respectively the amplitude and the curvature scale
of the Gaussian peak:

$\displaystyle\mu$ $\displaystyle=-\left.\zeta_{\rm G}\right|_{r=0}\,,$ (5)
$\displaystyle k_{*}^{2}$ $\displaystyle=\frac{\left.\triangle\zeta_{\rm
G}\right|_{r=0}}{\mu}\,.$ (6)

The mean $\bar{\zeta}_{\rm G}(r)$ and the variance $\sigma_{\zeta}$ are
respectively given as

$\displaystyle-\frac{\bar{\zeta}_{\rm G}(r)}{\sigma_{0}}$
$\displaystyle=\frac{\mu/\sigma_{0}}{1-\gamma^{2}}\left(\psi+\frac{1}{3}R_{*}^{2}\triangle\psi\right)-\frac{\mu
k_{*}^{2}/\sigma_{2}}{\gamma(1-\gamma^{2})}\left(\gamma^{2}\psi+\frac{1}{3}R_{*}^{2}\triangle\psi\right)\,,$
(7) $\displaystyle\frac{\sigma_{\zeta}^{2}}{\sigma_{0}^{2}}$
$\displaystyle=1-\frac{1}{1-\gamma^{2}}\psi^{2}-\frac{1}{\gamma^{2}(1-\gamma^{2})}\left(2\gamma^{2}\psi+\frac{1}{3}R_{*}^{2}\triangle\psi\right)\frac{1}{3}R_{*}^{2}\triangle\psi$
$\displaystyle\quad-\frac{5}{\gamma^{2}}\left(\frac{\psi^{\prime}}{r}-\frac{1}{3}\triangle\psi\right)^{2}R_{*}^{2}-\frac{1}{\gamma^{2}}R_{*}^{2}{\psi^{\prime}}^{2}\,,$
(8)

with $\gamma=\sigma_{1}^{2}/(\sigma_{0}\sigma_{2})$,
$R_{*}=\sqrt{3}\sigma_{1}/\sigma_{2}$ and

$\psi(r)=\frac{1}{\sigma_{0}^{2}}\int\frac{dk}{k}\frac{\sin(kr)}{kr}{\cal
P}(k)\,.$ (9)

In the above, $\sigma_{n}$ is a $n$-th order gradient moment for $\zeta_{\rm
G}$ defined by

$\sigma_{n}^{2}\equiv\int\frac{dk}{k}k^{2n}{\cal P}(k)\,.$ (10)

As in (2), here, the non-Gaussian curvature perturbation $\zeta$ is explicitly
given as a function of $\zeta_{\rm G}$. Then, by using (4), we can derive the
typical profile $\bar{\zeta}(r;\mu,k_{*})$ for $\zeta(r)$ as

$\displaystyle-\bar{\zeta}(r;\mu,k_{*})$ $\displaystyle=-\int\zeta
P(\zeta_{\rm G}|\mu,k_{*})d\zeta_{\rm G}$ (11)
$\displaystyle=-\int\left[\zeta_{\rm G}-\frac{3}{5}f_{\rm NL}\left(\zeta_{\rm
G}^{2}-\left\langle\zeta_{\rm G}^{2}\right\rangle\right)\right]P(\zeta_{\rm
G}|\mu,k_{*})d\zeta_{\rm G}$ (12) $\displaystyle=-\bar{\zeta}_{\rm
G}+\frac{3}{5}f_{\rm NL}\left({\bar{\zeta}_{\rm
G}}^{2}+\sigma_{\zeta}^{2}-\sigma_{0}^{2}\right)\,.$ (13)

##  3 General expression for the PBH fraction

Hereafter the calculations are mostly implicit and can be numerically done in
general, and the basic procedure of the calculations are shown in [33]. Thus,
we just give a sketch of the implicit calculations with some details being
given in Appendix B).

We consider that the criterion for PBH formation is given in terms of the so-
called compaction function ${\cal C}$, which is defined as the mass excess at
a certain region. As shown Appendix B, ${\cal C}$ can be written as a function
of $\zeta$. PBH would be formed when the maximum value of the compaction
function, ${\cal C}^{\rm max}$, is larger than a threshold ${\cal C}_{\rm
th}$. Evaluating (B.2) with respect to the typical profile $\bar{\zeta}$, we
may obtain the value of $r_{\rm m}$, at which the compaction function takes
the maximum value ${\mathcal{C}}^{\rm max}$, as a function of $\mu$ and
$k_{*}$. That is, $r_{\rm m}=\bar{r}_{\rm m}(\mu,k_{*})$ where an overbar
means it is evaluated by using $\bar{\zeta}$, which is a function of $\mu$ and
$k_{*}$. Then, ${\cal C}^{\rm max}$ can be also expressed as a function of
$\mu$ and $k_{*}$ as ${\cal C}^{\rm max}=\bar{\cal C}^{\rm max}(\mu,k_{*})$.
Note that $\bar{r}_{\rm m}$ and $\bar{\cal C}^{\rm max}$ depend on the non-
linear parameter $f_{\rm NL}$ through the non-Gaussian correction in the
typical profile $\bar{\zeta}$ as shown in (13). By identifying $\bar{\cal
C}^{\rm max}(\mu,k_{*})$ as the threshold value ${\cal C}_{\rm th}\approx
0.267$, we define the threshold $\mu_{\rm th}$ as a function of $k_{*}$:

$\mu_{\rm th}=\mu_{\rm th}(k_{*})\,.$ (14)

We may also regard the PBH mass $M$, which is defined by $r_{\rm m}$ [see
(B.4) in Appendix B], as a function of $\mu$ and $k_{*}$ by identifying
$r_{\rm m}=\bar{r}_{\rm m}(\mu,k_{*})$ and $\zeta=\bar{\zeta}(\bar{r}_{\rm
m};\mu,k_{*})$:

$M=\bar{M}(\mu,k_{*})\,.$ (15)

Then, eliminating $k_{*}$ from the above two equations, the threshold value of
$\mu_{\rm th}$ is given as a function of the PBH mass $M$:

$\mu_{\rm th}=\mu_{\rm th}(M)\,,$ (16)

and this would be modified by the non-Gaussian effect. As will be shown later,
in particular, for the case with the extended power spectrum we may have a
minimum value of $\mu$ for a fixed value of $M$. Let the function $\mu_{\rm
min}(M)$ denote this minimum value. Then, the relevant region of $\mu$ for PBH
formation with the mass $M$ is given by

$\mu>\mu_{\rm b}\equiv\max\left\\{\mu_{\rm min}(M),\mu_{\rm
th}^{(M)}(M)\right\\}\,.$ (17)

Since the parameters $\mu$ and $k_{*}$ are based on the Gaussian random
variable, we can use the standard expression for the peak number density:

$n_{\rm pk}(\mu,k_{*})d\mu dk_{*}=\frac{2\cdot 3^{3/2}}{(2\pi)^{3/2}}\mu
k_{*}\frac{\sigma_{2}^{2}}{\sigma_{0}\sigma_{1}^{3}}f\left(\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)P_{1}\left(\frac{\mu}{\sigma_{0}},\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)d\mu dk_{*}\,,$ (18)

where the function $f$ and the probability $P_{1}$ are given by

$\displaystyle f(x)$ $\displaystyle=\frac{1}{2}x(x^{2}-3)\left[{\rm
erf}\left(\frac{1}{2}\sqrt{\frac{5}{2}}x\right)+{\rm
erf}\left(\sqrt{\frac{5}{2}}x\right)\right]$ (19)
$\displaystyle\quad+\sqrt{\frac{2}{5\pi}}\left[\left(\frac{8}{5}+\frac{31}{4}x^{2}\right)e^{-5x^{2}/8}+\left(-\frac{8}{5}+\frac{1}{2}x^{2}\right)e^{-5x^{2}/2}\right]\,,$
(20) $\displaystyle P_{1}\left(\frac{\mu}{\sigma_{0}},\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)d\mu dk_{*}$ $\displaystyle=\frac{\mu
k_{*}}{\pi\sigma_{0}\sigma_{2}\sqrt{1-\gamma^{2}}}\exp\left\\{-\frac{\mu^{2}}{2}\left[\frac{1}{\sigma_{0}^{2}}+\frac{1}{\sigma_{2}^{2}(1-\gamma^{2})}\left(k_{*}^{2}-\gamma\frac{\sigma_{2}}{\sigma_{1}}\right)^{2}\right]\right\\}d\mu
dk_{*}\,.$ (21)

Further, since the direct observable is not $k_{*}$ but the PBH mass $M$, we
change the variable from $k_{*}$ to $M$ as follows:

$\displaystyle n_{\rm pk}(\mu,M)d\mu dM$
$\displaystyle=\frac{3^{3/2}}{(2\pi)^{3/2}}\frac{\sigma_{2}^{2}}{\sigma_{0}\sigma_{1}^{3}}\mu
k_{*}f\left(\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)P_{1}\left(\frac{\mu}{\sigma_{0}},\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)$ (22)
$\displaystyle\quad\times\left|\left[\frac{1}{\bar{r}_{\rm
m}}-\bar{\zeta}^{\prime}(\bar{r}_{\rm m})\right]\frac{d\bar{r}_{\rm
m}}{dk_{*}}-\left(\frac{d\bar{\zeta}}{dk_{*}}\right)_{r=\bar{r}_{\rm
m}}\right|^{-1}d\mu d\log M\,,$ (23)

where $k_{*}$ should be regarded as a function of $\mu$ and $M$ as in (15). By
counting the number of peaks whose $\mu$ is larger than the threshold value
$\mu_{\rm b}$, we can evaluate the number density of PBHs as

$n_{\rm PBH}d\log M=\left[\int^{\infty}_{\mu_{\rm b}}d\mu~{}n_{\rm
pk}(\mu,M)\right]Md\log M\,.$ (24)

We also note that the scale factor $a$ is a function of $M$ as shown in (B.4).
Then, the fraction of PBHs to the total density $\beta_{0}d\log M$ can be
given by

$\displaystyle\beta_{0}d\log M$ $\displaystyle=\frac{Mn_{\rm PBH}}{\rho
a^{3}}d\log M$ (25) $\displaystyle=\frac{4\pi}{3}n_{\rm PBH}k_{\rm
eq}^{-3}\left(\frac{M}{M_{\rm eq}}\right)^{3/2}d\log M$ (26)
$\displaystyle=\frac{2\cdot 3^{1/2}k_{\rm
eq}^{-3}}{(2\pi)^{1/2}}\frac{\sigma_{2}^{2}}{\sigma_{0}\sigma_{1}^{3}}\left(\frac{M}{M_{\rm
eq}}\right)^{3/2}\Bigg{\\{}\int^{\infty}_{\mu_{\rm b}}d\mu\mu
k_{*}f\left(\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)P_{1}\left(\frac{\mu}{\sigma_{0}},\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)$ (27) $\displaystyle\hskip
110.00017pt\times\left|\left[\frac{1}{\bar{r}_{\rm
m}}-\bar{\zeta}^{\prime}(\bar{r}_{\rm m})\right]\frac{d\bar{r}_{\rm
m}}{dk_{*}}-\left(\frac{d\bar{\zeta}}{dk_{*}}\right)_{r=\bar{r}_{\rm
m}}\right|^{-1}\Bigg{\\}}d\log M\,,$ (28)

where $M_{\rm eq}$ and $k_{\rm eq}$ are respectively corresponding to the
formed PBH mass and the comoving wave number reentering the horizon, at the
matter radiation equality.

##  4 Effect of local non-Gaussianity on PBH abundance

Given the general expression for the PBH abundance (28), now we employ
explicit examples with specific power spectra and investigate the effect of
the local non-Gaussianity on the abundance of PBHs.

###  4.1 Monochromatic power spectrum

Let us consider the monochromatic power spectrum given by

${\cal P}(k)=\sigma_{0}^{2}k_{0}\delta(k-k_{0})\,.$ (29)

For this case, gradient moments are explicitly calculated as
$\sigma_{n}^{2}=\sigma_{0}^{2}k_{0}^{2n}$. Then, replacing $k_{*}$ by $k_{0}$
in (7), we find

$-\bar{\zeta}_{\rm G}=\mu\psi(r)=\mu\frac{\sin(k_{0}r)}{k_{0}r}\,,$ (30)

where $\psi$ is calculated from (9). Since the value of $k_{*}$ is fixed to be
$k_{*}=k_{0}$, we do not have $k_{*}$ dependence for $\bar{r}_{\rm m}$,
$\bar{\cal C}^{\rm max}$, $\mu_{\rm th}$ and $\bar{M}$. The value of $\mu_{\rm
th}$ is converted into the threshold value of $M_{\rm th}$ through the
function $\bar{M}(\mu)$, i.e. $M_{\rm th}=\bar{M}(\mu_{\rm th})$. Through the
non-Gaussian correction in the typical profile $\bar{\zeta}$, as mentioned in
the previous section, $\bar{r}_{\rm m}$, $\bar{\cal C}^{\rm max}$ and also
$\bar{M}$ are dependent on $f_{\rm NL}$. Figure 1 shows the non-Gaussian
effects on $\bar{r}_{\rm m}$, $\bar{\cal C}^{\rm max}$ and $\bar{M}$ as
functions of the peak amplitude $\mu$. Here, we take $\sigma_{0}=0.06$ and in
each panel the blue, red, and green lines are respectively for the Gaussian
case, i.e. $f_{\rm NL}=0$, $3f_{\rm NL}/5=1$ and $3f_{\rm NL}/5=-1/4$. In the
panel for $\bar{\cal C}^{\rm max}$, the yellow horizontal line represents the
threshold for PBH formation. The dashed lines in the panel for $\bar{M}$ are
also the thresholds with respect to $M$.

Figure 1: (Left panel) $\bar{r}_{\rm m}(\mu)$, (middle panel) $\bar{\cal
C}^{\rm max}(\mu)$ and (right panel) $\bar{M}(\mu)$ for $\sigma_{0}=0.06$.
Here, $F_{\rm NL}\equiv 3f_{\rm NL}/5$ for notational simplicity.

Since the probability distribution $P_{1}$ is given by

$P_{1}\left(\frac{\mu}{\sigma_{0}},\frac{\mu
k_{*}^{2}}{\sigma_{2}}\right)\underset{\gamma\rightarrow
1}{\longrightarrow}\frac{1}{2\sqrt{2\pi}}\frac{\sigma_{2}}{\mu
k_{*}}\delta(k_{*}-k_{0})\exp\left(-\frac{\mu^{2}}{2\sigma_{0}^{2}}\right)\,,$
(31)

the peak number density can be written as

$n_{\rm
pk}(\mu)d\mu=\frac{3^{3/2}}{(2\pi)^{2}}\frac{1}{\sigma_{0}}k_{0}^{3}f\left(\frac{\mu}{\sigma_{0}}\right)\exp\left(-\frac{\mu^{2}}{2\sigma_{0}^{2}}\right)d\mu\,.$
(32)

Changing the variable from $\mu$ to $M$ through the function $M=\bar{M}(\mu)$,
and taking the threshold value $M_{\rm th}$ into account, we obtain the
following expression for the PBH number density:

$n_{\rm PBH}d\log
M=\frac{3^{3/2}}{(2\pi)^{2}}\frac{1}{\sigma_{0}}k_{0}^{3}f\left(\frac{\mu}{\sigma_{0}}\right)\exp\left(-\frac{\mu^{2}}{2\sigma_{0}^{2}}\right)M\frac{d\mu}{dM}\Theta(M-M_{\rm
th})dM\,,$ (33)

where $\mu$ should be regarded as a function of $M$ through the relation
$M=\bar{M}(\mu)$. Finally, we obtain

$\beta_{0}d\log M=\frac{3^{1/2}}{\pi}\frac{1}{\sigma_{0}}\left(\frac{M}{M_{\rm
eq}}\right)^{3/2}\left(\frac{k_{0}}{k_{\rm
eq}}\right)^{3}f\left(\frac{\mu}{\sigma_{0}}\right)\exp\left(-\frac{\mu^{2}}{2\sigma_{0}^{2}}\right)M\frac{d\mu}{dM}\Theta(M-M_{\rm
th})dM\,.$ (34)

Note that the functional form of this expression is completely the same as
(80) in [33]. The non-Gaussian correction in $\beta_{0}$ appears just through
the modifications of $M_{\rm th}$ and the value of $\mu$ at $M_{\rm th}$. The
PBH fraction $\beta_{0}$ for $3f_{\rm NL}/5=-1/4$, $0$ and $1$ is depicted as
a function of $M$ with $\sigma_{0}=0.06$ and $k_{0}=10^{5}k_{\rm eq}$ in
Figure 2. As can be seen, for larger $f_{\rm NL}$ the cut-off scale for PBH
mass is shifted to smaller value. This is due to the smaller $M_{\rm th}$ for
larger $f_{\rm NL}$ as shown in the right panel in Figure 1. Furthermore, the
peak value of the PBH fraction $\beta_{0}$ increases with increasing the value
of $f_{\rm NL}$. This is due to the fact that the value of $\mu$ at $M=M_{\rm
th}$ becomes smaller for larger value of $f_{\rm NL}$, as can be seen in the
right panel in Figure 1.

Figure 2:  $\beta_{0}$ as a function of $M$ with $\sigma_{0}=0.06$ and
$k_{0}=10^{5}k_{\rm eq}$. The PBH fraction increases with increasing the value
of $f_{\rm NL}$.

###  4.2 Extended power spectrum

Let us consider the simple extended power spectrum given by

${\cal
P}(k)=3\sqrt{\frac{6}{\pi}}\sigma^{2}\left(\frac{k}{k_{0}}\right)^{3}\exp\left(-\frac{3}{2}\frac{k^{2}}{k_{0}^{2}}\right)\,.$
(35)

For this case, the functional form of $\psi(r)$ is given by

$\psi(r)=\exp\left(-\frac{k_{0}^{2}r^{2}}{6}\right)\,.$ (36)

The dependence of $r_{\rm m}(\mu,k_{*})$, $\bar{\cal C}^{\rm max}(\mu,k_{*})$
and $\bar{M}(\mu,k_{*})$ on $k_{*}$ and $\mu$ is shown in Figure 3 for
$\sigma=0.1$ and $3f_{\rm NL}/5=1$.

Figure 3: (Left panel) $r_{\rm m}(\mu,k_{*})$, (middle panel) $\bar{\cal
C}^{\rm max}(\mu,k_{*})$, ${\cal C}_{\rm th}$ and $\mu^{(k_{*})}_{\rm th}$,
and (right panel) $\bar{M}(\mu,k_{*})$ for $\sigma=0.1$ and $3f_{\rm NL}/5=1$.

As we have mentioned, for the extended power spectrum, we find that the
maximum value of $M$ is realized at $k_{*}=0$ for each $\mu$. Therefore, for a
given value of $M$, the minimum value of $\mu_{\rm min}$ is given by

$\mu_{\rm min}(M)=\mu(M,0)\,.$ (37)

Thus, in the estimation of the PBH fraction, we should take into account the
relevant region of $\mu$ with the mass $M$ given by (17). The dependence of
$\mu_{b}(M)$ on different values of $f_{\rm NL}$ is shown in the left panel in
Figure 4. The resultant PBH fraction $\beta_{0}$ for the extended power
spectrum is depicted as a function of $M$ in the right panel of Figure 4.

Figure 4:  (Left panel) $\mu_{\rm th}(M)$ and $\mu_{\rm min}(M)$ for
$\sigma=0.1$, and (right panel) $\beta_{0}$ as a function of $M$. The PBH
abundance increases with a larger value of $f_{\rm NL}$.

##  5 Summary

We have discussed the effect of the local type non-Gaussianity on the PBH
abundance. Our procedure is based on the peak theory for the Gaussian variable
$\zeta_{\rm G}$, which is related to the non-Gaussian curvature perturbation
$\zeta$ via (2), which is the simplest type of non-Gaussianity with the non-
linear parameter $f_{\rm NL}$. Summary of our procedure is schematically shown
in Figure 5. The value of $|f_{\rm NL}|\sim 1$ induces a similar effect to a
few factors of difference in the amplitude of the power spectrum. The
negative/positive value of $|f_{\rm NL}|$ of ${\cal O}(1)$ may reduce/enhance
the PBH abundance in several orders of magnitude. The calculation procedure
for estimating the PBH abundance shown in this article can be also applied to
more general types of primordial non-Gaussianity straightforwardly, once the
primordial curvature perturbation is given as a function of the Gaussian
variable $\zeta_{\rm G}$.

Figure 5: A flow chart of our procedure.

Note added: While preparing the manuscript, we found a paper by Atal, Garriga
and Marcos-Caballero [60], which has some overlap with our work.

### Acknowledgements

This work is supported in part by the Japan Society for the Promotion of
Science (JSPS) and the National Research Foundation of Korea (NRF) under the
Japan-Korea Basic Scientific Cooperation Program (NRF-2018K2A9A2A08000127). CY
is supported in part by JSPS KAKENHI Grant (JP19H01895). JG is supported in
part by the Basic Science Research Program through NRF Research Grant
(2016R1D1A1B03930408). SY is supported by MEXT KAKENHI Grant Number 15H05888
and 18H04356.

##  Appendix A Peaks of $\zeta_{\rm G}$ and $\zeta$

First, let us introduce the following Taylor expansion of the fields
$\zeta(x_{i})$ an $\zeta_{\rm G}(x_{i})$:

$\displaystyle\zeta$
$\displaystyle=\zeta_{0}+\zeta^{i}_{1}x_{i}+\frac{1}{2}\zeta_{2}^{ij}x_{i}x_{j}+{\cal
O}(x^{3})\,,$ (A.1) $\displaystyle\zeta$ $\displaystyle=\zeta_{\rm
G0}+\zeta_{\rm G1}^{i}x_{i}+\frac{1}{2}\zeta_{\rm G2}^{ij}x_{i}x_{j}+{\cal
O}(x^{3})\,.$ (A.2)

With (2), the coefficients satisfy the following relations:

$\displaystyle\zeta_{0}$ $\displaystyle=\zeta_{\rm G0}-\frac{3}{5}f_{\rm
NL}\left(\zeta_{\rm G0}^{2}-\sigma_{0}^{2}\right)\,,$ (A.3)
$\displaystyle\zeta_{1}^{i}$ $\displaystyle=\zeta_{\rm
G1}^{i}\left(1-\frac{6}{5}f_{\rm NL}\zeta_{\rm G0}\right)\,,$ (A.4)
$\displaystyle\zeta_{2}$ $\displaystyle=\zeta_{\rm
G2}\left(1-\frac{6}{5}f_{\rm NL}\zeta_{\rm G0}\right)-\frac{6}{5}f_{\rm
NL}\sum_{i}\left(\zeta_{\rm G1}^{i}\right)^{2}\,,$ (A.5)

where $\zeta_{2}=\zeta_{2}^{11}+\zeta_{2}^{22}+\zeta_{2}^{33}$ and $\zeta_{\rm
G2}=\zeta_{\rm G2}^{11}+\zeta_{\rm G2}^{22}+\zeta_{\rm G2}^{33}$. Here, we
adopt a simple criterion for PBH formation $\delta>\delta_{\rm th}$ at the
horizon entry to perform the order-of-magnitude estimation.

The value of $\delta$ can be expressed in terms of the Taylor expansion around
a certain point:

$\delta=\frac{4}{9}\frac{1}{a^{2}H^{2}}\left\\{e^{2\zeta_{0}}\left[\zeta_{2}-\frac{1}{2}\sum_{i}\left(\zeta_{1}^{i}\right)^{2}\right]\right\\}+{\cal
O}(y^{2})\,,$ (A.6)

where we have introduced a new spatial coordinate $\bm{y}$ to emphasize the
difference from the expansion around the peak of $\zeta$. Focusing on the
region which will collapse into a PBH, from the inequality $\delta>\delta_{\rm
th}$, we obtain the following inequality:

$\left.\frac{\zeta_{2}}{\bar{a}^{2}H^{2}}\right|_{\text{horizon
entry}}>\frac{9}{4}\delta_{\rm th}\sim 1\,,$ (A.7)

where we have defined the renormalized local scale factor
$\bar{a}=ae^{-\zeta_{0}}$.

Let us apply the inequality (A.7) to a peak of $\zeta$. We introduce the
following two variables to characterize the profile around a peak of $\zeta$:

$\displaystyle\mu_{\rm NG}$ $\displaystyle=-\zeta|_{\rm peak}\,,$ (A.8)
$\displaystyle\kappa_{\rm NG}$ $\displaystyle=\frac{\triangle\zeta|_{\rm
peak}}{\mu_{\rm NG}}\,.$ (A.9)

We note that, while we introduce $k_{*}^{2}$ for $\zeta_{\rm G}$ in the main
text because we are interested in the case $\triangle\zeta_{\rm G}>0$,
$\kappa_{\rm NG}$ can be negative in general. Adopting
$\bar{a}^{2}H^{2}\sim|\kappa|$ as the horizon entry condition, (A.7) can be
rewritten as

$\mu_{\rm NG}>\frac{9}{4}\delta_{\rm th}\sim 1\,.$ (A.10)

(A.3), (A.4) and (A.5) are then written as

$\displaystyle\mu_{\rm NG}$ $\displaystyle=\mu+\frac{3}{5}f_{\rm
NL}\left(\mu^{2}-\sigma_{0}^{2}\right)\,,$ (A.11) $\displaystyle 0$
$\displaystyle=\zeta_{\rm G1}^{i}\left(1+\frac{6}{5}f_{\rm NL}\mu\right)\,,$
(A.12) $\displaystyle\mu_{\rm NG}\kappa_{\rm NG}$
$\displaystyle=\mu\kappa\left(1+\frac{6}{5}f_{\rm
NL}\mu\right)+\frac{6}{5}f_{\rm NL}\sum_{i}\left(\zeta_{\rm
G1}^{i}\right)^{2}\,.$ (A.13)

If $1+6f_{\rm NL}\mu/5$ is negative, the sign of $\mu_{\rm NG}\kappa_{\rm NG}$
can be different from $\mu\kappa$. That is, a peak of $\zeta$ can be a minimum
of $\zeta_{\rm G}$. In order to avoid this situation, we simply assume
$3f_{\rm NL}/5\gtrsim-1/4$.

##  Appendix B PBH formation criteria

Here, we briefly review the estimate on the threshold for the formation of
PBHs following [31], introducing the basic perturbation variables used in the
main text.

Assuming the spherical symmetry around high peaks, we introduce a compaction
function, which represents an excess of the Misner-Sharp mass $\delta{M}$ in
the spherical region with the radius $r$, given by (see [31, 33] for details)

${\cal
C}(r)\equiv\frac{\delta{M}}{R}=\frac{1}{3}\left[1-\left(1-r\zeta^{\prime}\right)^{2}\right]\,,$
(B.1)

where $R\equiv are^{-\zeta}$ is the areal radius at the radius $r$, and the
second equality holds on the comoving slice. Assuming that ${\cal C}$ is a
smooth function of $r$, we may find a radius $r=r_{\rm m}$ at which ${\cal C}$
takes the maximum value ${\cal C}^{\rm max}={\cal C}(r_{\rm m})$, and ${\cal
C}^{\prime}(r_{\rm m})=0$ should be satisfied. From (B.1), ${\cal
C}^{\prime}(r_{\rm m})=0$ gives the condition for $\zeta$ at $r=r_{\rm m}$ as

$\left.\left(\zeta^{\prime}+r\zeta^{\prime\prime}\right)\right|_{r=r_{\rm
m}}=0\,.$ (B.2)

Then, as a criterion for PBH formation, we consider that PBH is formed when
the maximum value of the compaction function is larger than a threshold ${\cal
C}_{\rm th}$:

${\cal C}^{\rm max}>{\cal C}_{\rm th}\,.$ (B.3)

In this article, we use ${\cal C}_{\rm th}\approx 0.267$ as a reference value
[31, 33]. The specific value ${\cal C}\approx 0.267$ is imported from the
results of the numerical simulations for PBH formation performed in [31]. We
note that, as is shown in [31], at the moment of the horizon entry
$are^{-\zeta}=1/H$, the value of the compaction function is equivalent to the
half of the volume averaged density perturbation. Therefore the inequality
(B.3) can be expressed as $\bar{\delta}_{\rm max}>\delta_{\rm th}=2{\cal
C}_{\rm th}$ with $\bar{\delta}_{\rm max}$ being the maximum value of the
averaged density perturbation.

Finally, let us give the estimate of the PBH mass. Since the PBH mass is given
by $M=\alpha H^{-1}/2$ with $\alpha$ being a numerical factor, from the
horizon entry condition $are^{-\zeta}=1/H$, the PBH mass $M$ can be expressed
as follows:

$M=\frac{1}{2}\alpha H^{-1}=\frac{1}{2}\alpha ar_{\rm m}e^{-\zeta(r_{\rm
m})}=M_{\rm eq}k_{\rm eq}^{2}r_{\rm m}^{2}e^{-2\zeta(r_{\rm m})}\,,$ (B.4)

where we have used $H\propto a^{-2}$ during radiation dominated era, and
$M_{\rm eq}=\alpha H_{\rm eq}^{-1}/2$ and $k_{\rm eq}=a_{\rm eq}H_{\rm eq}$ at
the matter-radiation equality. The value of the numerical factor $\alpha$ is
rather ambiguous, and we set $\alpha=1$ as a fiducial value in the main text.

## References

  * [1] Y. B. Zel’dovich and I. D. Novikov, Soviet Ast. 19, 602 (1967).
  * [2] S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * [3] G. F. Chapline, Nature 253, 251 (1975).
  * [4] B. J. Carr, Astrophys. J. 201, 1 (1975).
  * [5] J. Garcia-Bellido, A. D. Linde and D. Wands, Phys. Rev. D 54, 6040 (1996) [astro-ph/9605094].
  * [6] K. Jedamzik and J. C. Niemeyer, Phys. Rev. D 59, 124014 (1999) [astro-ph/9901293].
  * [7] P. H. Frampton, M. Kawasaki, F. Takahashi and T. T. Yanagida, JCAP 1004, 023 (2010) [arXiv:1001.2308 [hep-ph]].
  * [8] M. Kawasaki, N. Kitajima and T. T. Yanagida, Phys. Rev. D 87, no. 6, 063519 (2013) [arXiv:1207.2550 [hep-ph]].
  * [9] K. Kohri, C. M. Lin and T. Matsuda, Phys. Rev. D 87, no. 10, 103527 (2013) [arXiv:1211.2371 [hep-ph]].
  * [10] B. Carr, F. Kuhnel and M. Sandstad, Phys. Rev. D 94, no. 8, 083504 (2016) [arXiv:1607.06077 [astro-ph.CO]].
  * [11] G. Chapline and P. H. Frampton, JCAP 1611, no. 11, 042 (2016) [arXiv:1608.04297 [gr-qc]].
  * [12] J. M. Ezquiaga, J. Garcia-Bellido and E. Ruiz Morales, Phys. Lett. B 776, 345 (2018) [arXiv:1705.04861 [astro-ph.CO]].
  * [13] S. Clesse and J. Garcia-Bellido, Phys. Dark Univ. 22, 137 (2018) [arXiv:1711.10458 [astro-ph.CO]].
  * [14] K. Kohri and T. Terada, Class. Quant. Grav. 35, no. 23, 235017 (2018) [arXiv:1802.06785 [astro-ph.CO]].
  * [15] B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, Phys. Rev. D 81, 104019 (2010) [arXiv:0912.5297 [astro-ph.CO]].
  * [16] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. 35, no. 6, 063001 (2018) [arXiv:1801.05235 [astro-ph.CO]].
  * [17] W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974).
  * [18] D. K. Nadezhin, I. D. Novikov and A. G. Polnarev, Soviet Ast. 22, 129 (1978).
  * [19] I. D. Novikov and A. G. Polnarev, Soviet Ast. 24, 147 (1980).
  * [20] J. C. Niemeyer and K. Jedamzik, Phys. Rev. D 59, 124013 (1999) [astro-ph/9901292].
  * [21] M. Shibata and M. Sasaki, Phys. Rev. D 60, 084002 (1999) [gr-qc/9905064].
  * [22] I. Musco, J. C. Miller and L. Rezzolla, Class. Quant. Grav. 22, 1405 (2005) [gr-qc/0412063].
  * [23] A. G. Polnarev and I. Musco, Class. Quant. Grav. 24, 1405 (2007) [gr-qc/0605122].
  * [24] I. Musco, J. C. Miller and A. G. Polnarev, Class. Quant. Grav. 26, 235001 (2009) [arXiv:0811.1452 [gr-qc]].
  * [25] I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013) [arXiv:1201.2379 [gr-qc]].
  * [26] A. G. Polnarev, T. Nakama and J. Yokoyama, JCAP 1209, 027 (2012) [arXiv:1204.6601 [gr-qc]].
  * [27] T. Harada, C. M. Yoo and K. Kohri, Phys. Rev. D 88, no. 8, 084051 (2013) Erratum: [Phys. Rev. D 89, no. 2, 029903 (2014)] [arXiv:1309.4201 [astro-ph.CO]].
  * [28] T. Nakama, T. Harada, A. G. Polnarev and J. Yokoyama, JCAP 1401, 037 (2014) [arXiv:1310.3007 [gr-qc]].
  * [29] S. Young, C. T. Byrnes and M. Sasaki, JCAP 1407, 045 (2014) [arXiv:1405.7023 [gr-qc]].
  * [30] T. Nakama, JCAP 1410, no. 10, 040 (2014) [arXiv:1408.0955 [gr-qc]].
  * [31] T. Harada, C. M. Yoo, T. Nakama and Y. Koga, Phys. Rev. D 91, no. 8, 084057 (2015) [arXiv:1503.03934 [gr-qc]].
  * [32] I. Musco, arXiv:1809.02127 [gr-qc].
  * [33] C. M. Yoo, T. Harada, J. Garriga and K. Kohri, PTEP 2018, no. 12, 123E01 (2018) [arXiv:1805.03946 [astro-ph.CO]].
  * [34] A. G. Doroshkevich, Astrofizika 6, 581 (1970).
  * [35] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, Astrophys. J. 304, 15 (1986).
  * [36] M. Kawasaki and H. Nakatsuka, Phys. Rev. D 99, no. 12, 123501 (2019) [arXiv:1903.02994 [astro-ph.CO]].
  * [37] V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto and C. Unal, arXiv:1904.00970 [astro-ph.CO].
  * [38] S. Young, I. Musco and C. T. Byrnes, arXiv:1904.00984 [astro-ph.CO].
  * [39] J. S. Bullock and J. R. Primack, Phys. Rev. D 55, 7423 (1997) [astro-ph/9611106].
  * [40] P. Ivanov, Phys. Rev. D 57, 7145 (1998) [astro-ph/9708224].
  * [41] J. Yokoyama, Phys. Rev. D 58, 083510 (1998) [astro-ph/9802357].
  * [42] J. C. Hidalgo, arXiv:0708.3875 [astro-ph].
  * [43] R. Saito, J. Yokoyama and R. Nagata, JCAP 0806, 024 (2008) [arXiv:0804.3470 [astro-ph]].
  * [44] E. Bugaev and P. Klimai, Phys. Rev. D 85, 103504 (2012) [arXiv:1112.5601 [astro-ph.CO]].
  * [45] C. T. Byrnes, E. J. Copeland and A. M. Green, Phys. Rev. D 86, 043512 (2012) [arXiv:1206.4188 [astro-ph.CO]].
  * [46] A. Linde, S. Mooij and E. Pajer, Phys. Rev. D 87, no. 10, 103506 (2013) [arXiv:1212.1693 [hep-th]].
  * [47] E. V. Bugaev and P. A. Klimai, Int. J. Mod. Phys. D 22, 1350034 (2013) [arXiv:1303.3146 [astro-ph.CO]].
  * [48] S. Young and C. T. Byrnes, JCAP 1308, 052 (2013) [arXiv:1307.4995 [astro-ph.CO]].
  * [49] Y. Tada and S. Yokoyama, Phys. Rev. D 91, no. 12, 123534 (2015) [arXiv:1502.01124 [astro-ph.CO]].
  * [50] S. Young and C. T. Byrnes, JCAP 1504, no. 04, 034 (2015) [arXiv:1503.01505 [astro-ph.CO]].
  * [51] S. Young, D. Regan and C. T. Byrnes, JCAP 1602, no. 02, 029 (2016) [arXiv:1512.07224 [astro-ph.CO]].
  * [52] K. Ando, K. Inomata, M. Kawasaki, K. Mukaida and T. T. Yanagida, Phys. Rev. D 97, no. 12, 123512 (2018) [arXiv:1711.08956 [astro-ph.CO]].
  * [53] G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, JCAP 1803, no. 03, 016 (2018) [arXiv:1801.09415 [astro-ph.CO]].
  * [54] J. M. Ezquiaga and J. García-Bellido, JCAP 1808, 018 (2018) [arXiv:1805.06731 [astro-ph.CO]].
  * [55] K. Ando, M. Kawasaki and H. Nakatsuka, Phys. Rev. D 98, no. 8, 083508 (2018) [arXiv:1805.07757 [astro-ph.CO]].
  * [56] V. Atal and C. Germani, Phys. Dark Univ. , 100275 [arXiv:1811.07857 [astro-ph.CO]].
  * [57] S. Passaglia, W. Hu and H. Motohashi, Phys. Rev. D 99, no. 4, 043536 (2019) [arXiv:1812.08243 [astro-ph.CO]].
  * [58] D. H. Lyth, K. A. Malik and M. Sasaki, JCAP 0505, 004 (2005) [astro-ph/0411220].
  * [59] Y. Akrami et al. [Planck Collaboration], arXiv:1905.05697 [astro-ph.CO].
  * [60] V. Atal, J. Garriga and A. Marcos-Caballero, arXiv:1905.13202 [astro-ph.CO].