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

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)

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].

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.

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.

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)

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)$.