Pulsar Timing Array Constraints on the Induced Gravitational Waves

TL;DR

This paper investigates how primordial black holes and their associated gravitational waves, detectable by pulsar timing arrays, can constrain early universe models, especially considering non-Gaussian curvature perturbations.

Contribution

It demonstrates that non-Gaussianity in curvature perturbations can alleviate tensions between induced gravitational wave predictions and PTA constraints.

Findings

Induced GWs from primordial black holes are in mild tension with current PTA constraints.

Non-Gaussianity with $f_{NL} extgreater 10$ can relax this tension.

Future observations by the Square Kilometer Array will be able to detect these GWs within a decade.

Abstract

If the black holes detected by LIGO/VIRGO are primordial black holes (PBHs) sourcing from a large primordial curvature perturbation on small scales, the corresponding induced gravitational waves (GWs) would peak at nanohertz that is detectable by the current and future observations of pulsar timing array (PTA). In this paper we show that with the mass function estimated from the merger rate of LIGO O1 and O2 events, the induced GWs from such a curvature perturbation with a Gaussian narrow peak at some small scale would be in a seemingly mild tension with current constraints from PTA. However, if the curvature perturbation is of local-type non-Gaussianity with a non-linear parameter $f_\text{NL}\gtrsim\mathcal{O}(10)$, the tension could be relieved. Nevertheless, such an induced GWs must be detectable by the Square Kilometer Array in a decade or less.

Equations

$$f_{\text{peak}}\sim 6.7\times 10^{-9}{\left(\frac{M_{\text{PBH}}}{M_{\odot }}\right)}^{-1/2}\text{Hz}.$$

$${\beta }_{m_{\mathrm{H}}}=2{\int }_{{\Delta }_c}^{\infty }\frac{m}{m_{\mathrm{H}}}P(\Delta )d\Delta ,$$

$$P(\Delta )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{{\Delta }^2}{2{\sigma }^2}\right),$$

$${\sigma }^2={\int }_0^{\infty }\frac{dq}{q}{\tilde{W}}^2(R,q){\mathcal{P}}_{\Delta }(\eta ,q),$$

$$\Delta (t,q)=\frac{2(1+w)}{5+3w}{\left(\frac{q}{aH}\right)}^2\mathcal{R}(q),$$

$${\sigma }^2={\int }_0^{\infty }\frac{dq}{q}{\tilde{W}}^2(R,q){\left(\frac{4}{9}\right)}^2{\left(\frac{q}{aH}\right)}^4{\mathcal{P}}_{\mathcal{R}}(\eta ,q),$$

$${\sigma }^2(k)={\int }_0^{\infty }\frac{dq}{q}{\tilde{W}}^2(R=k^{-1},q){\left(\frac{4}{9}\right)}^2{\left(\frac{q}{k}\right)}^4{T}^2(\eta =k^{-1},q){\mathcal{P}}_{\mathcal{R}}(q),$$

$${\beta }_{\mathrm{eq}}=\frac{a_{\mathrm{eq}}}{a_{m_{\mathrm{H}}}}{\beta }_{m_{\mathrm{H}}},$$

$${\Omega }_{\mathrm{PBH},\mathrm{eq}}={\int }_{-\infty }^{\infty }d\ln m_{\mathrm{H}}\frac{a_{\mathrm{eq}}}{a_{m_{\mathrm{H}}}}{\beta }_{m_{\mathrm{H}}}.$$

$$f(m)\equiv \frac{1}{{\Omega }_{\mathrm{CDM}}}\frac{d{\Omega }_{\mathrm{PBH}}}{d\ln m}=\frac{{\Omega }_m}{{\Omega }_{\mathrm{CDM}}}\frac{d{\Omega }_{\mathrm{PBH},\mathrm{eq}}}{d\ln m},$$

$$f_{\mathrm{PBH}}\equiv \frac{{\Omega }_{\mathrm{PBH}}}{{\Omega }_{\mathrm{CDM}}}={\int }_{-\infty }^{\infty }f(m)d\ln m.$$

$$m=m_{\mathrm{H}}K(\Delta -{\Delta }_c{)}^{\gamma },$$

$$f(m)$$

$$\Delta ={\left(\frac{m}{K{m}_{\mathrm{H}}}\right)}^{1/\gamma }+{\Delta }_c.$$

$$d{s}^2=a^2(\eta )\left\{-(1+2\Phi )d{\eta }^2+\left[(1-2\Phi ){\delta }_{ij}+\frac{1}{2}{h}_{ij}\right]d{x}^{i}d{x}^j\right\},$$

$$h_{ij}(\eta ,\mathbf{x})=\int \frac{d^{3}k}{(2\pi {)}^{3/2}}\sum _{\lambda =+,\times }{e}_{ij}^{\lambda }(\hat{k})h_{\mathbf{k},\lambda }(\eta )e^{i\mathbf{k}\cdot \mathbf{x}},$$

$$h_{\mathbf{k}}^{\prime \prime }+2\mathcal{H}{h}_{\mathbf{k}}^{\prime }+k^2{h}_{\mathbf{k}}=2{\mathcal{P}}_{ij}^{lm}{e}^{ij}{T}_{lm}(\mathbf{k},\eta ),$$

$$T_{lm}=-2\Phi {\partial }_l{\partial }_m\Phi +{\partial }_l\left(\Phi +{\mathcal{H}}^{-1}{\Phi }^{\prime }\right){\partial }_m\left(\Phi +{\mathcal{H}}^{-1}{\Phi }^{\prime }\right).$$

$$h_{\mathbf{k}}(\eta )=\frac{2}{a(\eta )}{\int }_0^{\infty }d\tilde{\eta }a(\tilde{\eta })G_k(\eta ;\tilde{\eta }){\mathcal{P}}_{ij}^{lm}(\hat{k})e^{ij}{T}_{lm}(\tilde{\eta },\mathbf{k}).$$

$$⟨h_{\mathbf{k}}(\eta )h_{\mathbf{l}}(\eta )⟩={\delta }^{(3)}(\mathbf{k}+\mathbf{l})\frac{2{\pi }^2}{k^3}{\mathcal{P}}_h(\eta ,k).$$

$${\Omega }_{\mathrm{GW}}(\eta ,k)=\frac{1}{24}{\left(\frac{k}{\mathcal{H}(\eta )}\right)}^2\overline{{\mathcal{P}}_h(\eta ,k)},$$

$${\Omega }_{\mathrm{GW}}(k,\eta )=\frac{1}{54}{\int }_0^{\infty }\mathrm{d}v{\int }_{|1-v|}^{1+v}\mathrm{d}u\mathcal{T}(u,v){\mathcal{P}}_{\mathcal{R}}(ku){\mathcal{P}}_{\mathcal{R}}(kv),$$

$$\mathcal{T}(u,v)$$

$$\times \frac{1}{2}\left\{{\left[-4uv+(u^2+v^2-3)\ln \left|\frac{3-(u+v{)}^2}{3-(u-v{)}^2}\right|\right]}^2+{\left[\pi (u^2+v^2-3){\Theta }_{\frac{1}{2}}(u+v-\sqrt{3})\right]}^2\right\}.$$

$${\mathcal{P}}_{\mathcal{R}}(k)=\frac{{\mathcal{A}}_{\mathcal{R}}{k}_{\ast }}{\sqrt{2\pi }{\sigma }_{\ast }}\exp \left[-\frac{(k-k_{\ast }{)}^2}{2{\sigma }_{\ast }^2}\right],$$

$${\sigma }^2(k)=\frac{16{\mathcal{A}}_{\mathcal{R}}}{\sqrt{2\pi }{ϵ}_{\ast }}{\int }_0^{\infty }\frac{dz}{z}\exp (-z^2){\left[\frac{\sqrt{3}}{z}\sin (\frac{z}{\sqrt{3}})-\cos (\frac{z}{\sqrt{3}})\right]}^2\exp \left[-\frac{(xz-1{)}^2}{2{ϵ}_{\ast }^2}\right],$$

$${\sigma }^2(k)\approx$$

$$\displaystyle-\sqrt{2\pi}x\ \mathrm{e}^{\frac{x^{2}}{2x^{2}\epsilon_{*}^{2}+4\epsilon_{*}^{4}}}\left(x^{2}\left(3\epsilon_{*}^{2}+1\right)+6\epsilon_{*}^{4}\right)Q\left(-\frac{1}{2},\frac{x^{2}}{2x^{2}\epsilon_{*}^{2}+4\epsilon_{*}^{4}}\right)\Biggr{\}},$$

$$Q(a,z)\equiv \frac{\Gamma (a,z)}{\Gamma (a)},\Gamma (a,z)\equiv {\int }_z^{\infty }{t}^{a-1}{\mathrm{e}}^{-t}dt,\Gamma (a)\equiv {\int }_0^{\infty }{t}^{a-1}{\mathrm{e}}^{-t}dt.$$

$$m_{\mathrm{H}}=\frac{4\pi }{3}\rho H^{-3}.$$

$$g_{\ast s}(T)T^3{a}^3=\mathrm{constant},$$

Full text

Preprint: IPMU19-0095

1 Introduction

LIGO/VIRGO have detected gravitational waves (GWs) from mergers of binary black holes (BBHs) and binary neutron stars [1, 2, 3, 4, 5, 6]. It is possible that these black holes we observed are of the primordial origin [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], which are formed by gravitational collapse when the curvature perturbation is of $\mathcal{O}(1)$ at the horizon reentry [21, 22, 23, 24, 25]. During inflation, the primordial curvature perturbation $\mathcal{R}$ is generated by quantum fluctuations, which is then frozen when it is stretched outside the horizon by inflation, and finally seeds the temperature anisotropies and the large scale structure inhomogeneities on large scales observed today. It is already well constrained by the observations of the cosmic microwave background (CMB) that $\mathcal{R}$ is nearly scale invariant and Gaussian, and has an amplitude of order $10^{-5}$ on scales roughly larger than $1\text{Mpc}$ [26, 27]. However, the curvature perturbation on small scales is unknown, because the resolution of the CMB experiments is limited, and the other observations are not so accurate. Substantial primordial black holes (PBHs) can be the seeds for galaxy formation [28, 29, 30, 31], the dark matter candidate [32, 33, 34, 35, 36], or the sources of LIGO/VIRGO detection, depending on their masses and abundance at formation, which can be constrained by PBH remnants that survive Hawking radiation, star-capture processes, microlensing, CMB $\mu$-distortion, and so on.

Such a peak in the curvature perturbation is possible in some inflation models, which has attracted much attention recently [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55]. It is also well known that the large scalar perturbation can induce large tensor perturbations since they are coupled at the nonlinear order. If the curvature perturbation becomes large at small scales, the tensor perturbation it couples to at the nonlinear order gets also large, which can evolve to today as an isotropic stochastic bacground of GWs [56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88]. The peak frequency of the GW spectrum is solely determined by the peak scale of the curvature perturbation spectrum, which is connected to the PBH mass generated from this peak by

[TABLE]

If the black holes detected by LIGO/VIRGO are primordial, we must have such accompanied stochastic background of GWs induced by such scalar perturbations, which peaks around nHz. The peak amplitude of the spectrum for the stochastic GW is determined by the peak value of the scalar power spectrum that determines the PBH abundance, which in turn can be constrained by the black hole merger rate from LIGO/VIRGO events.

Millisecond pulsars have very stable rotation periods, thus the pulses we observe can be used as a tool to detect the spacetime metric between the observer and pulsars. By studying the time of arrivals from many pulsars in the sky, we can extract the information of the stochastic GWs from this pulsar timing array (PTA) [89, 90]. Currently, there are three PTA projects: the European Pulsar Timing Array (EPTA [91]), the Parkes Pulsar Timing Array (PPTA [92]), and the North American Observatory for Gravitational Waves (NANOGrav [93]), while the International Pulsar Timing Array (IPTA [94]) is the combination of them. The GW frequency range is set by the total observation time, which is roughly $f\sim 1/T_{\text{tot}}\sim 10^{-9}~{}\text{to}~{}10^{-7}~{}\mathrm{Hz}$, right in the range for the induced GWs accompanied with the solar-mass PBHs. For now there is no evidence of stochastic GW background. For instance, 2015 EPTA data sets an upper bound $\Omega_{\text{GW}}h^{2}<1.1\times 10^{-9}$ at $2.8~{}\text{nHz}$ [95] at $95\%$ confidence level. However, this is in conflict with what the LIGO detection has predicted, if the black holes detected by LIGO/VIRGO are primordial [96, 97].

Previous works assume that the scalar perturbation is Gaussian and has a narrow peak on a specific small scale. It is known that a broad peak can make the conflict even worse [98, 99]. However, introducing non-Gaussianity could be helpful. In the inflation models which predict a peak in the power spectrum of the curvature perturbation, it is quite natural to have non-negligible non-Gaussianity. The nonlinear relation between density contrast and the curvature perturbation provides another source of non-Gaussianity. As PBHs are produced at the large amplitude tail of the probability distribution of the curvature perturbation, positive non-negligible non-Gaussianity can enhance the abundance of the PBH formation. Equivalently, given the PBH abundance, a non-Gaussian scalar perturbation may induce smaller GWs, which may escape from the current PTA constraint. Ref. [100] has considered non-Gaussian scalar perturbations. By choosing a monochromatic mass function of PBH and an order estimation for the induced GWs, they found that the PTA constraint can be evaded if non-Gaussianity is large enough.

In this paper, we study this possibility carefully by reconstructing the mass function from the parameters of the LIGO events in O1 and O2, following the method established in [19] and calculating the GWs induced by the non-Gaussian scalar perturbations by the formula obtained previously in [75]. With the EPTA constraints, we found that it is still possible for the black holes detected by LIGO/VIRGO to be primordial if $f_{\text{NL}}\gtrsim\mathcal{O}(10)$.

This paper is organized as follows: in Section 2, we review the derivations of mass function of PBHs and energy-density spectrum of induced GWs; in Section 3, we calculate the mass function of PBHs which originate from Gaussian primordial curvature perturbations and find that it can be well parametrized; in Section 4, we make a maximum likelihood analysis to find the best-fit parameters by a fitting formula, and confirm the corresponding induced GW is incompatible with the constraints from PTA experiments; in Section 5, we discuss how to evade the conflict by considering the primordial curvature perturbations with local-type non-Gaussianity. We conclude in Section 6.

2 Mass function of PBHs and Induced GWs

In this section we review the PBH mass function and its associated induced GWs. Assuming the peak of the power spectrum for the curvature perturbation reenters the Hubble horizon in the radiation dominated era, PBH will form soon after this reentry time. The PBH mass fraction at formation $\beta_{m_{\mathrm{H}}}$ of horizon mass $m_{\mathrm{H}}$ is usually calculated by the Press-Schechter approach [101, 102, 103], which is given by222For an alternative approach based on the peak theory, see  [104, 105]

[TABLE]

where $\Delta=(\delta\rho-\rho)/\rho$ is the density contrast, $P(\Delta)$ is its probability distribution function, and $\Delta_{c}$ is the critical density contrast above which PBHs will form. Assuming $P(\Delta)$ is Gaussian,

[TABLE]

where $\sigma^{2}$ is the variance of density contrast given by integrating the density power spectrum $\mathcal{P}_{\Delta}(\eta,q)$,

[TABLE]

where $\tilde{W}(R,q)=\exp(-q^{2}R^{2}/2)$ is Fourier transform of a volume-normalised Gaussian window smoothing function, $R$ is horizon scale at a given conformal time $\eta$. For the subtlety of choosing the window function, see [106, 107]. Using the relation between curvature perturbation $\mathcal{R}$ and density contrast $\Delta$,

[TABLE]

where $w$ is the equation-of-state parameter, which is $1/3$ in the radiation dominated era, $\sigma^{2}$ can be written as

[TABLE]

where $\mathcal{P}_{\mathcal{R}}(\eta,q)$ is the power spectrum of curvature perturbations.

For a wavenumber $k$ corresponding to horizon reentry in the radiation dominated era, we have $k=\eta^{-1}=aH=\mathcal{H}$, $R=\mathcal{H}^{-1}=k^{-1}$, and then

[TABLE]

where $T(\eta,q)=3(\sin y-y\cos y)/y^{3}$ with $y\equiv q\eta/\sqrt{3}$ is the transfer function and $\mathcal{P}_{\mathcal{R}}(q)$ is the power spectrum of primordial curvature perturbations.

It is worth noting that the energy density of PBHs evolves like matter during radiation dominated era thus increases proportional to $a$ until equality, and $\beta_{m_{\mathrm{H}}}$ represents the fraction of PBHs when horizon mass is $m_{\mathrm{H}}$, thus at equality, the fraction $\beta_{\mathrm{eq}}$ is

[TABLE]

where $a_{\mathrm{eq}}$ is the scale factor at equality and $a_{m_{\mathrm{H}}}$ is the scale factor when horizon mass is $m_{\mathrm{H}}$. Integrating $\beta_{\mathrm{eq}}$ over all possible PBH formation time, we obtain the total abundance of PBHs at equality,

[TABLE]

The mass function of PBHs $f(m)$ is defined as the fraction of cold dark matter (CDM) made up of PBHs with mass $m$,

[TABLE]

where $\Omega_{\mathrm{CDM}}/\Omega_{m}\approx 0.84$ [27]. The total fraction of PBHs in CDM is

[TABLE]

The relation of the PBH mass and the density contrast is given by the critical collapse [108, 109]

[TABLE]

where the constants $K=3.3$, $\gamma=0.36$, $\Delta_{c}=0.45$ during radiation-dominated era are given by numerical studies. Therefore, the mass function of PBHs can be written as

[TABLE]

with

[TABLE]

Another important phenomenon generated by the peak of the power spectrum for the curvature perturbation is the induced GWs. We will shortly review it here for completeness following [81, 110]. The perturbed metric in the conformal Newton gauge reads,

[TABLE]

where $\eta$ is the conformal time, $\Phi$ is the Newton potential and $h_{ij}$ is the tensor mode of the metric perturbation in the transverse-traceless gauge. As we neglect the anisotropic stress tensor, the scalar perturbations in time and diagonal spatial components are identical. $h_{ij}(\eta,\mathbf{x})$ can be expanded as follows,

[TABLE]

where $+$ and $\times$ are the two polarizations, respectively. The equation of motion for $h_{ij}$ can be derived from the perturbed Einstein equation up to the second order, which in the momentum space can be written as

[TABLE]

for each polarization, where $\mathcal{H}$ is the conformal Hubble parameter, $\mathscr{P}_{ij}^{lm}$ is the projection operator to the transverse-traceless part, and a prime denotes the derivative with respect to $\eta$. The source term $T_{lm}(\mathbf{k},\eta)$ is of second order in scalar perturbations,

[TABLE]

This equation can be solved with Green’s function method in the momentum space, and its solution is

[TABLE]

The two-point function of $h_{\mathbf{k}}$ can be written as

[TABLE]

It is usual to describe the stochastic GWs by its energy density per logarithmic frequency interval normalized by the critical density,

[TABLE]

where the two polarization modes of GWs have been summed over, and the overline denotes average over a few wavelengths. Using the relation between scalar perturbations $\Phi$ and curvature perturbations $\mathcal{R}$ in the radiation-dominated era, $\Phi=-(2/3)\mathcal{R}$, we can calculate the GW spectrum induced by the curvature perturbation during radiation-dominated era as

[TABLE]

with

[TABLE]

3 A fitting formula for the PBH mass function

The PBHs will form when the rare peaks of the curvature perturbations reenter the Hubble horizon in the radiation dominated era. However, the amplitude of the curvature perturbations on CMB scales are highly constrained, which makes the PBH abundance completely negligible if the curvature perturbations on small scales are extropolated linearly by its red-tilted power spectrum on CMB scales. If the BBHs detected by LIGO are composed of PBHs, we need substantial formation of such PBHs, which requires a relatively large amplitude of the curvature perturbations of $\mathcal{O}(1)$ on the corresponding small scales. This is consistent with the power spetrum on CMB scales if we also consider the running of running from the current observational constraints [111]. For the sake of convenience, such a peaked power spectrum can be described by a narrow Gaussian peak with a characteristic wavenumber $k_{*}$ and a width $\sigma_{*}\ll k_{*}$,

[TABLE]

where the near-scale-invariant part from the CMB scales is neglected.

By using (2.6), we can calculate the variance of density contrast,

[TABLE]

where $\epsilon_{*}\equiv\sigma_{*}/k_{*}$, $x\equiv k/k_{*}$, $z\equiv q/k$. As the window function will greatly suppress the contribution from the superhorizon modes of $z\gtrsim 1$, the main contribution of the integral comes from inside the horizon, which means that we can expand the integrand in (3.2) in the small $z$ limit. This gives rise to

[TABLE]

where $Q(a,z)$ is the regularized incomplete Gamma function, defined as

[TABLE]

In the left panel of Fig. 1, the exact $\sigma^{2}$ from Eq.(3.2) and approximate $\sigma^{2}$ from Eq.(3.3) with $\epsilon_{*}=0,0.05,0.1$ are shown in solid and dashed lines, respectively, which justifies our approximation.

In order to calculate the mass function of PBHs by Eq.(2.12), one needs to know $a_{\mathrm{eq}}/a_{m_{\mathrm{H}}}$, and the relation between the wavenumber corresponding to horizon reentry $k$ and the horizon mass $m_{\mathrm{H}}$. By definition, the horizon mass is

[TABLE]

We can then use the

[TABLE]

and the Friedmann equation in radiation dominated era

[TABLE]

where $g_{*s}$ is the effective degrees of freedom for the entropy density and $g_{*r}$ is the effective degrees of freedom for relativistic particles. Then we get

[TABLE]

and

[TABLE]

where $T_{m_{\mathrm{H}}}$ is temperature when horizon mass is $m_{\mathrm{H}}$, and dependence of the effective degree of freedom on the temperature is determined by the thermal history of the universe [112, 27, 79, 106]. Therefore, the mass function of PBHs can be computed as

[TABLE]

In the right panel of figure 1, the mass functions of PBHs $f(m)$ with $\mathcal{A}_{\mathcal{R}}=0.05$, $k_{*}=10^{6}\mathrm{Mpc}^{-1}$, $\epsilon_{*}=0,0.05,0.1$ are shown from top to bottom, respectively, which tells us that the effect of $\epsilon_{*}$ to $f(m)$ is negligible once $\epsilon_{*}\ll 1$. Recalling that we study the case that $\sigma_{*}\ll k_{*}$ ($\epsilon_{*}\ll 1$), taking $\epsilon_{*}=0.05$ for concreteness, there are only two free parameters $\mathcal{A}_{\mathcal{R}}$ and $k_{*}$ in $f(m)$. In this case, we find a fitting formula for $f(m)$ as

[TABLE]

and thus

[TABLE]

The exact $f(m)$ (Eq.(3.10)) can be fitted by the corresponding parametrized $f(m)$ (Eq.(3.11)) with $\mathcal{A}_{\mathcal{R}}=0.05,\ k_{*}=10^{5.5},10^{6.0},10^{6.5}\mathrm{Mpc}^{-1}$ and $\mathcal{A}_{\mathcal{R}}=0.04,~{}0.05,~{}0.06$, $k_{*}=10^{6}\mathrm{Mpc}^{-1}$ with high precision as shown in figure 2. We can freely convert $\mathcal{A}_{\mathcal{R}}$ and $k_{*}$ into $\hat{A}$ and $\hat{m}$, and vice versa.

4 Best-fit parameters from the observed merger rate

In the previous discussion, we assumed that the peak of the power spectrum for the curvature perturbation reenters Hubble horizon in the radiation dominated era. After their formation, PBHs are dynamically coupled to cosmic expansion. However, as the background energy density decays faster than that of the PBHs, the local density may become greater than ambient energy density, thus the PBHs decouple from the background evolution of the universe and form a gravitational bound state. The closest PBHs fall toward to each other, yet with the torques from a nearby third PBH as well as some other matter inhomogeneities, the head-on collision may be replaced by the formation of a binary of PBHs [113].

Taking into account the torques caused by the surrounding PBHs and linear density perturbations, the merger rate of PBH binary reads [114, 115, 116]

[TABLE]

where $M\equiv m_{1}+m_{2}$, $\mu\equiv m_{1}m_{2}/M^{2}$. We also define the normalized mass function of PBHs

[TABLE]

and the suppression factor

[TABLE]

where $\sigma_{\mathrm{eq}}\approx 0.005$ is the variance of the density perturbations of the ambient fluid of the universe at equality.

Assuming all the BBH mergers observed by LIGO/Virgo during first and second observing runs [117] are PBHs which originate from (3.1), and by following the method in [115], we take a maximum likelihood analysis to find the best-fit parameters of our fitting formula of PBH mass function (3.11). The log-likelihood function is

[TABLE]

where

[TABLE]

is the probability density of having a BBH consisting of two black holes with mass $m_{1}$ and $m_{2}$ merging at redshift $z$, and $V_{c}(z)$ is the comoving volume. The experimental uncertainties of black hole mass and redshift are accounted by $p^{(j)}(m^{(j)}|m)$ and $p^{(j)}(z^{j}|z)$, which denote the probability of observing black hole mass $m^{(j)}$ and redshift $z^{(j)}$ given that the true mass is $m$ and the true redshift is $z$. $p^{(j)}$ is taken to be a Gaussian distribution with mean and variance given in [19]. The detectability of detectors based on the signal-to-noise ratio (SNR) $\rho$ is implemented by step function $\theta$, and the detectability threshold is taken to be $\rho_{c}=8$. The results of log-likelihood $\mathcal{L}$ and normalized likelihood $L$ are shown in figure 3, and $\hat{m}=17.1_{-3.9}^{+5.1}$ at $99\%$ confidence level.

Since the maximum likelihood analysis can only find the best fit of $\hat{m}$, we use the number of BBH merger events that LIGO observed to constraint $\hat{A}$. The number of BBH merger events that LIGO should observe during period $\Delta t$ is

[TABLE]

During the first and second runs $\Delta t\simeq 165$ days, LIGO has observed $N_{\mathrm{obs}}=10$ BBH merger events [117], which can be converted into the constraint of $\hat{A}$. The results are shown in figure 4. The best-fit parameters are $\hat{m}=17.1M_{\odot}$ and $\hat{A}=4.59\times 10^{-9}$, which can be converted into $\mathcal{A}_{\mathcal{R}}=0.0404$ and $k_{*}=1.06\times 10^{6}\mathrm{Mpc}^{-1}$, and are compatible with the observational constraints to the fraction of PBHs in dark matter as is shown in figure LABEL:fig:fpbh.

Given the best-fit parameters $\mathcal{A}_{\mathcal{R}}=0.0404$ and $k_{*}=1.06\times 10^{6}\mathrm{Mpc}^{-1}$, we can calculate the spectrum of induced GWs by Eq.(2.21) until equality. After the source term re-enters the horizon, it will decay rapidly, leaving the induced GWs evolving as free propagating radiation decaying as $a^{-4}$. The GW spectrum today can be connected to the equality by

[TABLE]

We plot $\Omega_{\mathrm{GW}}(k,\eta_{0})h^{2}$ in figure 6 labelled by $F_{\text{NL}}=0$. It is easy to see that the best-fit configuration from LIGO detections is ruled out by the EPTA limits on stochastic GW background. But we must emphasize that there are only 10 samples for the detection, which are not enough for statistics. If more heavier BBHs are observed by LIGO in future, it is possible for the best-fit $k_{*}$ to be smaller, thus the peak of $\Omega_{\mathrm{GW}}(\eta_{0},k)h^{2}$ may be shifted to the left and be compatible with EPTA. In this paper, we will work under the assumption that, even if future detections of more LIGO/VIRGO events would not significantly change the currently observed merger rate, one could still be able to resolve this conflict simply by introducing some non-Gaussianity in the scalar perturbation, as we will show in the next section.

5 Non-Gaussian curvature perturbation

The primordial non-Gaussianity of the curvature perturbation may come from the self coupling of the inflaton, the coupling of inflaton to other fields, or a non-Bunch-Davies initial condition. The non-Gaussianity may enhance or suppress the formation of PBHs, depending on the signature of the nonlinear parameter $f_{\text{NL}}$. Besides, the curvature perturbation is not the quantity we can observe directly today. The density contrast $\Delta$ that some authors use in the calculation of PBH formation has a nonlinear relation with the curvature perturbation $\mathcal{R}$ on scales comparable to the Hubble horizon, which becomes crucial at the horizon reentry for the PBH formation [126, 127, 84, 128]. All of these introduce non-Gaussianity in the curvature perturbation, which can have characterisitc impacts on the induced GWs  [100, 129, 75, 76, 81, 130]. It can be either a suppression or an enhancement on the peak value of the GW spectrum, and the former case can help to solve the conflict we encountered in the last section.

Following [131, 132, 133, 134, 135, 136, 137], we consider a local-type non-Gaussian curvature perturbation $\mathcal{R}(\mathbf{x})$,

[TABLE]

where $F_{\text{NL}}$ is the nonlinear parameter for the curvature perturbation, relating to the nonlinear parameter $f_{\text{NL}}$ for the Newton potential $\Phi$ by $F_{\text{NL}}=(3/5)f_{\text{NL}}$. We will see that only positive $F_{\text{NL}}$ is useful for suppressing the GW spectrum, because it will increase the PBH abundance given the same variance. $\mathcal{R}_{\mathrm{G}}(\mathbf{x})$ is the Gaussian curvature perturbation, which can be solved from (5.1) as [138, 139]

[TABLE]

The probability distribution function of this full curvature perturbation $\mathcal{R}$ can be derived by requiring $\mathcal{R}_{\mathrm{G}}$ obeying the Gaussian distribution,

[TABLE]

with the Jacobian

[TABLE]

The initial PBH mass fraction can be written as

[TABLE]

where $\mathcal{R}_{c}\simeq\frac{9}{2\sqrt{2}}\Delta_{c}\simeq 1.4$ is the critical curvature perturbation for PBH formation [101], and $\langle\mathcal{R}_{\mathrm{G}}^{2}\rangle\simeq\mathcal{A}_{\mathcal{R}}$ for $\sigma_{*}\ll k_{*}$.

In figure 5, $\mathcal{A}_{\mathcal{R}}$ with different $\beta$ as a function of $F_{\text{NL}}$ are shown. Given a fixed $\beta$ from observations, we see that $F_{\text{NL}}\mathcal{A}_{\mathcal{R}}$ is constant while $F_{\text{NL}}\rightarrow\infty$,

[TABLE]

[TABLE]

where $\mathrm{Ierfc}(x)$ is the inverse function of $\mathrm{erfc}(x)$. For $\beta\ll 1$, we can approximate

[TABLE]

The GW spectrum from non-Gaussian curvature perturbations is given by [81, 75]

[TABLE]

where

[TABLE]

for $\sigma_{*}\ll k_{*}$. Roughly speaking, Eq.(5.9) can be estimated as

[TABLE]

where $\Omega^{(0)}(k)$, $\Omega^{(2)}(k)$, $\Omega^{(4)}(k)$ are the corresponding integral terms, and they are all of order $\mathcal{O}(1)$ at the peak. Therefore, fixing the PBH abundance $\beta$ that is determined by our fit, the GW spectrum induced by the Gaussian part of the curvature perturbation is $\Omega_{\text{GW}}^{(0)}\sim A_{0}^{2}$, while the extremely non-Gaussian part of GW spectrum is $\Omega_{\mathrm{GW}}^{f_{\text{NL}}\rightarrow\infty}\sim\mathcal{O}(A_{\infty}^{4})$. This gives

[TABLE]

with $A_{0}=0.0404$ obtained from the best-fit parameter in section 4. Therefore, we can see that the peak value of the GW spectrum is greatly suppressed, which makes it possible to evade the constraint from EPTA. Accurate calculations for $\Omega_{\text{GW}}$ from (5.9) is done and depicted in Fig. 6 for a fixed PBH abundance $\beta=3.4\times 10^{-12}$ obtained from Eq.(2.1), for $F_{\text{NL}}=0,~{}10,~{}100$ and $F_{\text{NL}}\rightarrow\infty$. It is shown explicitly that all black holes observed by LIGO can be PBHs, if the curvature perturbation is non-Gaussian with $F_{\text{NL}}\gtrsim 10$.

Here we would like to mention that the future radio telescope project Square Kilometer Array (SKA) [143] can extend the detecting accuracy to $\Omega_{\mathrm{GW}}\lesssim 10^{-14}$ at around $4~{}\text{nHz}$. We depict the expected SKA sensitivity curve in Fig. 6. It is easy to see that whether the LIGO/VIRGO detection events are mostly consists of PBHs can be easily checked by SKA. We can also see that as the non-Gaussian peak has a higher frequency, we do not need 20 years to see that peak. 5 to 10 years of observation will be enough to see whether the non-Gaussian peak in GW spectrum exists or not.

6 Conclusion and Discussions

The possibility that LIGO-detected black holes may be PBHs has attracted much attention since the PBH merger rate fits that of the LIGO estimation [8, 7]. These PBHs can be distinguished with black holes formed by astrophysical processes by studying their spin distribution [144, 145, 146, 147, 148], redshift dependence of the merger rate [149, 150], mass dependence of the merger rate [151, 152], and eccentricity [153, 154]. Currently there is no clear evidence for the primordial or astrophysical nature of the BBHs, yet there is some tension in the PTA constraints on stochastic GWs accompanied with the PBHs.

Solar-mass PBHs form substantially if the primordial curvature perturbation has a peak on comoving scale of $0.1~{}\text{pc}$, which also induces a stochastic background of GWs with its peak frequency within the reaches of PTA observations. Therefore, currently null detection of such GW background put a conservative bound on the peak amplitude of curvature perturbation and hence on the mass function of PBHs, which is in a seemingly mild tension with the best-fit parameters obtained by the maximum likelihood analysis of merger rate from the LIGO/VIRGO O1 and O2 events. To resolve this potential tension, we consider the GWs induced by curvature perturbations with local-type non-Gaussianity, and find that current constraints on the stochastic background of GWs from PTA observations could be evaded, if the nonlinear parameter is positive and larger than $\mathcal{O}(10)$.

We must emphasize that the statistic error may be large as the current number of samples given by the LIGO/VIRGO O1 and O2 events is limited. Also, it may be possible that some of the events are not primordial thus the induced GWs are overestimated. More merger observations are needed to determine the nature of the black holes and their mass function, which will be provided by the third run of LIGO/VIRGO and the upcoming KAGRA [155] and LIGO-India [156], and the future space-based experiments like LISA [157, 158, 159], Taiji [160], Tianqin [161], BBO [162, 163] and DECIGO [164, 165]. Our paper provide a possible solution to the tension if it still exists in the future statistics with more events, and we expect SKA can verify or falsify the primordial nature of the binary black holes, which is independent of the other criteria.

Acknowledgement

We thank Misao Sasaki and Ying-li Zhang for useful discussions. RGC and XYY are supported by the National Natural Science Foundation of China Grants Nos.11435006, 11647601, 11690022, 11821505, 11851302, and by the Strategic Priority Research Program of CAS Grant NO. XDB23010500 and No.XDB23030100, and by the Key Research Program of Frontier Sciences of CAS. SP is supported by the MEXT/JSPS KAKENHI No. 15H05888 and the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. SJW is supported by the postdoctoral scholarship of Tufts University.