# Non-Gaussianity effects on the primordial black hole abundance for sharply-
peaked primordial spectrum

Takahiko Matsubara  and Misao Sasaki

###### Abstract

We perturbatively study the effect of non-Gaussianities on the mass fraction
of primordial black holes (PBHs) at the time of formation by systematically
taking its effect into account in the one-point probability distribution
function of the primordial curvature perturbation. We focus on the bispectrum
and trispectrum and derive formulas that describe their effects on the
skewness and kurtosis of the distribution function. Then considering the case
of narrowly peaked spectra, we obtain simple formulas that concisely express
the effect of the bi- and trispectra. In particular, together with the $g_{\rm
NL}$ and $\tau_{\rm NL}$ parameters of the trispectrum, we find that non-
Gaussianity parameters for various types of the bispectrum are linearly
combined to give an effective parameter, $f_{\rm NL}^{\rm eff}$, that
determines the PBH mass fraction in the narrow spectral shape limit.

##  1 Introduction

The possibility that black holes may be formed in the very early universe was
suggested about half a century ago [1, 2, 3]. Since then it has been a topic
of constant interest in cosmology, but it has never been explored in depth.
However, thanks to the rapidly growing interest in gravitational wave
cosmology in recent years, as well as to the theoretical progress and
technical developments, the primordial black holes (PBHs) have become one of
the hottest topics in cosmology today (for recent reviews, see e.g., [4, 5]).

One of the most studied PBH formation mechanisms is the collapse of a region
with a sufficiently large curvature perturbation, presumably produced from
inflation, during the radiation-dominated early universe. While such high
peaks in the curvature perturbation hardly exist in the conventional models of
inflation, since there are virtually no stringent observational constraint on
cosmologically very small scales that are relevant for the PBH formation,
various models that can produce sufficiently abundant PBHs have been discussed
in the literature.

The abundance of PBHs produced in the early Universe plays a crucial role in
the present Universe. Depending on the typical mass of PBHs, various effects
which are observationally detectable are anticipated if the abundance is
sufficiently large. There is no compelling evidence for the existence of PBHs
at any mass scale up to the present time. Hence, what we have so far is the
upper limit in the PBH abundance in broad ranges of mass scale (For a recent
review, see Ref. [6]). If PBHs are formed from high peaks of the curvature
perturbation, an upper bound on the PBH abundance places an upper bound on the
amplitude of the primordial power spectrum. The observed amplitude of the
power spectrum on scales of the cosmic microwave background (CMB) and the
large-scale structure (LSS) is too small to produce PBHs. However, scales
relevant to the PBH formation are typically much smaller than those scales.
Thus, various mechanisms that give rise to the amplitude of the power spectrum
on small scales much larger than that on the CMB scale have been proposed
(see, e.g., Ref. [7]). In many of those models, however, the non-Gaussianities
in the curvature perturbation are negligible or simply ignored. Nevertheless,
even if they are small, they may considerably affect the PBH formation, as it
is acutely sensitive to non-Gaussian features in the tail of the probability
distribution function [8].

The effects of non-Gaussianity upon the abundance of PBHs have been
investigated by many authors using various methods [9, 10, 11, 12, 13, 14, 15,
16, 19, 17, 18, 20, 21, 22]. However, in most of the previous work, only the
local-type of non-Gaussianity is considered. Even for the bispectrum, which is
the lowest non-trivial order in perturbation where the non-Gaussianity
appears, only a limited number of authors consider the other types, such as
equilateral- and orthogonal-types of non-Gaussianity [12, 14]. The effect of
the trispectrum of non-Gaussianity on the abundance have not been discussed
much, while the clustering of PBHs in the presence of local-type trispectra
are considered in, e.g., Refs. [23, 24, 25]. Thus, it is desirable to
understand the effect of various types of non-Gaussianity on the PBH formation
more systematically.

In this paper, we analytically study the effect of non-Gaussianity up through
the trispectrum in a model-independent way as much as possible. To
analytically address the problem, we use a set of approximations that are
valid in many situations. In the presence of skewness and kurtosis in the
primordial fluctuations with hierarchical orders, we derive a formula for the
non-Gaussian corrections to the abundance of PBHs in a high-peaks limit,
generalizing the known formula for the threshold. We also derive integral
formulas to calculate the values of skewness and kurtosis in popular models of
non-Gaussianity, including local-, equilateral-, folded- and orthogonal-type
models for bispectrum, and generalized local-type models for trispectrum.
These integral formulas are further reduced to asymptotic formulas of analytic
forms by taking the limit of a sharply-peaked shape of the primordial
spectrum. The derived formulas provide useful apparatus for predicting the
abundance of PBHs in a variety of models.

This paper is organized as follows: In section 2, after the basic notations
and definitions we use in this paper are given, high-peaks formulas for the
PBH formation are introduced. In section 3, we derive integral formulas for
skewness and kurtosis parameters for a class of non-Gaussian models mentioned
above, for an arbitrary shape of the primordial spectrum. In section 4, the
derived integral formulas are reduced to analytic forms by taking the sharply-
peaked limit of the spectrum. In section 5, we numerically evaluate the
behavior of derived formulas and demonstrate it in several cases.

##  2  The abundance of PBHs from non-Gaussian initial conditions

###  2.1 Definitions

First, we define the fundamental quantities used in this paper. The PBHs are
considered to be formed from large positive perturbations of the 3-dimensional
curvature $\delta R^{(3)}$ at an early stage of the Universe, well before the
astrophysical structure formation takes place. The 3-dimensional curvature
perturbation on comoving slices is characterized by the curvature perturbation
in Fourier space, denoted by ${\cal R}(\bm{k})$. In linear order, we have
$\delta R^{(3)}=4(k^{2}/a^{2}){\cal R}$, where $a$ is the scale factor. Since
the Einstein equations imply $\delta R^{(3)}\approx 6H^{2}\varDelta$ on or
above Hubble scales at linear order, where $H=\dot{a}/a$ is the Hubble
parameter and $\varDelta$ is the energy density contrast on comoving slices,
large positive 3-dimensional curvature perturbations are equivalent to large
positive energy density fluctuations.

If its probability distribution is Gaussian, the statistical properties of the
comoving curvature perturbation are completely characterized by the power
spectrum $P_{\cal R}(k)$. The power spectrum is the two-point correlation of
the fluctuations in Fourier space, and is defined by

$\langle{\cal R}(\bm{k}_{1}){\cal
R}(\bm{k}_{2})\rangle_{\mathrm{c}}=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2})P_{\cal
R}(k_{1}),$ (2.1)

where $\delta_{\mathrm{D}}^{3}(\bm{k})$ is the Dirac’s delta function and
$\langle\cdots\rangle_{\mathrm{c}}$ denotes the cumulants or connected part of
the statistical average. For the two-point cumulant above, it can be replaced
by a simple average, assuming the mean value of curvature fluctuations is
zero, $\langle{\cal R}\rangle=0$.

The presence of non-Gaussianity gives rise to higher-order correlations, such
as the bispectrum, trispectrum, and so forth. The bispectrum $B_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ and the trispectrum $T_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$ are three- and four-point
functions in Fourier space, defined by

$\displaystyle\langle{\cal R}(\bm{k}_{1}){\cal R}(\bm{k}_{2}){\cal
R}(\bm{k}_{3})\rangle_{\mathrm{c}}$
$\displaystyle=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})B_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}),$ (2.2) $\displaystyle\langle{\cal
R}(\bm{k}_{1}){\cal R}(\bm{k}_{2}){\cal R}(\bm{k}_{3}){\cal
R}(\bm{k}_{4})\rangle_{\mathrm{c}}$
$\displaystyle=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}+\bm{k}_{4})T_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4}).$ (2.3)

The three-point cumulant in Eq. (2.2) again can be replaced by simple average
because the mean value is zero. However, the cumulant in Eq. (2.3) cannot be
replaced by simple average, and we have $\langle{\cal R}_{1}{\cal R}_{2}{\cal
R}_{3}{\cal R}_{4}\rangle_{\mathrm{c}}=\langle{\cal R}_{1}{\cal R}_{2}{\cal
R}_{3}{\cal R}_{4}\rangle-\langle{\cal R}_{1}{\cal R}_{2}\rangle\langle{\cal
R}_{3}{\cal R}_{4}\rangle-\langle{\cal R}_{1}{\cal R}_{3}\rangle\langle{\cal
R}_{2}{\cal R}_{4}\rangle-\langle{\cal R}_{1}{\cal R}_{4}\rangle\langle{\cal
R}_{2}{\cal R}_{3}\rangle$. The appearance of Dirac’s delta functions in the
above definition is due to the assumed statistical homogeneity in 3-space with
translational symmetry.

In general, in linear order, the relation between the comoving curvature
perturbation ${\cal R}$ and the density contrast $\varDelta$ on comoving
slices is given by [26, 27, 28]

$\varDelta(\bm{k};t)=\mathcal{M}(k){\cal R}(\bm{k})\,,$ (2.4)

where

$\mathcal{M}(k;t)=\frac{2+2w}{5+3w}\left(\frac{k}{aH}\right)^{2},$ (2.5)

and $w=p/\rho$ is the equation of state parameter. The PBH formation criteria
are typically described by smoothed density fields with a smoothing radius of
the horizon scale, $R=(aH)^{-1}$. In the following, we assume the formation of
PBHs takes place in a radiation-dominated epoch, $w=1/3$, and thus we use

$\mathcal{M}(k)=\frac{4}{9}k^{2}R^{2}\,,$ (2.6)

for the coefficient of Eq. (2.4).

The smoothed density field in configuration space is given by

$\varDelta_{R}(\bm{x})=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\bm{k}\cdot\bm{x}}\varDelta(\bm{k})W(kR)$
(2.7)

where $W(kR)$ is the window function for the smoothing. In this paper, we
adopt the Gaussian window function,

$W(kR)=\exp\left(-\frac{k^{2}R^{2}}{2}\right)\,.$ (2.8)

###  2.2 The abundance of PBHs

The abundance of PBHs is frequently modeled by the initial mass fraction
$\beta$ of the universe that turns into PBHs at the time of formation. One of
the simplest and most commonly used criteria for the PBH formation is to set a
threshold density contrast $\varDelta_{\mathrm{c}}$. Then $\beta$ may be
estimated by computing the fraction of space with
$\varDelta_{R}(\bm{x})>\varDelta_{\mathrm{c}}$ in the initial density field.
This gives

$\beta_{\mathrm{th}}=\int_{\varDelta_{\mathrm{c}}}^{\infty}P(\varDelta_{R};R)\,d\varDelta_{R},$
(2.9)

where $P(\varDelta_{R};R)$ is the probability distribution function of the
initial density field smoothed over a horizon scale $R$. In the analogy to the
Press-Schechter theory of structure formation [29], a fudge factor 2 is
sometimes put in front of the right-hand side (RHS) of the above equation. Our
discussion below does not depend on whether this fudge factor is present or
not. Another criteria of the PBH formation is to use the peak theory [30, 27].
Given the number density of peaks $n_{\mathrm{pk}}(\varDelta_{\mathrm{c}};R)$
of the smoothed density field above the threshold, $\varDelta_{\mathrm{c}}$,
the mass fraction is given by

$\beta_{\mathrm{pk}}=(2\pi)^{3/2}R^{3}n_{\mathrm{pk}}(\varDelta_{\mathrm{c}};R),$
(2.10)

where the prefactor $(2\pi)^{3/2}R^{3}$ corresponds to the effective volume of
the Gaussian filter. There is yet another formation for the threshold based on
the so-called Compaction function [31, 32, 22]. But as its relation to the
probability distribution function seems rather non-trivial, we leave this case
for future studies.

When the initial density field is given by a random Gaussian field, the one-
point probability distribution function is given by
$P=(2\pi{\sigma}^{2})^{-1/2}e^{-{\varDelta_{R}}^{2}/(2{\sigma}^{2})}$, where
${\sigma}^{2}=\langle{\varDelta_{R}}^{2}\rangle$ is the variance of the
smoothed density field. In this case, Eq. (2.9) is given by

$\beta_{\mathrm{th}}^{\mathrm{G}}=\frac{1}{2}\mathrm{erfc}\left(\frac{\nu}{\sqrt{2}}\right)\approx\frac{1}{\sqrt{2\pi}}\frac{e^{-\nu^{2}/2}}{\nu},$
(2.11)

where $\nu\equiv\varDelta_{\mathrm{c}}/\sigma$ is the normalized threshold.
The last expression of the above equation is an asymptotic form for a high
threshold of $\nu\gg 1$. The number density of peaks in a random Gaussian
field is calculated from the joint probability distribution of field
derivatives and is analytically given by an integral form [30]. An asymptotic
form of the result for high peaks ($\nu\gg 1$) is given by

$n_{\mathrm{pk}}=\frac{1}{(2\pi)^{2}}\left(\frac{{\sigma_{1}}^{2}}{3{\sigma}^{2}}\right)^{3/2}(\nu^{2}-1)e^{-\nu^{2}/2},$
(2.12)

where

${\sigma_{j}}^{2}=\int\frac{k^{2}dk}{2\pi^{2}}k^{2j}P_{\varDelta}(k)W^{2}(kR)$
(2.13)

and $P_{\varDelta}(k)$ is the power spectrum of the density field $\varDelta$.
The corresponding expression for the mass fraction is given by

$\beta_{\mathrm{pk}}^{\mathrm{G}}=\frac{1}{\sqrt{2\pi}}\left(\frac{R\sigma_{1}}{\sqrt{3}\sigma}\right)^{3}(\nu^{2}-1)e^{-\nu^{2}/2}.$
(2.14)

There are differences in the predictions of the two different models, Eqs.
(2.11) and (2.14) with the Gaussian initial conditions. However, the overall
shape is relatively close to each other if we take a higher threshold value
$\nu$ for the peak theory and the threshold theory for reasonable ranges in
the PBH mass [27, 28].

###  2.3 Non-Gaussian distributions

The one-point distribution function of non-Gaussian fields are characterized
by higher-order cumulants $\langle{\varDelta_{R}}^{n}\rangle_{\mathrm{c}}$. We
assume a hierarchical scaling of higher-order cumulants,
$\langle{\varDelta_{R}}^{n}\rangle_{\mathrm{c}}\propto{\sigma}^{2n-2}$, and
define reduced cumulants

$S_{n}\equiv\frac{\langle{\varDelta_{R}}^{n}\rangle_{\mathrm{c}}}{{\sigma}^{2n-2}},$
(2.15)

which are considered to be of order unity or less. In general, the probability
distribution function in the integrand of Eq. (2.9) can be expanded in an
Edgeworth series,

$P(\varDelta_{R};R)=\frac{e^{-\nu^{2}/2}}{\sqrt{2\pi}\sigma}\Biggl{\\{}1+\frac{S_{3}}{3!}H_{3}(\nu)\sigma+\left[\frac{S_{4}}{4!}H_{4}(\nu)+\frac{1}{2!}\left(\frac{S_{3}}{3!}\right)^{2}H_{6}(\nu)\right]\sigma^{2}\\\
+\left[\frac{S_{5}}{5!}H_{4}(\nu)+\frac{S_{3}S_{4}}{3!4!}H_{7}(\nu)+\frac{1}{3!}\left(\frac{S_{3}}{3!}\right)^{3}H_{9}(\nu)\right]\sigma^{3}+\cdots\Biggr{\\}},$
(2.16)

where $H_{n}(\nu)=e^{\nu^{2}/2}(-d/d\nu)^{n}e^{-\nu^{2}/2}$ is the Hermite
polynomial. In practice, one can truncate the series in some order when the
inequality $\sigma\ll\nu^{-3}$ is satisfied. The Edgeworth expansion has been
used to investigate the abundance of PBHs [12]. For high peaks, $\nu\gg 1$,
the truncated Edgeworth expansion is only applicable for sufficiently small
$\sigma$. Integrating the Edgeworth series above, the mass fraction in the
threshold model, Eq. (2.9) reduces to

$\beta_{\mathrm{th}}=\beta_{\mathrm{th}}^{\mathrm{G}}+\frac{e^{-\nu^{2}/2}}{\sqrt{2\pi}}\Biggl{\\{}\frac{S_{3}}{3!}H_{2}(\nu)\sigma+\left[\frac{S_{4}}{4!}H_{3}(\nu)+\frac{1}{2!}\left(\frac{S_{3}}{3!}\right)^{2}H_{5}(\nu)\right]\sigma^{2}\\\
+\left[\frac{S_{5}}{5!}H_{3}(\nu)+\frac{S_{3}S_{4}}{3!4!}H_{6}(\nu)+\frac{1}{3!}\left(\frac{S_{3}}{3!}\right)^{3}H_{8}(\nu)\right]\sigma^{3}+\cdots\Biggr{\\}}.$
(2.17)

In the high-peaks limit, the mass fraction in the threshold model has been
derived in the context of biased structure formation, and the result is given
by [33, 34]

$\beta_{\mathrm{th}}\approx\frac{1}{\sqrt{2\pi}}\frac{e^{-\nu^{2}/2}}{\nu}\exp\left(\sum_{n=3}^{\infty}\frac{\nu^{n}}{n!}\frac{\langle{\varDelta_{R}}^{n}\rangle}{\sigma^{n}}\right)=\frac{\nu^{-1}}{\sqrt{2\pi}}\exp\left[-\frac{\nu^{2}}{2}\left(1-2\sum_{n=3}^{\infty}\frac{{\varDelta_{\mathrm{c}}}^{n-2}}{n!}S_{n}\right)\right].$
(2.18)

This result can also be obtained by summing up all the infinite series of the
leading contributions in the high-peaks limit $H_{n}(\nu)\rightarrow\nu^{n}$
in the Eq. (2.17). In this way, the abundance of PBHs for non-Gaussian initial
conditions in the high-peaks limit is characterized by the series of reduced
cumulants $S_{n}$.

The above formula for the non-Gaussianity in the high-peaks limit can be
generalized to the peaks model of Eq. (2.14). The details of the derivation
are given in Appendix A. One finds that the non-Gaussian contributions are the
same as those in the high-peaks limit of the threshold model. Namely,

$\beta_{\mathrm{pk}}\approx\frac{1}{\sqrt{2\pi}}\left(\frac{R\sigma_{1}}{\sqrt{3}\sigma}\right)^{3}(\nu^{2}-1)\exp\left[-\frac{\nu^{2}}{2}\left(1-2\sum_{n=3}^{\infty}\frac{{\varDelta_{\mathrm{c}}}^{n-2}}{n!}S_{n}\right)\right].$
(2.19)

Therefore, irrespective of the formation models of PBHs, the effect of non-
Gaussianity in the high-peaks limit may be expressed in the generic form,

$\beta\approx\beta^{\mathrm{G}}\exp\left(\nu^{2}\sum_{n=3}^{\infty}\frac{{\varDelta_{\mathrm{c}}}^{n-2}}{n!}S_{n}\right)=A(\nu)\exp\left[-\frac{\nu^{2}}{2}\left(1-2\sum_{n=3}^{\infty}\frac{{\varDelta_{\mathrm{c}}}^{n-2}}{n!}S_{n}\right)\right],$
(2.20)

where

$\beta^{\mathrm{G}}=A(\nu)e^{-\nu^{2}/2}$ (2.21)

is the mass fraction for the Gaussian initial condition, and $A(\nu)$ is the
prefactor that depends on the formation models of PBHs, i.e.,

$A(\nu)\equiv\frac{1}{\sqrt{2\pi}}\times\begin{cases}\displaystyle\nu^{-1}&\mathrm{(threshold\
model)},\\\
\displaystyle\left(\frac{R\sigma_{1}}{\sqrt{3}\sigma}\right)^{3}(\nu^{2}-1)&\mathrm{(peaks\
model)}.\end{cases}$ (2.22)

We note that the asymptotic formula (2.20) is consistent only when

$2\sum_{n=3}^{\infty}\frac{{\varDelta_{\mathrm{c}}}^{n-2}}{n!}S_{n}<1\,,$
(2.23)

otherwise, the mass fraction $\beta$ exceeds unity in the limit $\nu\gg 1$. In
this paper, we assume the above condition is satisfied.

When the values of reduced cumulants $S_{n}$ are of order unity, the above
condition is safely satisfied. For example, when all the reduced cumulants
$S_{n}$ have the same value, the condition (2.23) with
$\varDelta_{\mathrm{c}}\simeq 1/3$ implies $S_{n}\lesssim 8$. Thus it is not
too restrictive. On the other hand, if the condition (2.23) is not met, the
asymptotic formula (2.20) breaks down. This is because the resummation of
leading contributions in the expansion in Eq. (2.17) for the threshold model
is not justified as $S_{n}$ are no longer of $\mathcal{O}(1)$. The same
applies to the resummation in the peaks model, Eq. (A.16). For example, when
higher-order cumulants of a non-Gaussian model satisfy
$\langle{\varDelta_{R}}^{n}\rangle_{\mathrm{c}}\sim\mathcal{O}(\sigma^{n})$
[instead of $\sim\mathcal{O}(\sigma^{2n-2})$, c.f., Eq. (2.15)], the
resummation does not work [18].

Comparing Eqs. (2.20) and (2.21), and noting
$\nu=\varDelta_{\mathrm{c}}/\sigma$, we find that the effects of non-
Gaussianity may be conveniently taken into account by replacing the threshold
value for the PBH formation in the exponent of the Gaussian prediction, Eq.
(2.21), as

$\varDelta_{\mathrm{c}}\rightarrow\varDelta_{\mathrm{c}}^{\mathrm{eff}}\equiv\varDelta_{\mathrm{c}}\sqrt{1-S}\,;\quad
S\equiv 2\sum_{n=3}^{\infty}\frac{{\varDelta_{\mathrm{c}}}^{n-2}}{n!}S_{n}\,,$
(2.24)

apart from the prefactor $A(\nu)$. When $S$ is positive, the tail of the
distribution function increases, which reduces the effective threshold in
comparison with the Gaussian prediction, and the number of PBHs increases. On
the contrary, if $S$ is negative, the number of PBHs decreases. Thus, the
expected number density of PBHs is exponentially sensitive to non-Gaussianity.
We note that the effective threshold above,
$\varDelta_{\mathrm{c}}^{\mathrm{eff}}$, does not apply to the prefactor
$A(\nu)$. Hence the effects of non-Gaussianity are not completely degenerate
with the Gaussian case, though the change in the exponent dominates the effect
on the number density of PBHs.

When the observational constraint is given by $\beta<\beta_{0}$, where
$\beta_{0}\ll 1$ and thus $\ln(1/\beta_{0})>0$, Eq. (2.20) implies

$\sigma^{2}<\frac{{\varDelta_{\mathrm{c}}}^{2}}{2\ln(1/\beta_{0})}\left(1-S\right),$
(2.25)

where the logarithm of the prefactor $\ln A(\nu)$ is ignored, assuming
$\ln(1/\beta_{0})\gg\ln A(\nu)$. Therefore, an observational constraint on the
upper limit of the amplitude of primordial spectrum $\sqrt{\mathcal{P}_{\cal
R}}$ is tighter for $S>0$. If we only keep the leading term in $S$, $S\propto
S_{3}$, this is in qualitative agreement with the results in Refs. [8, 14] at
linear order in $S_{3}$ (or in $f_{\rm NL}$ as discussed in the next section).
Nonlinear behaviors discussed in these references might be explained by the
effect of ${f_{\mathrm{NL}}}^{2}$ in the variance,
$\sigma^{2}\sim{\sigma_{\mathrm{G}}^{2}}+\mathrm{(const.)}\times{f_{\mathrm{NL}}}^{2}{\sigma_{\mathrm{G}}^{4}}$.
However, for a narrowly peaked power spectrum, the prefactor (const.) of this
relation is suppressed by the width of the spectrum, $\epsilon=\Delta
k/k_{0}\ll 1$, where $k_{0}$ is the peak of the spectrum.

As seen from Eq. (2.20), the parameters of reduced cumulants, $S_{n}$, with
arbitrary higher orders may equally contribute if there are of the same order
of magnitude. In the following, however, we mainly focus on the effects of
bispectrum and trispectrum of the primordial fluctuations, which are
responsible to the skewness parameter $S_{3}$ and the kurtosis parameter
$S_{4}$. When the higher-order cumulants $S_{n}$ with $n\geq 5$ are absent,
Eq. (2.20) reduces to

$\beta\approx\beta^{\mathrm{G}}\exp\left[\nu^{2}\varDelta_{\mathrm{c}}\left(\frac{S_{3}}{6}+\frac{{\varDelta_{\mathrm{c}}}S_{4}}{24}\right)\right].$
(2.26)

The condition of Eq. (2.23) in this case is given by

$S_{3}+\frac{S_{4}}{12}\lesssim 9,$ (2.27)

for $\varDelta_{\mathrm{c}}\simeq 1/3$. The skewness and kurtosis parameters,
$S_{3}$ and $S_{4}$, are determined by non-Gaussian initial conditions. In the
next section, we derive useful formulas to calculate these parameters in
typical models of primordial non-Gaussianity.

##  3  Skewness and kurtosis in models of primordial non-Gaussianity

###  3.1 Definitions

The reduced cumulants of third and fourth orders, $S_{3}$ and $S_{4}$, are
called skewness and kurtosis parameters, respectively. They are related to the
bispectrum $B(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ and the trispectrum
$T(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$ of density field $\varDelta$
by

$\displaystyle S_{3}$
$\displaystyle=\frac{1}{\sigma^{4}}\int_{\bm{k}_{123}=\bm{0}}B(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})W(k_{1}R)W(k_{2}R)W(k_{3}R),$
(3.1) $\displaystyle S_{4}$
$\displaystyle=\frac{1}{\sigma^{6}}\int_{\bm{k}_{1234}=\bm{0}}T(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})W(k_{1}R)W(k_{2}R)W(k_{3}R)W(k_{4}R),$
(3.2)

where we use an abbreviated notation, $\bm{k}_{1\cdots
n}\equiv\bm{k}_{1}+\cdots+\bm{k}_{n}$, and

$\int_{\bm{k}_{1\cdots
n}=\bm{0}}\cdots\equiv\int\frac{d^{3}k_{1}}{(2\pi)^{3}}\cdots\frac{d^{3}k_{n}}{(2\pi)^{3}}(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}_{1}+\cdots+\bm{k}_{n})\cdots.$
(3.3)

The bispectrum and trispectrum of the density field are related to those of
the curvature perturbation by

$\displaystyle B(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$
$\displaystyle=\mathcal{M}(k_{1})\mathcal{M}(k_{2})\mathcal{M}(k_{3})B_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}),$ (3.4) $\displaystyle
T(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$
$\displaystyle=\mathcal{M}(k_{1})\mathcal{M}(k_{2})\mathcal{M}(k_{3})\mathcal{M}(k_{4})T_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4}).$ (3.5)

Substituting Eqs. (2.6), (2.8), (3.4) and (3.5) into Eqs. (3.1) and (3.2), we
have

$\displaystyle S_{3}$
$\displaystyle=\frac{1}{\sigma^{4}}\left(\frac{4}{9}\right)^{3}\int_{\bm{k}_{123}=\bm{0}}(k_{1}R)^{2}(k_{2}R)^{2}(k_{3}R)^{2}e^{-({k_{1}}^{2}+{k_{2}}^{2}+{k_{3}}^{2})R^{2}/2}B_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}),$ (3.6) $\displaystyle S_{4}$
$\displaystyle=\frac{1}{\sigma^{6}}\left(\frac{4}{9}\right)^{4}\int_{\bm{k}_{1234}=\bm{0}}(k_{1}R)^{2}(k_{2}R)^{2}(k_{3}R)^{2}(k_{4}R)^{2}e^{-({k_{1}}^{2}+{k_{2}}^{2}+{k_{3}}^{2}+{k_{4}}^{2})R^{2}/2}T_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4}).$ (3.7)

The variance $\sigma^{2}$ is the same as ${\sigma_{0}}^{2}$ defined by Eq.
(2.13). That is

$\sigma^{2}=\left(\frac{4}{9}\right)^{2}\int\frac{k^{2}dk}{2\pi^{2}}(kR)^{4}e^{-k^{2}R^{2}}P_{\cal
R}(k)\,,$ (3.8)

where $P_{\cal R}(k)$ is the power spectrum of the comoving curvature
perturbation.

###  3.2 Skewness

There are many models of primordial non-Gaussianity. One of the most commonly
assumed model for the bispectrum is the local-type [35, 36, 37], which is
given by

$B_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\frac{6}{5}f_{\mathrm{NL}}\left[P_{\cal
R}(k_{1})P_{\cal R}(k_{2})+\mathrm{cyc.}\right],$ (3.9)

where $f_{\mathrm{NL}}$ is a parameter of non-Gaussianity amplitude,
$+\,\mathrm{cyc.}$ represents the two terms obtained by cyclic permutations of
the preceding term with respect to $k_{1}$, $k_{2}$, $k_{3}$. Among other
types of non-Gaussianity for the bispectrum, popular alternatives are the
equilateral [38], folded [39], and orthogonal [40] types. These may be
constructed from the following elements of the bispectrum:

$B^{\mathrm{I}}_{123}\equiv P_{1}P_{2}+\mathrm{cyc.},\quad
B^{\mathrm{II}}_{123}\equiv(P_{1}P_{2}P_{3})^{2/3},\quad
B^{\mathrm{III}}_{123}\equiv{P_{1}}^{1/3}{P_{2}}^{2/3}P_{3}+\mathrm{5\
perm.},$ (3.10)

where we denote $P_{1}=P_{\cal R}(k_{1})$, $P_{2}=P_{\cal R}(k_{2})$,
$P_{3}=P_{\cal R}(k_{3})$ for simplicity, and $+\,\mathrm{5\ perm.}$
represents the 5 terms obtained by permutations of the preceding term. In
terms of these elements, the bispectrum for each type of non-Gaussianity is
given by

$\displaystyle B^{\mathrm{loc}}_{\cal R}$
$\displaystyle=\frac{6}{5}f_{\mathrm{NL}}^{\mathrm{loc}}B^{\mathrm{I}}_{123}\,,$
(3.11) $\displaystyle B^{\mathrm{eql}}_{\cal R}$
$\displaystyle=\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{eql}}\left(-B^{\mathrm{I}}_{123}-2B^{\mathrm{II}}_{123}+B^{\mathrm{III}}_{123}\right),$
(3.12) $\displaystyle B^{\mathrm{fol}}_{\cal R}$
$\displaystyle=\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{fol}}\left(B^{\mathrm{I}}_{123}+3B^{\mathrm{II}}_{123}-B^{\mathrm{III}}_{123}\right),$
(3.13) $\displaystyle B^{\mathrm{ort}}_{\cal R}$
$\displaystyle=\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{ort}}\left(-3B^{\mathrm{I}}_{123}-8B^{\mathrm{II}}_{123}+3B^{\mathrm{III}}_{123}\right).$
(3.14)

We introduce the skewness parameter elements corresponding to the above three
elements $B^{A}_{123}$ ($A={\rm I},{\rm II},{\rm III}$) of the bispectrum as

$S^{A}_{3}\equiv\frac{1}{\sigma^{4}}\left(\frac{4}{9}\right)^{3}\int_{\bm{k}_{123}=\bm{0}}(k_{1}R)^{2}(k_{2}R)^{2}(k_{3}R)^{2}W(k_{1}R)W(k_{2}R)W(k_{3}R)B^{A}_{123}\,.$
(3.15)

The skewness parameters for local, equilateral, folded, and orthogonal types,
which we denote by $S^{\mathrm{loc}}_{3}$, $S^{\mathrm{eql}}_{3}$,
$S^{\mathrm{fol}}_{3}$, $S^{\mathrm{ort}}_{3}$, respectively, are given by
linear superpositions of $S^{A}_{3}$ ($A={\rm I},{\rm II},{\rm III}$) in
exactly the same forms as Eqs. (3.11) – (3.14),

$\displaystyle S^{\mathrm{loc}}_{3}$
$\displaystyle=\frac{6}{5}f_{\mathrm{NL}}^{\mathrm{loc}}S^{\mathrm{I}}_{3},$
(3.16) $\displaystyle S^{\mathrm{eql}}_{3}$
$\displaystyle=\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{eql}}\left(-S^{\mathrm{I}}_{3}-2S^{\mathrm{II}}_{3}+S^{\mathrm{III}}_{3}\right),$
(3.17) $\displaystyle S^{\mathrm{fol}}_{3}$
$\displaystyle=\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{fol}}\left(S^{\mathrm{I}}_{3}+3S^{\mathrm{II}}_{3}-S^{\mathrm{III}}_{3}\right),$
(3.18) $\displaystyle S^{\mathrm{ort}}_{3}$
$\displaystyle=\frac{18}{5}f_{\mathrm{NL}}^{\mathrm{ort}}\left(-3S^{\mathrm{I}}_{3}-8S^{\mathrm{II}}_{3}+3S^{\mathrm{III}}_{3}\right).$
(3.19)

We note that, excluding the signs of $f_{\rm NL}^{X}$
($X=\mathrm{loc},\mathrm{eql},\mathrm{fol},\mathrm{ort}$) in the coefficients,
$S^{\mathrm{\rm loc}}_{3}$ is positive definite, while the sign of the other
two is indeterminate. Given the values of $f_{\mathrm{NL}}^{X}$, these
skewness parameters are uniquely determined once the primordial power spectrum
$P_{\cal R}(k)$ is known. In the high-peaks limit, the effects of skewness in
each non-Gaussian type on the abundance of PBHs are given by substituting the
results into Eq. (2.26).

The skewness parameter elements are given by substituting Eqs. (3.10) into Eq.
(3.15). To represent the results in a convenient, compact form, we introduce
the dimensionless curvature perturbation power spectrum,

$\mathcal{P}_{\cal R}(k)\equiv\frac{k^{3}P_{\cal R}(k)}{2\pi^{2}}\,.$ (3.20)

Changing the integration variables as $\bm{p}=\bm{k}_{1}R$,
$\bm{q}=\bm{k}_{2}R$ and $r=|\bm{p}+\bm{q}|$, where the variable $r$ describes
the angular degrees of freedom,
$\mu=\bm{p}\cdot\bm{q}/(pq)=(r^{2}-p^{2}-q^{2})/(2pq)$, some of the angular
integrations may be analytically calculated to give

$\displaystyle S^{\mathrm{I}}_{3}$
$\displaystyle=\frac{3}{\sigma^{4}}\left(\frac{4}{9}\right)^{3}\int_{0}^{\infty}dp\,dq\,e^{-p^{2}-q^{2}}pq\left[(p^{2}+q^{2}+2)\frac{\sinh(pq)}{pq}-2\cosh(pq)\right]\mathcal{P}_{\cal
R}\left(\frac{p}{R}\right)\mathcal{P}_{\cal R}\left(\frac{q}{R}\right),$
(3.21) $\displaystyle S^{\mathrm{II}}_{3}$
$\displaystyle=\frac{1}{\sigma^{4}}\left(\frac{4}{9}\right)^{3}\int_{0}^{\infty}dp\,dq\,e^{-(p^{2}+q^{2})/2}pq\left[\mathcal{P}_{\cal
R}\left(\frac{p}{R}\right)\mathcal{P}_{\cal
R}\left(\frac{q}{R}\right)\right]^{2/3}\int_{|p-q|}^{p+q}\frac{dr}{2}\,r\,e^{-r^{2}/2}\left[\mathcal{P}_{\cal
R}\left(\frac{r}{R}\right)\right]^{2/3},$ (3.22) $\displaystyle
S^{\mathrm{III}}_{3}$
$\displaystyle=\frac{6}{\sigma^{4}}\left(\frac{4}{9}\right)^{3}\int_{0}^{\infty}dp\,dq\,p^{2}q\,e^{-(p^{2}+q^{2})/2}\left[\mathcal{P}_{\cal
R}\left(\frac{p}{R}\right)\right]^{1/3}\left[\mathcal{P}_{\cal
R}\left(\frac{q}{R}\right)\right]^{2/3}\int_{|p-q|}^{p+q}\frac{dr}{2}\,e^{-r^{2}/2}\,\mathcal{P}_{\cal
R}\left(\frac{r}{R}\right).$ (3.23)

Thus, all the skewness parameters in the three types of non-Gaussianity can be
computed from the above equations for an arbitrary power spectrum $P_{\cal
R}(k)$. The variance $\sigma^{2}$ of Eq. (3.8) is similarly given by a one-
dimensional integral:

$\sigma^{2}=\left(\frac{4}{9}\right)^{2}\int
dp\,p^{3}e^{-p^{2}}\mathcal{P}_{\cal R}\left(\frac{p}{R}\right).$ (3.24)

###  3.3 Kurtosis

In contrast to the case of the bispectrum, not so many variations in the types
of the primordial trispectrum have been proposed. Lacking considerations on
the general types of the trispectrum, here we focus on the most popular type,
that is the generalized local-type [41],

$T_{\cal
R}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})=\frac{54}{25}g_{\mathrm{NL}}\left[P_{1}P_{2}P_{3}+\mathrm{3\
perm.}\right]+\tau_{\mathrm{NL}}\left[P_{1}P_{2}P_{23}+\mathrm{11\
perm.}\right],$ (3.25)

where $P_{23}=P_{\cal R}(|\bm{k}_{2}+\bm{k}_{3}|)$, with $g_{\mathrm{NL}}$ and
$\tau_{\mathrm{NL}}$ being the parameters. If the primordial perturbations
emerge from the quantum fluctuations of a single scalar field, the last
parameter is related to the parameter of the local-type bispectrum by
$\tau_{\mathrm{NL}}=(36/25){f_{\mathrm{NL}}}^{2}$ [42]. If multiple scalar
fields are involved, there is an inequality,
$\tau_{\mathrm{NL}}>(36/25){f_{\mathrm{NL}}}^{2}$ [43]. In this model, we
define elements of the trispectrum by

$T^{\mathrm{I}}_{1234}\equiv P_{1}P_{2}P_{3}+\mathrm{3\ perm.},\quad
T^{\mathrm{II}}_{1234}\equiv P_{1}P_{2}P_{23}+\mathrm{11\ perm.}.$ (3.26)

The corresponding elements of kurtosis are given by

$S^{A}_{4}\equiv\frac{1}{\sigma^{6}}\left(\frac{4}{9}\right)^{4}\int_{\bm{k}_{1234}=\bm{0}}(k_{1}R)^{2}(k_{2}R)^{2}(k_{3}R)^{2}(k_{4}R)^{2}W(k_{1}R)W(k_{2}R)W(k_{3}R)W(k_{4}R)T^{A}_{1234},$
(3.27)

where $A=\mathrm{I},\mathrm{II}$. We note that both elements are positive
definite, as clear from their definitions. The resulting kurtosis for the
local-type is given by

$S_{4}=\frac{54}{25}g_{\mathrm{NL}}S^{\mathrm{I}}_{4}+\tau_{\mathrm{NL}}S^{\mathrm{II}}_{4}.$
(3.28)

Once the primordial power spectrum is given, one can evaluate the kurtosis.
The effect on the abundance of PBHs in the high-peaks limit is given by
substituting the result into Eq. (2.26).

The kurtosis parameter elements are given by substituting Eqs. (3.26) into Eq.
(3.27). Changing integration variables as $\bm{p}=\bm{k}_{1}R$,
$\bm{q}=\bm{k}_{2}R$, $\bm{r}=(\bm{k}_{2}+\bm{k}_{3})R$,
$\mu=\bm{p}\cdot\bm{r}/(pr)$, $\mu^{\prime}=-\bm{q}\cdot\bm{r}/(qr)$, and
expressing $\mu^{\prime}$ in terms of $s=|\bm{q}-\bm{r}|$, some of the angular
integrations can be analytically calculated (see Ref. [44] for the same type
of calculation). The results are

$\displaystyle S^{\mathrm{I}}_{4}$
$\displaystyle=\frac{4}{\sigma^{6}}\left(\frac{4}{9}\right)^{4}\int_{0}^{\infty}dp\,dq\,dr\,e^{-p^{2}}e^{-(q^{2}+r^{2})/2}pr$
$\displaystyle\hskip
42.0pt\times\left[(p^{2}+r^{2}+2)\frac{\sinh(pr)}{pr}-2\cosh(pr)\right]\mathcal{P}_{\cal
R}\left(\frac{p}{R}\right)\mathcal{P}_{\cal
R}\left(\frac{q}{R}\right)\int_{|q-r|}^{q+r}\frac{ds}{2}e^{-s^{2}/2}\mathcal{P}_{\cal
R}\left(\frac{s}{R}\right),$ (3.29) $\displaystyle S^{\mathrm{II}}_{4}$
$\displaystyle=\frac{12}{\sigma^{6}}\left(\frac{4}{9}\right)^{4}\int_{0}^{\infty}dp\,dq\,dr\,e^{-p^{2}-q^{2}-r^{2}}\frac{pq}{r}\left[(p^{2}+r^{2}+2)\frac{\sinh(pr)}{pr}-2\cosh(pr)\right]$
$\displaystyle\hskip
108.0pt\times\left[(q^{2}+r^{2}+2)\frac{\sinh(qr)}{qr}-2\cosh(qr)\right]\mathcal{P}_{\cal
R}\left(\frac{p}{R}\right)\mathcal{P}_{\cal
R}\left(\frac{q}{R}\right)\mathcal{P}_{\cal R}\left(\frac{r}{R}\right).$
(3.30)

The above formulas enable us to compute all the kurtosis parameters of the
local type for an arbitrary primordial power spectrum $P_{\cal R}(k)$.

##  4  Power spectrum with a narrow peak

In the previous section, we derived formulas for the skewness and kurtosis
that can be used to evaluate their effects on the abundance of PBHs for a
general primordial power spectrum. In this section, we focus on the case of a
narrowly peaked spectrum, which would lead to a nearly monochromatic mass
function of PBHs.

We consider the case that the primordial power spectrum peaked at a wavenumber
$k_{0}$. When the sharpness of the peak is extreme, we can approximately
substitute

$f(k)\mathcal{P}_{\cal R}(k)\rightarrow f(k_{0})\mathcal{P}_{\cal R}(k),$
(4.1)

where $f(k)$ is an arbitrary function of $k$. Applying this substitution, the
variance of Eq. (3.24) reduces to

$\sigma^{2}\simeq\left(\frac{4}{9}\right)^{2}e^{-{k_{0}}^{2}R^{2}}(k_{0}R)^{4}\int\frac{dk}{k}\,\mathcal{P}_{\cal
R}(k)\,,$ (4.2)

and the skewness parameter elements in Eqs. (3.21) – (3.23) reduce to

$\displaystyle S^{\mathrm{I}}_{3}$
$\displaystyle\simeq\frac{27}{2}\,\frac{1}{(k_{0}R)^{4}}\left[({k_{0}}^{2}R^{2}+1)\frac{\sinh({k_{0}}^{2}R^{2})}{{k_{0}}^{2}R^{2}}-\cosh({k_{0}}^{2}R^{2})\right]\,,$
(4.3) $\displaystyle S^{\mathrm{II}}_{3}$
$\displaystyle\simeq\frac{9}{8}\,\frac{e^{{k_{0}}^{2}R^{2}/2}}{{k_{0}}^{2}R^{2}}\frac{{C_{2/3}}^{3}}{{C_{1}}^{2}},$
(4.4) $\displaystyle S^{\mathrm{III}}_{3}$
$\displaystyle\simeq\frac{27}{4}\,\frac{e^{{k_{0}}^{2}R^{2}/2}}{{k_{0}}^{2}R^{2}}\frac{C_{1/3}C_{2/3}}{C_{1}},$
(4.5)

where we have introduced

$C_{\alpha}\equiv\int\frac{dk}{k}\left[\mathcal{P}_{\cal
R}(k)\right]^{\alpha}.$ (4.6)

The integral range of the above is localized in the vicinity of $k_{0}$. We
note that $S^{\mathrm{I}}_{3}$ takes the minimum value
$S^{\mathrm{I}}_{3}|_{\rm min}\approx 4.93$ at $k_{0}R\approx 1.7$. We also
note that $S^{\mathrm{II}}_{3}$ and $S^{\mathrm{III}}_{3}$ depend on the
detailed shape of the peak in the power spectrum through the integrals
$C_{\alpha}$.

Due to the assumption of a narrow peak, the values of integrals $C_{\alpha}$
are small. If we characterize the narrowness of the power spectrum by
$\epsilon$ in the space of wavenumber, we have
$C_{\alpha}\sim\mathcal{O}(\epsilon)$ for $\alpha>0$. Thus the parameters
$S^{\mathrm{II}}_{3}$ and $S^{\mathrm{III}}_{3}$ are of order
$\sim\epsilon\,S^{\mathrm{I}}_{3}$,

$S^{\mathrm{I}}_{3}\gg S^{\mathrm{II}}_{3},S^{\mathrm{III}}_{3}.$ (4.7)

This is a remarkable property of a narrow power spectrum that simplifies the
bispectrum in Eqs. (3.11) – (3.14). For example, if the narrow shape of the
power spectrum is given by a rectangular function of width $k_{0}\epsilon$,

$\mathcal{P}_{\cal R}(k)=\begin{cases}A_{0},&\displaystyle\mathrm{if\ \
}1-\frac{\epsilon}{2}\leq\frac{k}{k_{0}}\leq 1+\frac{\epsilon}{2},\\\
0&\mathrm{otherwise},\end{cases}$ (4.8)

where $A_{0}$ is a normalization constant, then we have
$C_{\alpha}={A_{0}}^{\alpha}\epsilon$, which gives

$S^{\mathrm{II}}_{3}\simeq\frac{9\epsilon}{8}\frac{e^{{k_{0}}^{2}R^{2}/2}}{(k_{0}R)^{2}},\qquad
S^{\mathrm{III}}_{3}\simeq\frac{27\epsilon}{4}\frac{e^{{k_{0}}^{2}R^{2}/2}}{(k_{0}R)^{2}}.$
(4.9)

If the narrow shape is given by a Gaussian function of the full width at half
maximum (FWHM)111FWHM for a normal distribution of standard deviation $\sigma$
is given by $2\sqrt{2\ln 2}\,\sigma\simeq 2.35482\sigma$. $k_{0}\epsilon$,

$\mathcal{P}_{\cal R}(k)=A_{0}\exp\left[-\frac{4\ln
2}{k_{0}^{2}\epsilon^{2}}(k-k_{0})^{2}\right],$ (4.10)

then we have $C_{\alpha}=\sqrt{\pi/\ln
2}\,{\alpha}^{-1/2}{A_{0}}^{\alpha}\epsilon$, which gives

$S^{\mathrm{II}}_{3}\simeq\frac{81\pi^{1/2}\epsilon}{32\sqrt{6\ln
2}}\frac{e^{{k_{0}}^{2}R^{2}/2}}{(k_{0}R)^{2}},\qquad
S^{\mathrm{III}}_{3}\simeq\frac{81\pi^{1/2}\epsilon}{8\sqrt{2\ln
2}}\frac{e^{{k_{0}}^{2}R^{2}/2}}{(k_{0}R)^{2}}.$ (4.11)

In any case, these parameters are much smaller than $S^{\mathrm{I}}_{3}$ by a
factor of $\epsilon$.

Due to this property, the skewness parameters in Eqs. (3.16) – (3.19) are
solely determined by a single element $S^{\mathrm{I}}_{3}$. Taking account of
contributions from all types, the skewness of the density field is given by

$S_{3}\simeq\frac{6}{5}\left(f_{\mathrm{NL}}^{\mathrm{loc}}-3f_{\mathrm{NL}}^{\mathrm{eql}}+3f_{\mathrm{NL}}^{\mathrm{fol}}-9f_{\mathrm{NL}}^{\mathrm{ort}}\right)S^{\mathrm{I}}_{3}.$
(4.12)

Thus different types of non-Gaussianity give the same form of skewness with
different amplitudes. This means that one cannot distinguish non-Gaussian
types only with the amplitude of the skewness. Apart from the numerical
coefficients and the signs in (4.12), this fact may be considered rather
trivial because the abundance of PBHs is just a single number. Different non-
Gaussian types cannot be distinguished by a single number.

Similarly, the substitution of Eq. (4.1) in Eqs. (3.29) and (3.30), the
kurtosis parameter elements are given by

$\displaystyle S^{\mathrm{I}}_{4}$
$\displaystyle\simeq\frac{243}{8}\frac{1}{(k_{0}R)^{8}}\Biggl{\\{}-e^{-{k_{0}}^{2}R^{2}}\sinh\left(2{k_{0}}^{2}R^{2}\right)+\frac{e^{3{k_{0}}^{2}R^{2}/2}}{4k_{0}R}\sqrt{2\pi}\left[2-3\,\mathrm{erfc}\left(\frac{k_{0}R}{\sqrt{2}}\right)+\mathrm{erfc}\left(\frac{3k_{0}R}{\sqrt{2}}\right)\right]\Biggr{\\}},$
(4.13) $\displaystyle S^{\mathrm{II}}_{4}$
$\displaystyle\simeq\frac{243}{(k_{0}R)^{8}}\left[({k_{0}}^{2}R^{2}+1)\frac{\sinh({k_{0}}^{2}R^{2})}{{k_{0}}^{2}R^{2}}-\cosh({k_{0}}^{2}R^{2})\right]^{2}\simeq\frac{4}{3}\left(S^{\mathrm{I}}_{3}\right)^{2}.$
(4.14)

The narrowness parameter $\epsilon$ discussed above cancels in this case, and
therefore $S^{\mathrm{I}}_{4}$ and $S^{\mathrm{II}}_{4}$ are equally of order
unity in $\epsilon$. Hence both parameters do not depend on the detailed shape
of the peak, as in the case of $S^{\mathrm{I}}_{3}$.

Taking into account all the leading contributions to the skewness and
kurtosis, the mass fraction with non-Gaussianity given by Eq. (2.26) in the
high peaks limit with a narrow spectrum may be expressed in a single formula,

$\beta\simeq\beta^{\mathrm{G}}\exp\left\\{\nu^{2}\varDelta_{\mathrm{c}}\left[\frac{1}{5}f_{\mathrm{NL}}^{\mathrm{eff}}S^{\mathrm{I}}_{3}+\frac{{\varDelta_{\mathrm{c}}}}{12}\left(\frac{27}{25}g_{\mathrm{NL}}S^{\mathrm{I}}_{4}+\frac{2}{3}\tau_{\mathrm{NL}}(S^{\mathrm{I}}_{3})^{2}\right)\right]\right\\},$
(4.15)

where

$f_{\mathrm{NL}}^{\mathrm{eff}}\equiv
f_{\mathrm{NL}}^{\mathrm{loc}}-3f_{\mathrm{NL}}^{\mathrm{eql}}+3f_{\mathrm{NL}}^{\mathrm{fol}}-9f_{\mathrm{NL}}^{\mathrm{ort}}\,,$
(4.16)

irrespectively of the detailed shape of the narrow peak. The evaluation of the
above formula requires only two parameters of skewness and kurtosis,
$S^{\mathrm{I}}_{3}$ and $S^{\mathrm{I}}_{4}$.

For the Gaussian part of the mass fraction, $\beta^{\mathrm{G}}$, the peaks
model, Eq. (2.22), contains a factor $R\sigma_{1}/\sigma$. Taking the limit of
a narrow spectrum in Eq. (2.13), we obtain $\sigma_{j}={k_{0}}^{j}\sigma$.
Therefore, we have

$\beta^{\mathrm{G}}=\frac{e^{-\nu^{2}/2}}{\sqrt{2\pi}}\times\begin{cases}\displaystyle\nu^{-1},&\mathrm{(threshold\
model)},\\\
\displaystyle\left(\frac{k_{0}R}{\sqrt{3}}\right)^{3}(\nu^{2}-1),&\mathrm{(peaks\
model)},\end{cases}$ (4.17)

in the limit of a narrow spectrum.

##  5  Numerical demonstrations

In this section, to obtain an intuitive sense of our results, and to find the
range of validity of the narrow-peak limit approximation of the spectrum, we
numerically evaluate the formulas we derived in the previous section.

Figure 1:  Left panel: Elements of skewness parameters as a function of
$k_{0}R$ in the narrow limit of spectral shapes. The elements
$S_{3}^{\mathrm{I}}$, $S_{3}^{\mathrm{II}}$ and $S_{3}^{\mathrm{III}}$ are
plotted. For the latter two, they depend on the shape of the narrow power
spectra. Rectangular shapes are indicated by “Rec” and Gaussian shapes are
indicated by “Gau” in the labels of the curves, and the narrowness parameter
is taken to be $\epsilon=0.03$ in the plot as an example. For another value of
the narrowness parameter, $S_{3}^{\mathrm{II}}$ and $S_{3}^{\mathrm{III}}$ are
simply proportional to $\epsilon$ in the plot. Right panel: Elements of
kurtosis parameters in the narrow limit of spectral shapes. The elements
$S_{4}^{\mathrm{I}}$ and $S_{4}^{\mathrm{II}}$ are plotted. These parameters
do not depend on the details of the shapes of the narrow spectra.

In the left panel of Fig. 1, the skewness parameter elements,
$S_{3}^{\mathrm{I}}$, $S_{3}^{\mathrm{II}}$ and $S_{3}^{\mathrm{III}}$
computed from Eqs. (4.3) – (4.5) in the narrow peak limit are shown as
functions of $k_{0}R$. The element $S_{3}^{\mathrm{I}}$ does not depend on the
precise shape of the spectrum in the narrow limit, while the other two
elements, $S_{3}^{\mathrm{II}}$ and $S_{3}^{\mathrm{III}}$, depend on the
shape of the narrow peak, hence the rectangular case, Eq. (4.8) and the
Gaussian case, Eq. (4.10) are presented. The width of the peak characterized
by the narrowness parameter $\epsilon$ is taken to be $\epsilon=0.03$ in the
plot just for illustration. We recall that the elements $S_{3}^{\mathrm{II}}$
and $S_{3}^{\mathrm{III}}$ are simply proportional to $\epsilon$. As seen from
Fig. 1, $S_{3}^{\mathrm{I}}$ dominates the others when $\epsilon$ is small.
For the value of $\epsilon=0.03$, $S_{3}^{\mathrm{I}}$ is approximately an
order of magnitude larger than the other two.

In the right panel of Fig. 1, the kurtosis parameter elements,
$S_{4}^{\mathrm{I}}$ and $S_{4}^{\mathrm{II}}$ computed by using Eqs. (3.29)
and (3.30) are shown in the narrow peak limit. Both of these two do not depend
on the precise shape of the narrowness parameter $\epsilon$, and both
contribute to the kurtosis parameter, irrespective of the narrowness of the
spectrum. Since the values of the two elements $S_{4}^{\mathrm{I}}$ and
$S_{4}^{\mathrm{II}}$ are roughly of $\mathcal{O}(10^{1})$ \-
$\mathcal{O}(10^{2})$, Eq. (2.27) requires that the parameters
$g_{\mathrm{NL}}$ and $\tau_{\mathrm{NL}}$ should be somewhat smaller than
unity to satisfy the condition for the consistency of the approximation.

Figure 2:  Left panel: Comparisons between exact integrations and
approximations of the narrow limit of primordial spectrum for three elements
of the skewness parameter, $S_{3}^{\mathrm{I}}$, $S_{3}^{\mathrm{II}}$ and
$S_{3}^{\mathrm{III}}$, as functions of the narrowness parameter $\epsilon$.
The spectral shape is assumed to have the rectangular shape, and
$k_{0}R=\sqrt{2}$. In the limits that the narrowness parameter $\epsilon$ is
small, the approximation of the narrow limit and the results of exact
integrations converge to the same lines to the left. The lines of the narrow
limit are given by straight lines as the narrow limit of $S_{3}^{\mathrm{I}}$
is constant and those of the other are linearly proportional to $\epsilon$.
Curved lines above the straight lines correspond to the results of exact
integrations. Right panel: Same as the left panel, but elements of skewness
parameter are compared.

In the left panel of Fig. 2, comparisons of the narrow limit approximation
with the exact numerical results are made for the skewness parameter elements,
$S_{3}^{\mathrm{I}}$, $S_{3}^{\mathrm{II}}$ and $S_{3}^{\mathrm{III}}$, and
for the kurtosis parameter elements, $S_{4}^{\mathrm{I}}$ and
$S_{4}^{\mathrm{II}}$. The right panel shows the comparison of the exact
integrations of Eqs. (3.21) – (3.23) with the narrow limit approximation for
the elements of skewness. The narrow limit of $S_{3}^{\mathrm{I}}$ is
independent of $\epsilon$, while those of $S_{3}^{\mathrm{II}}$ and
$S_{3}^{\mathrm{III}}$ are linearly proportional to $\epsilon$. The curved
lines represent the exact numerical results, and they converge to the
corresponding results in the narrow limit approximation as $\epsilon\to 0$. As
can be seen from the plot, the approximation of the narrow limit is accurate
for $\epsilon\lesssim 0.3$. For larger values of $\epsilon$, the
approximations are still fairly good, although the value of
$S_{3}^{\mathrm{III}}$ is non-negligible in comparison to that of
$S_{3}^{\mathrm{I}}$.

In the right panel of Fig. 2, the kurtosis parameter elements are similarly
compared. The narrow limits are given by the constant lines, while the exact
numerical results are given by the curved lines. The narrow limit
approximation is also valid for sufficiently small $\epsilon$, though the
range of validity seems slightly smaller, $\epsilon\lesssim 0.2$ in comparison
with the case of the skewness.

Figure 3:  Same as the Fig. 2, but the Gaussian shape of the power spectrum is
assumed.

In Fig. 3, the same comparisons as Fig. 2 are made, but for the Gaussian shape
of the spectrum. The narrow limits of $S_{3}^{\mathrm{I}}$,
$S_{4}^{\mathrm{I}}$ and $S_{4}^{\mathrm{II}}$ are the same as those in Fig. 2
because they do not depend on the precise shape of the narrow spectrum.
However, the narrow limits of the elements $S_{3}^{\mathrm{II}}$ and
$S_{3}^{\mathrm{III}}$, as well as the exact numerical results of all the
elements do depend on the precise shape of the spectrum. As seen from the
figures, the narrow limit approximation is slightly worse than that in the
rectangular case for the same value of $\epsilon$. The gradient of
$S_{3}^{\mathrm{III}}$ as a function of $\epsilon$ is larger in the Gaussian
case than in the rectangular case. This results in the fact that
$S_{3}^{\mathrm{III}}$ is non-negligible in comparison with
$S_{3}^{\mathrm{I}}$ already for mildly small values of $\epsilon$.

Figure 4:  Left panel: The mass fraction as a function of amplitude $\sigma$
of the primordial spectrum with the threshold model. The narrow limit of the
spectral shape is assumed. The lower (blue) curve corresponds to the Gaussian
initial condition and the upper (orange) curve corresponds to the non-Gaussian
initial condition given by a model of Eq. (4.15) with the narrow limit of
spectrum and parameters are given by $\varDelta_{\mathrm{c}}=1/3$,
$f_{\mathrm{NL}}^{\mathrm{eff}}=g_{\mathrm{NL}}=\tau_{\mathrm{NL}}=0.2$, and
$k_{0}R=\sqrt{2}$. Right panel: Same as the right panel but with the peaks
model.

In Fig. 4, we show the PBH mass fraction $\beta$ for the threshold model (left
panel) and for the peaks model (right panel) in the narrow spectral shape
limit, given by Eq. (4.15). In this limit, the mass fraction does not depend
on the precise shape of the spectrum, nor the precise value of the narrowness
parameter. The adopted values of the parameters in the plot are
$\varDelta_{\mathrm{c}}=1/3$,
$f_{\mathrm{NL}}^{\mathrm{eff}}=g_{\mathrm{NL}}=\tau_{\mathrm{NL}}=0.2$, and
$k_{0}R=\sqrt{2}$. As can be seen in the figure, the PBH mass fraction is
significantly enhanced for positive values of the non-Gaussian parameters
$f_{\mathrm{NL}}^{\mathrm{eff}}$, $g_{\mathrm{NL}}$, and $\tau_{\mathrm{NL}}$.

The peaks model is known to predict a larger mass fraction than the threshold
model for Gaussian initial conditions [45]. In the narrow spectral shape
limit, however, the non-Gaussian factor in the exponent of Eq. (4.15) does not
depend on whether the formation criterion is determined by the threshold or
peaks model, because it contains only $S_{3}^{\mathrm{I}}$ and
$S_{4}^{\mathrm{I}}$, given by Eqs. (4.3) and (4.13), respectively, which are
independent of formation criteria. Therefore, the enhancement factor for the
non-Gaussian case in comparison with the Gaussian case is the same for both
panels.

##  6 Conclusions

In this paper, we considered the effect of primordial non-Gaussianities on the
abundance of PBHs in the conventional PBH formation scenario in which a PBH is
formed from a rare, large positive curvature perturbation at a radiation-
dominated stage in the early universe. We first presented a general series
form of the one-point distribution function with arbitrary non-Gaussianities.
Then we derived an asymptotic formula for the PBH mass fraction at the time of
formation in the high-peaks limit. The asymptotic formula for the threshold
model for the PBH formation has been known. In this paper, we specifically
derived the asymptotic formula for the peaks model and found that the
enhancement factor is identical in both models if the non-Gaussianities are
the same.

Next, we focused on the effect of skewness and kurtosis on the abundance for a
few specific models of non-Gaussianity. The skewness and kurtosis are
calculated from the bispectrum and trispectrum. For the bispectrum, we
considered non-Gaussian models which are described by superpositions of
specific forms of the bispectrum given by Eq. (3.10), which include popular
models such as the local-type, equilateral-type, folded-type and orthogonal-
type non-Gaussianities. The integral formulas to calculate the skewness
parameter in this series of models were derived in Eq. (3.21) – (3.24). For
the trispectrum, we considered a local-type model given by Eq. (3.25). The
integral formulas of kurtosis parameter were derived in Eqs. (3.29) and
(3.30).

The integral formulas to obtain skewness and kurtosis parameters are
numerically not so difficult to evaluate. However, to obtain a clear,
intuitive understanding of the results, we considered the case when the
primordial curvature perturbation power spectrum is sharply peaked at a scale
$k_{0}$ with the width $\Delta k=\epsilon k_{0}$ where $\epsilon$ ($<1$) is
the narrowness parameter. As a result, the precise shape of the spectrum gives
only subdominant contributions to the resulting skewness and kurtosis
parameters, as shown in Eqs. (4.3), (4.7), (4.13) and (4.14). In the limit
$\epsilon\to 0$, we obtained an asymptotic formula for the PBH mass fraction,
Eq. (4.15) for various types of non-Gaussianities of the bispectrum and
trispectrum. In particular, we found that the non-Gaussianity parameters for
various types of the bispectrum are linearly combined to give an effective
$f_{\rm NL}^{\rm eff}$ that determines the PBH mass fraction.

In this paper, we simply assume a linear relation between the curvature
perturbation and the density contrast in Eq. (2.4), and also assume formation
models of PBHs with threshold and peaks criteria with a fixed value of the
threshold. These assumptions are employed in order to make it possible for us
to analytically treat the problem. It is not obvious that these simplistic
assumptions hold in realistic situations [18, 21, 49]. However, since our non-
Gaussian corrections in the high-peaks limit change the abundance only through
a multiplicative factor, as given by Eq. (2.20), whose form is common to the
two different formation models, namely, the threshold model and the peaks
model. Therefore, one may expect that this non-Gaussian factor in the high-
peaks limit is universal, not depending on the details of the formation models
that can be more complicated than those we assume in this paper.
Unfortunately, however, it is difficult to prove this expectation within the
context of this paper, as the derivation of our formula in the case of the
peaks model, for instance, is already complicated enough. Therefore, we leave
this issue for future studies.

We systematically studied the effect of non-Gaussianities on the mass fraction
of PBHs at the time of formation at the level of bispectrum and trispectrum.
We only considered the one-point distribution function of the curvature
perturbation, we have no clue about non-Gaussianity effects on the spatial
distribution of PBHs. The effects of bispectrum and trispectrum on the PBH
clustering are investigated and discussed in Refs. [23, 24, 25]. It seems
there has been no further systematic analysis of non-Gaussianity effects on
the spatial distribution of PBHs in general situations. This issue certainly
deserves further studies.

## Acknowledgments

We thank Chris Byrnes for useful conversations at the early stage of this
work. This work was supported by JSPS KAKENHI Grants Nos. JP19K03835 (TM),
21H03403 (TM), 19H01895 (MS), 20H04727(MS), and 20H05853(MS).

##  Appendix A High-peaks limit of the peaks model with non-Gaussianity

In this Appendix, we derive the formula of Eq. (2.19) for the mass fraction in
the high-peaks limit of the peaks model with primordial non-Gaussianity. We
rely on the fact that the number of peaks above a threshold $\nu$
asymptotically approach to the Euler characteristic of the three-dimensional
body of the regions where $\varDelta_{R}>\nu\sigma$ is satisfied. The expected
value of Euler characteristic in non-Gaussian fields above a threshold $\nu$
is formally given by [46, 47, 48]

$n_{\chi}(\nu)=\left\langle\exp\left(\sum_{n=3}^{\infty}\frac{1}{n!}\sum_{\mu_{1},\cdots,\mu_{n}}M^{(n)}_{\mu_{1}\cdots\mu_{n}}\frac{\partial^{n}}{\partial
A_{\mu_{1}}\cdots\partial
A_{\mu_{n}}}\right)F_{\chi}(\bm{A},\nu)\right\rangle_{\mathrm{G}}.$ (A.1)

Various quantities in this equation are defined below in order. The
10-dimensional vector of variables $\bm{A}$ is composed of normalized
variables,

$\bm{A}=\left(\alpha,y_{i},Z_{ij}\right),$ (A.2)

with $1\leq i\leq j\leq 3$, where

$\alpha=\frac{\varDelta_{R}}{\sigma_{0}},\quad
y_{i}=\sqrt{3}\,\frac{\partial_{i}\varDelta_{R}}{\sigma_{1}},\quad
Z_{ij}=\frac{3}{\gamma}\frac{\partial_{i}\partial_{j}\varDelta_{R}}{\sigma_{2}}+\delta_{ij}\frac{\varDelta_{R}}{\sigma_{0}},$
(A.3)

$\partial_{i}=\partial/\partial x_{i}$ is a coordinate derivative, and

$\gamma=\frac{{\sigma_{1}}^{2}}{\sigma_{0}\sigma_{2}}$ (A.4)

is a spectral parameter, and $\sigma_{j}$ is defined by Eq. (2.13). The
$n$-point cumulants of the variables $\bm{A}$ are denoted by

$M^{(n)}_{\mu_{1}\cdots\mu_{n}}=\left\langle A_{\mu_{1}}\cdots
A_{\mu_{n}}\right\rangle_{\mathrm{c}}.$ (A.5)

A local function for the number density of the Euler characteristic is given
by

$F_{\chi}(\bm{A},\nu)=\left(\frac{\sigma_{1}}{\sqrt{3}\sigma_{0}}\right)^{3}\Theta(\alpha-\nu)\delta_{\mathrm{D}}^{3}(\bm{y})\det(\alpha
I-Z),$ (A.6)

where $\Theta(x)$ is a step function and $\delta_{\mathrm{D}}^{3}(\bm{y})$ is
the 3-dimensional Dirac’s delta function, and $I$ is the $3\times 3$ unit
matrix. Finally, we define a Gaussian average by

$\langle\cdots\rangle_{\mathrm{G}}=\frac{1}{(2\pi)^{5}\sqrt{\det M}}\int
d^{10}A\cdots\exp\left(-\frac{1}{2}\bm{A}^{\mathrm{T}}{M}^{-1}\bm{A}\right),$
(A.7)

where $M=(M^{(2)}_{\mu_{1}\mu_{2}})$ is a $10\times 10$ matrix of two-point
cumulants. The two-point cumulants are explicitly given by

$\displaystyle\left\langle\alpha^{2}\right\rangle=1,\quad\left\langle\alpha\eta_{i}\right\rangle=0,\quad\left\langle\alpha
Z_{ij}\right\rangle=0,\quad\left\langle
y_{i}y_{j}\right\rangle=\delta_{ij},\quad\left\langle\eta_{i}Z_{jk}\right\rangle=0,$
$\displaystyle\left\langle
Z_{ij}Z_{kl}\right\rangle=\frac{3}{5\gamma^{2}}\left(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)-\delta_{ij}\delta_{kl}.$
(A.8)

Expanding the operators in the exponent of Eq. (A.1), we have the following
type of factor in the expansion:

$\left\langle\frac{\partial^{m_{0}}}{\partial\alpha^{m_{0}}}\frac{\partial^{m_{1}}}{\partial
y_{i_{1}^{\prime}}\cdots\partial
y_{i_{m_{1}}^{\prime}}}\frac{\partial^{m_{2}}}{\partial
Z_{i_{1}j_{1}}\cdots\partial
Z_{i_{m_{2}}j_{m_{2}}}}F_{\chi}\right\rangle_{\mathrm{G}}.$ (A.9)

This factor is equivalent to [48]

$\left(\frac{\sigma_{1}}{\sqrt{3}\sigma_{0}}\right)^{3}\int_{\nu}^{\infty}\\!d\nu\,\left(-\frac{d}{d\nu}\right)^{m_{0}}\left[\left\langle\delta_{\mathrm{D}}(\alpha-\nu)\right\rangle_{\mathrm{G}}\left\langle\frac{\partial^{m_{1}}\delta_{\mathrm{D}}^{3}(\bm{y})}{\partial
y_{i_{1}^{\prime}}\cdots\partial
y_{i_{m_{1}}^{\prime}}}\right\rangle_{\mathrm{G}}\left\langle\frac{\partial^{m_{2}}\det(\nu
I-Z)}{\partial Z_{i_{1}j_{1}}\cdots\partial
Z_{i_{m_{2}}j_{m_{2}}}}\right\rangle_{\mathrm{G}}\right],$ (A.10)

because the three kinds of variables $\alpha$, $y_{i}$ and $Z_{ij}$ are
mutually independent in the Gaussian averages with two-point cumulants. The
Gaussian distribution functions of the variables $\alpha$ and $y_{i}$ are
given by

$P^{(0)}_{\mathrm{G}}(\alpha)=\frac{e^{-\alpha^{2}/2}}{\sqrt{2\pi}},\quad
P^{(1)}_{\mathrm{G}}(\bm{y})=\frac{e^{-|\bm{y}|^{2}/2}}{(2\pi)^{3/2}},$ (A.11)

and the first two factors in the square bracket reduce to

$\left\langle\delta_{\mathrm{D}}(\alpha-\nu)\right\rangle_{\mathrm{G}}=\frac{e^{-\nu^{2}/2}}{\sqrt{2\pi}},\quad\left\langle\frac{\partial^{m_{1}}\delta_{\mathrm{D}}^{3}(\bm{y})}{\partial
y_{i_{1}^{\prime}}\cdots\partial
y_{i_{m_{1}}^{\prime}}}\right\rangle_{\mathrm{G}}=\begin{cases}\displaystyle\frac{(-1)^{m_{1}/2}(m_{1}-1)!!}{(2\pi)^{3/2}}\delta_{(i_{1}^{\prime}i_{2}^{\prime}}\cdots\delta_{i_{m_{1}-1}^{\prime}i_{m_{1}}^{\prime})},&(m_{1}:\mathrm{even}),\\\
0,&(m_{1}:\mathrm{odd}),\end{cases}$ (A.12)

where parentheses in the subscript of Kronecker’s delta indicate the
symmetrization of the indices inside them. For the last factor, we have an
identity for the $3\times 3$ matrices,

$\det(\nu
I-Z)=\nu^{3}-\nu^{2}\,\mathrm{tr}\,Z+\frac{\nu}{2}\left[\left(\mathrm{tr}\,Z\right)^{2}-\mathrm{tr}\,\left(Z^{2}\right)\right]-\det
Z.$ (A.13)

For a fixed number of derivatives in Eq. (A.9), $m\equiv m_{0}+m_{1}+m_{2}$,
one can see that the dominant contribution in the high-peaks limit,
$\nu\rightarrow\infty$, is given by a case of $m_{0}=m$, $m_{1}=m_{2}=0$. In
this case, we have $\langle\det(\nu I-Z)\rangle_{\mathrm{G}}=H_{3}(\nu)$ due
to Eqs. (A.8) and (A.13), and Eq. (A.10) reduces to

$\frac{1}{(2\pi)^{2}}\left(\frac{\sigma_{1}}{\sqrt{3}\sigma_{0}}\right)^{3}e^{-\nu^{2}/2}H_{m+2}(\nu)$
(A.14)

in the high-peaks limit. When the replacement $H_{m+2}(\nu)\rightarrow
H_{2}(\nu)\nu^{m}$ in the same limit is applied for each term of the expanded
series of Eq. (A.1), the resulting infinite series is resummed again. Thus we
have

$n_{\chi}(\nu)\approx\frac{1}{(2\pi)^{2}}\left(\frac{\sigma_{1}}{\sqrt{3}\sigma_{0}}\right)^{3}e^{-\nu^{2}/2}H_{2}(\nu)\exp\left(\sum_{n=3}^{\infty}\frac{\nu^{n}}{n!}\left\langle\alpha^{n}\right\rangle_{\mathrm{c}}\right).$
(A.15)

Identifying the number density of peaks with the Euler number density in the
high-peaks limit, the mass fraction of peaks is given by

$\beta_{\mathrm{pk}}\approx(2\pi)^{3/2}R^{3}n_{\chi}(\nu)\approx\frac{1}{\sqrt{2\pi}}\left(\frac{R\sigma_{1}}{\sqrt{3}\sigma}\right)^{3}(\nu^{2}-1)e^{-\nu^{2}/2}\exp\left(\sum_{n=3}^{\infty}\frac{\nu^{n}}{n!}\left\langle\alpha^{n}\right\rangle_{\mathrm{c}}\right).$
(A.16)

The cumulants of $\alpha$ are given by
$\langle\alpha^{n}\rangle_{\mathrm{c}}=\sigma^{n-2}S_{n}$ where $S_{n}$ is the
reduced cumulants defined by Eq. (2.15), and the above equation is equivalent
to Eq. (2.20).

## References

  * [1] S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971)
  * [2] B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399-415 (1974)
  * [3] B. J. Carr, Astrophys. J. 201, 1-19 (1975)
  * [4] B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355-394 (2020) [arXiv:2006.02838 [astro-ph.CO]].
  * [5] B. Carr and F. Kuhnel, SciPost Phys. Lect. Notes 48, 1 (2022) [arXiv:2110.02821 [astro-ph.CO]].
  * [6] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, Rept. Prog. Phys. 84, no.11, 116902 (2021) [arXiv:2002.12778 [astro-ph.CO]].
  * [7] A. M. Green, Fundam. Theor. Phys. 178, 129-149 (2015) [arXiv:1403.1198 [gr-qc]].
  * [8] C. T. Byrnes, E. J. Copeland and A. M. Green, Phys. Rev. D 86, 043512 (2012) [arXiv:1206.4188 [astro-ph.CO]].
  * [9] J. S. Bullock and J. R. Primack, Phys. Rev. D 55, 7423-7439 (1997) [arXiv:astro-ph/9611106 [astro-ph]].
  * [10] P. Pina Avelino, Phys. Rev. D 72, 124004 (2005) [arXiv:astro-ph/0510052 [astro-ph]].
  * [11] D. Seery and J. C. Hidalgo, JCAP 07, 008 (2006) [arXiv:astro-ph/0604579 [astro-ph]].
  * [12] S. Shandera, A. L. Erickcek, P. Scott and J. Y. Galarza, Phys. Rev. D 88, no.10, 103506 (2013) [arXiv:1211.7361 [astro-ph.CO]].
  * [13] S. Young and C. T. Byrnes, JCAP 08, 052 (2013) [arXiv:1307.4995 [astro-ph.CO]].
  * [14] S. Young, D. Regan and C. T. Byrnes, JCAP 02, 029 (2016) [arXiv:1512.07224 [astro-ph.CO]].
  * [15] V. Atal, J. Garriga and A. Marcos-Caballero, JCAP 09, 073 (2019) [arXiv:1905.13202 [astro-ph.CO]].
  * [16] C. M. Yoo, J. O. Gong and S. Yokoyama, JCAP 09, 033 (2019) [arXiv:1906.06790 [astro-ph.CO]].
  * [17] M. Taoso and A. Urbano, JCAP 08, 016 (2021) [arXiv:2102.03610 [astro-ph.CO]].
  * [18] F. Riccardi, M. Taoso and A. Urbano, JCAP 08, 060 (2021) [arXiv:2102.04084 [astro-ph.CO]].
  * [19] N. Kitajima, Y. Tada, S. Yokoyama and C. M. Yoo, JCAP 10, 053 (2021) [arXiv:2109.00791 [astro-ph.CO]].
  * [20] S. Pi and M. Sasaki, [arXiv:2112.12680 [astro-ph.CO]].
  * [21] S. Young, JCAP 05, no.05, 037 (2022) [arXiv:2201.13345 [astro-ph.CO]].
  * [22] A. Kehagias, I. Musco and A. Riotto, JCAP 12, 029 (2019)
  * [23] Y. Tada and S. Yokoyama, Phys. Rev. D 91, no.12, 123534 (2015)
  * [24] T. Suyama and S. Yokoyama, PTEP 2019, no.10, 103E02 (2019) [arXiv:1906.04958 [astro-ph.CO]].
  * [25] T. Matsubara, T. Terada, K. Kohri and S. Yokoyama, Phys. Rev. D 100, no.12, 123544 (2019) [arXiv:1909.04053 [astro-ph.CO]].
  * [26] A. R. Liddle and D. H. Lyth, “Cosmological inflation and large scale structure,” (Cambridge University Press, Cambridge, England, 2000).
  * [27] A. M. Green, A. R. Liddle, K. A. Malik and M. Sasaki, Phys. Rev. D 70, 041502 (2004) [arXiv:astro-ph/0403181 [astro-ph]].
  * [28] S. Young, C. T. Byrnes and M. Sasaki, JCAP 07, 045 (2014) [arXiv:1405.7023 [gr-qc]].
  * [29] W. H. Press and P. Schechter, Astrophys. J. 187, 425-438 (1974).
  * [30] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, Astrophys. J. 304, 15-61 (1986).
  * [31] I. Musco, V. De Luca, G. Franciolini and A. Riotto, Phys. Rev. D 103, no.6, 063538 (2021) [arXiv:2011.03014 [astro-ph.CO]].
  * [32] I. Musco, Phys. Rev. D 100, no.12, 123524 (2019) [arXiv:1809.02127 [gr-qc]].
  * [33] S. Matarrese, F. Lucchin and S. A. Bonometto, Astrophys. J. Lett. 310, L21-L26 (1986)
  * [34] G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, JCAP 03, 016 (2018) [arXiv:1801.09415 [astro-ph.CO]].
  * [35] A. Gangui, F. Lucchin, S. Matarrese and S. Mollerach, Astrophys. J. 430, 447-457 (1994) [arXiv:astro-ph/9312033 [astro-ph]].
  * [36] L. Verde, L. M. Wang, A. Heavens and M. Kamionkowski, Mon. Not. Roy. Astron. Soc. 313, L141-L147 (2000) [arXiv:astro-ph/9906301 [astro-ph]].
  * [37] E. Komatsu and D. N. Spergel, Phys. Rev. D 63, 063002 (2001) [arXiv:astro-ph/0005036 [astro-ph]].
  * [38] P. Creminelli, A. Nicolis, L. Senatore, M. Tegmark and M. Zaldarriaga, JCAP 05, 004 (2006) [arXiv:astro-ph/0509029 [astro-ph]].
  * [39] P. D. Meerburg, J. P. van der Schaar and P. S. Corasaniti, JCAP 05, 018 (2009) [arXiv:0901.4044 [hep-th]].
  * [40] L. Senatore, K. M. Smith and M. Zaldarriaga, JCAP 01, 028 (2010) [arXiv:0905.3746 [astro-ph.CO]].
  * [41] C. T. Byrnes, M. Sasaki and D. Wands, Phys. Rev. D 74, 123519 (2006) [arXiv:astro-ph/0611075 [astro-ph]].
  * [42] L. Boubekeur and D. H. Lyth, Phys. Rev. D 73, 021301 (2006) [arXiv:astro-ph/0504046 [astro-ph]].
  * [43] T. Suyama and M. Yamaguchi, Phys. Rev. D 77, 023505 (2008) [arXiv:0709.2545 [astro-ph]].
  * [44] T. Matsubara, C. Hikage and S. Kuriki, Phys. Rev. D 105, no.2, 023527 (2022) [arXiv:2012.00203 [astro-ph.CO]].
  * [45] C. M. Yoo, T. Harada, J. Garriga and K. Kohri, PTEP 2018, no.12, 123E01 (2018) [arXiv:1805.03946 [astro-ph.CO]].
  * [46] T. Matsubara, Astrophys. J. Lett. 434, L43 (1994) [arXiv:astro-ph/9405037 [astro-ph]].
  * [47] T. Matsubara, Astrophys. J. 457, 13-17 (1996) [arXiv:astro-ph/9501055 [astro-ph]].
  * [48] T. Matsubara and S. Kuriki, Phys. Rev. D 104, no.10, 103522 (2021) [arXiv:2011.04954 [astro-ph.CO]].
  * [49] I. Musco, J. C. Miller and A. G. Polnarev, Class. Quant. Grav. 26, 235001 (2009) [arXiv:0811.1452 [gr-qc]].