Effect of nonlinearity between density and curvature perturbations on the primordial black hole formation

TL;DR

This paper investigates how nonlinear relationships between density and curvature perturbations influence primordial black hole formation, revealing a slight suppression effect comparable to local non-Gaussianity with f_NL around -1.

Contribution

It introduces a detailed analysis of nonlinear effects on PBH formation, quantifies non-Gaussianity, and estimates the suppression of PBH abundance during the Radiation era.

Findings

Nonlinear effects slightly suppress PBH formation.

Larger curvature perturbations are needed due to nonlinearity.

Suppression is comparable to local non-Gaussianity with f_NL ~ -1.

Abstract

We study the effect of the nonlinear relation between density and curvature perturbations on the formation of PBHs. By calculating the variance and skewness of the density perturbation we derive the non-Gaussian property. As a criterion for PBH formation, the compaction function is used and it is found that larger curvature perturbations are required due to the nonlinear effect. We estimate the PBH abundance based on the Press-Schechter formalism with non-Gaussian probability density function during Radiation dominated era. It is found that the nonlinear effect slightly suppresses the PBH formation and the suppression is comparable to that expected if the primordial curvature perturbation would have the local form of non-Gaussianity with nonlinear parameter $f_\text{NL} \sim -1$.

Figures

Equations

$${\text{d}}^{2}s=-({\alpha }^2-{\psi }^2{a}^2{\beta }^2{r}^2)\text{d}{t}^2+2{\psi }^4{a}^2\beta r\text{d}r\text{d}t+{\psi }^4{a}^2(\text{d}{r}^2+r^2\text{d}{\Omega }^2),$$

$$\mathcal{C}(r)$$

$$\delta (r)=-\frac{4(1+w)}{5+3w}\frac{\Delta (e^{\zeta /2})}{e^{5\zeta /2}}{\left(\frac{1}{a{H}_b}\right)}^2+\mathcal{O}({ϵ}^3)$$

$$\bar{\delta }=\frac{4\pi {\int }_0^{r}dxR(x{)}^2{R}^{\prime }(x)\delta (x)}{4\pi R(r{)}^3/3},$$

$$\mathcal{C}(r)=$$

$${\delta }_k=\frac{4}{9}{\left(\frac{k}{a{H}_b}\right)}^2{\zeta }_k+\mathcal{O}({\zeta }^2).$$

$$\bar{\delta }(\mathbf{x},\rho )=\int {\text{d}}^{3}yW(|\mathbf{x}-\mathbf{y}|;\rho )\delta (\mathbf{y})=\int \frac{{\text{d}}^{3}k}{(2\pi {)}^3}\tilde{W}(k\rho )e^{i\mathbf{k}\cdot \mathbf{x}}{\delta }_k$$

$$W(x;\rho )$$

$$\tilde{W}(k\rho )$$

$$T(k\eta )=$$

$$P_k(\bar{\delta })$$

$${\sigma }_k^2$$

$$\beta (k)$$

$$M$$

$$\simeq 30M_{\odot }\frac{\gamma }{0.2}{\left(\frac{g_{\ast }}{10.75}\right)}^{-1/6}{\left(\frac{k}{3.43\times 10^5{\mathrm{Mpc}}^{-1}}\right)}^{-2}$$

$$\simeq 30M_{\odot }\frac{\gamma }{0.2}{\left(\frac{g_{\ast }}{10.75}\right)}^{-1/2}{\left(\frac{T}{31.6\mathrm{MeV}}\right)}^{-2}.$$

$$f(M)$$

$${\bar{\delta }}_{\text{HC}}(r_m)=$$

$$=$$

$$\times \left[\Delta \zeta (r)+(\zeta (r)-\zeta (r_m))\Delta \zeta (r)+\frac{1}{2}{\zeta }^{\prime }(r)(5{\zeta }^{\prime }(r)+2r{\zeta }^{\prime \prime }(r))\right]+\mathcal{O}({ϵ}^3)$$

$$=$$

$$\zeta (\mathbf{x})=$$

$${\zeta }_{lm}(k)=$$

$$⟨{\zeta }_{l_1{m}_1}(k_1){\zeta }_{l_2{m}_2}(k_2{)}^{\ast }⟩=$$

$$\sigma (r_m{)}^2$$

$${\mu }_3(r_m)$$

$${\bar{\delta }}_{\text{HC}}^{(1)}(r_m)$$

$$=\frac{4}{9}\int \text{d}\ln (r_{m}p)(r_{m}p{)}^4\tilde{W}(r_{m}p)\frac{r_m^{-2}{\zeta }_{00}(p)}{\sqrt{2}\pi },$$

$${\bar{\delta }}_{\text{HC}}^{(2)}(r_m)$$

$$=\frac{-4}{9}{\int }_0^{\infty }\text{d}(p_1{r}_m)\text{d}(p_2{r}_m)\frac{r_m^{-4}}{2{\pi }^2}{\zeta }_{00}(p_1){\zeta }_{00}(p_2)$$

$$\times \frac{3}{2}\left[\sin (p_1{r}_m)-(p_1{r}_m)\cos (p_1{r}_m)\right]\left[\sin (p_2{r}_m)-(p_2{r}_m)\cos (p_2{r}_m)\right],$$

$$\sigma (r_m{)}^2$$

$$={\left(\frac{4}{9}\right)}^2\int \text{d}(\log r_{m}p)|\tilde{W}(r_{m}p){|}^2(r_{m}p{)}^4|T(r_{m}p){|}^2{\mathcal{P}}_{\zeta }(p)+\mathcal{O}({\zeta }^4).$$

$${\mu }_3(r_m)$$

$$=-9{\sigma }^{-3}{\left(\frac{4}{9}\right)}^3\int \text{d}(p_1{r}_m)\text{d}(p_2{r}_m)\tilde{W}(p_1{r}_m)\tilde{W}(p_2{r}_m)|T(p_1{r}_m){|}^2{\mathcal{P}}_{\zeta }(p_1)|T(p_2{r}_m){|}^2{\mathcal{P}}_{\zeta }(p_2)$$

$$\times \left(\sin (p_1{r}_m)-(p_1{r}_m)\cos (p_1{r}_m)\right)\left(\sin (p_2{r}_m)-(p_2{r}_m)\cos (p_2{r}_m)\right)$$

$$=-\frac{9}{4}\sigma (r_m)$$

$${\bar{\delta }}_{\text{HC}}[{\delta }_g]={\delta }_g+\frac{{\mu }_3(r_m)}{6\sigma (r_m)}({\delta }_g^2-{\sigma }^2),$$

$$P_g({\delta }_g)=\frac{1}{\sqrt{2\pi }\sigma (r_m)}\exp \left(-\frac{1}{2}\frac{{\delta }_g^2}{\sigma (r_m{)}^2}\right).$$

Full text

††preprint: IPMU19-0029

1. I Introduction
2. II PBH formation criterion and Press-Schechter formalism
3. III Effect of nonlinearity
4. IV Numerical results
5. V Conclusion

I Introduction

Primordial black holes (PBHs) has been attracting much interest since the LIGO-VIRGO collaboration succeeded in direct detection of the gravitational waves generated by mergers of binary black holes with masses around 30 ${M_{\odot}}$. Up to now about ten merger events have been found, from which the event rate is estimated as $53.2^{+58.5}_{-28.8}\mathrm{Gpc}^{-3}\text{yr}^{-1}$ Abbott et al. (2018). It is pointed out that PBHs with the abundance of $\Omega_{\text{PBH}}/\Omega_{\text{DM}}\sim 10^{-3}$ can be a promissing candidate for the merger events Bird et al. (2016); Clesse and García-Bellido (2017); Sasaki et al. (2016). PBHs are also an interesting candidate for the dark matter (DM). Although the available mass range for PBH to be DM is stringently constrained by the observations of $\gamma$-ray background, micro-lensing, CMB and others Carr et al. (2010); Allsman et al. (2001); Tisserand et al. (2007); Wyrzykowski et al. (2011); Niikura et al. (2017); Ali-Haïmoud and Kamionkowski (2017); Graham et al. (2015), PBHs with mass $M_{\text{PBH}}\sim 10^{-13}{M_{\odot}}(\simeq 10^{20}\text{g})$ and $10^{-15}{M_{\odot}}(\simeq 10^{18}\text{g})$ can account for all DM of the universe becasue they avoid the micro(femto)-lensing constraint due to the wave effect and finite source size effect Katz et al. (2018); Niikura et al. (2017).

PBHs are also a unique tool to study the small scale perturbations and their generation in the early universe. The formation of PBHs requires large density fluctuations of $\mathcal{O}(0.1)$ compared with $\mathcal{O}(10^{-5})$ at CMB scales, which is hardly realized in simple single-field inflation models and forces one to consider some mechanism to enhance small-scale density fluctuations such as hybrid inflation Bugaev and Klimai (2012), double inflation Inomata et al. (2018) and curvaton models Kohri et al. (2013); Ando et al. (2017). PBHs are formed when overdense regions enter the horizon and gravitationally collapse Carr (1975). The criterion of over-density required for PBH production was studied by the numerical relativistic simulations Shibata and Sasaki (1999); Polnarev and Musco (2007). These simulations derived the threshold values on the curvature perturbation or the density perturbation using different gauge conditions. Harada et. al. Harada et al. (2015) compared those simulations and showed equivalence between their results. In Young et al. (2014); Harada et al. (2015) it was also pointed out that the threshold value of the curvature perturbation depends on environmental effects i.e. density profiles in the surrounding regions, which suggests that use of the averaged density perturbation is more suitable for the criterion than the curvature perturbation.

From point of view of model building, what we usually predict is the curvature perturbation power spectrum $\mathcal{P}_{\zeta}(k)$. To obtain the PBH abundance, one needs to evaluate the abundance of the overdense regions from the power spectrum. The Press–Schechter formalism Press and Schechter (1974) is mostly used to count the number of overdense regions. Other estimation based on the peak theory Bardeen et al. (1986) is also studied Green et al. (2004); Yoo et al. (2018). In the Press-Schechter formalism, we calculate variance $\sigma^{2}$ of the coarse-grained (averaged) density perturbation using a linear relation between the curvature and density perturbations, and then estimate the abundance of the overdense regions assuming that the coarse-grained density perturbation follows the Gaussian distribution with variance $\sigma^{2}$. However, the relation between the curvature and density perturbations is non-linear, which could produce non-Gaussianity in the density perturbations.

Although the non-Gaussianity is suppressed by the slow-roll parameters in the standard cosmology with the single field canonical inflation Maldacena (2003), some models like curvaton scenarios generate the perturbations with large non-Gaussianity. Since the non-Gaussianity changes the probability distribution of the perturbation, it can help or suppress the PBH formation. Thus, there exist extensive researches considering the effect of non-Gaussianity on the PBH formation Pina Avelino (2005); Byrnes et al. (2012); Young and Byrnes (2013); Sherkatghanad and Brandenberger (2015); Young et al. (2016); Franciolini et al. (2018).

In this paper, we revisit the relation between the curvature perturbation power spectrum and the PBH formation. In particular, we study the effects of the nonlinear relation between the curvature perturbation $\zeta$ and density perturbation $\delta$ on the PBH formation. Although the most previous calculations used the linear relation $\delta\propto\Delta\zeta$, the linear relation is not valid for large curvature perturbations. The nonlinear relation can be obtained using the gradient expansion [see Eq. (3)]. Thus, taking the nonlinearity into account, we estimate the threshold value on the average density perturbation from the result of the numerical simulations. Here as the criterion for the PBH formation, we use the Compaction function $\mathcal{C}(r)$ introduced in Ref. Shibata and Sasaki (1999). We also calculate the variance and skewness of the density perturbations and estimate the non-Gaussianity. We thus obtain the PBH abundance based on the Press-Schechter formalism.

In the next section II, we review the criterion of the PBH formation and estimate the formation rate using the Press–Schechter formalism and the linear density perturbation.. In Sec. III, the nonlinear effect on the PBH formation rate is analytically calculated. We use the spherical harmonics expansion of the curvature perturbation and evaluate the variance and skewness of averaged density perturbation. In Sec. IV, we numerically estimate the nonlinear effect for a simple power spectrum. Sec V is conclusion of this paper.

II PBH formation criterion and Press-Schechter formalism

When an overdensity region enters the horizon, a PBH is formed if the density perturbation grows gravitationally and overcomes the pressure. One can evaluate the threshold value of the overdensity for the PBH formation by performing numerical relativity simulations Shibata and Sasaki (1999); Polnarev and Musco (2007). In this section, we review the PBH formation criterion obtained by the results of the numerical relativity simulations following Ref. Harada et al. (2015).

In order to calculate the superhorizon dynamics of the metric and energy-momentum tensors beyond the linear approximation, we use the long wavelength approximation (the gradient expansion), where the Einstein equations are expanded by the space derivative. We introduce $\epsilon$ by $\partial_{i}\to\epsilon\partial_{i}$ and expand the equations in small $\epsilon$ on the flat background FLRW space. Since the superhorizon perturbations mainly contribute to the PBH formation, the gradient expansion is appropriate to the study of the formation criterion.

The well-known criterion of the PBH formation is that the averaged density fluctuation over the horizon exceeds a certain threshold value. Let us derived the threshold value following Shibata and Sasaki (1999); Harada et al. (2013, 2015); Yoo et al. (2018). We consider the spherically symmetric curved space-time in the conformally flat spatial coordinates. The metric with the polar coordinate $(r,\theta,\varphi)$ is written as

[TABLE]

where $a(t)$ is the scale factor, $\text{d}\Omega^{2}=\text{d}\theta^{2}+\sin^{2}\theta d\varphi^{2}$, $\alpha$ and $\beta$ are the lapse function and the radial component of the shift vector. In the comoving gauge, the comoving curvature perturbation $\zeta$ is given by $\psi=\exp(\zeta/2)$.111In our definition, the density perturbation $\delta$ and the curvature perturbation $\zeta$ has the same sign at linear order in momentum space. The areal radius is given by $R(r)=\psi(r)^{2}ar=e^{\zeta(r)}ar$. Shibata and Sasaki Shibata and Sasaki (1999) introduced the compaction function $\mathcal{C}(r)$ defined by

[TABLE]

where $\rho_{b}$ is the energy density of the background universe, $\delta$ is the density perturbation ($\delta=(\rho-\rho_{b})/\rho_{b}$) and ′ denotes $d/dr$. In this paper, the subscript ”$b$” denotes the background value. In the long-wave approximation the relation between the comoving density perturbation $\delta$ and the curvature perturbation $\zeta$ is written as Harada et al. (2015)

[TABLE]

where $w=P/\rho$ is the equation of state and $H_{b}$ is the background Hubble parameter. In this paper, we focus on the radiation dominated era ($w=1/3$).

Since $R\propto a$ and $\delta M\propto\rho_{b}R(aH_{b})^{-2}\propto a$, $\mathcal{C}(r)$ is independent of $a$ and $H_{b}$. $\mathcal{C}(r)$ can be calculated using the initial profile of $\zeta(r)$ since it is frozen in the super horizon. For a given profile, $\mathcal{C}(r)$ has the maximum value $\mathcal{C}_{m}$ at some point $r=r_{m}$. In the constant mean curvature slicing, the numerical simulation gives the threshold value of the PBH formation as $\mathcal{C}_{m;\text{CMC}}^{th}\simeq 0.4$ during the radiation dominated eta Shibata and Sasaki (1999). Once the initial profile gives $\mathcal{C}_{m}$ larger than the threshold value, we assume that the region collapses into the PBH. Some uncertainty on the threshold value comes from the different shapes of the initial profiles Shibata and Sasaki (1999). Detailed analysis about the shape of the profile was performed in Harada et al. (2015); Nakama et al. (2014) and it was shown that the sharp profiles tend to require large overdensities to collapse. In this paper, we only consider the condition on the averaged density, for simplicity. In the comoving gauge, the threshold value is given by $\mathcal{C}_{m;\text{Com}}^{th}=\frac{2}{3}\mathcal{C}_{m;\text{CMC}}^{th}\simeq 0.27$ Harada et al. (2015).

The averaged density fluctuation $\bar{\delta}$ is defined as

[TABLE]

which is related to the compaction function $\mathcal{C}(r)$ as Harada et al. (2015)

[TABLE]

Since the density perturbation depends on the cosmic time $t$, we should specify the time to evaluate $\bar{\delta}$. The usual choice is the horizon crossing time $H_{b}R(r_{m})=H_{b}ar_{m}e^{\zeta(r_{m})}=1$, which depends on the curvature profile itself. At the horizon crossing time, the threshold value of the density perturbation is given by $\bar{\delta}^{th}=2\mathcal{C}^{th}_{m;\text{Com}}\simeq 0.53$ which is adopted in our paper.

Now we estimate the PBH formation rate in the linear order calculation based on the Press–Schechter formalism Press and Schechter (1974). The density perturbation is related to $\zeta$ in momentum space as

[TABLE]

The averaged density fluctuation over the radius $\rho$ is given by

[TABLE]

where $W(r;\rho)$ and $\tilde{W}(k\rho)$ is the window functions in the real and momentum space. It is pointed out that the choice of the window function leads to large uncertainties Ando et al. (2018). When we use the average density fluctuation defined by Eq. (4), we should use the real space top-hat window,

[TABLE]

In the subhorizon, the perturbations start to oscillate and their amplitude decreases. The time evolution of the perturbations is described by the transfer function given by

[TABLE]

which is used together with the real space top hat window function because high momentum modes give a significant contribution to $\bar{\delta}(\mathbf{x},\rho)$.222For the Gaussian window function, which is often used in the literature, high momentum modes are exponentially suppressed and hence the effect of the transfer function is negligible. Assuming that $\zeta$ is Gaussian variable, the probability density function of $\bar{\delta}$ is given by

[TABLE]

In the linear order approximation, the horizon crossing condition is given by $aH_{b}=k$, which is slightly different from the nonlinear one. The PBH formation rate is given by

[TABLE]

The mass of PBHs is nearly equal to the horizon mass when the PBHs are formed Carr (1975). Thus, in the radiation dominated era, the PBH mass is written as

[TABLE]

where $g_{*}$ is relativistic degree of freedom, $\rho_{b}(t_{\text{HC}})$ is the background radiation energy density at horizon crossing and $\gamma$ is the ratio of PBH mass to the horizon mass. In this paper, we take the analytical estimation $\gamma=3^{-3/2}\simeq 0.2$ Carr (1975). The mass spectrum is then given by

[TABLE]

III Effect of nonlinearity

Now we consider the effect of the nonliear relation between $\bar{\delta}$ (or $\mathcal{C}$) and $\zeta$ found in Eqs. (3), (4) and (5). In the Press–Schechter formalism, we evaluate the probability density function of the averaged density for a given horizon scale $r_{m}$. When the averaged density exceeds $\bar{\delta}^{th}$, the region collapses into a PBH.333Although $r_{m}$ is defined for each initial profile in numerical simulations, we now take $r_{m}$ as horizon scale to perform the statistical estimation. This picture has the double counting problem when $\mathcal{C}(r)$ has many maximum points. The smooth window function or the transfer function helps to avoid this problem by suppressing the subhorizon modes. At horizon crossing $H_{b}R(r_{m})=1$, the averaged density perturbation is calculated as Harada et al. (2015)

[TABLE]

where $\bar{\delta}_{\text{HC}}^{(n)}$ denotes the $\mathcal{O}(\zeta^{n})$ term. Thus, the averaged density is a nonlinear function of $\zeta$. When we substitute the threshold value $\bar{\delta}^{th}=0.53$, the linear order solution is given by $r_{m}\zeta^{\prime}(r_{m})\simeq-0.40$ and the exact one is $r_{m}\zeta^{\prime}(r_{m})\simeq-0.55$, which are differ by a factor about $1.4$. Fig.1 shows the contribution of $\bar{\delta}_{\text{HC}}^{(1)}(r_{m})$ and $\bar{\delta}_{\text{HC}}^{(2)}(r_{m})$ for $\zeta(r)=\zeta_{0}e^{-r^{2}/2}$ at $r_{m}=\sqrt{2}$. The non-linearity slightly suppresses $\bar{\delta}_{\text{HC}}$ compared to the linear order calculation $\bar{\delta}_{\text{HC}}^{(1)}(r_{m})$. This suggests that PBHs are difficult to form and hence the formation rate decreases due to non-linearity. We will confirm it using a non-Gaussian distribution function.

In the following section, we treat $\zeta$ and $\bar{\delta}_{\text{HC}}$ as stochastic variables and estimate the probability density function of $\bar{\delta}_{\text{HC}}$ to calculate the PBH formation rate Eq.(13). Since the criterion Eq.(16) is obtained assuming spherical symmetry, we use the spherical harmonics expansion Yoo and Desjacques (2013) as444Since PBH formation requires the highly curved space, the scale $k$ may depend on the curvature of each region. Here we simply assume that the Fourier transformation on the flat background space.

[TABLE]

where $j_{l}(kr)$ is the spherical Bessel function of the first kind, $Y_{lm}(\hat{\bm{x}})$ is the spherical harmonics and $\int\text{d}\Omega_{k}$ is the integration over the direction of $\bm{k}$. Due to spherically symmetry, only $\zeta_{00}(k)$ has a non zero value. In this paper, we assume that $\zeta_{00}(k)$ is a Gaussian variable, which is justified for slow-roll inflation Maldacena (2003). On the other hand, since $\bar{\delta}_{\text{HC}}(r_{m})$ is the nonlinear function of $\zeta$, $\bar{\delta}_{\text{HC}}$ follows the non-Gaussian distribution. Following the treatment of non-Gaussianity of the PBH formations Young et al. (2016); Franciolini et al. (2018), we evaluate the probability density function using the variance $\sigma$ and the skewness $\mu_{3}$ defined by Ezquiaga and García-Bellido (2018)

[TABLE]

At first, we simplify $\bar{\delta}_{\text{HC}}^{(1)}(r_{m})$ and $\bar{\delta}_{\text{HC}}^{(2)}(r_{m})$ using Eq.(20) as

[TABLE]

where $\tilde{W}(z)(=3(\sin z-z\cos z)/z^{3})$ is the Fourier transformation of the real space top-hat window function. Although the real space top-hat window function is derived from the definition of the compaction function, it may have some troubles. (1) Our calculation is only valid in small $k$ because of the long wavelength approximation. On the other hand, the averaged density with top-hat window function tends to include a significant contribution from large $k$ modes compared to the Gaussian window function or the momentum space top-hat window function. (2) As pointed out in Yoo et al. (2018) the contribution of large $k$ in the top-hat window function comes from the sharp edge in the real space. Thus, we need to suppress the high momentum modes, for example, multiplying the transfer function or using the Gaussian or the momentum space top-hat window function instead of the real space top-hat window function. In this paper, since we concentrate on the non-linearity of the compaction function, we continue to use the real space top-hat window function and multiply it by the transfer function, $\mathcal{P}_{\zeta}(p)\to|T(r_{m}p)|^{2}\mathcal{P}_{\zeta}(p)$ in an ad-hoc way.555Using monochromatic spectrum $\zeta(r)=A\sin(kr)/(kr)$, $\mathcal{C}(r)$ shows approximately periodic behavior and has many peaks. The transfer function helps to pick up only the first peak with the smallest $r_{m}$. Instead, we can use the step function to pick up the first peak, which results in the larger PBH abundance.

The variance Eq.(23) up to $\mathcal{O}(\zeta^{2})$ is written as

[TABLE]

Because we only include $\braket{\bar{\delta}_{\text{HC}}^{(1)2}(r_{m})}$, $\sigma^{2}$ reproduces the linear order calculation in Eq.(12). As for the skewness we calculate $\mu_{3}(r_{m})$ up to $\mathcal{O}(\zeta^{4})$,

[TABLE]

Since the linear and nonlinear terms both depend on $r_{m}\zeta^{\prime}(r_{m})$ in Eq.(18), $\mu_{3}$ is expressed by $\sigma$. Using the variance $\sigma^{2}(r_{m})$ and the skewness $\mu_{3}(r_{m})$, we approximately obtain the probability density function of $\bar{\delta}_{\text{HC}}(r_{m})$. When the skewness is regarded as the small perturbation, $\bar{\delta}_{\text{HC}}$ is wriiten with use of the Gaussian variable $\delta_{g}$ as

[TABLE]

where $\delta_{g}$ follows the Gaussian probability density function $P_{g}(\delta_{g})$,

[TABLE]

Therefore, the probability density function of $\bar{\delta}_{\text{HC}}$ is obtained as Young and Byrnes (2013); Byrnes et al. (2012); Young et al. (2016)

[TABLE]

where $\delta_{g;\pm}(\bar{\delta}_{\text{HC}})$ are two solutions of $\bar{\delta}_{\text{HC}}=\bar{\delta}_{\text{HC}}[\delta_{g}]$ in Eq.(31) given by

[TABLE]

We then rewrite the PBH formation rate Eq.(13) as

[TABLE]

with a error function $\text{Erf}(x)=\frac{2}{\sqrt{\pi}}\int^{z}_{0}\text{d}t~{}e^{-t^{2}}$.

To calculate the mass spectrum, we need to relate the PBH mass and the horizon scale $r_{m}$. The PBH mass is more directly related to areal radius $R(r_{m})$, not $r_{m}$ because the horizon scale is given by $H_{b}^{-1}=R(r_{m})$. The difference between two scales may change the PBH mass by $e^{2\zeta(r_{m})}=\mathcal{O}(1)$. However, other uncertainties also change masses of PBHs. For example, the critical phenomena of PBH formation Niemeyer and Jedamzik (1999) and the choice of the initial profiles leads to uncertainties of $\mathcal{O}(1)$. In this paper, thus, we obtain PBH mass using Eq.(II) with $k_{m}=r_{m}^{-1}$ as an approximation.

IV Numerical results

To estimate the nonlinear effect, we calculate PBH mass spectrum for the following simple curvature power spectrum:

[TABLE]

Inserting $P_{\zeta}(k)$ into Eq.(28) and Eq.(30), the variance and the skewness are written as

[TABLE]

The negative $\mu_{3}$ suppresses the PBH formation rate, as expected. We show the scale dependence of the variance and the skewness with $A=1$ in Fig.2. The black solid line represents the variance of $\bar{\delta}_{\text{HC}}$, which has the maximum value at $k_{*}r_{m}\sim 2.5$. The dotted red line shows $-\mu_{3}(r_{m})$.

Next, we calculate the mass spectrum of PBHs. We consider $30{M_{\odot}}$ PBHs with $\Omega_{\text{PBH}}/\Omega_{\text{DM}}\simeq 10^{-3}$ which account for the gravitational wave events discovered by the LIGO-Virgo collaboration Sasaki et al. (2016) as a demonstration. First, we consider the linear order calculation in Eq.(11) with the top-hat window function and the transfer function, which is equivalent to Eq.(35) with $\mu_{3}=0$. Inserting $\sigma$(37) into the Gaussian probability distribution function (11), we calculate the PBH mass spectrum (15). We use the power spectrum (36) with $A_{c}=0.00605$ and $k_{*}\sim 8.6\times 10^{5}\mathrm{Mpc}^{-1}$ to achieve $30{M_{\odot}}$ PBHs with $\Omega_{\text{PBH}}/\Omega_{\text{DM}}\simeq 10^{-3}$. The red line in Fig.3 shows the results of the linear order calculation.

Second, we include the skewness $\mu_{3}$ (IV) into the PBH formation rate (35). The solid black line in Fig.3 represent the mass spectrum on the nonlinear calculation on the power spectrum with $A=A_{c}$ and $k=k_{*}$. As expected, the negative $\mu_{3}$ suppress the PBH formation rate and abundance. Thus, the larger power spectrum is needed to reproduce the result of the linear calculation as shown by the black dashed line in Fig. 3. It is seen that the nonlinear calculation well reproduce the linear one with the additional factor $\simeq 1.4^{2}$.

At last, we compare the non-Gaussianity from nonlinear relation between the density and curvature perturbations (16) with the case where there would be primordial non-Gaussianity of the curvature perturbation. When the curvature perturbation itself follows a non-Gaussian distribution, $\bar{\delta}_{\text{HC}}$ also has non-Gaussianity. As an example of non-Gaussianity of $\zeta$, we focus on the local form whose effect on the PBH formation was studied in Young et al. (2016). The bispectrum $B_{\zeta}$ of the curvature perturbation is written as

[TABLE]

where $f_{NL}$ is the non-linear parameter. Using linear relation between $\zeta$ and $\bar{\delta}$ [Eqs.(6) and (7)], we calculate the contribution of $B_{\zeta}^{\text{local}}$ to the skewness of $\bar{\delta}$ as

[TABLE]

where $\cos\theta=\bm{p}\cdot\bm{q}/(|\bm{p}||\bm{q}|)$ and $\sigma$ is obtained from (12). Multiplying $\mathcal{P}_{\zeta}$ by the transfer function $T$, $\mu_{3:f_{\text{NL}}}$ for the power spectrum (36) is given by

[TABLE]

We calculate $\mu_{3:f_{\text{NL}}=1}(r_{m})$ with the top-hat window function and show the result by the blue dashed line in Fig.2. It is seen that the effect of the nonlinear $\bar{\delta}-\zeta$ relation is comparable to the primordial non-Gaussianity with $f_{\text{NL}}\sim-1$.

V Conclusion

Using the long-wave approximation one obtains the density perturbation as a function of the curvature perturbation and it nonlinearly depends on the curvature perturbation. We have studied the effect of the nonlinear relation on formation of PBHs. As for the PBH formation criterion, we have adopted the compaction function which is also nonlinearly dependent on the curvature perturbation. It is found that required curvature perturbation for the PBH formation increases due to the nonlinearity.

Since the criterion of PBH formation is obtained for spherically symmetric density profiles in the numerical simulation, we have used the spherical harmonics expansion of the curvature perturbation. We then calculated the variance and the skewness of the density perturbation and derived the non-Gaussian property. We estimated the PBH abundance based on the Press–Schechter formalism with non-Gaussian probability density function. As a result, we found that the nonlinear term slightly suppresses the PBH formation. This suppression is comparable to that expected if the primordial curvature perturbation would have the local form of non-Gaussianity with non-linear parameter $f_{\text{NL}}\sim-1$. Comparing the linear order calculation, this suppression requires $\sim 1.4^{2}$ larger power spectrum.

Acknowledgements

This work was supported by JSPS KAKENHI Grant Nos. 17H01131 (M.K.) and 17K05434 (M.K.), MEXT KAKENHI Grant No. 15H05889 (M.K.), World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan (M.K., H.N.), and Advanced Leading Graduate Course for Photon Science (H.N.).