Clustering of primordial black holes with non-Gaussian initial fluctuations

TL;DR

This paper investigates how primordial black hole clustering depends on the nature of initial fluctuations, showing Gaussian fluctuations do not induce clustering on large scales, but non-Gaussianity can enhance it.

Contribution

It formulates the two-point correlation function of PBHs at formation using a functional integration approach, highlighting the impact of non-Gaussianity on clustering.

Findings

Gaussian fluctuations do not induce super-Hubble scale clustering.

Non-Gaussianity via local-type trispectrum enhances PBH clustering.

Clustering is sensitive to the primordial curvature perturbation statistics.

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.

Equations

$${\theta }_{\mathrm{local}}(\mathbf{x}):=\int d^{3}y{W}_{\mathrm{local}}(\mathbf{x}-\mathbf{y})\theta (\mathbf{y}),$$

$${\theta }_{\mathrm{local}}\ge {\theta }_{\mathrm{th}}.$$

$$P_1(\mathbf{x})=\int [D\theta ]P[\theta ]{\int }_{{\theta }_{\mathrm{th}}}^{\infty }d\alpha {\delta }_D({\theta }_{\mathrm{local}}(\mathbf{x})-\alpha ),$$

$$P_2({\mathbf{x}}_1,{\mathbf{x}}_2)=\int [D\theta ]P[\theta ]{\int }_{{\theta }_{\mathrm{th}}}^{\infty }d{\alpha }_1{\delta }_D({\theta }_{\mathrm{local}}({\mathbf{x}}_1)-{\alpha }_1){\int }_{{\theta }_{\mathrm{th}}}^{\infty }d{\alpha }_2{\delta }_D({\theta }_{\mathrm{local}}({\mathbf{x}}_2)-{\alpha }_2).$$

$${\delta }_D(x)=\int \frac{d\varphi }{2\pi }{e}^{i\varphi x},$$

$$P_1(\mathbf{x})=\int [D\theta ]P[\theta ]{\int }_{{\theta }_{\mathrm{th}}}^{\infty }d\alpha {\int }_{-\infty }^{\infty }\frac{d\varphi }{2\pi }\exp \left[i\varphi \int d^{3}y{W}_{\mathrm{local}}(\mathbf{x}-\mathbf{y})\theta (\mathbf{y})-i\varphi \alpha \right]$$

$$P_2({\mathbf{x}}_1,{\mathbf{x}}_2)$$

$$Z[J]:=\int [D\theta ]P[\theta ]\exp \left[i\int d^{3}yJ(\mathbf{y})\theta (\mathbf{y})\right]=⟨\exp \left[i\int d^{3}yJ(\mathbf{y})\theta (\mathbf{y})\right]⟩,$$

$$\xi_{\theta(c)}({\bm{x}}_{1},{\bm{x}}_{2},...,{\bm{x}}_{n}):=\frac{1}{i^{n}}\frac{\delta^{n}\log Z[J]}{\delta J({\bm{x}}_{1})\delta J({\bm{x}}_{2})\cdots\delta J({\bm{x}}_{n})}\biggr{|}_{J=0}\,.$$

$$\log Z[J]=\sum _{n=1}^{\infty }\frac{i^n}{n!}\int d^3{y}_1{d}^3{y}_2\cdots d^3{y}_n{\xi }_{\theta (c)}({\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_n)J({\mathbf{y}}_1)J({\mathbf{y}}_2)\cdots J({\mathbf{y}}_n).$$

$$P_1(\mathbf{x})$$

$$=$$

$${\xi }_{\mathrm{local}(c)}^{(n)}:=\int d^3{y}_1{d}^3{y}_2\cdots d^3{y}_n{\xi }_{\theta (c)}({\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_n)\prod _{r=1}^n{W}_{\mathrm{local}}(\mathbf{x}-{\mathbf{y}}_r).$$

$$P_2({\mathbf{x}}_1,{\mathbf{x}}_2)$$

$${\xi }_{\mathrm{local}(c),m}^{(n)}$$

$$=$$

$${\int }_{-\infty }^{\infty }\frac{d\varphi }{2\pi }{e}^{-i\varphi \alpha }\exp \left[\sum _{n=2}^{\infty }\frac{i^n}{n!}{\varphi }^n{\xi }_{\mathrm{local}(c)}^{(n)}\right]$$

$$P_1=\frac{1}{\sqrt{2\pi }}{\int }_{\nu }^{\infty }dw\exp \left[\sum _{n=3}^{\infty }\frac{(-1{)}^n}{n!}\frac{{\xi }_{\mathrm{local}(c)}^{(n)}}{{\sigma }_{\mathrm{local}}^n}\frac{d^n}{d{w}^n}\right]e^{-\frac{w^2}{2}},$$

$$P_2({\mathbf{x}}_1,{\mathbf{x}}_2)$$

$$P_1$$

$$P_2$$

$$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 }_{\mathrm{local}(c)}^{(n)}}{{\sigma }_{\mathrm{local}}^n}{H}_n\left(\frac{\nu }{\sqrt{2}}\right)\right],$$

$$P_2({\mathbf{x}}_1,{\mathbf{x}}_2)$$

$${\xi }_{\mathrm{PBH}}({\mathbf{x}}_1,{\mathbf{x}}_2)$$

$$\frac{{\xi }_{\mathrm{local}(c),m}^{(n)}}{{\sigma }_{\mathrm{local}}^n}\ll \frac{2^{n/2}m!(n-m)!}{H_m(\nu /\sqrt{2})H_{n-m}(\nu /\sqrt{2})},$$

$$\frac{{\xi }_{\mathrm{local}(c)}^{(2)}({\mathbf{x}}_1,{\mathbf{x}}_2)}{{\sigma }_{\mathrm{local}}^2}\ll \frac{2}{H_1(\nu /\sqrt{2}{)}^2}.$$

$${\xi }_{\mathrm{PBH}}({\mathbf{x}}_1,{\mathbf{x}}_2)$$

$$\delta (\mathbf{x})=-\frac{4(1+w)}{3w+5}{e}^{-5{\mathcal{R}}_c(\mathbf{x})/2}\frac{△}{a^2{H}^2}(e^{{\mathcal{R}}_c(\mathbf{x})/2}),$$

$$\delta (\mathbf{x})=-\frac{4}{9}\frac{1}{a^2{H}^2}△{\mathcal{R}}_c(\mathbf{x}).$$

$$⟨\delta (\mathbf{x})\delta (\mathbf{y})⟩={\left(\frac{4}{9a^2{H}^2}\right)}^2{△}^2⟨{\mathcal{R}}_c(\mathbf{x}){\mathcal{R}}_c(\mathbf{y})⟩.$$

$$\frac{⟨{\delta }^2(\mathbf{x}){\delta }^2(\mathbf{y})⟩}{{⟨{\delta }^2(\mathbf{x})⟩}^2}-1=2{⟨\delta (\mathbf{x})\delta (\mathbf{y})⟩}^2,$$

Full text

**Clustering of primordial black holes with non-Gaussian initial fluctuations **

Teruaki Suyamaa, and Shuichiro Yokoyamab,c

a* Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan

b* Kobayashi Maskawa Institute, Nagoya University, Aichi 464-8602, Japan

c* Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

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 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 have assumed that the distribution of PBHs is spatially uniform. As discussed in Ref. [3] (see also [4, 5, 6]),#1#1#1 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 the binary formation but also on the formation of the supermassive BHs. 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.#2#2#2There have been several works about late-time clustering of PBHs, e.g. in dark matter halos (see, e.g., Refs. [9, 10]). 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 also Refs. [15, 16] estimated an expected constraint on it by making use of future 21cm 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, Ali-Haimoud:2018dauDesjacques:2018wuu] 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 would be 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, 23, 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}({\bm{x}})$ as

[TABLE]

where the window function $W_{\rm local}({\bm{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 the formation. In this section, to maintain generality we do not specify a particular gauge for defining the primordial fluctuations $\theta({\bm{x}})$, and assume the criterion for PBH formation is given by#3#3#3Strictly 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.

[TABLE]

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

2.1 Probability of PBH formation

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

[TABLE]

where $P[\theta]$ is a probability distribution function for the primordial fluctuations, $\theta({\bm{x}})$; also, the probability that two points ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$ are PBHs can be given by

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Choosing $J({\bm{y}}):=\phi~{}W_{\rm local}({\bm{x}}-{\bm{y}})$, the one-point probability of PBH formation, $P_{1}$, can be written as

[TABLE]

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

[TABLE]

In the same way, the two-point probability of PBH formation, $P_{2}$, can be expressed as

[TABLE]

where

[TABLE]

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

[TABLE]

and then the expression for $P_{1}$ can be reduced to

[TABLE]

where $w:=\alpha/\sigma_{\rm local}$, $\nu:=\theta_{\rm th}/\sigma_{\rm local}$, and $\sigma_{\rm local}:=\left(\xi^{(2)}_{{\rm local}(c)}\right)^{1/2}$. We can also obtain a reduced form for $P_{2}$:

[TABLE]

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

[TABLE]

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)}({\bm{x}}_{1},{\bm{x}}_{2})/\sigma_{\rm local}^{2}\ll 1$, we can also obtain

[TABLE]

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 the PBH formation, we can perform the integration approximately in terms of $w$ and obtain approximate expressions for $P_{1}$ and $P_{2}$ as

[TABLE]

and

[TABLE]

Then, by using these expressions for $P_{1}$ and $P_{2}$, the two-point correlation function of the PBHs can be evaluated as

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

This expression corresponds to the result obtained in Refs. [23, 24] in the context of the halo bias. Under the assumptions given by Eqs. (45) and (46), this equation is valid for $\xi_{\rm PBH}({\bm{x}}_{1},{\bm{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({\bm{x}})$, in the above formulation, we use the primordial curvature perturbations on the comoving slice, often denoted by ${\mathcal{R}}_{c}({\bm{x}})$. In the long-wavelength approximation, comoving density fluctuations can be given in terms of ${\mathcal{R}}_{c}({\bm{x}})$ as  [27]

[TABLE]

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({\bm{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}({\bm{x}})}\to a$ [25], at the linear order, the above expression can be reduced to

[TABLE]

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

[TABLE]

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

Due to the locality, at the leading order in $\delta$ the PBH abundance at point ${\bm{x}}$ would be determined by the local variance $\langle\delta^{2}({\bm{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

[TABLE]

and it should be suppressed on super-Hubble scales, as can be seen by Eq. (56). 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

[TABLE]

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

[TABLE]

For super-Hubble distance $|{\bm{x}}-{\bm{y}}|\gg{(aH)}^{-1}$, we obtain

[TABLE]

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 a 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

[TABLE]

and the right-hand side of Eq. (60) is obtained by ignoring terms proportional to the two-point correlation function of ${\mathcal{R}}_{c}$, $\langle{\cal R}_{c}({\bm{x}}){\cal R}_{c}({\bm{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. (60) becomes

[TABLE]

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. (LABEL:eq:PBHcorr), 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#4#4#4We use the Fourier transform expression given by

$f({\bm{x}})=\int\frac{d^{3}k}{(2\pi)^{3}}F({\bm{k}})e^{i{\bm{k}}\cdot{\bm{x}}}.$

(63)

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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where ${\bm{r}}:={\bm{x}}_{1}-{\bm{x}}_{2}$. Substituting Eq. (LABEL:eq:PBHcorr) into this equation and noting that we apply the comoving density fluctuations, $\delta({\bm{x}})$, to $\theta_{\rm local}({\bm{x}})$ used in the previous section, we obtain

[TABLE]

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.

For primordial local-type non-Gaussianity as defined by Eq. (70), $P_{\delta},B_{\delta},~{}$and $T_{\delta}$ are respectively given in terms of the non-linearity parameters and the power spectrum of the primordial curvature perturbations as

[TABLE]

Substituting Eq. (94) into Eq. (LABEL:eq:Ppbh1), we have

[TABLE]

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 (LABEL:eq:PBHcorr). 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. (LABEL:eq:PBHcorr) (or Eqs. (38) and (41)). 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 $|{\bm{k}}+{\bm{p}}|\to p$, the above expression becomes

[TABLE]

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

[TABLE]

where a constant, $C$, corresponds to the $k$-independent terms in Eq. (130). 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. (70) is valid for CMB scales which are super-Hubble ones at the PBH formation. Although the slow-roll condition would be violated when 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. (70) 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.

Acknowledgments

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 JSPS Grant-in-Aid for Young Scientists (B) No.15K17632, and by the MEXT Grant-in-Aid for Scientific Research on Innovative Areas Nos.15H05888, 17H06359, 18H04338, and 19K03864.