Can we distinguish astrophysical from primordial black holes via the stochastic gravitational wave background?

How are black holes forming in the Universe? There are two channels to form black holes, namely the astrophysical channel from the deaths of stars and the formation of black holes in the early Universe from primordial fluctuations (Zel’dovich & Novikov, 1967; Hawking, 1971; Carr, 1975; Khlopov & Polnarev, 1980; Khlopov et al., 1985; Carr, 2005; Sasaki et al., 2018). The discovery of primordial black holes (PBHs) would change our fundamental understanding of the Universe and address one of the long-standing questions in the standard model cosmology about the nature of dark matter (Bird et al., 2016; Carr et al., 2016). Even though PBHs are not candidates for the bulk of the dark matter, at least in the solar mass range, there remains an interesting question about whether a population of PBHs exists in the Universe that contributes to the observed gravitational wave events. One of the primary challenges is to decide whether one may be able to distinguish between astrophysical black holes (ABHs) forming from stellar deaths and PBHs arising from an early universe origin.

The discovery of gravitational waves (Abbott et al., 2016b, d) has opened a new window for detection of PBHs by exploring the properties of the GW merger rates and associated source populations, such as the masses of individual sources and the spins of the GW sources (Bird et al., 2016; Clesse & García-Bellido, 2017; Sasaki et al., 2016; Gow et al., 2020; Jedamzik, 2020, 2021; De Luca et al., 2020b). Several recent studies have considered the properties of the mass distributions of GW sources from individual events in order to distinguish between ABHs and PBHs (Hall et al., 2020; Wong et al., 2021; De Luca et al., 2021; Hütsi et al., 2021; Franciolini et al., 2021). However such studies are vulnerable to freedom in higher generation merging models that enable the upper mass gap, a powerful indicator of first-generation ABHs, to be bridged (Rodriguez et al., 2019; Kimball et al., 2020; Hamers et al., 2021).

Perhaps the ultimate signature that can potentially distinguish between different formation channels of black holes is through the evolution of the merger rate with redshift (Raidal et al., 2017; Raidal et al., 2019; Vaskonen & Veermäe, 2020; Atal et al., 2020; De Luca et al., 2020a). GW sources produced from the deaths of stars form after the formation of the first stars, whereas PBHs exist in large numbers at a very high redshift before any stars formed. This is one of the key differences that can be used to distinguish between a population of ABHs and PBHs. For the ABHs, even though one would expect the merger rate of black holes to be related to the star formation rate of the Universe, one of the major sources of uncertainty in the merger rate is due to the time delay between the formation and merger of GW sources (Dominik et al., 2012, 2015; Lamberts et al., 2016; Dvorkin et al., 2016; Eldridge et al., 2019; Vitale et al., 2019; Safarzadeh et al., 2020; Santoliquido et al., 2021). For GW sources with small-time delays ($t^{eff}_{d}<100$ Myr), the peak of the mergers of GW sources is around the peak of the star formation rate ($z\approx 2$), whereas for the scenarios with larger time delays, the peak of the merger of the ABHs can be shifted towards lower redshifts. However for PBHs, the merger rate is always an increasing function of redshift. The number of mergers at high redshift for GW sources of primordial origin is always going to surpass the ABH merger rate (Ali-Haïmoud et al., 2017; Raidal et al., 2017; Raidal et al., 2019; Vaskonen & Veermäe, 2020; Atal et al., 2020; De Luca et al., 2020a). The formation of PBHs can also give additional sources of stochastic GW background (Kohri & Terada, 2018; Espinosa et al., 2018; Wang et al., 2019).

This distinction may be largely academic since GW sources at high redshifts cannot be detected as individual events (apart from rare lensed events (Nakamura, 1998; Wang et al., 1996; Dai et al., 2017; Broadhurst et al., 2018, 2019; Oguri, 2019; Mukherjee et al., 2021)). However the high merger rate of the GW sources at high redshift leads to a stochastic gravitational wave background due to the contribution from unresolved sources (Allen, 1996; Phinney, 2001; Regimbau & Chauvineau, 2007; Wu et al., 2012; Romano & Cornish, 2017; Abbott et al., 2016c, 2018b, 2019, 2021). We show that PBH mergers contribute to a potentially detectable stochastic GW background signal (Wang et al., 2018; Mandic et al., 2016).

Measurement of the stochastic GW background (even in the absence of a detection) can probe the high redshift merger rate and its evolution with redshift (Mukherjee & Silk, 2019; Boco et al., 2019; Callister et al., 2020). Using data from the third observational run, LIGO/Virgo has estimated the stochastic GW background power spectrum (Abbott et al., 2021). Though the measurement has not detected the stochastic GW background, it has provided an upper bound on the stochastic GW background power spectrum(Abbott et al., 2021).

Here we construct a hybrid merger rate model of ABHs and PBHs, by taking into account the time delay between formation and mergers of the astrophysical sources and incorporating a general redshift dependence of the PBH merger rate. Our hybrid model is driven by the use of the Madau-Dickinson star formation rate history (SFR) to derive the ABH merger rate and a power-law model for the PBH merger rate. The free parameters of this model are the time delay parameter, the local merger rate, the index of the power-law model of PBH merger rate, the characteristic mass-scale of PBHs, and the fraction of PBHs over ABHs. This five-parameter model makes it possible to perform a relatively fast MCMC search of GW sources to jointly probe the parameter space of ABHs and PBHs. In this paper, we show the current constraints on the parameters of this hybrid model of the merger rate from the LIGO/Virgo third observing run of the stochastic GW background (Abbott et al., 2021) and also using the bounds on the local merger rate from the individual events of the first half of the third observation run of O3a (Abbott et al., 2020a, b). We also provide forecasts for the O5 observation run of the LIGO/Virgo detectors in its design sensitivity (Aasi et al., 2015; Acernese et al., 2015) and for the enhanced A+ sensitivity(Abbott et al., 2018a; Barsotti et al., 2020). We do not consider the possible contributions from population-II/population-III astrophysical sources at high redshift, as there is no detection of the stochastic GW background from O3 observations. The method proposed here can easily be extended to include the contribution from such sources (Inayoshi et al., 2021).

We consider a parametric hybrid model of GW merger rates as a function of redshift to search for ABHs and PBHs. The model is composed of these two components with corresponding probability distributions of masses $P_{ABH}(m_{i})$ and $P_{PBH}(m_{i})$ written as

where $\mathcal{N}$ is the normalization factor such that the local merger rate integrated over the mass distribution agrees with the value inferred from the individual detected events. $\mathcal{R}_{ABH}(z_{m},M)$ is the merger rate for the ABHs that are expected to follow the cosmic star formation rate history. For the low redshift universe, we assume that the Madau-Dickinson fitting form (Madau & Dickinson, 2014) is a proxy for the star formation rate history. We write the model of the ABHs as

where $R_{SFR}(z)$ is the star-formation rate motivated by the Madau-Dickinson relation (Madau & Dickinson, 2014)

In Eq. 2, the time-delay distribution between the formation and mergers of the ABHs is denoted by $P(t^{eff}_{d},M)$, where $t^{eff}_{d}$ denotes the time-delay between formation and mergers. The redshift evolution of the ABH merger rate denoted in Eq. (2) is solely decided by the time-delay model. Depending on the probability distribution of the time delay, the peak of the merger rate, the slope of the merger rate at low redshift, and the slope of the merger rate at high redshift are specified. The only degree of freedom in this model is considered to be the adopted time-delay distribution. Currently, we have very little idea of the value of the time-delay parameter. Stellar population synthesis models suggest that the mergers of the different kinds of GW sources can be delayed by as much as a few hundreds of Myr up to about the age of the Universe (Dominik et al., 2012, 2015; Lamberts et al., 2016; Eldridge et al., 2019; Santoliquido et al., 2021).

The second term in Eq. (1) captures the redshift evolution of the PBHs. The merger rate of PBHs can be modeled depending on whether the binaries are dominated by Poisson statistics or whether there is clustering. In the presence of strong clustering, the current merger rate is exponentially decreased, even for a large fraction of PBHs as dark matter. The redshift evolution of the merger rate is one of the key signatures for distinguishing between the clustered and Poissonian scenarios. The general behavior of the PBH merger rate is that there is an increase in the merger rate with increasing redshift, and so we model the PBH merger rate as

where $\alpha$ is a positive index in the power-law model and $R_{PBH}(0,m_{1},m_{2})$ denotes the local merger rate of the GW sources of primordial origin. We consider $\alpha$ as a free parameter to search for PBHs in a model-independent way from the stochastic GW background. However, for most of the known scenarios of black formation, the value of $\alpha\sim 1.3$ for the Poisson distribution (Raidal et al., 2017; Sasaki et al., 2018; Raidal et al., 2019). For the forecast studies in the latter part of the paper, we consider the fiducial value of $\alpha=1.3$. The merger rate of the GW sources is degenerate between the clustering signal and the PBH fraction (Raidal et al., 2017; Raidal et al., 2019; Young & Byrnes, 2020; Vaskonen & Veermäe, 2020; Atal et al., 2020; De Luca et al., 2020a). If the spatial clustering $\xi_{PBH}>>1$ is very large, then the local merger rate can be exponentially suppressed, even if the fraction of PBH in dark matter $f_{PBH}=1$. We can express the local merger rate for the extremely clustered scenario as

As a result, the relatively low observed merger rate of GW sources does not necessarily imply that $f_{PBH}<10^{-2}$ (Raidal et al., 2017; Vaskonen & Veermäe, 2020; Atal et al., 2020). However, the merger rate at high redshift is going to be large (Raidal et al., 2017; Atal et al., 2020) making it a key signature of PBH as dark matter. The exact dependence depends on the model for generating the PBHs (Raidal et al., 2017; Vaskonen & Veermäe, 2020; Atal et al., 2020). For the Poisson distribution of PBHs (without clustering), the merger rate of the PBHs can be written as (Raidal et al., 2019; Clesse & Garcia-Bellido, 2020)

where $\eta\equiv m_{1}m_{2}/(m_{1}+m_{2})^{2}$ is the symmetric mass ratio, $M=m_{1}+m_{2}$ is the total mass of the binaries, and $f_{\text{sup}}$ is the suppression factor which depends on the effect from the surrounding matter distribution and also on the effects from other PBHs. One can explore both clustering and Poisson scenarios to explore the PBH population by using the evolution of the merger rate and its dependence on the population of black holes.

The form of the PBH merger rate as a function of redshift is fairly model-independent and can probe different populations of PBH sources. The form differs from the ABH source population, governed by the Madau-Dickinson SFR. A simple power-law model with two free parameters, namely the power-law index $\alpha$ and the fraction of PBHs over ABHs, defined as

can cover a broad range of PBH generation scenarios (Raidal et al., 2017; Raidal et al., 2019; Vaskonen & Veermäe, 2020; Atal et al., 2020; De Luca et al., 2020a). We will see that any positive value of the power-law index $\alpha$ and ratio $f_{PBH/ABH}$ can rule out the contribution of PBHs as dark matter candidates for the mass ranges accessible from the LIGO/Virgo detector network. The power-law functional form given in Eq. (7) is also capable of capturing models beyond the PBH scenario such as the contribution from population-II/population-III sources below redshift $z=6$ (Inayoshi et al., 2021). In future work, we will apply this technique to the population-II/population-III sources to distinguish these sources from the low-redshift ABHs and the PBH population.

With these two components, the hybrid model, including both ABH and PBH parts, captures both the low redshift GW sources from the stellar origin, which has a bump at a redshift $z_{p}$ depending on the value of time-delay, and also a monotonically increasing function which captures the redshift evolution of the PBH merger rate.

We show in Fig. 1 a few examples of the hybrid merger rate for a power-law distribution of the time-delays $(t^{eff}_{d})^{-1}$ with a minimum value of the time-delay parameter of $10$ Myr. The PBH component is shown for $f_{PBH/ABH}=1$ (by dotted lines) and $f_{PBH/ABH}=0.1$ (solid lines). The value of the local merger rate is taken to be $R^{GW}_{0}=30$ Gpc^-3 yr^-1. The bump around redshift $z\sim 2$, due to the peak in the star formation rate, is evident for values of $\alpha<1.28$. The power-law index of the PBH merger rate and the fraction of PBH over ABH $f_{PBH/ABH}$ determines the relative strengths of the ABH and PBH contributions if the local merger rate is fixed. So, for any large values of the parameter $\alpha>1.28$, and higher values of $f_{PBH/ABH}$, the bump due to the SFR gets obscured, whereas when the value of $\alpha$ and $f_{PBH/ABH}$ is small, the SFR bump is more prominent. The bump moves towards a lower value of redshift if the time delay is large. The models shown here are for different values of $\alpha$, but different scenarios of PBH formation are described by the same merger rate. Larger values of $\alpha$, with the local merger rate in agreement with the LIGO/Virgo GWTC-2 (Abbott et al., 2020a, b), mimic the cases where the PBHs have large spatial clustering and the fraction of PBH as dark matter is close to unity (Atal et al., 2020). Similarly the cases with small values of $\alpha$ denote captures models with less spatial clustering and fraction of PBH in dark matter less than $0.1$ (Sasaki et al., 2016; Raidal et al., 2017; Raidal et al., 2019; Young & Byrnes, 2020; Atal et al., 2020; De Luca et al., 2020b).

Mergers of the GW sources at high redshift contribute to the stochastic GW background which can be written as (Allen, 1996; Phinney, 2001)

where $d_{L}$ is the luminosity distance, $\theta$ denotes the source properties such as masses of the binary black holes, their spins, inclination angle, and $\frac{dE_{GW}}{df_{r}}(\theta)$ is the energy emission per frequency bin in the source frame, written in terms of the source properties and chirp masses $\mathcal{M}_{c}$ of the gravitational wave sources as

Here $\mathcal{G}(f_{r})$ captures the frequency dependence during the inspiral, merger, and ringdown phases of the gravitational wave signal (Ajith et al., 2008)

where $f_{x}=c^{3}(a_{1}\eta^{2}+a_{2}\eta+a_{3})/\pi GM$ written in terms of total mass $M=m_{1}+m_{2}$ and symmetric mass ratio $\eta=m_{1}m_{2}/M^{2}$. A GW binary will be emitting gravitational waves in the inspiral part up to frequency $f_{merg}$, followed by the ringdown part up to frequency $f_{ring}$, and will stop emitting a gravitational wave signal after $f_{cut}$. $f_{w}$ denotes the width of the Lorentzian function. The values of the parameters $a_{1},a_{2}$, and $a_{3}$ are given in table 1 (Ajith et al., 2008). In this analysis, we have ignored the ffect spin of the GW sources which does not make a significant difference in the GW spectrum at low frequencies Zhu et al. (2011).

We consider two different scenarios for the probability distribution of the black hole masses. For the ABHs, we consider the mass distribution of individual sources to be power-law with $m_{i}^{-2.3}$ for the heavier mass and flat in log-space for the lighter mass, following the previous stochastic analysis by the LIGO/Virgo collaboration (Abbott et al., 2019). One can also consider other mass distributions for estimating the stochastic background. We will show later, that the bounds obtained from the current data are not susceptible to changes in the mass distribution. For the PBH case, we consider two different mass distributions, (i) power-law profile as for the ABH case, (ii) a log-normal mass distribution with a characteristic mass-scale $M_{c}$ and standard deviation $\sigma$ as

which is motivated by the small-scale density fluctuations (Dolgov & Silk, 1993; Carr et al., 2017). We consider two different values of the characteristic mass-scale $M_{c}$: $30M_{\odot}$ and $1M_{\odot}$, and the corresponding value of the parameter $\sigma$ as $0.5$. We show in Fig. 2 the power spectrum of the stochastic GW background for different values of the time-delay parameter for ABHs ($t_{d}^{eff}=100\,\text{Myr},t_{d}^{eff}=1\,\text{Gyr},t_{d}^{eff}=10\,\text{Gyr}$), and also for PBHs with characteristic mass ($M_{c}=30\,M_{\odot}$ and $M_{c}=1\,M_{\odot}$) and power-law index ($\alpha=1.28,0.5$). The amplitude and shape of the stochastic GW background power spectrum varies with the changes in the model parameters. This leads to a distinguishable signature which can be exploited to help identify ABH and PBH sources. We also show the detector noise curve for the O5 observation run of LIGO/Virgo (Aasi et al., 2015; Acernese et al., 2015) and for A+ sensitivity (Abbott et al., 2018a; Barsotti et al., 2020) in Fig. 2.

We have set up a Bayesian framework (using the Bayes theorem (Price, 1763)) to estimate the parameters of the hybrid merger rate. Using the Bayes theorem, we can estimate the posterior of the power-law index $\alpha$, time-delay parameter $t^{eff}_{d}$, local merger-rate $R^{GW}_{0}$, and characteristic mass scale $M_{c}$ as

where $\Pi(\alpha)$, $\Pi{(t^{eff}_{d})}$, and $\Pi(R^{GW}_{0})$ denote the priors on the parameters $\vec{\theta}\equiv\{\alpha,t^{eff}_{d},R^{GW}_{0}\}$ and $\mathcal{L}(\hat{\Omega}_{GW}|\vec{\theta})$ denotes the likelihood, taken as Gaussian

where $\Omega^{GW}(f)$ is the model of the stochastic GW background signal which depends on the parameters $\alpha,t^{eff}_{d},$ and $R^{GW}_{0}$. The measured cross-correlation signal $\hat{\Omega}^{GW}_{IJ}(f)=20\pi^{2}f^{3}\text{Re}[d_{I}(f)^{*}d_{J}(f)]/3H_{0}^{2}T\gamma_{IJ}(f)$ between the data $d_{I,J}(f)$ from the detectors $I$ and $J$ in the Fourier domain. In this expression, $\gamma_{IJ}(f)$ denotes the overlap reduction function , $T$ denotes the observation time duration, and $H_{0}$ denotes the Hubble constant. The corresponding variance can be written in terms of the one-sided noise power spectrum of the detectors $P_{I,J}$ and frequency resolution $\Delta f$ as $\Sigma_{IJ}(f)=50\pi^{4}f^{6}P_{I}(f)P_{J}(f)/9H_{0}^{4}T\Delta f\gamma^{2}_{IJ}(f)$.

For the analysis of the third observation run of the LIGO/Virgo (Abbott et al., 2021), we have taken flat priors¹¹1$\mathcal{U}[a,b]$ denotes uniform distribution for the values in the range a and b.: on the $\alpha$ parameter from $\mathcal{U}[0,10]$, and on the time-delay parameter $t^{eff}_{d}$ from $\mathcal{U}[0.01,13]$ Gyr. We have taken a flat prior on the local merger rate parameter as $\mathcal{U}[15,38]$ Gpc^-3 yr^-1 according to GWTC-2 for GW compact objects with masses heavier than $5$ $M_{\odot}$ (Abbott et al., 2020a, b), and for masses below $5$ $M_{\odot}$, we have taken a flat prior on the local merger rate as $\mathcal{U}[15,710]$ Gpc^-3 yr^-1. The contributions to the local merger rate can arise from both ABHs and PBHs, and we adopt an agnostic view with no assumptions about the separate populations. We consider that the joint contribution of the local merger rate agrees with the merger rate inferred from GWTC-2 (Abbott et al., 2020a, b). In this analysis, we consider three different models of the mass distribution, (i) log-normal distribution of the model of the PBH masses for a fixed value of the characteristic mass $M_{c}=30\,M_{\odot}$, (ii) log-normal distribution of the model of the PBH masses for a fixed value of the characteristic mass $M_{c}=1\,M_{\odot}$, and (iii) power-law mass distribution as for the ABH distribution with a mass range of $[5,50]M_{\odot}$. Though this mass cutoff at the PISN mass-scale is not motivated for any PBH mass distribution (and only possible for ABHs), we consider this to show whether the current stochastic GW observation from O3 is able to distinguish anything between the log-normal mass distribution and astrophysical mass-distribution of black holes. For all three cases, we have considered the probability distribution of the ABH masses as power-law.

By using the data of stochastic GW background from O3 (Abbott et al., 2021) and the bounds on the local merger rate from the O3a observation (Abbott et al., 2020a, b), we obtain constraints on the log-normal model of the PBH mass distribution for the value of $M_{c}=30\,M_{\odot}$ and ABH mass distribution with minimum mass $M_{min}=5\,M_{\odot}$, and maximum mass $M_{max}=50\,M_{\odot}$ for a power-law distribution, shown in Fig. 3(a). The value of $f_{PBH/ABH}=1$ is kept fixed in this analysis. As we currently do not have any measurement of the stochastic GW signal from O3 observations, these bounds show the limits even in a scenario when $f_{PBH/ABH}=1$, i.e. $50\%$ of the total detected black holes are of primordial origin. The bounds on the parameter $\alpha$ will get weaker if the value of $f_{PBH/ABH}<1$. The three upper panels of each plot show the 1-D posteriors on the parameters $\alpha$, $t_{d}^{eff}$, and $R_{0}^{GW}$ along with the 2-D joint posteriors between the parameters in the lower panels, which is obtained using Eq. (12). The plot also shows the $68\%$ and $95\%$ contours on the parameters. Although we do not obtain any constraint on the time-delay parameter from the non-detection of the stochastic GW background from the O3-run, there is a cut-off in the power-law index for $\alpha>6$. This happens because, for large values of $\alpha$, the merger rates at high redshift are larger, and are constrained by the non-detection of the stochastic GW background. The constraints on the time-delay parameter are weak because the peak of the GW mergers shifts to a lower redshift in the presence of time-delay, and the relative strength at the peak is constrained by the local merger rate, as shown in Fig. 1. For the log-normal case with the PBH mass distribution having $M_{c}=1\,M_{\odot}$, we show the possible constraints from the O3 observation in Fig. 3(b) which is similar in nature to the case $M_{c}=30$ $M_{\odot}$. The constraints on the power-law index $\alpha$ are stronger for the higher values of the local merger rate, as can be seen from Fig. 3(b). The merger rate parameters for the power-law model of the mass distribution exhibits similar constraints as for the log-normal distribution (see Fig. 3(c)). Even though for a power-law black hole mass distribution with a possible mass-cut off at the PISN mass scale (which is a possible scenario only for the ABHs) and differs from the log-normal distribution of PBH masses, the current bounds are not at all susceptible to this difference. From these results, we can conclude that the bounds which are obtained from the non-detection of the stochastic GW background signal from the O3 observation are nearly independent of the mass model used for PBHs. The upper bound from the third observation run on the power-law index of the PBH merger rate are $\sim 2.56_{-1.76}^{+1.64}$ and the time delay parameter is $\sim 6.7_{-4.74}^{+4.22}$ Gyr at $68\%$ C.L. The bounds on the time-delay parameter are driven by the choice of prior (flat prior $[0.01,13]$ Gyr) and there is no constraining power in limiting the value of time-delay from the bound on the stochastic GW background from O3 data. The bound on the time-delay parameter from individual events (Fishbach & Kalogera, 2021) is in agreement with our result.

Although the constraints from the current LIGO/Virgo observations are weak, in the future, the stochastic GW background will become a powerful probe for distinguishing between the populations of ABHs and PBHs. To study the feasibility of the joint estimation of the ABH and PBH merger rates from future stochastic GW background data, we simulate the stochastic GW background signal with the hybrid model of the merger rate adopting the parameters time-delay $t^{eff}_{d}=100$ Myr, $\alpha=1.3$, and local merger rate $R^{GW}_{0}=30$ Gpc^-3 yr^-1 (Abbott et al., 2020b). We assume the relative fraction of sources in ABHs and PBHs are the same ($f_{PBH/ABH}=1$). The mass distribution of the ABHs is taken to be power-law with minimum mass $5\,M_{\odot}$, and maximum mass $50\,M_{\odot}$. The mass distribution for the PBHs is taken to be log-normal with characteristic mass scale $M_{c}=30\,M_{\odot}$ and corresponding standard deviation $\sigma=0.5$. We consider the noise power spectrum of O5²²2https://dcc.ligo.org/LIGO-P1500222-v29/public and A+ ³³3https://dcc.ligo.org/public/0149/T1800042/004/T1800042-v4.pdf for this forecast (Aasi et al., 2015; Acernese et al., 2015; Abbott et al., 2018a; Barsotti et al., 2020). The LIGO/Virgo observation at its design sensitivity for the fifth observation run with $50\%$ duty cycle is denoted by O5. The A+ sensitivity is a future upgrade of the advanced LIGO detectors almost by a factor of two (Barsotti et al., 2020). We consider the A+ noise curve integrated for two years with a duty cycle of $50\%$.

We consider a four-parameter model, namely the time-delay parameter $t^{eff}_{d}$, and the spectral index of the PBH merger rate $\alpha$, the characteristic mass scale $M_{c}$, and the fraction of PBH and ABH $f_{PBH/ABH}$ as the free parameters with posterior given by

where $\Pi(\alpha)$, $\Pi{(t^{eff}_{d})}$, $\Pi(M_{c})$ and $\Pi(f_{PBH/ABH})$ denote the priors on the parameters $\vec{\theta}\equiv\{\alpha,t^{eff}_{d},M_{c},f_{PBH/ABH}\}$. For the forecast, we have taken a flat prior on the power-law index $\alpha$ from $\mathcal{U}[0,10]$, time-delay parameter $t_{d}^{eff}$ from $\mathcal{U}[0.01,10]$ Gyr, characteristic mass from $\mathcal{U}[1,50]\,M_{\odot}$, and fraction of PBHs over ABHs $f_{PBH/ABH}$ from $\mathcal{U}[0,3]$. The local merger rate is kept fixed at $R_{0}^{GW}=30$ Gpc^-3 yr^-1, assuming that it can be constrained from the individual detected events until the O5 and A+ runs. If there is a population of GW sources that has a merger rate that increases with redshift, then it can be isolated from the population of sources which is governed by the Madau-Dickinson law, even if the choice of $f_{PBH/ABH}$ is kept fixed, for fast estimation of the parameters.

The joint estimations are shown in Fig. 4 for the O5 sensitivity in green and for the A+ sensitivity in magenta. The four top panels show the 1-D posteriors on the parameters $\alpha$, $t_{d}^{eff}$, $M_{c}$, and $f_{PBH/ABH}$. The joint 2-D posterior distributions between the parameters are shown in the lower panels in Fig. 4 which are obtained using Eq. (14). We also show the $68\%$ and $95\%$ contours on the parameters in Fig. 4. The results show that we can measure the power-law index $\alpha=1.3$ of the PBH merger rate from future observations with O5 and A+ sensitivities. The lower values of the characteristic mass scale $M_{c}$ of the PBHs can be limited from the stochastic GW background observations for characteristic mass scale $M_{c}=30\,M_{\odot}$. However, the time-delay parameter $t^{eff}_{d}>100$ Myr cannot be inferred very well from the stochastic GW observations of O5 and only a partial improvement is possible from A+. For the fiducial case with a minimum time delay of $100$ Myr, from the stochastic GW background data from A+, we can expect to weakly limit the values of the time-delay parameter greater than about $9$ Gyr. Any interesting bounds from O5 on the time delay parameter are unlikely. The parameter related to the fractional PBH and ABH $f_{PBH/ABH}$ ratio can be constrained using the stochastic GW background data for both O5 and A+ sensitivities. This indicates that joint estimation of the hybrid merger rates will provide an interesting avenue for distinguishing between ABHs and PBHs. We show the $16^{th},\,50^{th},\,84^{th}$ percentiles of the posterior distribution in Table 2. Our forecast shows that for non-zero injected values, the parameters $\alpha,M_{c},$ and $f_{PBH/ABH}$ can be related to the PBH fraction by these observation run. The measurement of non-zero values of the parameters $f_{PBH/ABH}$ and $\alpha$ from future data on the stochastic GW background would confirm the existence of a population of black holes which are different from the population of sources that follow the Madau-Dickinson star formation rate (Madau & Dickinson, 2014) in a model-independent way. Individual PBH production scenarios can be tested by using this technique. With a longer duration of observation time $t$, the constraints on the parameters will improve by $\sqrt{t}$.

In this paper, we show how one may distinguish between different populations of GW sources using the stochastic GW background. The merger rate of the ABHs is likely to follow the star formation rate which can be modeled by the Madau Dickinson relation (Madau & Dickinson, 2014), whereas, for the PBHs (or equally for population-II/population-III sources), the merger rate is going to be different at high redshift (Raidal et al., 2017; Raidal et al., 2019; Vaskonen & Veermäe, 2020; Atal et al., 2020; De Luca et al., 2020a). The merger rate of PBHs and their redshift evolution also going to depend on whether they are spatially clustered or having a Poisson distribution. As a result, different populations of GW sources present at early epochs can be distinguished by using the stochastic GW background to explore the high redshift Universe. In the future a joint multi-messenger study of the stochastic GW background and different probes of star formation rate using an electromagnetic signal will further improve the capability to distinguish between the population of ABHs and PBHs.

We construct a hybrid model of the GW merger rates and source populations that is characterized by four parameters, namely the local merger rate of GW sources $R_{0}^{GW}$, the astrophysical time delay between formation and merger for ABHs $t^{eff}_{d}$, the power-law index of the merger rate for PBHs $\alpha$, and the characteristic mass scale $M_{c}$ of PBHs. We obtain constraints on $R_{0}^{GW},\,\alpha,\,t^{eff}_{d}$ from the stochastic GW background data of third observing run of the LIGO/Virgo collaboration O3 (Abbott et al., 2021) and the bounds on the local merger rate from individual events of the O3a (Abbott et al., 2020a, b) for three different mass choices (see Fig. 3). Current data can only rule out very large values of the parameter $\alpha$ and can impose no constraints on the time-delay parameter. However, in the near future from O5 and A+ sensitivities, the stochastic GW background should be a powerful probe that is capable of distinguishing between different populations of the GW sources. We show a forecast in Fig. 4, which indicates that bounds on the parameter $\alpha=1.28$, and weak constraints on the time-delay parameter $t_{d}^{eff}=100$ Myr, characteristic mass scale of PBHs $M_{c}=30$ M_⊙, and the fraction of PBH/ABH $f_{PBH/ABH}=1$ should be feasible with O5 (Abbott et al., 2016a) and A+ (Abbott et al., 2018a; Barsotti et al., 2020) sensitivities. Due to the correlations between the parameters, our study does not give a very large gain in constraining these parameters. But interesting conclusions on whether there exists a population of black holes with merger rate $(1+z)^{\alpha}$ and the possible mass-scale of these sources should be achievable from A+, if not feasible from the O5 sensitivity. By using the posteriors on the PBH population parameters such as $\alpha$, $M_{c}$ and $f_{PBH/ABH}$, one can explore different models of PBH formation scenarios and can constrain the parameter space of those models (Zel’dovich & Novikov, 1967; Hawking, 1971; Carr, 1975; Khlopov & Polnarev, 1980; Khlopov et al., 1985; Carr, 2005; Clesse & García-Bellido, 2017; Ali-Haïmoud et al., 2017; Raidal et al., 2017; Raidal et al., 2019; Sasaki et al., 2018; Jenkins & Sakellariadou, 2020; Atal et al., 2020).

In future work, we will explore the measurability of the stochastic GW signal from space-based GW detectors such as Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al., 2017) and the next generation ground-based detectors such as the Einstein Telescope (Punturo et al., 2010), and the Cosmic Explorer (Reitze et al., 2019). For future GW detectors, we will be able to identify individual sources up to higher redshift ($z\sim 50$) than the detector horizon of current generation networks of detectors ($z\sim 1$). As a result, sources that are present only at a very high redshift ($z\gtrsim 50$) will contribute to the stochastic GW background signal for the next generation GW detectors. We should be able to unveil the population of high redshift GW sources. It should also be possible to explore the possible time dependence of the stochastic GW background (Mukherjee & Silk, 2019, 2020). This could provide an independent means of distinguishing between different kinds of contributing sources.

The authors are very thankful to Jishnu Suresh for reviewing the manuscript and providing very useful comments. S. M. also thanks to the group members of the LSC stochastic group for stimulating discussion on this topic. This work is part of the Delta ITP consortium, a program of the Netherlands Organisation for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture, and Science (OCW). This analysis is carried out at the Horizon cluster hosted by Institut d’Astrophysique de Paris. We thank Stephane Rouberol for smoothly running the Horizon cluster. We acknowledge the use of following packages in this work: Astropy (Astropy Collaboration et al., 2013, 2018), Corner (Foreman-Mackey, 2016), emcee: The MCMC Hammer (Foreman-Mackey et al., 2013), IPython (Pérez & Granger, 2007), Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), and SciPy (Jones et al., 01). The authors would like to thank the LIGO/Virgo scientific collaboration for providing the noise curves. 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. This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation.

The data underlying this article will be shared at request to the corresponding author.