Clustering of primordial black holes with non-Gaussian initial fluctuations

Abstract

We formulate the two-point correlation function of primordial black holes (PBHs) at their formation time, based on the functional integration approach which has often been used in the context of halo clustering. We find that PBH clustering on super-Hubble scales could never be induced in the case where the initial primordial fluctuations are Gaussian, while it can be enhanced by the so-called local-type trispectrum (four-point correlation function) of the primordial curvature perturbations.

1. Introduction

Thanks to the recent detections of gravitational waves from binary black holes by the LIGO/VIRGO collaboration, primordial black holes (PBHs) are attracting attention as a candidate for such binary black hole systems. In Ref. [1] (see also Ref. [2]) we estimated the merger rate of PBH binary systems and found that if the primordial black holes account for |$0.1 \% \sim 1.0\%$| of the dark matter in the Universe, the PBH scenario can be consistent with the first LIGO gravitational wave (GW) event, GW150914. In the analysis in Ref. [1], we assumed that the distribution of PBHs is spatially uniform. As discussed in Ref. [3] (see also Refs. [4–6]),¹ if PBHs spatially clustered at the formation, which can be characterized by the two-point correlation function, it would affect the probability of PBH binary formation and the estimation of the merger rate of PBH binaries.

There are several works about the initial clustering of PBHs.² As a pioneer work, Ref. [11] (and also Ref. [12]) discussed galaxy formation due to the spatial fluctuations in PBH number density. More detailed discussion was given in Ref. [13]. The PBH two-point correlation function and the power spectrum were estimated in Ref. [13] by making use of the peak formalism developed in the context of the clustering of galaxies/halos. As a result, the power spectrum of PBH distribution on large scales should be dominated by the Poisson noise and it behaves as matter isocurvature perturbations. Based on such a PBH isocurvature due to the Poisson noise, Ref. [14] put a constraint on the abundance of PBHs from Ly-|$\alpha$| forest observations, and Refs. [15,16] estimated an expected constraint on it by making use of future 21|$\,$|cm observations. Recently, Ref. [17] studied the initial clustering of PBHs in more detail and found that even on much smaller scales PBHs should not be clustered beyond Poisson. Furthermore, Ref. [18] investigated the dependence of the clustering feature on the shape of the primordial curvature power spectrum. While Refs. [13,17,18] assumed Gaussian initial fluctuations, Refs. [19–20] investigated the impact of primordial local-type non-Gaussianity on the super-Hubble density fluctuations of PBHs, based on the peak-background split picture. Although they found that local-type non-Gaussianity could induce the super-Hubble correlations of the PBH density fluctuations, they focused on a specific type of non-Gaussianity and did not obtain a formula for a PBH two-point correlation function which is applicable for more general types of non-Gaussianity. Recently, Ref. [21] also investigated the effect of primordial non-Gaussianities, not only on the abundance but also on the clustering property of PBHs. In Ref. [21], the threshold of PBH formation is supposed to be given in terms of primordial curvature perturbations. Their result indicates that PBHs are clustered beyond Poisson on super-Hubble scales even for Gaussian primordial curvature perturbations, which is apparently inconsistent with the results obtained in previous works.

In this paper we estimate the two-point correlation function of PBHs by making use of a functional integration approach. This method is powerful in studying correlation functions of biased objects since it allows us to systematically include the effect of non-Gaussian properties of the underling density fluctuations [22–24]. Actually, this approach has also been used in Ref. [21]. In the radiation-dominated era, PBHs are actually formed soon after horizon reentry if the amplitude of primordial fluctuations is greater than a certain threshold. The primordial fluctuation often used to study a criterion for PBH formation is the density contrast in a comoving slice. This quantity represents a local three-curvature and is in good accordance with the physical argument that PBH formation should be determined by local dynamics (see, e.g., Ref. [25]). Thus, contrary to Ref. [21], we employ the density contrast in a comoving slice as a critical quantity for PBH formation.

This paper is organized as follows. In the next section, based on the functional integration approach (or path integral method; see, e.g., Ref. [22]), we formulate a two-point correlation function for PBHs which is applicable to non-Gaussian primordial perturbations. In Sect. 3 we investigate the possibility that PBHs are clustered on large scales due to primordial non-Gaussianities. In order to discuss how PBHs are clustered, we simply assume scale-independent local-type non-Gaussianity of the primordial curvature perturbations and show the relation between the PBH two-point correlation on large scales and the non-linearity parameter. Section 4 is devoted to the conclusion.

2. Formulation

Since during the radiation-dominated era PBHs are formed in the overdense (or positive spatial curvature) region with the Hubble horizon size, the criterion for PBH formation is supposed to be locally determined independently of super-Hubble scale fluctuations. Thus we introduce local smoothed primordial fluctuations |$\theta_{\rm local}({\boldsymbol x})$| as

$$\begin{equation} \theta_{\rm local} ({\boldsymbol x}) := \int d^3 y\, W_{\rm local} ({\boldsymbol x} - {\boldsymbol y}) \, \theta ({\boldsymbol y}), \label{eq:local} \end{equation}$$

(1)

where the window function |$W_{\rm local}({\boldsymbol x})$| is a smoothing function with scale |$R$| that also removes wavelength modes longer than the scale |$R$|⁠. For PBH formation, the scale |$R$| is roughly matched with the Hubble horizon size at formation. In this section, to maintain generality we do not specify a particular gauge for defining the primordial fluctuations |$\theta ({\boldsymbol x})$|⁠, and assume the criterion for PBH formation is given by³

$$\begin{equation} \theta_{\rm local} \geq \theta_{\rm th}. \end{equation}$$

(2)

In the next section we will choose the density contrast on the comoving slice as |$\theta_{\rm local}({\boldsymbol x})$|⁠.

2.1. Probability of PBH formation

The probability that a point |${\boldsymbol x}$| becomes a PBH can be given by

$$\begin{equation} P_1({\boldsymbol x}) = \int [D\theta] P[\theta] \int_{\theta_{\rm th}}^{\infty} d\alpha\, \delta_D (\theta_{\rm local}({\boldsymbol x}) - \alpha), \end{equation}$$

(3)

where |$P[\theta]$| is a probability distribution function for the primordial fluctuations, |$\theta ({\boldsymbol x})$|⁠; also, the probability that two points |${\boldsymbol x}_1$| and |${\boldsymbol x}_2$| are PBHs can be given by

$$\begin{equation} P_2 ({\boldsymbol x}_1,\,{\boldsymbol x}_2) = \int [D\theta] P[\theta] \int_{\theta_{\rm th}}^{\infty} d\alpha_1\, \delta_D (\theta_{\rm local}({\boldsymbol x}_1) - \alpha_1) \int_{\theta_{\rm th}}^{\infty} d\alpha_2\, \delta_D (\theta_{\rm local}({\boldsymbol x}_2) - \alpha_2). \end{equation}$$

(4)

By using the expression for the one-dimensional Dirac delta function given as

$$\begin{equation} \delta_D (x) = \int \frac{d\phi}{2 \pi} e^{i \phi x}, \end{equation}$$

(5)

and the expression for the local smoothed fluctuations given by Eq. (1), we have

$$\begin{equation} P_1({\boldsymbol x}) = \int [D\theta] P[\theta] \int_{\theta_{\rm th}}^{\infty} d\alpha\,\int^\infty_{-\infty} \frac{d \phi}{2\pi} \exp \left[ i \phi \int d^3 y \,W_{\rm local} ({\boldsymbol x} - {\boldsymbol y}) \,\theta ({\boldsymbol y}) - i \phi \, \alpha \right] \end{equation}$$

(6)

and

$$\begin{eqnarray} P_2 ({\boldsymbol x}_1,\,{\boldsymbol x}_2) &=& \int [D\theta] P[\theta] \int_{\theta_{\rm th}}^{\infty} d\alpha_1\,\int_{\theta_{\rm th}}^{\infty} d\alpha_2\, \int^\infty_{-\infty} \frac{d \phi_1}{2\pi} \int^\infty_{-\infty} \frac{d \phi_2}{2\pi} \exp \left[ - i \phi_1\,\alpha_1 - i \phi_2 \, \alpha_2 \right] \cr\cr &&\times \exp \left\{ i \int d^3 y \left[ \phi_1 W_{\rm local} ({\boldsymbol x}_1 - {\boldsymbol y}) + \phi_2 W_{\rm local} ({\boldsymbol x}_2 - {\boldsymbol y}) \right] \theta ({\boldsymbol y})\right\}\!. \end{eqnarray}$$

(7)

2.2. |$P_1$| and |$P_2$| in terms of the |$n$|-point correlators of the primordial fluctuations

Let us introduce a generating function given by

$$\begin{equation} Z[J] := \int [D \theta] P[\theta] \exp \left[ i \int d^3 y J({\boldsymbol y}) \theta ({\boldsymbol y}) \right] = \left \langle \exp \left[ i \int d^3 y J({\boldsymbol y}) \theta ({\boldsymbol y}) \right] \right \rangle, \end{equation}$$

(8)

with |$Z[0] = 1$|⁠. The “connected” |$n$|-point correlation functions of |$\theta ({\boldsymbol x})$| can be given in terms of the generating function as

$$\begin{equation} \xi_{\theta (c)} ({\boldsymbol x}_1, {\boldsymbol x}_2, \ldots , {\boldsymbol x}_n) := \frac{1}{i^n} \frac{ \delta^n \log Z[J]}{\delta J({\boldsymbol x}_1) \delta J({\boldsymbol x}_2) \cdots \delta J ({\boldsymbol x}_n)} \biggr|_{J = 0} . \end{equation}$$

(9)

Inversely, we can obtain the expression for |$\log Z[J]$| in terms of the “connected” |$n$|-point correlation functions as

$$\begin{equation} \log Z[J] = \sum^\infty_{n=1} \frac{i^n}{n!} \int d^3 y_1 d^3 y_2 \cdots d^3 y_n \, \xi_{\theta (c)} ({\boldsymbol y}_1,{\boldsymbol y}_2, \ldots, {\boldsymbol y}_n) \, J({\boldsymbol y}_1) J({\boldsymbol y}_2) \cdots J({\boldsymbol y}_n) . \end{equation}$$

(10)

Choosing |$J({\boldsymbol y}) : = \phi \, W_{\rm local} ({\boldsymbol x} - {\boldsymbol y})$|⁠, the one-point probability of PBH formation, |$P_1$|⁠, can be written as

$$\begin{eqnarray} P_1 ({\boldsymbol x}) &=& \int_{\theta_{\rm th}}^{\infty} d\alpha\,\int^\infty_{-\infty} \frac{d \phi}{2\pi} e^{-i\phi\,\alpha} \,Z\left[\phi\, W_{\rm local} ({\boldsymbol x} - {\boldsymbol y})\right] \cr\cr &=& \int_{\theta_{\rm th}}^{\infty} d\alpha\,\int^\infty_{-\infty} \frac{d \phi}{2\pi} \exp \left[ - i \phi \, \alpha \right] \exp\left[ \sum^\infty_{n=2} \frac{i^n}{n!} \phi^n \, \xi^{(n)}_{{\rm local} (c)} \right] , \end{eqnarray}$$

(11)

where |$\xi^{(n)}_{{\rm local} (c)}$| correspond to the moments of |$\theta_{\rm local}$|⁠:

$$\begin{equation} \xi_{{\rm local} (c)}^{(n)} := \int d^3 y_1 d^3 y_2 \cdots d^3 y_n \, \xi_{\theta (c)} ({\boldsymbol y}_1,{\boldsymbol y}_2, \ldots, {\boldsymbol y}_n) \prod^n_{r=1} W_{\rm local} ({\boldsymbol x} - {\boldsymbol y}_r) . \end{equation}$$

(12)

In the same way, the two-point probability of PBH formation, |$P_2$|⁠, can be expressed as

$$\begin{eqnarray} P_2 ({\boldsymbol x}_1,\,{\boldsymbol x}_2) &=& \int_{\theta_{\rm th}}^{\infty} d\alpha_1\,\int_{\theta_{\rm th}}^{\infty} d\alpha_2\, \int^\infty_{-\infty} \frac{d \phi_1}{2\pi} \int^\infty_{-\infty} \frac{d \phi_2}{2\pi} \exp \left[ - i \phi_1\,\alpha_1 - i \phi_2 \, \alpha_2 \right] \cr\cr && \qquad \times \exp \left[ \sum^\infty_{n=2} i^n \sum_{m = 0}^n \frac{\phi_1^m \phi_2^{n-m}}{m! (n-m)!} \xi_{{\rm local}(c),m}^{(n)}\right] , \end{eqnarray}$$

(13)

where

$$\begin{eqnarray} \xi_{{\rm local}(c),m}^{(n)} & := &\xi_{{\rm local}(c)}^{(n)}( \underbrace{{\boldsymbol x}_1, {\boldsymbol x}_1, \ldots, {\boldsymbol x}_1}_{{\rm total}~m},\underbrace{{\boldsymbol x}_2, {\boldsymbol x}_2, \ldots {\boldsymbol x}_2}_{{\rm total}~n-m}) \cr\cr & = & \int d^3 y_1 d^3 y_2 \cdots d^3 y_n \, \xi_{\theta (c)} ({\boldsymbol y}_1,{\boldsymbol y}_2, \ldots, {\boldsymbol y}_n) \cr\cr && \qquad\qquad \times \prod^m_{r_1 = 1} W_{\rm local} ({\boldsymbol x}_1 - {\boldsymbol y_{r_1}})\prod^n_{r_2 = m+1} W_{\rm local}({\boldsymbol x}_2 - {\boldsymbol y}_{r_2}) . \end{eqnarray}$$

(14)

This can be regarded as a cross correlation between the |$m$|th and |$(n-m)$|th moments of the local smoothed primordial fluctuations.

In the above expression for |$P_1$| we can perform the integration with respect to |$\phi$| as

$$\begin{align} \int^\infty_{-\infty} \frac{d \phi}{2\pi} e^{- i \phi \, \alpha } \exp\left[ \sum^\infty_{n=2} \frac{i^n}{n!} \phi^n \, \xi^{(n)}_{{\rm local} (c)} \right] & = \int^\infty_{-\infty} \frac{d \phi}{2\pi} \exp\left[ \sum^\infty_{n=3} \frac{i^n}{n!} \phi^n \, \xi^{(n)}_{{\rm local} (c)} \right] e^{ - i \phi \, \alpha - \frac{1}{2} \phi^2 \sigma_{\rm local}^2} \cr & = \int^\infty_{-\infty} \frac{d \phi}{2\pi} \exp\left[ \sum^\infty_{n=3} \frac{i^n}{n!} \frac{ \xi^{(n)}_{{\rm local} (c)}\partial^n}{\partial (-i\alpha)^n} \right] e^{ - i \phi \, \alpha - \frac{1}{2} \phi^2 \sigma_{\rm local}^2} \nonumber\\ & = \exp\left[ \sum^\infty_{n=3} \frac{(-1)^n}{n!} \xi^{(n)}_{{\rm local} (c)} \frac{d^n}{d \alpha^n} \right] \frac{e^{-\frac{\alpha^2}{2 \sigma_{\rm local}^2}}}{ \sqrt{2 \pi \sigma_{\rm local}^2}} , \end{align}$$

(15)

and then the expression for |$P_1$| can be reduced to

$$\begin{eqnarray} \label{eq:exP1_1} P_1 = \frac{1}{\sqrt{2 \pi}}\int^\infty_{\nu} dw \exp\left[ \sum^\infty_{n=3} \frac{(-1)^n}{n!}\frac{\xi^{(n)}_{{\rm local} (c)}}{\sigma_{\rm local}^n} \frac{d^n}{d w^n} \right] e^{-\frac{w^2}{2}} , \end{eqnarray}$$

(16)

where |$w := \alpha / \sigma_{\rm local}$|⁠, |$\nu := \theta_{\rm th} / \sigma_{\rm local}$|⁠, and |$\sigma_{\rm local}:= \big(\xi^{(2)}_{{\rm local}(c)}\big)^{1/2}$|⁠.

We can also obtain a reduced form for |$P_2$|⁠:

$$\begin{align} \label{eq:exP2_1} P_2 ({\boldsymbol x}_1,\,{\boldsymbol x}_2) & = \frac{1}{2 \pi} \int_{\nu}^{\infty} dw_1\,\int_{\nu}^{\infty} dw_2 \exp\left[ \frac{\xi_{{\rm local}(c)}^{(2)}({\boldsymbol x}_1, {\boldsymbol x}_2)}{ \sigma_{\rm local}^2} \frac{\partial^2}{\partial w_1 \partial w_2} \right. \cr & \qquad\qquad\qquad \left. + \sum^\infty_{n = 3} (-1)^n\sum_{m=0}^n \frac{\xi_{{\rm local}(c),m}^{(n)}/\sigma_{\rm local}^n}{m! (n-m)!} \frac{\partial^n}{\partial w_1^m \partial w_2^{n-m}}\right] e^{- \frac{w_1^2 + w_2^2}{2}} . \end{align}$$

(17)

2.3. PBH two-point correlation function

Let us obtain an approximate formula for the PBH two-point correlation function by performing the integration with respect to |$w$| in the above expressions. To do so we employ two approximations: (a) a weak non-Gaussian limit, and (b) a high peak limit.

In the weak non-Gaussian limit, we assume that |$\xi^{(n)}_{{\rm local}(c)}/\sigma_{\rm local}^n \ll 1$| for |$n \geq 3$|⁠. Under this assumption, up to the linear order in |$\xi^{(n)}_{{\rm local}(c)}/\sigma_{\rm local}^n$| we have

$$\begin{eqnarray} P_1 &=& \frac{1}{\sqrt{2 \pi}} \int^\infty_\nu dw \left[1 + \sum_{n=3}^\infty \frac{(-1)^n}{n!}\frac{\xi^{(n)}_{{\rm local} (c)}}{\sigma_{\rm local}^n} \frac{d^n}{d w^n} \right] e^{-\frac{w^2}{2}} + O\left( (\xi^{(n)}_{{\rm local} (c)}/\sigma_{\rm local}^n)^2 \right)\!. \end{eqnarray}$$

(18)

With the assumptions |$\xi^{(n)}_{{\rm local}(c),m}/\sigma_{\rm local}^n \ll 1$| for |$n \geq 3$| and |$\xi_{{\rm local}(c)}^{(2)}({\boldsymbol x}_1, {\boldsymbol x}_2) / \sigma_{\rm local}^2 \ll 1$|⁠, we can also obtain

$$\begin{eqnarray} P_2 &=& \frac{1}{2 \pi} \int^\infty_\nu dw_1 \int^\infty_\nu dw_2 \left[1 + \frac{\xi_{{\rm local}(c)}^{(2)}({\boldsymbol x}_1, {\boldsymbol x}_2)}{ \sigma_{\rm local}^2} \frac{\partial^2}{\partial w_1 \partial w_2} \right.\nonumber\\[4pt] &&\qquad\qquad \qquad \qquad \qquad \left. + \sum^\infty_{n = 3} (-1)^n\sum_{m=0}^n \frac{\xi_{{\rm local}(c),m}^{(n)}/\sigma_{\rm local}^n}{m! (n-m)!} \frac{\partial^n}{\partial w_1^m \partial w_2^{n-m}}\right] e^{- \frac{w_1^2 + w_2^2}{2}}\nonumber\\[4pt] && \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad +\, O\left( (\xi^{(n)}_{{\rm local} (c)}/\sigma_{\rm local}^n)^2 \right)\!. \end{eqnarray}$$

(19)

Furthermore, by making use of the Hermite polynomials, |$H_n(x)$|⁠, given as |$H_n(x) := (-1)^n e^{x^2} (d/dx)^n e^{-x^2}$|⁠, we can replace the derivatives in terms of |$w$| with the Hermite polynomials. Then, finally, by employing the high peak approximation, |$\nu \gg 1$|⁠, which might naturally be valid for PBH formation, we can perform the integration approximately in terms of |$w$| and obtain approximate expressions for |$P_1$| and |$P_2$| as

$$\begin{equation} \label{eq:exP1_a} P_1 \approx \frac{e^{-\nu^2/2}}{\sqrt{2 \pi} \nu} \left[1 + \sum_{n=3}^\infty \frac{1}{2^{n/2}n!}\frac{\xi^{(n)}_{{\rm local} (c)}}{\sigma_{\rm local}^n} H_n\left(\frac{\nu}{\sqrt{2}} \right) \right] \end{equation}$$

(20)

and

$$\begin{eqnarray}\label{eq:exP2_a} P_2 ({\boldsymbol x}_1, {\boldsymbol x}_2) &\approx& \frac{e^{-\nu^2}}{2 \pi \nu^2} \left[ 1 + \frac{\xi_{{\rm local}(c)}^{(2)} ({\boldsymbol x}_1, {\boldsymbol x}_2)}{2 \sigma_{\rm local}^2} H_1 \left( \frac{\nu}{\sqrt{2}} \right)^2 \right.\nonumber\\ && \qquad\qquad \left. + \sum_{n=3}^{\infty}\sum_{m=0}^{n} \frac{\xi_{{\rm local}(c),m}^{(n)}/\sigma_{\rm local}^n}{2^{n/2}m! (n-m)!}H_m\left(\frac{\nu}{\sqrt{2}} \right)H_{n-m}\left(\frac{\nu}{\sqrt{2}} \right) \right]\!. \end{eqnarray}$$

(21)

Then, by using these expressions for |$P_1$| and |$P_2$|⁠, the two-point correlation function of the PBHs can be evaluated as

$$\begin{eqnarray} \xi_{\rm PBH} ({\boldsymbol x}_1, {\boldsymbol x}_2) &:=& \frac{P_2({\boldsymbol x}_1,{\boldsymbol x}_2) }{P_1^2} - 1 \cr\cr &\approx& \frac{H_1(\nu/\sqrt{2})^2}{2 \sigma_{\rm local}^2} \xi_{{\rm local}(c)}^{(2)} ({\boldsymbol x}_1, {\boldsymbol x}_2) + \sum_{n=3}^{\infty}\sum_{m=1}^{n-1} \frac{H_m(\nu/\sqrt{2})H_{n-m}(\nu/\sqrt{2})}{2^{n/2}m! (n-m)!\sigma_{\rm local}^n}\xi_{{\rm local}(c),m}^{(n)} . \nonumber\\ \end{eqnarray}$$

(22)

Note that, strictly speaking, in order for the above expression to be valid we should require stronger assumptions given as

$$\begin{eqnarray} \frac{\xi^{(n)}_{{\rm local}(c),m}}{\sigma_{\rm local}^n} \ll \frac{2^{n/2} m!(n-m)!}{H_m (\nu/\sqrt{2}) H_{n-m}(\nu/\sqrt{2})} \label{eq:assum_1} \end{eqnarray}$$

(23)

and

$$\begin{equation} \frac{\xi^{(2)}_{{\rm local}(c)}({\boldsymbol x}_1,{\boldsymbol x}_2)}{\sigma_{\rm local}^2} \ll \frac{2}{H_1 (\nu/\sqrt{2})^{2}}. \label{eq:assum_2} \end{equation}$$

(24)

Up to the four-point correlation function of the local smoothed primordial fluctuations, by using an approximate form of the Hermite polynomials, |$H_n (x) \sim 2^n x^n$| for |$x \gg 1$|⁠, we have

$$\begin{eqnarray} \xi_{\rm PBH} ({\boldsymbol x}_1, {\boldsymbol x}_2) &\sim& \frac{\nu^2}{\sigma_{\rm local}^2} \xi_{{\rm local}(c)}^{(2)}({\boldsymbol x}_1, {\boldsymbol x}_2) \cr\cr && + \frac{1}{2} \frac{\nu^3}{\sigma_{\rm local}^3} \left( \xi^{(3)}_{{\rm local}(c)} ({\boldsymbol x}_1, {\boldsymbol x}_1, {\boldsymbol x}_2) ({\boldsymbol x}_1 \leftrightarrow {\boldsymbol x}_2) \right) \cr\cr && \frac{1}{4} \frac{\nu^4}{\sigma_{\rm local}^4} \xi^{(4)}_{{\rm local}(c)} ({\boldsymbol x}_1, {\boldsymbol x}_1, {\boldsymbol x}_2,{\boldsymbol x}_2) \cr\cr && + \frac{1}{6} \frac{\nu^4}{\sigma_{\rm local}^4} \left( \xi^{(4)}_{{\rm local}(c)} ({\boldsymbol x}_1, {\boldsymbol x}_2, {\boldsymbol x}_2,{\boldsymbol x}_2) + ({\boldsymbol x}_1 \leftrightarrow {\boldsymbol x}_2)\right)\!. \label{eq:PBHcorr} \end{eqnarray}$$

(25)

This expression corresponds to the result obtained in Refs. [23,24] in the context of the halo bias. Under the assumptions given by Eqs. (23) and (24), this equation is valid for |$\xi_{\rm PBH}({\boldsymbol x}_1, {\boldsymbol x}_2) \ll 1$|⁠.

3. PBH clustering with local-type non-Gaussianity

Before a quantitative discussion for the amplitude of the PBH two-point correlation function, let us first provide an intuitive idea of why we take into account up to the four-point correlation function of the primordial fluctuations in the above formulation.

Hereafter, as the primordial fluctuations, |$\theta ({\boldsymbol x})$|⁠, in the above formulation we use the primordial curvature perturbations on the comoving slice, often denoted by |${\mathcal R}_c({\boldsymbol x})$|⁠. In the long-wavelength approximation, comoving density fluctuations can be given in terms of |${\mathcal R}_c({\boldsymbol x})$| as [27]

$$\begin{eqnarray} \delta ({\boldsymbol x}) = -\frac{4 (1 + w)}{3 w + 5} e^{- 5 {\mathcal R}_c({\boldsymbol x})/2} \frac{\triangle}{a^2 H^2} ( e^{{\mathcal R}_c({\boldsymbol x})/2}), \label{eq:delta} \end{eqnarray}$$

(26)

where |$w$|⁠, |$a$|⁠, and |$H$| are respectively an equation of state of the Universe, a scale factor, and the Hubble parameter. As can be seen in the above expression, if we use the primordial curvature perturbations on the comoving slice as |$\theta ({\boldsymbol x})$| in the previous section, a natural variable for the local primordial fluctuations |$\theta_{\rm local}$| would be the comoving density fluctuations. In fact, this quantity represents a local three-curvature and is in good accordance with a physical argument that the criterion for PBH formation should be determined by local dynamics (i.e. within the Hubble horizon) and be free from the addition of super-Hubble modes. By introducing a local scale factor |$a e^{{\mathcal R}_c ({\boldsymbol x})} \to a$| [25] at the linear order, the above expression can be reduced to

$$\begin{eqnarray} \delta ({\boldsymbol x}) = -\frac{4}{9} \frac{1}{a^2 H^2} \triangle {\mathcal R}_c({\boldsymbol x}) . \label{eq:pertDelta} \end{eqnarray}$$

(27)

Here, we take |$w = 1/3$| in the radiation-dominated era. The two-point correlation function of |$\delta$| is given by

$$\begin{equation} \langle \delta ({\boldsymbol x}) \delta ({\boldsymbol y}) \rangle={\left( \frac{4}{9a^2 H^2} \right)}^2 \triangle^2 \langle {\cal R}_c ({\boldsymbol x}) {\cal R}_c ({\boldsymbol y}) \rangle. \label{dxdy} \end{equation}$$

(28)

Because of the two Laplacians, |$\langle \delta ({\boldsymbol x}) \delta ({\boldsymbol y}) \rangle$| rapidly approaches zero for |$| {\boldsymbol x}-{\boldsymbol y}| \gg {(a H)}^{-1}$| unless |${\cal R}_c$| is extremely red-tilted. Thus, in general, |$\langle \delta ({\boldsymbol x}) \delta ({\boldsymbol y}) \rangle$| is suppressed on super-Hubble scales.

Due to the locality, at the leading order in |$\delta$| the PBH abundance at point |${\boldsymbol x}$| would be determined by the local variance |$\langle \delta^2 ({\boldsymbol x}) \rangle$|⁠, and then the PBH two-point correlation function is given by the its correlation. If |$\delta$| is Gaussian, the correlation of the local variance is given as

$$\begin{equation} \frac{\langle \delta^2 ({\boldsymbol x}) \delta^2 ({\boldsymbol y}) \rangle}{{\langle \delta^2 ({\boldsymbol x}) \rangle}^2}-1 =2 {\langle \delta ({\boldsymbol x}) \delta ({\boldsymbol y}) \rangle}^2 , \end{equation}$$

(29)

and it should be suppressed on super-Hubble scales, as can be seen by Eq. (28). Thus, PBHs are produced by the same amount in every super-Hubble size region. In other words, PBHs are not clustered on super-Hubble scales. If, on the other hand, |$\delta$| is non-Gaussian, the correlation of the local variance may remain on super-Hubble scales. In order to see this explicitly, let us focus on the simple case where |${\cal R}_c$| consists of two uncorrelated Gaussian fields |$\phi, \chi$| as

$$\begin{equation} {\cal R}_c ({\boldsymbol x})=(1+\alpha \chi ({\boldsymbol x})) \phi ({\boldsymbol x}). \label{p-tnl} \end{equation}$$

(30)

Here it is assumed that |$\chi$| has super-Hubble scale correlation and |$\phi$| gives a dominant contribution to PBH formation. For such a case, the density contrast on the comoving slice can be given as

$$\begin{equation} \delta ({\boldsymbol x}) = -(1 + \alpha \chi({\boldsymbol x}) ) \frac{4}{9} \frac{1}{a^2 H^2} \triangle \phi ({\boldsymbol x}). \end{equation}$$

(31)

For super-Hubble distance |$| {\boldsymbol x}-{\boldsymbol y}| \gg {(a H)}^{-1}$|⁠, we obtain

$$\begin{equation} \frac{\langle \delta^2 ({\boldsymbol x}) \delta^2 ({\boldsymbol y}) \rangle}{{\langle \delta^2 ({\boldsymbol x}) \rangle}^2}-1 \approx 4\alpha^2 \langle \chi ({\boldsymbol x}) \chi ({\boldsymbol y}) \rangle+{\cal O}(\alpha^4). \label{scl} \end{equation}$$

(32)

There are two remarks regarding this result. First, on super-Hubble scales the correlation of the local variance is directly proportional to the correlation function of |$\chi$|⁠. This result simply reflects our naive intuition that a local quantity can possess correlation over super-Hubble distance only when the quantity is sourced by another quantity having correlation over super-Hubble distance. Secondly, the correlation of the local variance is proportional to part of the connected part of the four-point function of |${\cal R}_c$|⁠. More explicitly, the connected part of the four-point function of |${\cal R}_c$| is given by

$$\begin{equation} \frac{\langle {\cal R}_c^2 ({\boldsymbol x}) {\cal R}_c^2 ({\boldsymbol y}) \rangle_c}{{\langle {\cal R}_c^2 ({\boldsymbol x}) \rangle}^2} =4\alpha^2 \left[ \langle \chi ({\boldsymbol x}) \chi ({\boldsymbol y}) \rangle+ \frac{{\langle {\cal R}_c ({\boldsymbol x}) {\cal R}_c ({\boldsymbol y}) \rangle}^2}{{\langle {\cal R}_c^2 ({\boldsymbol x}) \rangle}^2} \left( \langle \chi^2 ({\boldsymbol x}) \rangle + \langle \chi ({\boldsymbol x}) \chi ({\boldsymbol y}) \rangle \right) \right] +{\cal O}(\alpha^4),\quad \end{equation}$$

(33)

and the right-hand side of Eq. (32) is obtained by ignoring terms proportional to the two-point correlation function of |${\mathcal R}_c$|⁠, |$\langle {\cal R}_c ({\boldsymbol x}) {\cal R}_c ({\boldsymbol y}) \rangle$|⁠, in the above equation (i.e. only the first term).

In particular, if the power spectrum of |$\chi$| is the same as that of |$\phi$|⁠, the parameter |$\alpha$| is related to the local-type trispectrum parameter |$\tau_{\rm NL}$| as |$\tau_{\rm NL}=\alpha^2$| [28], and Eq. (32) becomes

$$\begin{equation} \frac{\langle \delta^2 ({\boldsymbol x}) \delta^2 ({\boldsymbol y}) \rangle}{{\langle \delta^2 ({\boldsymbol x}) \rangle}^2}-1 \approx 4\tau_{\rm NL} \langle {\cal R}_c ({\boldsymbol x}) {\cal R}_c ({\boldsymbol y}) \rangle+{\cal O}(\alpha^4). \end{equation}$$

(34)

Thus, the correlation of the local variance is proportional to |$\tau_{\rm NL}$| and the correlation function of the curvature perturbation. Notice that the bispectrum parameter is |$f_{\rm NL}=0$| in the present case, and it is actually the trispectrum (not the bispectrum) that determines the clustering of PBHs over the super-Hubble distance.

On the other hand, if |$f_{\rm NL}\neq 0$|⁠, |$\tau_{\rm NL}$| is also non-zero with the lower bound |$\tau_{\rm NL} \geq \frac{36}{25} f_{\rm NL}^2$| [29]. Thus, PBHs are necessarily clustered on super-Hubble scales in this case, which is consistent with Refs. [19,20] which showed that the clustering is characterized by |$f_{\rm NL}$|⁠.

In the following discussion we actually employ the local-type ansatz for the non-Gaussianities of primordial curvature perturbations. Based on our result in Eq. (25), we explicitly show that the PBH two-point correlation function is proportional to |$\tau_{\rm NL}$| and the two-point correlation function of the primordial curvature perturbations.

3.1. Primordial local-type non-Gaussianity

The local-type primordial non-Gaussianity in Fourier space⁴ has conventionally been characterized by introducing three constant parameters, so-called non-linearity parameters, |$f_{\rm NL}$|⁠, |$g_{\rm NL}$|⁠, and |$\tau_{\rm NL}$|⁠, which respectively represent the amplitudes of the bispectrum and trispectrum of the primordial curvature perturbations as [28]

$$\begin{align} \langle {\mathcal R}_c ({\boldsymbol k}_1) {\mathcal R}_c ({\boldsymbol k}_2) {\mathcal R}_c ({\boldsymbol k}_3)\rangle & := (2 \pi)^3 \delta^{(3)}({\boldsymbol k}_1 + {\boldsymbol k}_2 + {\boldsymbol k}_3)\, \frac{6}{5} f_{\rm NL} [ P_{{\mathcal R}_c}(k_1) P_{{\mathcal R}_c}(k_2) + 2~{\rm perms} ] , \cr \langle {\mathcal R}_c ({\boldsymbol k}_1) {\mathcal R}_c ({\boldsymbol k}_2) {\mathcal R}_c ({\boldsymbol k}_3){\mathcal R}_c ({\boldsymbol k}_4)\rangle & := (2 \pi)^3 \delta^{(3)}({\boldsymbol k}_1 + {\boldsymbol k}_2 + {\boldsymbol k}_3 + {\boldsymbol k}_4)\nonumber\\ & \quad \times \left\{ \frac{54}{25}g_{\rm NL} \left[ P_{{\mathcal R}_c}(k_1) P_{{\mathcal R}_c}(k_2) P_{{\mathcal R}_c}(k_3) 3~{\rm perms} \right] \right. \cr & \qquad \quad \left. + \tau_{\rm NL} [ P_{{\mathcal R}_c}(k_1) P_{{\mathcal R}_c}(k_2) P_{{\mathcal R}_c}(|{\boldsymbol k}_1 + {\boldsymbol k}_3|) 11~{\rm perms} ] \vphantom{\frac{54}{25}}\right\}\!,\nonumber\\ \label{eq:localnG} \end{align}$$

(35)

where |$P_{{\mathcal R}_c}(k)$| is the power spectrum of the primordial curvature perturbations given as

$$\begin{equation} \langle{\mathcal R}_c ({\boldsymbol k}) {\mathcal R}_c ({\boldsymbol k}') \rangle := (2 \pi)^3 \delta^{(3)}({\boldsymbol k} + {\boldsymbol k}') P_{{\mathcal R}_c}(k) . \end{equation}$$

(36)

3.2. PBH power spectrum with non-Gaussian primordial fluctuations

In the previous section we discussed the two-point correlation function of the spatial distribution of PBHs. The two-point correlation function can be expressed by the PBH power spectrum as

$$\begin{eqnarray} \xi_{\rm PBH}({\boldsymbol x}_1, {\boldsymbol x}_2) := \langle \delta_{\rm PBH}({\boldsymbol x}_1) \delta_{\rm PBH} ({\boldsymbol x}_2) \rangle = \int \frac{d^3k}{(2 \pi)^3} P_{\rm PBH}(k)\, e^{i {\boldsymbol k}\cdot ({\boldsymbol x}_1 - {\boldsymbol x}_2)} , \end{eqnarray}$$

(37)

where |$\delta_{\rm PBH}({\boldsymbol x})$| is the number density field of PBHs, and |$P_{\rm PBH} (k)$| is the PBH power spectrum:

$$\begin{equation} \langle\delta_{\rm PBH} ({\boldsymbol k}) \delta_{\rm PBH} ({\boldsymbol k}') \rangle := (2 \pi)^3 \delta^{(3)}({\boldsymbol k} + {\boldsymbol k}') P_{{\rm PBH}}(k) , \end{equation}$$

(38)

with |$\delta_{\rm PBH}({\boldsymbol k})$| being a Fourier transform of the PBH number density field. Assuming statistical isotropy, the PBH power spectrum can be inversely given by

$$\begin{equation} P_{\rm PBH}(k) = \int d^3 r \, \xi_{\rm PBH} ({\boldsymbol r}) \,e^{-i{\boldsymbol k}\cdot {\boldsymbol r}} , \end{equation}$$

(39)

where |${\boldsymbol r} := {\boldsymbol x}_1 - {\boldsymbol x}_2$|⁠. Substituting Eq. (25) into this equation and noting that we apply the comoving density fluctuations |$\delta ({\boldsymbol x})$| to |$\theta_{\rm local} ({\boldsymbol x})$| used in the previous section, we obtain

$$\begin{eqnarray} P_{\rm PBH}(k) &\simeq& \frac{\nu^2}{\sigma^2_R}W_R(k)^2 P_{\delta} (k) \cr\cr && + \frac{\nu^3}{\sigma_R^3}\int \frac{d^3 p}{(2 \pi)^3} W_R(p) W_R(|{\boldsymbol k} - {\boldsymbol p}|) W_R(k) B_\delta ({\boldsymbol p},-{\boldsymbol k},{\boldsymbol k} - {\boldsymbol p}) \cr\cr && + \frac{1}{3}\frac{\nu^4}{\sigma_R^4}\int \frac{d^3 p_1 d^3 p_2}{(2 \pi)^6}\cr\cr && \qquad \times W_R(p_1)W_R(p_2) W_R(|{\boldsymbol k} - {\boldsymbol p}_1 - {\boldsymbol p}_2|) W_R(k) T_\delta ({\boldsymbol p}_1,{\boldsymbol p}_2, -{\boldsymbol k}, {\boldsymbol k} - {\boldsymbol p}_1 - {\boldsymbol p}_2) \cr\cr && + \frac{1}{4}\frac{\nu^4}{\sigma_R^4} \int \frac{d^3 p_1 d^3 p_2}{(2 \pi)^6} \cr\cr && \qquad \times W_R(p_1)W_R(p_2) W_R(|{\boldsymbol k} - {\boldsymbol p}_1|)W_R(|{\boldsymbol p}_2 + {\boldsymbol k}|) T_\delta ({\boldsymbol p}_1,{\boldsymbol p}_2, {\boldsymbol k}-{\boldsymbol p}_1, -{\boldsymbol k} - {\boldsymbol p}_2) , \nonumber\\ \label{eq:Ppbh1} \end{eqnarray}$$

(40)

where |$P_\delta$|⁠, |$B_\delta$|⁠, and |$T_\delta$| are respectively power, bi-, and trispectra of the comoving density fluctuations, |$W_R(k)$| is a window function smoothed with the comoving scale |$R = (aH)^{-1}$| at the PBH formation in Fourier space.

$$\begin{eqnarray} P_\delta (k) &=& \left(\frac{4}{9}\right)^2 (kR)^4 P_{{\mathcal R}_c}(k), \cr\cr B_\delta (k_1, k_2, k_3) &=& \left(\frac{4}{9}\right)^3 (k_1R)^2 (k_2R)^2 (k_3R)^2 \frac{6}{5}f_{\rm NL} \left[ P_{{\mathcal R}_c}(k_1) P_{{\mathcal R}_c}(k_2) + 2~{\rm perms} \right], \cr\cr T_\delta ({\boldsymbol k}_1, {\boldsymbol k}_2, {\boldsymbol k}_3, {\boldsymbol k}_4) &=& \left(\frac{4}{9}\right)^4 (k_1R)^2 (k_2R)^2 (k_3R)^2(k_4R)^2 \cr\cr && \times \left\{ \frac{54}{25}g_{\rm NL} \left[ P_{{\mathcal R}_c}(k_1) P_{{\mathcal R}_c}(k_2) P_{{\mathcal R}_c}(k_3) 3~{\rm perms} \right] \right. \cr\cr &&\qquad \left. + \tau_{\rm NL} [ P_{{\mathcal R}_c}(k_1) P_{{\mathcal R}_c}(k_2) P_{{\mathcal R}_c}(|{\boldsymbol k}_1 + {\boldsymbol k}_3|) 11~{\rm perms} ] \vphantom{\frac{54}{25}} \right\} . \label{eq:Deltaspec} \end{eqnarray}$$

(41)

Substituting Eq. (41) into Eq. (40), we have

$$\begin{eqnarray} P_{\rm PBH}(k) &\simeq& \left(\frac{4\nu}{9\sigma_R}\right)^2 W_{\rm local}(k)^2 P_{{\mathcal R}_c} (k) \cr\cr && + \frac{6}{5}f_{\rm NL} \left(\frac{4\nu}{9\sigma_R} \right)^3 W_{\rm local}(k) \cr\cr && \times \int \frac{d^3 p}{(2 \pi)^3} W_{\rm local}(p) W_{\rm local}(|{\boldsymbol k} - {\boldsymbol p}|) \cr\cr && \quad \times \left[ 2 P_{{\cal R}_c}(p) P_{{\cal R}_c}(k) + P_{{\cal R}_c}(p) P_{{\cal R}_c}(|{\boldsymbol k} - {\boldsymbol p}|) \right] \cr\cr && + \frac{18}{25}g_{\rm NL}\left(\frac{4\nu}{9\sigma_R} \right)^4 W_{\rm local}(k) \cr\cr && \times \int \frac{d^3 p_1 d^3 p_2}{(2 \pi)^6}W_{\rm local}(p_1)W_{\rm local}(p_2) W_{\rm local}(|{\boldsymbol k} - {\boldsymbol p}_1 - {\boldsymbol p}_2|) \cr\cr && \quad\quad \times \left[ 3 P_{{\mathcal R}_c}(k)P_{{\mathcal R}_c}(p_1)P_{{\mathcal R}_c}(p_2) + P_{{\mathcal R}_c}(p_1)P_{{\mathcal R}_c}(p_2)P_{{\mathcal R}_c}(|{\boldsymbol k} - {\boldsymbol p}_1 - {\boldsymbol p}_2|) \right] \cr\cr && + \frac{\tau_{\rm NL}}{3}\left(\frac{4\nu}{9\sigma_R} \right)^4 W_{\rm local}(k) \cr\cr && \times \int \frac{d^3 p_1 d^3 p_2}{(2 \pi)^6}W_{\rm local}(p_1)W_{\rm local}(p_2) W_{\rm local}(|{\boldsymbol k} - {\boldsymbol p}_1 - {\boldsymbol p}_2|) \cr\cr && \quad\times \left[ 6 P_{{\mathcal R}_c}(k)P_{{\mathcal R}_c}(p_1)P_{{\mathcal R}_c}(|{\boldsymbol p}_1 + {\boldsymbol p}_2|) + 6 P_{{\mathcal R}_c}(p_1)P_{{\mathcal R}_c}(p_2)P_{{\mathcal R}_c}(|{\boldsymbol k} - {\boldsymbol p}_1 |) \right] \cr\cr && + \frac{54}{25}g_{\rm NL}\left(\frac{4\nu}{9\sigma_R}\right)^4 \cr\cr && \times \int \frac{d^3 p_1 d^3 p_2}{(2 \pi)^6}W_{\rm local}(p_1)W_{\rm local}(p_2) W_{\rm local}(|{\boldsymbol k} + {\boldsymbol p}_1|)W_{\rm local}(|{\boldsymbol p}_2 - {\boldsymbol k}|) \cr\cr && \qquad \qquad \qquad \times P_{{\mathcal R}_c}(p_1) P_{{\mathcal R}_c}(p_2) P_{{\mathcal R}_c}(|{\boldsymbol k}+{\boldsymbol p}_1|) \cr\cr && + \frac{\tau_{\rm NL}}{4}\left(\frac{4\nu}{9\sigma_R}\right)^4 \cr\cr && \times \int \frac{d^3 p_1 d^3 p_2}{(2 \pi)^6}W_{\rm local}(p_1)W_{\rm local}(p_2) W_{\rm local}(|{\boldsymbol k} + {\boldsymbol p}_1|)W_{\rm local}(|{\boldsymbol p}_2 - {\boldsymbol k}|) \cr\cr && \times \left[ 4 P_{{\cal R}_c}(p_1)P_{{\cal R}_c}(p_2) P_{{\cal R}_c}(k) \right.\cr\cr && \quad\left.+\,4 P_{{\cal R}_c}(p_1) \left( P_{{\cal R}_c}(p_2) P_{{\cal R}_c}(|{\boldsymbol k} + {\boldsymbol p}_1 - {\boldsymbol p}_2|) + P_{{\cal R}_c} (|{\boldsymbol k} + {\boldsymbol p}_1|)P_{{\cal R}_c}(|{\boldsymbol p}_1 + {\boldsymbol p}_2|) \right) \right]\!, \nonumber\\ \label{eq:Ppbh2} \end{eqnarray}$$

(42)

where |$W_{\rm local}(k) := (kR)^2 \,W_R(k)$|⁠. As can be seen in the above equation, the contributions from non-zero |$g_{\rm NL}$| and |$\tau_{\rm NL}$| are of the order of |$P_{{\mathcal R}_c}^3$|⁠. Note that this equation is derived based on the approximate expression in Eq. (25). With respect to the order of |$P_{{\mathcal R}_c}$|⁠, |$(\xi^{(2)})^2$|⁠, |$(\xi^{(2)})^3$|⁠, and |$\xi^{(2)} \xi^{(3)}$| terms should also be of the order of |$P_{{\mathcal R}_c}^3$| and hence they should be taken into account in the expression in Eq. (25) (or Eqs. (20) and (21)). However, one can show that these quadratic and cubic terms are included in the |$k$|-independent contribution as in the following expression. Thus, for the super-Hubble correlation of PBHs, these terms can be neglected.

For the |$kR \ll 1$| limit, which we focus on in this paper, noting that we can take |$W_{\rm local}(k) \to 0$| and |$|{\boldsymbol k}+ {\boldsymbol p}| \to p$|⁠, the above expression becomes

$$\begin{eqnarray} P_{\rm PBH}(k) &\simeq & \tau_{\rm NL} \left(\frac{4\nu}{9\sigma_R}\right)^4 \left[\int \frac{d^3p}{(2\pi)^3} W_{\rm local}(p)^2 P_{{\cal R}_c}(p) \right]^2 P_{{\cal R}_c}(k) + (k \mbox{-} {\rm independent}~{\rm terms}) \cr\cr &=& \tau_{\rm NL} \nu^4 P_{{\cal R}_c}(k) + (k \mbox{-} {\rm independent}~{\rm terms}) . \label{eq:powerPBH} \end{eqnarray}$$

(43)

Thus, in the case where the primordial curvature perturbations have local-type non-Gaussianity parameterized by |$\tau_{\rm NL}$|⁠, the PBH power spectrum does not decay even on super-Hubble scales and is proportional to the power spectrum of the primordial curvature perturbations. Inversely, the two-point correlation function of the PBHs is obtained as

$$\begin{eqnarray} \xi^{(2)}_{\rm PBH}({\boldsymbol r}) = \int \frac{d^3 k}{(2 \pi)^3} P_{\rm PBH}(k) e^{i{\boldsymbol k}\cdot{\boldsymbol r}} = \tau_{\rm NL} \nu^4 \xi^{(2)}_{{\mathcal R}_c} ({\boldsymbol r}) +C\,\delta^{(3)}({\boldsymbol r}), \end{eqnarray}$$

(44)

where a constant, |$C$|⁠, corresponds to the |$k$|-independent terms in Eq. (43). A typical value of the enhancement factor is |$\tau_{\rm NL} \nu^4 \sim O(10^6) \times (\tau_{\rm NL}/10^3)$|⁠. For super-Hubble scales only the first term on the right-hand side is relevant, and thus, if |$\tau_{\rm NL}$|-type non-Gaussianity exists, it would give a large effect on the clustering behavior of PBHs even on super-Hubble scales (⁠|$r \gg R$|⁠) at formation.

As discussed in Refs. [19,20], if PBHs contribute to the dark matter component of the Universe, their spatial distribution behaves as dark matter isocurvature perturbations. Recent cosmic microwave background (CMB) observations give a tight constraint on the fraction of dark matter isocurvature perturbations [30]. By introducing a parameter representing the mass fraction of PBHs in the total dark matter, |$f_{\rm PBH}$|⁠, the observational constraint roughly means |$f_{\rm PBH}^{2} \tau_{\rm NL} \nu^4 \lesssim O(10^{-2})$|⁠. If PBHs comprise a dominant component of dark matter, adopting |$\nu=10$| as an approximate value necessary for PBH formation gives an upper limit on |$\tau_{\rm NL}$| of |$\tau_{\rm NL} \lesssim 10^{-6}$|⁠.

Note that this constraint for |$\tau_{\rm NL}$| is obtained by assuming that the simple expression for the primordial trispectrum given by Eq. (35) is valid for CMB scales, which are super-Hubble ones at PBH formation. Although the slow-roll condition would be violated when the PBH formation scale exits the horizon, in general the CMB scales are considered to exit the horizon in the slow-roll phase. As discussed in, e.g., Ref. [31], it might be difficult to realize simple expressions for the primordial non-Gaussianity given by Eq. (35) in single-field slow-roll inflation. Thus, the above constraint for |$\tau_{\rm NL}$| cannot be simply applied in single-field inflation and it does not tightly constrain the PBH formation scenario in such a single-field class.

4. Conclusion

We have investigated the clustering behavior of PBHs at formation during the radiation-dominated era. Since during the radiation-dominated era PBHs would be formed from the direct collapse of the overdense region with Hubble scales, for PBH clustering, super-Hubble fluctuations might be important. We formulated the two-point correlation function of PBHs by making use of a functional integration approach which takes into account the non-Gaussian property of the primordial fluctuations. Our result shows that PBHs are never clustered at formation as long as the primordial fluctuations obey Gaussian statistics and the super-Hubble two-point correlation function is induced by the connected part of the four-point correlation function of the primordial fluctuations. In order to evaluate the super-horizon two-point correlation function of PBHs quantitatively, we considered non-Gaussian primordial perturbations of the local type up to the four-point function (trispectrum) parametrized by two constant parameters, |$g_{\rm NL}$| and |$\tau_{\rm NL}$|⁠. We found that the |$\tau_{\rm NL}$|-type non-Gaussianity determines the super-Hubble two-point correlation function. Thus, to estimate the clustering behavior of PBHs we should carefully investigate the non-Gaussian property of the primordial fluctuations, which strongly depends on the generation mechanism, such as the inflation model.

Acknowledgements

We would like to thank Misao Sasaki, Tomohiro Harada, and Chul-Moon Yoo for useful discussions. SY would also like to thank Jinn-Ouk Gong and Takahiko Matsubara for useful discussions about the functional integration approach formulation. SY is supported by MEXT KAKENHI Grant Numbers 15H05888 and 18H04356. TS is supported by Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Young Scientists (B) No. 15K17632, and by MEXT Grants-in-Aid for Scientific Research on Innovative Areas Nos. 15H05888, 17H06359, 18H04338, and 19K03864.

Footnotes

¹Clustering of primordial black holes has also been discussed in Refs. [7,8]. These papers have investigated the observational effects of PBH clustering, not only on binary formation but also on the formation of supermassive BHs.

²There have been several works about late-time clustering of PBHs, e.g. in dark matter halos (see, e.g., Refs. [9,10]).

³Strictly speaking, the threshold depends on the perturbation profile (see, e.g., Ref. [26]). In this paper we ignore this dependence and assume that the threshold is the same for all profiles. This assumption does not affect the main result of this paper.

⁴We use the Fourier transform expression given by

$$\begin{equation*} f ({\boldsymbol x}) = \int \frac{d^3 k}{(2 \pi)^3} F ({\boldsymbol k}) e^{i {\boldsymbol k}\cdot {\boldsymbol x}} . \end{equation*}$$

References

[1]

Sasaki

Suyama

Tanaka

, and

Yokoyama

Phys. Rev. Lett.

117

061101

(

2016

); [erratum]. (

https://doi.org/10.1103/PhysRevLett.117.061101

)

121

059901

(

2018

) (

https://doi.org/10.1103/PhysRevLett.121.059901

)

Crossref

Search ADS

PubMed

[2]

Ali-Haïmoud,

Kovetz,

E. D.

and

Kamionkowski,

Phys. Rev. D

123523

(

2017

). (

https://doi.org/10.1103/PhysRevD.96.123523

)

Crossref

Search ADS

[3]

Raidal

Vaskonen,

and

Veermäe,

J. Cosmol. Astropart. Phys.

1709

037

(

2017

). (

https://doi.org/10.1088/1475-7516/2017/09/037

)

Crossref

Search ADS

[4]

Ballesteros

Serpico,

P. D.

and

Taoso,

J. Cosmol. Astropart. Phys.

1810

043

(

2018

). (

https://doi.org/10.1088/1475-7516/2018/10/043

)

Crossref

Search ADS

[5]

Bringmann

Depta

P. F.

Domcke

, and

Schmidt-Hoberg

Phys. Rev. D

063532

(

2019

). (

https://doi.org/10.1103/PhysRevD.99.063532

)

Crossref

Search ADS

[6]

Ding

Nakama

Silk

, and

Wang

, arXiv:1903.07337 [astro-ph.CO] [Search inSPIRE].

[7]

Dokuchaev,

V. I.

Eroshenko,

Yu. N.

and

Rubin,

S. G.

Astron. Rep.

779

(

2008

). (

https://doi.org/10.1134/S1063772908100016

)

Crossref

Search ADS

[8]

Belotsky

K. M.

Dokuchaev,

V. I.

Eroshenko,

Y. N.

Esipova,

E. A.

Khlopov,

M. Yu.

Khromykh,

L. A.

Kirillov,

A. A.

Nikulin,

V. V.

Rubin,

S. G.

and

Svadkovsky,

I. V.

Eur. Phys. J. C

246

(

2019

). (

https://doi.org/10.1140/epjc/s10052-019-6741-4

)

Crossref

Search ADS

[9]

Clesse

and

García-Bellido,

Phys. Dark Univ.

142

(

2017

). (

https://doi.org/10.1016/j.dark.2016.10.002

)

Crossref

Search ADS

[10]

García-Bellido

and

Clesse,

Phys. Dark Univ.

144

(

2018

). (

https://doi.org/10.1016/j.dark.2018.01.001

)

Crossref

Search ADS

[11]

Meszaros,

Astron. Astrophys.

(

1975

). (

http://adsabs.harvard.edu/full/1975A%26A....38....5M

)

[12]

Carr

B. J.

Astrophys. J.

201

(

1975

). (

https://doi.org/10.1086/153853

)

Crossref

Search ADS

[13]

Chisholm

J. R.

Phys. Rev. D

083504

(

2006

). (

https://doi.org/10.1103/PhysRevD.73.083504

)

Crossref

Search ADS

[14]

Afshordi,

McDonald,

and

Spergel,

D. N.

Astrophys. J.

594

L71

(

2003

). (

https://doi.org/10.1086/378763

)

Crossref

Search ADS

[15]

Gong

J.-O.

and

Kitajima

J. Cosmol. Astropart. Phys.

1708

017

(

2017

). (

https://doi.org/10.1088/1475-7516/2017/08/017

)

Crossref

Search ADS

[16]

Gong

J.-O.

and

Kitajima

J. Cosmol. Astropart. Phys.

1811

041

(

2018

). (

https://doi.org/10.1088/1475-7516/2018/11/041

)

Crossref

Search ADS

[17]

Ali-Haïmoud,

Phys. Rev. Lett.

121

081304

(

2018

). (

https://doi.org/10.1103/PhysRevLett.121.081304

)

Crossref

Search ADS

PubMed

[18]

Desjacques

and

Riotto

Phys. Rev. D

123533

(

2018

). (

https://doi.org/10.1103/PhysRevD.98.123533

)

Crossref

Search ADS

[19]

Tada

and

Yokoyama

Phys. Rev. D

123534

(

2015

). (

https://doi.org/10.1103/PhysRevD.91.123534

)

Crossref

Search ADS

[20]

Young

and

Byrnes

C. T.

J. Cosmol. Astropart. Phys.

1504

034

(

2015

). (

https://doi.org/10.1088/1475-7516/2015/04/034

)

Crossref

Search ADS

[21]

Franciolini,

Kehagias,

Matarrese,

and

Riotto,

J. Cosmol. Astropart. Phys.

1803

016

(

2018

). (

https://doi.org/10.1088/1475-7516/2018/03/016

)

Crossref

Search ADS

[22]

Matarrese

Lucchin

, and

Bonometto

S. A.

Astrophys. J.

310

L21

(

1986

). (

https://doi.org/10.1086/184774

)

Crossref

Search ADS

[23]

Matarrese

and

Verde

Astrophys. J.

677

L77

(

2008

). (

https://doi.org/10.1086/587840

)

Crossref

Search ADS

[24]

Gong

J.-O.

and

Yokoyama

Mon. Not. Roy. Astron. Soc.

417

L79

(

2011

). (

https://doi.org/10.1111/j.1745-3933.2011.01124.x

)

Crossref

Search ADS

[25]

Young

Byrnes

C. T.

, and

Sasaki

J. Cosmol. Astropart. Phys.

1407

045

(

2014

). (

https://doi.org/10.1088/1475-7516/2014/07/045

)

Crossref

Search ADS

[26]

Nakama

Harada,

Polnarev,

A. G.

and

Yokoyama,

J. Cosmol. Astropart. Phys.

1401

037

(

2014

). (

https://doi.org/10.1088/1475-7516/2014/01/037

)

Crossref

Search ADS

[27]

Harada

Yoo

C.-M.

Nakama

, and

Koga

Phys. Rev. D

084057

(

2015

). (

https://doi.org/10.1103/PhysRevD.91.084057

)

Crossref

Search ADS

[28]

Byrnes

C. T.

Sasaki

, and

Wands

Phys. Rev. D

123519

(

2006

). (

https://doi.org/10.1103/PhysRevD.74.123519

)

Crossref

Search ADS

[29]

Suyama

and

Yamaguchi

Phys. Rev. D

023505

(

2008

). (

https://doi.org/10.1103/PhysRevD.77.023505

)

Crossref

Search ADS

[30]

Akrami

et al. [Planck Collaboration], arXiv:1807.06211 [astro-ph.CO] [Search INSPIRE].

[31]

Cabass

Pajer

, and

Schmidt

J. Cosmol. Astropart. Phys.

1701

003

(

2017

). (

https://doi.org/10.1088/1475-7516/2017/01/003

)

Crossref

Search ADS

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Download all slides

Month:	Total Views:
October 2019	169
November 2019	62
December 2019	43
January 2020	31
February 2020	19
March 2020	7
April 2020	11
May 2020	14
June 2020	10
July 2020	25
August 2020	13
September 2020	9
October 2020	8
November 2020	8
December 2020	13
January 2021	7
February 2021	7
March 2021	14
April 2021	15
May 2021	16
June 2021	8
July 2021	14
August 2021	13
September 2021	6
October 2021	12
November 2021	9
December 2021	2
January 2022	7
February 2022	16
March 2022	7
April 2022	8
May 2022	7
June 2022	6
July 2022	10
August 2022	12
September 2022	13
October 2022	7
November 2022	39
December 2022	53
January 2023	33
February 2023	24
March 2023	3
April 2023	14
May 2023	21
June 2023	3
July 2023	10
August 2023	11
September 2023	27
October 2023	8
November 2023	15
December 2023	13
January 2024	11
February 2024	6
March 2024	14
April 2024	10
May 2024	17
June 2024	9
July 2024	13
August 2024	5
September 2024	5
October 2024	12
November 2024	15
December 2024	2
January 2025	10
February 2025	4
March 2025	7
April 2025	7
May 2025	1

Article Contents

Clustering of primordial black holes with non-Gaussian initial fluctuations

Abstract

1. Introduction

2. Formulation

2.1. Probability of PBH formation

2.2. |$P_1$| and |$P_2$| in terms of the |$n$|-point correlators of the primordial fluctuations

2.3. PBH two-point correlation function

3. PBH clustering with local-type non-Gaussianity

3.1. Primordial local-type non-Gaussianity

3.2. PBH power spectrum with non-Gaussian primordial fluctuations

4. Conclusion

Acknowledgements

Footnotes

References

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Clustering of primordial black holes with non-Gaussian initial fluctuations

Abstract

1. Introduction

2. Formulation

2.1. Probability of PBH formation

2.2. |$P_1$| and |$P_2$| in terms of the |$n$|-point correlators of the primordial fluctuations

2.3. PBH two-point correlation function

3. PBH clustering with local-type non-Gaussianity

3.1. Primordial local-type non-Gaussianity

3.2. PBH power spectrum with non-Gaussian primordial fluctuations

4. Conclusion

Acknowledgements

Footnotes

References

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only