Distinguishing Primordial Black Holes from Astrophysical Black Holes by Einstein Telescope and Cosmic Explorer

TL;DR

This paper explores how next-generation gravitational-wave detectors like Einstein Telescope and Cosmic Explorer can differentiate primordial black holes from astrophysical black holes through sub-solar mass detection and merger rate analysis.

Contribution

It proposes methods to identify primordial black holes using sub-solar mass detections and merger rate evolution, advancing black hole classification techniques.

Findings

Detectable limits for sub-solar mass PBHs in dark matter.

Forecasted event rate distributions for PBH and ABH mergers.

Potential to distinguish super-solar mass PBHs from ABHs.

Abstract

We investigate how the next generation gravitational-wave (GW) detectors, such as Einstein Telescope (ET) and Cosmic Explorer (CE), can be used to distinguish primordial black holes (PBHs) from astrophysical black holes (ABHs). Since a direct detection of sub-solar mass black holes can be taken as the smoking gun for PBHs, we figure out the detectable limits of the abundance of sub-solar mass PBHs in cold dark matter by the targeted search for sub-solar mass PBH binaries and binaries containing a sub-solar mass PBH and a super-solar mass PBH, respectively. On the other hand, according to the different redshift evolutions of merger rate for PBH binaries and ABH binaries, we forecast the detectable event rate distributions for the PBH binaries and ABH binaries by ET and CE respectively, which can serve as a method to distinguish super-solar mass PBHs from ABHs.

Equations

$$R(z)=3.9\times 10^6\times {\left(\frac{t(z)}{t_0}\right)}^{-\frac{34}{37}}{m}^{-32/37}{f}^2{\left(f^2+{\sigma }_{\mathrm{eq}}^2\right)}^{-21/74},$$

$$N_{\mathrm{obs}}=\int R(z)\frac{\mathrm{d}VT}{\mathrm{d}z}\mathrm{d}z,$$

$$\frac{\mathrm{d}VT}{\mathrm{d}z}=\frac{\mathrm{d}{V}_c}{\mathrm{d}z}\frac{T_{\mathrm{obs}}}{1+z}f(z|\xi ),$$

$$R_{90}=\frac{2.303}{VT},$$

$$VT=\int \frac{\mathrm{d}VT}{\mathrm{d}z}\mathrm{d}z.$$

$$N_{\mathrm{temp}}\propto {\left(M_{\min }{f}_{\min }\right)}^{-8/3}.$$

$$P(m)=\{\begin{array}{ll}{P}_0,& 0.2M_{\odot }\le m<1M_{\odot }\\ P_1,& 1M_{\odot }\le m<30M_{\odot }\\ P_2,& 30M_{\odot }\le m<60M_{\odot }\\ P_3,& 60M_{\odot }\le m\le 100M_{\odot }\end{array}$$

$$\int P(m)\mathrm{d}m=0.8P_0+29P_1+30P_2+40P_3=1.$$

$${\mathcal{R}}_{12}$$

$$\mathcal{R}(t|\vec{\theta })=\int {\mathcal{R}}_{12}(t|\vec{\theta })\mathrm{d}{m}_1\mathrm{d}{m}_2.$$

$${\mathcal{R}}_{12}(t_0|\vec{\theta })=Rp(m_1,m_2|\vec{\theta }),$$

$$p(\vec{d}|\vec{\theta },R)\propto R^N{e}^{-R\beta (\vec{\theta })}\prod _i^N\int \mathrm{d}\vec{\lambda }p(d_i|\vec{\lambda })p(\vec{\lambda }|\vec{\theta }),$$

$$\beta (\vec{\theta })\equiv \int \mathrm{d}\vec{\lambda }VT(\vec{\lambda })p(\vec{\lambda }|\vec{\theta }),$$

$$p(\vec{\theta },R|\vec{d})\propto p(\vec{d}|\vec{\theta },R)p(\vec{\theta },R),$$

$$p(\vec{\theta },R)\propto \frac{1}{R}.$$

$$p(\vec{\theta }|\vec{d})\propto {\left[\beta (\vec{\theta })\right]}^{-N}\prod _i^N\int \mathrm{d}\vec{\lambda }p(d_i|\vec{\lambda })p(\vec{\lambda }|\vec{\theta }),$$

$$f_{PBH0}\equiv f_{\mathrm{PBH}}{P}_0\Delta m_0,$$

$${\mathcal{R}}_{12}(z)={\int }_{t_{\text{min}}}^{t_{\max }}{R}_{\mathrm{birth}}(t(z)-t_d,m_1)\times P_d(t_d)d{t}_d,$$

$$R_{\mathrm{birth}}(t,m_{\mathrm{bh}})=\int \psi [t-\tau (m)]\varphi (m)\delta (m-g_{\mathrm{bh}}^{-1}(m_{\mathrm{bh}}))dm,$$

$$\frac{m_{\mathrm{bh}}}{m}=A{\left(\frac{m}{40M_{\odot }}\right)}^{\beta }\frac{1}{{\left(\frac{Z(z)}{0.01Z_{\odot }}\right)}^{\gamma }+1},$$

$$\psi (z)=k\frac{a\exp [b(z-z_m)]}{a-b+b\exp [a(z-z_m)]},$$

$$\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{d}z}=\int \mathrm{d}{m}_1\mathrm{d}{m}_2{\mathcal{R}}_{12}(z)\frac{\mathrm{d}VT}{\mathrm{d}z}.$$

$$N_{\mathrm{obs}}=\int \mathrm{d}z\frac{\mathrm{d}{N}_{\mathrm{obs}}}{\mathrm{d}z}.$$

$${\Omega }_{\mathrm{GW}}(\nu )=\frac{\nu }{{\rho }_c}\frac{\mathrm{d}{\rho }_{\mathrm{GW}}}{\mathrm{d}\nu },$$

$${\Omega }_{\mathrm{GW}}(\nu )=\frac{\nu }{{\rho }_c{H}_0}{\int }_0^{z_{\max }}\mathrm{d}z\frac{R(z)}{(1+z)E(z)}\frac{\mathrm{d}{E}_{\mathrm{GW}}}{\mathrm{d}{\nu }_s},$$

$$\frac{\mathrm{d}{E}_{\mathrm{GW}}}{\mathrm{d}{\nu }_s}=\frac{{\pi }^{2/3}{M}^{5/3}\eta }{3}\{\begin{array}{ll}{\nu }_s^{-1/3},& {\nu }_s<{\nu }_1,\\ \frac{{\nu }_s}{{\nu }_1}{\nu }^{-1/3},& {\nu }_1\le {\nu }_s<{\nu }_2,\\ \frac{{\nu }_s^2}{{\nu }_1{\nu }_2^{4/3}}\frac{{\nu }_4^4}{{\left(4{\left({\nu }_s-{\nu }_2\right)}^2+{\nu }_4^2\right)}^2},& {\nu }_2\le {\nu }_s<{\nu }_3,\end{array}$$

$$\rho =\sqrt{2T_{\mathrm{obs}}}{\left[\int \mathrm{d}\nu \sum _{I=1}^n\sum _{J>I}^n\frac{{\Gamma }_{IJ}^2(\nu )S_h^2(\nu )}{P_{nI}(\nu )P_{nJ}(\nu )}\right]}^{1/2},$$

$$S_h(\nu )=\frac{3H_0^2}{2{\pi }^2}\frac{{\Omega }_{\mathrm{GW}}(\nu )}{{\nu }^3}.$$

Full text

I Introduction

Ten binary black hole (BBH) mergers were detected during LIGO/Virgo O1 and O2 observing runs Abbott et al. (2016a, b, c, 2017a, 2017b, 2017c, 2018a). Understanding the origin of these BBHs is an essential scientific goal, which is still under intensive investigation (see e.g. Bird et al. (2016); Sasaki et al. (2016); Chen and Huang (2018); Clesse and García-Bellido (2018); Fishbach et al. (2017); Clesse and García-Bellido (2017); Antonini and Rasio (2016); Inayoshi et al. (2017); Ali-Haïmoud et al. (2017); Perna et al. (2019); Kavanagh et al. (2018); Rodriguez et al. (2015, 2016); Park et al. (2017); Wu (2020); Belczynski et al. (2014, 2016); Woosley (2016); Rodriguez and Loeb (2018); Choksi et al. (2018); de Mink et al. (2010); de Mink and Mandel (2016)). The fact that the component masses of these BBHs observed by gravitational waves (GWs) exhibit a much heavier mass distribution than the one inferred from X-ray observations Wiktorowicz et al. (2013); Casares and Jonker (2014); Corral-Santana et al. (2013, 2016) has triggered the interest in the community to speculate that the observed BBH mergers might come from the stellar mass primordial black holes (PBHs) Bird et al. (2016); Sasaki et al. (2016); Chen and Huang (2018); Clesse and García-Bellido (2018).

PBHs are the black holes (BHs) that form in the early universe by gravitational collapse of primordial density perturbations Hawking (1971); Carr and Hawking (1974); Khlopov (2010); Sasaki et al. (2018) and undergo quite different evolutionary histories than the astrophysical black holes (ABHs) which originate from the demise of massive stars. PBHs may contribute to a fraction of cold dark matter (CDM), and the abundance of PBHs in CDM, $f_{\mathrm{PBH}}$, has been constrained by a variety of experiments, *e.g. *extra-galactic gamma-ray Carr et al. (2010), femtolensing of gamma-ray bursts Barnacka et al. (2012), the existence of white dwarfs in our local galaxy Graham et al. (2015), Subaru/HSC microlensing Niikura et al. (2017), Kepler milli/microlensing Griest et al. (2013), OGLE microlensing Niikura et al. (2019), EROS/MACHO microlensing Tisserand et al. (2007); Calcino et al. (2018), dynamical heating of ultra-faint dwarf galaxies Brandt (2016), X-ray/radio constraints Gaggero et al. (2017), cosmic microwave background (CMB) spectrum Ali-Haïmoud and Kamionkowski (2017); Aloni et al. (2017); Horowitz (2016); Chen et al. (2016); Poulin et al. (2017), and GWs Wang et al. (2018); Abbott et al. (2018b, 2019a); Magee et al. (2018); Wang et al. (2019); Chen et al. (2019a); Yuan et al. (2019); Chen et al. (2019b).

An alternative way to explain LIGO/Virgo BBHs is through ABH models. The formation and merger of ABHs are guided by evolutionary environments. There exist three main channels in the literature. The first one is the dynamical formation channel, in which BHs are formed through the evolution of massive stars and segregated to the cluster core to pair as BBHs Rodriguez et al. (2015, 2016); Park et al. (2017). The second one is the classical isolated binary evolution channel, in which the BBHs are formed through highly non-conservative mass transfer or common envelope ejection Belczynski et al. (2014, 2016); Woosley (2016); Rodriguez and Loeb (2018); Choksi et al. (2018). The third one is the chemically homogeneous evolution, in which stars evolve almost chemically homogeneously to form BHs because of the mixing of helium produced in the center throughout the envelope de Mink et al. (2010); de Mink and Mandel (2016). Properties of BBHs, such as the spin Farr et al. (2017); Tiwari et al. (2018); Ng et al. (2018); Stevenson et al. (2017); Bogomazov et al. (2018); Lopez et al. (2018); Arca Sedda and Benacquista (2019); Farr et al. (2018), redshift Fishbach et al. (2018); Emami and Loeb (2018); Bai et al. (2018), and eccentricity distributions Samsing et al. (2014); Samsing (2018); Samsing et al. (2018); Lower et al. (2018) have been proposed to discriminate different channels of astrophysical origin BBH (AOBBH) models.

In this paper, we will explore and forecast the possibility of distinguishing PBHs from ABHs by using GW observations, especially by the third generation ground-based GW detectors like Einstein Telescope (ET) Punturo et al. (2010) and Cosmic Explorer (CE) Abbott et al. (2017d), which are expected to detect many more BBHs than current LIGO/Virgo, at an order of $\mathcal{O}(10^{5})$ events per year Regimbau et al. (2017); Vitale and Farr (2018). The rest of this paper is organized as follows. In Sec. II, we focus on the sub-solar mass ($\lesssim 1M_{\odot}$) BBHs. Because ABHs are expected to be heavier than the Chandrasekhar mass limit $\sim 1.4M_{\odot}$ Chandrasekhar and Milne (1931); Chandrasekhar (1931), direct detection of sub-solar mass BHs can be evidence of PBHs. In Sec. II.1, assuming PBHs have a monochromatic mass distribution, we estimate the detectable limits on the abundance of PBHs from the targeted search by ET and CE, respectively. In Sec. II.2, considering that PBHs have a broad mass distribution and all the BBH events detected by LIGO/Virgo originate from PBHs, we adopt a model-independent approach to constrain the abundance of PBHs with super-solar mass ($\gtrsim 1M_{\odot}$) from LIGO/Virgo events, and then explore the detectable limits on the abundance of sub-solar mass PBHs by searching for the BBHs containing a sub-solar mass PBH and a super-solar mass PBH. Sec. III is dedicated to the super-solar mass BBHs. The redshift evolution of the merger rate for primordial origin BBHs (POBBHs) and AOBBHs can be quite different, which results in different redshift distributions of the expected number of observable BBHs. We estimate and forecast the event number distributions of the PBH and ABH models for ET and CE respectively, which can serve as a complementary tool to distinguish PBHs from ABHs. Finally, we summarize and discuss our results in Sec. IV

II Distinguish PBHs from ABHs by Sub-solar Mass BHs

Direct detection of sub-solar mass BHs can be taken as a smoking gun for PBHs. Nonetheless, the event rate relies both on the merger rate of POBBHs and the sensitivity of GW detectors. In the following two subsections, we will explore the abilities to detect sub-solar mass BHs for different GW detectors by considering the cases when PBHs have a monochromatic and a general mass functions, respectively.

II.1 Monochromatic mass function

In this subsection, we assume all the PBHs have the same mass and estimate the detectable limits of $f_{\mathrm{PBH}}$ by the targeted search of POBBHs.

The redshift $z$ evolution of the local merger rate $R(z)$ in the comoving frame for the monochromatic mass function, by accounting the angular momentum exerted both by all PBHs and the background inhomogeneity, is given by Ali-Haïmoud et al. (2017); Chen and Huang (2018)

[TABLE]

where $mM_{\odot}$ is the component mass of BBHs measured in the source frame, and $\sigma_{\mathrm{eq}}$ is the variance of density perturbations of the rest DM on scale of order $\mathcal{O}(10^{0}\sim 10^{3})M_{\odot}$ at radiation-matter equality. Following Ali-Haïmoud et al. (2017); Chen and Huang (2018), we choose $\sigma_{\mathrm{eq}}\approx 0.005$. Here $f_{\mathrm{PBH}}\equiv\Omega_{\mathrm{PBH}}/\Omega_{\mathrm{CDM}}$ is the energy density fraction of PBHs in CDM and is related to the total abundance of PBHs in non-relativistic matter, $f$, by $f_{\mathrm{PBH}}\approx f/0.85$. Besides, $t(z)$ is the cosmic time at the redshift of $z$ and $t_{0}\equiv t(0)$ is the age of our universe. Note that throughout this paper, we adopt the units in which the Newtonian constant $G$ and the speed of light in vacuum $c$ are set to unity. Here we ignore the effects of early formed PBH clusters which may alter the merger rate of PBH binaries and present a rich phenomenology Vaskonen and Veermäe (2020); Trashorras et al. (2020); Jedamzik (2020). We hope to come back to this issue in the future.

The expected number of detection, $N_{\mathrm{obs}}$, then follows Chen and Huang (2018); Kavanagh et al. (2018)

[TABLE]

where $\mathrm{d}VT/\mathrm{d}z$ is the spacetime sensitivity of a GW detector as a function of redshift and accounts for the selection effects of that detector. Generally, $\mathrm{d}VT/\mathrm{d}z$ depends on the properties $\xi$ (*e.g. *masses and spin) of a binary and is defined as Abbott et al. (2016d, e)

[TABLE]

where $V_{c}$ is the comoving volume Hogg (1999), $T_{\mathrm{obs}}$ is the observing time, and the denominator $1+z$ accounts for the converting of cosmic time from source frame to detector frame due to the cosmic redshift. Here $0<f(z|\xi)<1$ is the probability of detecting a BBH with the given parameters $\xi$ at redshift $z$ O’Shaughnessy et al. (2010). The $90\%$ confidential upper limit on the binary merger rate can then be obtained by using the loudest event statistic formalism Biswas et al. (2009),

[TABLE]

where

[TABLE]

We adopt the semi-analytical approximation from Abbott et al. (2016d, e) to calculate $VT$ by neglecting the effect of spins for BHs and using the IMRPhenomPv2 waveform to simulate the BBH templates. Furthermore, we set a single-detector signal-to-noise ratio (SNR) threshold $\rho_{\mathrm{th}}=8$ as a criterion of detection, which roughly corresponds to a network threshold of $12$. Fig. 1 shows the estimated $90\%$ upper limits on the merger rate of equal-mass BBHs, $R_{90}$, which are estimated from Eq. (4).

The detectable limit of $f_{\mathrm{PBH}}$ by the targeted search from LIGO O1 & O2, LIGO Design, ET and CE are shown in Fig. 2. The effective observing time of LIGO O1 & O2 is set to their total running time of $165.6$ days Abbott et al. (2016c, 2017e). Meanwhile LIGO Design, ET and CE are supposed to operate at $1$ year with full duty. The upper limit for $f_{\mathrm{PBH}}$ of sub-solar mass PBHs in the mass range $\left[0.2,1\right]M_{\odot}$ has been reported in Abbott et al. (2018b); Magee et al. (2018) by the null targeted search result of BBHs in that mass range. We extrapolate the results of Abbott et al. (2018b); Magee et al. (2018) in several aspects. Firstly, we adopt the merger rate presented in Ali-Haïmoud et al. (2017), which takes a more careful examination of the dynamical evolution of the binary systems than the one used by Abbott et al. (2018b); Magee et al. (2018) from the results of Sasaki et al. (2016). Secondly, we also estimate the detectable limits of $f_{\mathrm{PBH}}$ by the proposed third generation GW detectors like CE and ET. Lastly, we do not limit to the mass range of $\left[0.2,1\right]M_{\odot}$ but extend to the range constrained by the detectors automatically. In particular, since whether the super-solar mass BBHs observed so far are of POBBHs or not is still under debate, we use the dashed lines to illustrate the detectable limits of $f_{\mathrm{PBH}}$ for the super-solar mass PBHs. Furthermore, as the masses go down below $0.2M_{\odot}$, the search difficulty arises because the number of templates required in the template bank, $N_{\mathrm{temp}}$, scales both as the minimum mass ${M_{\mathrm{min}}}$ and the starting frequency ${f_{\mathrm{min}}}$ Magee et al. (2018)

[TABLE]

The dramatic increment of computational resource limits the current GW search pipeline that it cannot deal with the BBHs with component masses far below $0.2M_{\odot}$ efficiently. However, in the future, this issue may be addressed by the improvement of search algorithms or computational technologies. Note also that the upper limits of $f_{\mathrm{PBH}}$ in mass range of $[10,300]M_{\odot}$ has already been derived in Ali-Haïmoud et al. (2017); Kavanagh et al. (2018) by using LIGO/Virgo O1 data. We updated their constraints by using both O1 and O2 data. And therefore our upper limits are more stringent than those in Ali-Haïmoud et al. (2017); Kavanagh et al. (2018).

Besides, we update the constraint on $f_{\mathrm{PBH}}$ from the null search of stochastic gravitational-wave background (SGWB) from LIGO O1 in Wang et al. (2018) to both LIGO O1 and O2 runs Abbott et al. (2017f); Renzini and Contaldi (2019); Abbott et al. (2019b, c). The method adopted in this paper is described in the Appendix. Usually one may expect that the null detection in the targeted search would give a tighter constraint on the abundance of PBHs than that from SGWBs. However, the GW signals from light BBHs are so weak to be resolved individually by GW observations, and hence the null detection of SGWBs provides a more stringent constraint on the abundance of light PBHs because those weak signals can superpose to form a detectable SGWB. See a cross of the red and black curves around $0.1M_{\odot}$ in Fig. 2. Actually, it is also true for other detections, such as LIGO design, ET and CE.

II.2 General mass function – a model independent approach

A null search result of the sub-solar mass PBHs in LIGO’s O1 and O2 data has been reported in Abbott et al. (2018b); Magee et al. (2018); Abbott et al. (2019a). However, the search was only targeted at the POBBHs with component masses between $0.2M_{\odot}$ and $1M_{\odot}$. In this subsection, we propose to search for the POBBHs with a component of sub-solar mass and the other of super-solar mass, which are expected to emit stronger GWs.

Here we extend the discussion in the former subsection to the case in which PBHs exhibit a broad mass distribution, and assume that all the ten BBHs observed by LIGO/Virgo so far are POBBHs. Contrary to the previous works Chen et al. (2019b); Raidal et al. (2017); Kocsis et al. (2018); Raidal et al. (2018) by choosing some specific mass functions, *e.g. *a power-law or lognormal ones, which are pertinent to some specific formation models of PBHs, we take a model-independent approach by binning the mass function $P(m)$ from $0.2M_{\odot}$ to $100M_{\odot}$ as

[TABLE]

where $P_{i}=\{P_{0},P_{1},P_{2},P_{3}\}$ are four constants satisfying the normalization condition

[TABLE]

Here, only three out of the four $P_{i}$s are independent and we will choose $\vec{\theta}=\{P_{1},P_{2},P_{3}\}$ as free parameters. $\vec{\theta}$ will be fitted by the ten BBHs from LIGO’s O1 and O2 runs. In this subsection, we are only interested in the PBHs with masses in the range of $\left[{m_{\mathrm{min}}},{m_{\mathrm{max}}}\right]=\left[0.2M_{\odot},100M_{\odot}\right]$. Here ${m_{\mathrm{min}}}=0.2M_{\odot}$ corresponds to the lower mass bound of LIGO’s O1 search for the sub-solar mass ultracompact binaries in Abbott et al. (2018b).

The merger rate for POBBHs with a general mass function is given in Chen and Huang (2018).111In Raidal et al. (2017); Kocsis et al. (2018); Liu et al. (2019, 2018); Raidal et al. (2018), the authors ignored the relative distribution of PBHs, and cannot guarantee that PBH binaries be formed from the closest neighboring PBHs. The time-dependent comoving merger rate density for a general normalized mass function, $P(m|\vec{\theta})$, takes the form of

[TABLE]

in units of $\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$, where the component masses $m_{1}$ and $m_{2}$ are in units of $M_{\odot}$. The time (or redshift)-dependent merger rate can be obtained by integrating over the component masses

[TABLE]

The local merger rate density distribution then reads Chen et al. (2019b)

[TABLE]

where the local merger rate $R\equiv\mathcal{R}(t_{0}|\vec{\theta})$ is chosen such that the population distribution of BBH mergers, $p(m_{1},m_{2}|\vec{\theta})$, is normalized.

To extract the population parameters $\{\vec{\theta},R\}$ from the merger events observed by LIGO/Virgo, it is necessary to perform the hierarchical Bayesian inference on the BBHs’ mass distribution Abbott et al. (2016d, e, c); Wysocki et al. (2018); Fishbach et al. (2018); Mandel et al. (2018); Thrane and Talbot (2018). If we have the data of $N$ BBH detections, $\vec{d}=(d_{1},\dots,d_{N})$, then the likelihood for an inhomogeneous Poisson process is Wysocki et al. (2018); Fishbach et al. (2018); Mandel et al. (2018); Thrane and Talbot (2018)

[TABLE]

where $\vec{\lambda}\equiv\{m_{1},m_{2}\}$. The likelihood of an individual event $p(d_{i}|\vec{\lambda})$ is proportional to the posterior of that event $p(\vec{\lambda}|d_{i})$, as the standard priors on masses for each event in LIGO/Virgo analysis are taken to be uniform. We will use the publicly available posterior samples of ten BBHs Abbott et al. (2016c, 2018a) from LIGO/Virgo observations to evaluate the integral in Eq. (12). Meanwhile, $\beta(\vec{\theta})$ is defined as

[TABLE]

in which $VT(\vec{\lambda})$ is given by Eq. (5). The posterior probability distribution $p(\vec{\theta},R|\vec{d})$ can be directly estimated by

[TABLE]

where as usual the prior distribution $p(\vec{\theta},R)$ (Abbott et al., 2016d, 2017a) is chosen to be uniform for $\vec{\theta}$ parameters and log-uniform for local merger rate $R$,

[TABLE]

Integrating over $R$ in Eq. (14), it is then easy to obtain the marginalized posterior

[TABLE]

which has been widely used in previous population inferences (Abbott et al., 2016d, 2017a, c, e; Fishbach and Holz, 2017; Chen et al., 2019b).

Using ten BBH events from LIGO’s O1 and O2 runs, we find the median value and $90\%$ equal-tailed credible intervals for the parameters $\{\vec{\theta},R\}$ to be $P_{1}=2.1^{+0.7}_{-0.8}\times 10^{-2}$, $P_{2}=5.4^{+4.7}_{-3.1}\times 10^{-3}$, $P_{3}=5.1^{+15.2}_{-4.6}\times 10^{-4}$, and $R=308^{+193}_{-135}\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$, from which we also infer the fraction of PBHs in CDM to be $f_{\mathrm{PBH}}=3.3\,^{+2.3}_{-1.8}\times 10^{-3}$. Such an abundance of PBHs is consistent with previous estimations that $10^{-3}\lesssim f_{\mathrm{PBH}}\lesssim 10^{-2}$, confirming that the dominant fraction of CDM should not originate from PBHs in the mass range of $[0.2,100]M_{\odot}$ (Sasaki et al., 2016; Ali-Haïmoud et al., 2017; Raidal et al., 2017; Kocsis et al., 2018; Chen and Huang, 2018; Chen et al., 2019b). From now on we will investigate the possibility of detecting sub-solar mass BBHs. Here we denote the abundance of PBHs in the mass range of $[0.2,1]M_{\odot}$ as

[TABLE]

where $\Delta m_{0}=(1-0.2)M_{\odot}=0.8M_{\odot}$. As a consequence of the above analysis, it is then straightforward to infer the upper bound of $f_{\mathrm{PBH0}}$ to be $f_{\mathrm{PBH0}}\leq 1.8\times 10^{-3}$ by LIGO’s O1 and O2 runs. In the future, if third-generation ground-based GW detectors are in operation, the detection ability will be greatly enhanced and we have more chance to detect the sub-solar mass BBHs if they do exist. In addition to search for the BBHs with two sub-solar mass components as LIGO/Virgo have done, we also propose to search for the BBHs with one sub-solar mass component (with mass lying in $[0.2,1]M_{\odot}$) and another super-solar mass one (with mass lying in $[1,100]M_{\odot}$). Using the loudest event statistic formalism [see Eq. (4)] and the values of $\vec{\theta}$ inferred from LIGO’s O1 and O2 runs, $f_{\mathrm{PBH0}}$ can be constrained to an unprecedented level. Assuming no such BBHs will be detected, ET implies $f_{\mathrm{PBH0}}\leq 4.1\times 10^{-7}$ while CE implies $f_{\mathrm{PBH0}}\leq 4.5\times 10^{-8}$. The results of the constraints on $f_{\mathrm{PBH}}$ (and $f_{\mathrm{PBH0}}$) when PBHs have a broad mass distribution are shown in Fig. 3. Note that the red line in Fig. 3 shows the upper limit of $f_{\mathrm{PBH0}}\leq 1.6\times 10^{-2}$ from the null targeted search of BBHs with two sub-solar mass components, assuming that PBHs take a flat distribution in the mass range $[0.2,1]M_{\odot}$.

It is worthy to note that the targeted search of BBHs with a sub-solar mass and a super-solar mass components will improve the detectable limit of $f_{\mathrm{PBH}}$ by an order of $\mathcal{O}(10^{2}\sim 10^{3})$, comparing to the targeted search for BBHs with two sub-solar mass BHs as shown in Fig. 2.

III Distinguish PBHs from ABHs by Super-Solar Mass BHs

Besides the method of using sub-solar mass BBHs to distinguish PBHs from ABHs, there is another way by exploring the redshift evolution of the event rate of super-solar mass BBHs. In Chen and Huang (2018), we found that the merger rate of POBBHs increases as a function of redshift $z$, namely $\mathcal{R}(z)\propto t(z)^{-34/37}$, which is independent of both the abundance and mass function of PBHs. However, the merger rate predicted by the AOBBHs will firstly increase with $z$, then peaks at some low redshift, and lastly rapidly decreases with $z$. Fig. 4 shows the merger rate of POBBHs and AOBBHs as a function of redshift $z$. The difference between the merger rate of these two models increases at the higher redshift. Currently, LIGO can only observe BBHs at low redshifts (with $z<1$), but future GW detectors such as CE and ET will be able to probe much higher redshifts (with $z\geq 10$). In this section, we will demonstrate how well the third-generation ground-based detectors like CE and ET can be used to distinguish PBHs from ABHs according to the quite different event rate distributions at high redshifts.

To comply with the previous studies (Abbott et al., 2017a, 2018c), we restrict the component masses of BBHs to the range ${m_{\mathrm{min}}}\leq m_{2}\leq m_{1}$ and $m_{1}+m_{2}\leq{m_{\mathrm{max}}}$, with ${m_{\mathrm{min}}}=5M_{\odot}$ and ${m_{\mathrm{max}}}=100M_{\odot}$. To calculate the observable events rate, we first need to know the merger rate distribution. For PBHs, we assume they have a broad mass distribution of Eq. (7) and adopt the best-fits from Sec. II.2. For ABHs, the merger rate is a convolution of the birthrate of ABHs $R_{\mathrm{birth}}(z,m)$ with the distribution of the time delays $P_{d}(t_{d})$ between the formation and merger of AOBBHs (see e.g. Dvorkin et al. (2016))

[TABLE]

where $t_{d}$ is the time delay, and $t(z)$ is the age of the Universe at merger. The birthrate $R_{\mathrm{birth}}$ can be estimated by (Dvorkin et al., 2016)

[TABLE]

where $m_{\mathrm{bh}}$ is the mass of the remnant BH, $\tau(m)$ is the lifetime of a progenitor star of mass $m$ and $\phi(m)\propto m^{-2.35}$ is the initial mass function (IMF) Salpeter (1955). We consider the WWp model (Woosley and Weaver, 1995) of BH formation in which $m$ and $m_{\mathrm{bh}}$ are related by

[TABLE]

where $Z(z)$ is the metallicity and an explicit functional form can be found in Belczynski et al. (2016). The parameter values in Eq. (20) are $A=0.3$, $\beta=0.8$ and $\gamma=0.2$ (Dvorkin et al., 2016). Solving Eq. (20) yields the functional form $m=g_{\mathrm{bh}}^{-1}(m_{\mathrm{bh}})$. The star formation rate (SFR) $\psi(t)$ in Eq. (19) is given by (Nagamine et al., 2004)

[TABLE]

where z is the redshift and the parameter values of $\{k,a,b,z_{m}\}$ are dependent on the specific AOBBH formation scenarios. In this paper, we will consider $4$ different AOBBH models.

The first one is the Fiducial model in which the SFR is a fit to the observations of the luminous galaxies and the fit parameters are $k=0.178\,M_{\odot}\,\mathrm{yr}^{-1}\,\mathrm{Mpc}^{-3}$, $z_{m}=2.00$, $a=2.37$, $b=1.80$ as given by Vangioni et al. (2015), and the time delay distribution has the form of $P_{d}\propto t_{d}^{-1}$ with $t_{\mathrm{min}}<t_{d}<t_{\mathrm{max}}$, where the minimum delay time of a binary system to evolve until coalescence is set to be $t_{\mathrm{min}}=50$ Myr, and the maximum delay time $t_{\mathrm{max}}$ is set to the Hubble time (Abbott et al., 2016f). This model corresponds to the classical isolated binary evolution of AOBBH model Dvorkin et al. (2016).

The second one is the GRB-based model which differs from the fiducial model by using the SFR calibrated from the gamma ray bursts (GRB) at high redshift, and the fit parameters are given by $k=0.146\,M_{\odot}\,\mathrm{yr}^{-1}\,\mathrm{Mpc}^{-3}$, $z_{m}=1.72$, $a=2.80$, $b=2.46$ as given by Vangioni et al. (2015).

The third one is the LongDelay model which is identical to the Fiducial model but assumes a significantly longer minimum time delay $t_{\mathrm{min}}=5$ Gyr. The model is potentially consistent with the chemically homogeneous evolution of rapidly rotating massive stars in very tight binaries Mandel and de Mink (2016).

The last one is the FlatDelay model assuming a flat time delay distribution with $t_{\mathrm{min}}=50$ Myr and $t_{\mathrm{max}}=1$ Gyr. This model is inspired by the supposition that dynamical formation of the most massive binaries is likely to happen fairly early in the history of the host environment Abbott et al. (2016f).

Fig. 4 compares the redshift distribution of the normalized merger rate, $R(z)/R(0)$, for the POBBH and different AOBBH models, respectively. It is shown that the merger rate of LongDelay model decreases quickly as redshift increases and GRB-based model has a relative high merger rate at high redshift comparing to the Fiducial and the FlatDelay models. Moreover, the merger rate from all the AOBBH models decreases at high redshift while the merger rate from the POBBH model increases at high redshift. The redshift dependent observable events number density of a GW detector can be calculated by

[TABLE]

Integrating over the redshift $z$ results in the total number of observable events, $N_{\mathrm{obs}}$,

[TABLE]

Note that Eq. (2) is a special case of Eq. (23) when $m_{1}=m_{2}=m$. Fig. 5 shows the expected number of observable BBHs, $\mathrm{d}N_{\mathrm{obs}}/\mathrm{d}z$, as a function of redshift for CE and ET, respectively. Fig. 6 shows the total number of observable BBHs, $N_{\mathrm{obs}}$, as a function of redshift for CE and ET, respectively. The third-generation GW detectors like CE and ET are expected to detect $\mathcal{O}(10^{5})$ BBH mergers each year and dig much deeper at redshit, and the fact that the redshift distribution of $\mathrm{d}N_{\mathrm{obs}}/\mathrm{d}z$ for POBBHs and AOBBHs are quite different from each other can be taken as a complementary tool to distinguish these two formation models of BBHs. In particular, the contribution to the expected number of observable BBHs from high redshift (with $z>5$) for all the $4$ AOBBH models can be negligible, and hence the total number of observable events, $N_{\mathrm{obs}}$, approaches a constant number at $z>5$. However, for the POBBHs, the contribution from the higher redshift cannot be ignored, and therefore $N_{\mathrm{obs}}$ still increases when $z>5$. Therefore, an “excess” of total number of observable events after redshift $z=5$ could possibly point to a population of POBBHs.

IV Summary and Discussion

Even though several BBH coalescences have been detected by LIGO/Virgo, the origin of these BHs is still unknown. In this paper, we explore how well the next generation detectors, such as ET and CT, can be used to distinguish PBHs from ABHs.

Firstly, we investigate the possibility of direct detection of sub-solar mass BBHs, hence validating the existence of PBHs. For PBHs with a monochromatic mass function, we estimate and forecast the detectable limit of $f_{\mathrm{PBH}}$ from the targeted search of BBHs by LIGO, ET, and CE, respectively. Furthermore, in order to get better sensitivity, we propose to search for the BBHs containing a sub-solar mass PBH and a super-solar mass PBH. We predict that the abundance of PBHs in the mass range of $[0.2,1]M_{\odot}$ can be constrained to an order of $\mathcal{O}(10^{-7})$ and $\mathcal{O}(10^{-8})$ if no such BBHs are to be detected by ET and CE, respectively.

Secondly, we explore the possibility of utilizing the redshift evolution of the merger rate of super-solar mass BBHs to distinguish PBHs from ABHs. We estimate and forecast the redshift distribution of the expected number of observable BBHs for the PBH and ABH models, respectively. When the third generation ground-based GW detectors like CE and ET are in operation, it is expected to detect $\mathcal{O}(10^{5})$ BBH mergers each year and reach much deeper redshift ($z\gtrsim 10$), and the redshift distribution of detectable BBH events can serve as an alternative means to distinguish PBHs from ABHs.

Throughout this paper, we assume that all the LIGO/Virgo BBHs originate from the same formation channel. However, this assumption can be too oversimplified because the observed BBHs might be a mixing of POBBHs and AOBBHs. To identify each BBH as a POBBH or AOBBH will be quite difficult in this scenario. In addition to the mass and redshift distribution of BBHs, other information, *e.g. *spin distribution, will also be invaluable in order to find out the progenitors of the LIGO/Virgo BBHs. For instance, it is expected that PBHs formed in the early universe have negligible spins Chiba and Yokoyama (2017); Mirbabayi et al. (2019); De Luca et al. (2019), while the ABHs which originated from the Population III star binaries, favor a relatively high spin distribution Kinugawa et al. (2016).

Acknowledgements.

We would like to thank Lu Chen, Yun Fang, Fan Huang, Jun Li, Lang Liu, Shi Pi, You Wu, Yu Sang, Sai Wang, Hao Wei, Chen Yuan, and Xue Zhang for useful conversations. We acknowledge the use of HPC Cluster of ITP-CAS. This work is supported by grants from NSFC (grant No. 11690021, 11975019, 11947302, 11991053), the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDB23000000, XDA15020701), and Key Research Program of Frontier Sciences, CAS, Grant NO. ZDBS-LY-7009. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center Vallisneri et al. (2015) (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.

Appendix A Constrain $f_{\mathrm{PBH}}$ by SGWB

Another way to constrain $f_{\mathrm{PBH}}$ is through the SGWB Wang et al. (2018) which is a superposition of the energy spectra emitted by the BBHs that are unlikely to be resolved individually. This method differs from the one by targeted search mainly due to it also utilizes the redshift dependence of the merger rate (see Eq. (25) below). Notice that the targeted search is more sensitive to the local merger rate. The energy-density spectrum of a SGWB is characterized by the dimensionless quantity Allen and Romano (1999)

[TABLE]

where $\rho_{c}=3H_{0}^{2}/(8\pi)$ is the critical energy density of our universe, $d\rho_{\mathrm{GW}}$ is the energy density in the frequency interval $\left[\nu,\,\nu+d\nu\right]$, and we take the value of Hubble constant $H_{0}=67.74\,\mathrm{km}\sec^{-1}\mathrm{Mpc}^{-1}$ from Planck (Ade et al., 2016). For the binary mergers, the magnitude of a SGWB can be calculated via Phinney (2001); Regimbau and Mandic (2008); Zhu et al. (2011, 2013)

[TABLE]

where $\nu_{s}=(1+z)\,\nu$ is the frequency in source-frame, $E(z)=\sqrt{\Omega_{r}\left(1+z\right)^{4}+\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}}$ accounts for the evolution of our universe, and the factor $(1+z)$ in the denominator of Eq. (25) converts the merger rate, $R(z)$, from source frame to detector frame. We adopt the best-fit results from Planck (Ade et al., 2016) that $\Omega_{r}=9.15\times 10^{-5}$, $\Omega_{m}=0.3089$, and $\Omega_{\Lambda}=1-\Omega_{m}-\Omega_{r}$. The cutoff redshift is chosen to be $z_{\mathrm{max}}=\nu_{3}/\nu-1$ Wang et al. (2018), in which $\nu_{3}$ is given by Eq. (26) below. Furthermore, the energy spectrum $\mathrm{d}E_{\mathrm{GW}}/\mathrm{d}\nu_{s}$, emitted by an individual BBH with equal component masses $m_{1}=m_{2}=m$, is approximated by Cutler et al. (1993); Chernoff and Finn (1993); Zhu et al. (2011)

[TABLE]

where $M=m_{1}+m_{2}=2\,m$ is the total mass of the BBH, $\eta=m_{1}m_{2}/M^{2}=1/4$, and $\nu_{i}=\left(a_{i}\eta^{2}+b_{i}\eta+c_{i}\right)/\left(\pi M\right)$ with $i=\{1,2,3\}$. The coefficients $a_{i}$, $b_{i}$ and $c_{i}$ are presented in Table I of Ajith et al. (2008). Again, all masses are measured in the source frame and units of $M_{\odot}$.

For a network of $n$ individual detectors, the SNR $\rho$ for measuring the SGWB with an observation time $T_{\mathrm{obs}}$ is given by (Allen and Romano, 1999; Thrane and Romano, 2013)

[TABLE]

where $P_{nI}(\nu)$ is the auto power spectral density for the noise in detector $I$ and $\Gamma_{IJ}(\nu)$ is the overlap reduction function Christensen (1992); Flanagan (1993). In Eq. (27), $S_{h}$ is the strain power density spectrum of a SGWB, which is related to $\Omega_{\mathrm{GW}}$ through Thrane and Romano (2013)

[TABLE]

Here we set $\rho=1$, which corresponds to $1\sigma$ confidential level, as the criterion for the detection of SGWBs.