Possible Discrimination of Black Hole Origins from the Lensing Rate of DECIGO and B-DECIGO Sources

Up to now, ground-based detectors have reported about 90 gravitational wave (GW) events and most of them are from the coalescence of binary black holes (BBHs). These detections are very helpful for exploring various aspects of astrophysics such as testing general relativity (Abbott et al. 2021b, 2019a, 2019b) and constraining the Hubble constant (Abbott et al. 2021a, 2017). However, the origin of these binary black holes is still an open issue. Since the first detection of GW emitted from the coalescence of BBHs by the Advanced LIGO (Abbott et al. 2016), the possibility that this signal was emitted from binary primordial black holes (PBHs) coalescence (Bird et al. 2016; Sasaki et al. 2016; Clesse & García-Bellido 2017; Khalouei et al. 2021) instead of astrophysical black holes (ABHs) was widely investigated. PBHs could be the results of the gravitational collapse of primordial fluctuations in the early universe (Hawking 1971; Carr & Hawking 1974), or other mechanisms such as the quantum fluctuation (Clesse & Garcia-Bellido 2015; Chen & Cai 2019; Fu et al. 2019; Motohashi et al. 2020), the bubble collisions (Hawking et al. 1982), the cosmic string (Hogan 1984; Hawking 1989), and the domain wall (Caldwell et al. 1996). The diversity of formation mechanism means that PBHs might have a different evolutionary history with respect to ABHs.

The successful detection of GW with current ground-based detectors has significantly boosted this field to the new era. Several next-generation detectors, including both ground and space missions, will be put into operation in the near future with higher sensitivity. Meanwhile, covering wider frequency ranges by different facilities is also an important goal. The DECi-hertz Interferometer Gravitational wave Observatory (DECIGO; Seto et al. 2001; Kawamura et al. 2019) is a future plan of a Japanese space mission for observing GWs around f ∼ 0.1 Hz and a small-scale version called B-DECIGO is planned to be launched within this decade. DECIGO will be able to detect the GW signals of inspiralling double compact objects (DCOs) and will receive these signals earlier than ground-based detectors due to the lower detection frequency range (mHz–100 Hz). It will also has a higher sensitivity than the next-generation ground-based detectors. Thus, a large number of GW events are expected to be detected by this detector. It is shown that the third-generation ground-based detector, i.e., the Einstein Telescope, will detect 10⁵–10⁶ DCO GW events per year (Piórkowska et al. 2013; Biesiada et al. 2014; Ding et al. 2015). With a higher sensitivity, DECIGO is expected to detect the same order or even more GW events.

Given the large number of detections and further reach of the detection range, we can expect some GW events will be gravitationally lensed. LIGO-Virgo collaboration has tried to search the lensed GWs from current available observations, and concluded a null search result in the first half of O3 data (The LIGO Scientific Collaboration et al. 2021). Meanwhile, it was argued that, with certain priors, the lensed GW has already been detected (Diego et al. 2021). In any case, we will be able to accumulate numbers of lensed GW signals in the future with the help of next-generation detectors. It was shown that the lensed BH-BH systems could be observed at the rate of 50 yr⁻¹ in DECIGO and a few in B-DECIGO (Piorkowska-Kurpas et al. 2021).

In this paper, we calculated the predicted event rates of strongly lensed GW signals from PBHs and ABHs detected by DECIGO and B-DECIGO, and explored the differences in their event rates and redshift distributions. This can serve as an alternative method to distinguish PBHs and ABHs. This paper is organized as follows. In Section 2, we introduce three theoretical mass distributions of PBH and two different stellar formation models of ABH used in this paper. In Section 3, we calculate the detection rate of unlensed events. We discuss lensed events in Section 4. Finally, we summarize our results and draw conclusions in Section 5. Throughout this paper, we adopt a flat ΛCDM cosmological model with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3.

To calculate the observable events for ABH models, we take the merger rate of a convolution of the birthrate of ABHs as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{\mathrm{ABH}}(z,{m}_{1}) & = & {\int }_{{t}_{\min }}^{{t}_{\max }}{R}_{{\rm{b}}}(t(z)-{t}_{{\rm{d}}},{m}_{1})\\ & & \times P({t}_{{\rm{d}}}){{dt}}_{{\rm{d}}},\end{array}\end{eqnarray} \tag{ 1 }$$

where the minimum time delay is set as ${t}_{\min }$ = 50 Myr, ${t}_{\max }$ is set as the Hubble time, t_d is the time delay, and t(z) is the age of the universe at the moment of merger. In addition, R_b(t, m₁) can be estimated from (Dvorkin et al. 2016)

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{\rm{b}}}(t,{m}_{\mathrm{bh}}) & = & {\int }_{0}^{\infty }\psi [t-\tau (m^{\prime} )]\phi (m^{\prime} )\\ & & \times \delta (m^{\prime} -{g}_{{\rm{b}}}^{-1}({m}_{\mathrm{bh}})){dm}^{\prime} ,\end{array}\end{eqnarray} \tag{ 2 }$$

where m_bh is the mass of the remnant black hole, τ(m) is the lifetime of a progenitor star, and ϕ(m) ∝ m^−2.35 is the initial mass function. In order to get ${g}_{{\rm{b}}}^{-1}({m}_{\mathrm{bh}})$, we should consider the relation as

$$\begin{eqnarray}&&\displaystyle \frac{{m}_{\mathrm{bh}}}{m}=A{\left(\displaystyle \frac{m}{40\,{M}_{\odot }}\right)}^{\beta }\displaystyle \frac{1}{{\left(\tfrac{Z}{0.01\,{Z}_{\odot }}\right)}^{\gamma }+1},\end{eqnarray} \tag{ 3 }$$

where Z is the metallicity Belczynski et al. (2016). The values of parameters [A, β, γ] are taken as [0.3, 0.8, 0.2] (Dvorkin et al. 2016). The star formation rate ψ(t) in R_b(t, m₁) is given by

$$\begin{eqnarray}&&\psi (z)=k\displaystyle \frac{a\,\exp [b(z-{z}_{m})]}{a-b+b\,\exp [a(z-{z}_{m})]},\end{eqnarray} \tag{ 4 }$$

where the parameters [k, a, b, z_m ] are dependent on the specific ABH formation scenarios. Here, we consider two ABH models. The first one is the Fiducial model, which corresponds to the classical isolated binary evolution of ABH. The parameters {k, a, b, z_m } for the Fiducial model are from Vangioni et al. (2015), i.e., [0.178 M_⊙ yr⁻¹ Mpc⁻³, 2.37, 1.80, 2.00]. The second one is the Gamma-ray Burst (GRB)-based model, which uses the star formation rate calibrated from the GRB at high redshift. The values for these parameters are set as [0.146 M_⊙ yr⁻¹ Mpc⁻³, 2.80, 2.46, 1.72] (Vangioni et al. 2015). For local merger rate ${R}_{0,\mathrm{ABH}}\equiv {\int }_{\geqslant 1\,{M}_{\odot }}{{ \mathcal R }}_{\mathrm{ABH}}(z,{m}_{1}){{dm}}_{1}$, we take it as ${23.9}_{-8.6}^{+14.3}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$ in Abbott et al. (2021c).

In order to obtain the GW event rates of PBH systems, their intrinsic merger rate must be calculated first as well. For a general mass function P(m), the time-dependent comoving merger rate density is (Chen & Huang 2018)

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal R }}_{\mathrm{PBH}}({m}_{1},{m}_{2},z| \theta ,{f}_{\mathrm{PBH}})=2.8\times {10}^{6}{\left(\tfrac{t(z)}{{t}_{0}}\right)}^{-\tfrac{34}{37}}{f}_{\mathrm{PBH}}^{2}\\ \quad \times \,{\left(0.7{f}_{\mathrm{PBH}}^{2}+{\sigma }_{\mathrm{eq}}^{2}\right)}^{-\tfrac{21}{74}}\times \min \left(\tfrac{P({m}_{1}| \theta )}{{m}_{1}},\tfrac{P({m}_{2}| \theta )}{{m}_{2}}\right)\\ \quad \times \,{\left({m}_{1}{m}_{2}\right)}^{\tfrac{3}{37}}{\left({m}_{1}+{m}_{2}\right)}^{\tfrac{36}{37}}\left(\tfrac{P({m}_{1}| \theta )}{{m}_{1}}+\tfrac{P({m}_{2}| \theta )}{{m}_{2}}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

in units of ${M}_{\odot }^{-2}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$, where θ represents the parameter vector for each mass distribution model. ${f}_{\mathrm{PBH}}=\tfrac{{{\rm{\Omega }}}_{\mathrm{PBH}}}{{{\rm{\Omega }}}_{\mathrm{CDM}}}$ is the abundance of PHBs in the dark matter at the present universe. σ_ep is the variance of the rest dark matter density perturbations at scale ${ \mathcal O }({10}^{0}\mbox{--}{10}^{3}){M}_{\odot }$ and we usually take the value of it as σ_eq ≈ 0.005 (Ali-Haïmoud et al. 2017; Chen & Huang 2018). By integrating over the mass, we can get the merger rate of all PBHs at cosmic time t:

$$\begin{eqnarray}&&{R}_{\mathrm{PBH}}(z| \theta ,{f}_{\mathrm{PBH}})={\int }_{0}^{\infty }{{ \mathcal R }}_{\mathrm{PBH}}({m}_{1},{m}_{2},z| \theta ,{f}_{\mathrm{PBH}}){{dm}}_{1}{{dm}}_{2}\end{eqnarray} \tag{ 6 }$$

and the local merger rate is R_0,PBH ≡ R_PBH(z = 0∣θ, f_PBH). Here we adopt three mass functions of PBHs. These are the power-law mass function, the lognormal mass function, and the critical collapse mass function. The power-law mass function can be written as

$$\begin{eqnarray}&&P\left(m,\alpha ,{M}_{\min }\right)=\displaystyle \frac{\alpha -1}{{M}_{\min }}{\left(\displaystyle \frac{m}{{M}_{\min }}\right)}^{-\alpha }\end{eqnarray} \tag{ 7 }$$

with $m\geqslant {M}_{\min }$ and α > 1. This mass function is a natural prediction when PBHs are generated by scale-invariant density fluctuations or from the collapse of cosmic strings. The lognormal mass function is expressed as (Carr et al. 2017; Bellomo et al. 2018; Hütsi et al. 2021)

$$\begin{eqnarray}&&P\left(m,\sigma ,{m}_{{\rm{c}}}\right)=\displaystyle \frac{1}{\sqrt{2\pi }\sigma m}\exp \left(-\displaystyle \frac{{[\mathrm{ln}\left(m/{m}_{{\rm{c}}}\right)]}^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 8 }$$

where m_c and σ denote the peak mass of mP(m) and the width of mass spectrum, respectively. This mass function is often a good approximation if the PBHs are produced from a smooth symmetric peak in the inflationary power spectrum. And the last mass function—the critical collapse distribution—is (Carr et al. 2016, 2017; Hütsi et al. 2021)

$$\begin{eqnarray}&&P\left(m,{M}_{f}\right)=\displaystyle \frac{{\gamma }^{2}}{{M}_{f}^{1+\gamma }{\rm{\Gamma }}(1/\gamma )}{m}^{\gamma }\exp \left(-{\left(\displaystyle \frac{m}{{M}_{f}}\right)}^{\gamma }\right),\end{eqnarray} \tag{ 9 }$$

where M_f is a mass-scale that corresponds to the horizon mass at the collapse epoch, γ is a universal exponent that is related to the critical collapse of radiation, and Γ(x) denotes the gamma function. This mass function is supposed to be closely related to PBHs originating from density fluctuations with a δ-function power spectrum.

To extract the population parameters {θ, f_PBH} for PBH models from N_obs detections of GW events $d=\{{d}_{1},\ldots {d}_{{N}_{\mathrm{obs}}}\}$, we can perform the hierarchical Bayesian inference (Chen et al. 2019; Mandel et al. 2019; Chen & Huang 2020; Wu 2020; Hütsi et al. 2021). First, the likelihood for N_obs BBH events is

$$\begin{eqnarray}\begin{array}{rcl}p(d| \theta ,{f}_{\mathrm{PBH}}) & = & {\left(N(\theta ,{f}_{\mathrm{PBH}})\right)}^{{N}_{\mathrm{obs}}}{e}^{-N(\theta ,{f}_{\mathrm{PBH}})\xi (\theta ,{f}_{\mathrm{PBH}})}\\ & & \times \displaystyle \prod _{i}^{{N}_{\mathrm{obs}}}\int d\lambda L({d}_{i}| \lambda ){p}_{\mathrm{pop}}(\lambda | \theta ,{f}_{\mathrm{PBH}}),\end{array}\end{eqnarray} \tag{ 10 }$$

where λ ≡ {m₁, m₂, z}, the likelihood of one BBH event L(d_i ∣λ) is proportional to the posterior p(λ∣d_i ), and N(θ, f_PBH) is the total number of events in the model characterized by the set of population parameters {θ, f_PBH} as

$$\begin{eqnarray}&&N(\theta ,{f}_{\mathrm{PBH}})=\int d\lambda {T}_{\mathrm{obs}}{{ \mathcal R }}_{\mathrm{PBH}}(\lambda | \theta ,{f}_{\mathrm{PBH}})\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dz}},\end{eqnarray} \tag{ 11 }$$

where dV_c/dz is the differential comoving volume, the factor 1/(1 + z) accounts for the cosmological time dilation from source frame to the detector frame, and T_obs is effective observing time of LIGO O1-O3. p_pop(λ∣θ, f_PBH) is the distribution of BH masses and redshifts in coalescing binaries as

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{pop}}(\lambda | \theta ,{f}_{\mathrm{PBH}})=\displaystyle \frac{1}{N(\theta ,{f}_{\mathrm{PBH}})}\\ \,\times \left[{T}_{\mathrm{obs}}{{ \mathcal R }}_{\mathrm{PBH}}(\lambda | \theta ,{f}_{\mathrm{PBH}})\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dz}}\right].\end{array}\end{eqnarray} \tag{ 12 }$$

Meanwhile, ξ(θ, f_PBH) can be calculated by

$$\begin{eqnarray}&&\xi (\theta ,{f}_{\mathrm{PBH}})\approx \displaystyle \frac{1}{{N}_{\mathrm{inj}}}\displaystyle \sum _{k=1}^{{N}_{\det }}\displaystyle \frac{{p}_{\mathrm{pop}}({\lambda }_{k}| \theta ,{f}_{\mathrm{PBH}})}{{p}_{\mathrm{draw}}({\lambda }_{k})},\end{eqnarray} \tag{ 13 }$$

where N_inj is the total number of injections, ${N}_{\det }$ is the number of injections that are successfully detected, and p_draw is the probability distribution from which the injections are drawn. Then the posterior probability distribution p(θ, f_PBH∣d) can be calculated from

$$\begin{eqnarray}&&p(\theta ,{f}_{\mathrm{PBH}}| d)=p(d| \theta ,{f}_{\mathrm{PBH}})p(\theta ,{f}_{\mathrm{PBH}}),\end{eqnarray} \tag{ 14 }$$

where p(θ, f_PBH) is the prior distribution for θ and f_PBH. We use the flat and log-uniform prior for the θ and f_PBH, respectively.

We adopt the available BBHs data from GWTC-3 (Abbott et al. 2021); these GW data satisfy the following criteria: black hole masses (m₁ and m₂) larger than 3 M_⊙, effective spin (∣χ_eff∣) lower than 0.3, network-matched filter signal-to-noise ratio (S/N) larger than 9 and inverse false alarm rate (ifar) higher than 1 yr. Totally, there are 58 events from LIGO's O1, O2, O3a, and O3b runs. Based on these data, we incorporate the PBH population distribution into the ICAROGW to estimate the likelihood function and use Bilby to search over the parameter space (Ashton et al. 2019; Mastrogiovanni et al. 2021). Finally, the results of population parameters {θ, f_PBH} of different PBH models are shown in Table 1. Although, there are still many discussions about the origin mechanism of GWs of BBH. If the range of fraction of PBH in dark matter (DM) at the present universe can be from 10⁻³ to 10⁻², the merging model of PBH binaries (Equation (5)) can provide enough BBH events to explain the case theoretically (Chen & Huang 2018). Furthermore, possibilities of the above merging model of PBH binaries have been explored by combining the PBH population and GW events detected by LIGO-Virgo (Chen et al. 2019; Chen & Huang 2020; Wu 2020). In addition, some studies suggest that the current GW data prefer the binaries model of mixed PBH and ABH populations (De Luca et al. 2021; Hütsi et al. 2021; Franciolini et al. 2022). Therefore, this may require more GW data to discriminate in the future.

Table 1. The Results for Population Parameters {θ, f_PBH} of Different PBH Models and Best-fit Value of Local Merger Rate of Three PBH Models

Mass Function	Parameter 1	Parameter 2 (M⊙)	${\mathrm{log}}_{10}{f}_{\mathrm{PBH}}$	R0,PBH (Gpc−3 yr−1)
Power law	$\alpha ={1.91}_{-0.11}^{+0.10}$	${M}_{\min }={6.43}_{-0.45}^{+0.30}$	$-{2.74}_{-0.03}^{+0.03}$	22
Lognormal	$\sigma ={0.76}_{-0.06}^{+0.08}$	${m}_{c}={16.78}_{-1.54}^{+1.57}$	$-{2.81}_{-0.04}^{+0.04}$	20
Critical collapse	$\gamma ={0.86}_{-0.12}^{+0.12}$	${M}_{f}={8.13}_{-2.32}^{+2.32}$	$-{2.80}_{-0.04}^{+0.04}$	24

Download table as:  ASCIITypeset image

As shown in Figure 1, the merger rates of binary PBHs and ABHs evolve differently with respect to redshift, especially at high redshift. By studying the redshift distribution of GW events, we can distinguish PBHs and ABHs. We first calculated the unlensed event rates of these sources. The yearly detection rate of the DCOs at redshift z_s with S/N exceeding the detector's threshold is

$$\begin{eqnarray}&&\dot{N}(\gt {\rho }_{0}| {z}_{s})={\int }_{0}^{{z}_{s}}\displaystyle \frac{d\dot{N}(\gt {\rho }_{0})}{{dz}}{dz},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\dot{N}(\gt {\rho }_{0})}{{{dz}}_{s}} & = & 4\pi {\left(\displaystyle \frac{c}{{H}_{0}}\right)}^{3}\displaystyle \frac{{R}_{0}({z}_{s})}{1+{z}_{s}}\displaystyle \frac{{\tilde{r}}^{2}({z}_{s})}{E({z}_{s})}\\ & & \times {C}_{{\rm{\Theta }}}(x({z}_{s},{\rho }_{0})).\end{array}\end{eqnarray} \tag{ 16 }$$

Here we set the threshold ρ₀ = 8. R₀(z_s ) denotes the intrinsic merger rate of ABHs or PBHs at redshift z_s as mentioned above. ${\tilde{r}}^{2}({z}_{s})={\int }_{0}^{z}\tfrac{{z}^{{\prime} }}{E({z}^{{\prime} })}$ is the nondimensional comoving distance with E(z) = H(z)/H(0). ${C}_{{\rm{\Theta }}}(x)={\int }_{x}^{\infty }{P}_{{\rm{\Theta }}}d{\rm{\Theta }}$ with the probability distribution of Θ, ${P}_{{\rm{\Theta }}}({\rm{\Theta }})=5{\rm{\Theta }}{\left(4-{\rm{\Theta }}\right)}^{3}/256$ for 0 < Θ < 4. $x(z,\rho )=\tfrac{\rho }{8}{\left(1+z\right)}^{1/6}\tfrac{c}{{H}_{0}}\tfrac{\tilde{r}(z)}{{r}_{0}}{\left(\tfrac{1.2{M}_{\odot }}{{{ \mathcal M }}_{0}}\right)}^{5/6}$ which is correlated with the detector's characteristic distance parameter r₀. This parameter depends on the noise spectrum of the detector and can be estimated as

$$\begin{eqnarray}&&{r}_{0}^{2}=\displaystyle \frac{5}{192{\pi }^{4/3}}{\left(\displaystyle \frac{3}{20}\right)}^{5/3}\displaystyle \frac{{\left({{GM}}_{\odot }\right)}^{5/3}}{{c}^{3}}{f}_{7/3}.\end{eqnarray} \tag{ 17 }$$

where ${f}_{7/3}={\int }_{0}^{\infty }{[{f}^{7/3}{S}_{h}(f)]}^{-1}{df}$. The noise spectrum S_h (f) of DECIGO is (Kawamura et al. 2011)

$$\begin{eqnarray}\begin{array}{rcl}{S}_{h}(f) & = & {10}^{-48}\times \left[7.05\left(1+\displaystyle \frac{{f}^{2}}{{f}_{p}^{2}}\right)+4.8\times {10}^{-3}\right.\\ & & \left.\times \displaystyle \frac{{f}^{-4}}{1+\tfrac{{f}^{2}}{{f}_{p}^{2}}}+5.33\times {10}^{-4}{f}^{-4}\right]{\mathrm{Hz}}^{-1},\end{array}\end{eqnarray} \tag{ 18 }$$

where f_p = 7.36 Hz. And the noise spectrum of B-DECIGO is

$$\begin{eqnarray}\begin{array}{rcl}{S}_{h}(f) & = & {10}^{-46}\times \left[4.040+6.399\times {10}^{-2}{f}^{-4}\right.\\ & & \left.+6.399\times {10}^{-3}{f}^{2}\right]{\mathrm{Hz}}^{-1}.\end{array}\end{eqnarray} \tag{ 19 }$$

The confusion noise from unresolved DCO systems will influence the detection results, so we include the noise spectrum caused by these systems as well. The sky-averaged spectrum of these GW backgrounds can be expressed as

$$\begin{eqnarray}&&{S}_{h}^{\mathrm{GW},\mathrm{DCO}}=\displaystyle \frac{4}{\pi }{f}^{-3}{\rho }_{\mathrm{cr}}{{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{DCO}},\end{eqnarray} \tag{ 20 }$$

where ρ_cr is the critical density of the universe and ${{\rm{\Omega }}}_{\mathrm{GW}}^{\mathrm{DCO}}$ is the energy density parameter corresponding to the unresolved signals from the DCO systems (refer to Piorkowska-Kurpas et al. 2021 for details). The total spectrum is ${S}_{h}(f)+{S}_{h}^{\mathrm{GW},\mathrm{DCO}}(f)$. Substituting the total spectrum into Equation (17), the final characteristic distance r₀ is different in different BBH scenarios. Table 2 listed the results of DECIGO and B-DECIGO with the BBH systems used in this paper.

Zoom In Zoom Out Reset image size

Combining the detector characteristic distance and intrinsic merger rate for a given model, we can calculate the redshift distributions of observed BBH events for both DECIGO and B-DECIGO. The results are shown in Figures 2 and 3. As shown in Figure 2 ABHs and PBHs follow different redshift distributions. That is, the ABH events assemble at lower redshifts while the PBH events have a more flattened shape and have more detectable events at higher redshifts. The total yearly event rates were also estimated from Equation (15). The results are shown in Table 3. It suggests that there are significant differences between the event rates of ABH models and the ones of PBH models. We will expect to observe 1 order of magnitude more PBH events than ABH. For B-DECIGO, due to the smaller characteristic distance, fewer high-redshift merger events could be detected than DECIGO. Though the differences are not as significant as DECIGO, we can still distinguish the distributions of ABHs and PBHs in the redshift range z < 8 as shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With numerous DCO GW events supposed to be detected, we naturally expect some of such events to be gravitationally lensed by intervening galaxies. The lensing optical depth can be expressed as

$$\begin{eqnarray}&&{\tau }_{\pm }=\displaystyle \frac{16}{30}{\pi }^{3}{\left(\displaystyle \frac{c}{{H}_{0}}\right)}^{3}\tilde{r}{\left({z}_{s}\right)}^{3}{\left(\displaystyle \frac{{\sigma }_{* }}{c}\right)}^{4}{n}_{* }\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{4+\alpha }{\beta }\right)}{{\rm{\Gamma }}\left(\tfrac{\alpha }{\beta }\right)}{y}_{\pm ,\max }^{2}.\end{eqnarray} \tag{ 21 }$$

We assume conservatively that the lensing galaxies are all elliptical galaxies and can be characterized by the singular isothermal spheres model. The parameters of velocity dispersion distribution function are then taken as σ_* = 161 ±5 km s⁻¹, n_* = 8.0 × 10⁻³ h³ Mpc⁻³, α = 2.32 ± 0.10, and β = 2.67 ± 0.07 (Choi et al. 2007). τ₊ and τ₋ denote the optical depth of brighter image and fainter image, respectively. y is the misalignment of the source with respect to the optical axis of the lens and is related to the magnification factors by ${\mu }_{\pm }=\tfrac{1}{y}\pm 1$. Assuming the detection threshold ρ₀ = 8, the misalignment should not exceed the limiting value ${y}_{\pm ,\max }={\left[{\left(\displaystyle \frac{8}{{\rho }_{\mathrm{int}r}}\right)}^{2}\mp 1\right]}^{-1}$ (see Piórkowska et al. 2013; Biesiada et al. 2014; Ding et al. 2015 for details).

The intrinsic low S/N (ρ_int < 8) signal will be detectable due to the lensing magnification effect, this will increase the detection rate to

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\dot{N}}{\partial {z}_{s}\partial \rho }=4\pi {\left(\displaystyle \frac{c}{{H}_{0}}\right)}^{3}\displaystyle \frac{{R}_{0}({z}_{s})}{1+{z}_{s}}\displaystyle \frac{{\tilde{r}}^{2}({z}_{s})}{E({z}_{s})}{P}_{{\rm{\Theta }}}(x({z}_{s},\rho ))\displaystyle \frac{x\left({z}_{s},\rho \right)}{\rho }.\end{eqnarray} \tag{ 22 }$$

Combining this detection rate with the lensing depth, the rate of lensed GW events are

$$\begin{eqnarray}&&{\dot{N}}_{\mathrm{lensed}}={\int }_{0}^{{z}_{s}}{{dz}}_{s}{\int }_{0}^{{\rho }_{0}}\tau ({z}_{s},\rho )\displaystyle \frac{{\partial }^{2}\dot{N}}{\partial {z}_{s}\partial \rho }d\rho .\end{eqnarray} \tag{ 23 }$$

The total detection rate includes both the case with weaker intrinsic signal ρ_int < 8 (but exceeds 8 after lensing magnification) and the one with stronger signal ρ_int > 8. The above calculation considers the detection as a continuous survey. In practice, detectors have a finite mission duration, and we will miss some lensing signals due to the time delay. Therefore, with this effect taken into account, the optical depth should be corrected as

$$\begin{eqnarray}&&{\tau }_{{\rm{\Delta }},\pm }=\tau \left[1-\displaystyle \frac{1}{7}\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{\alpha +8}{\beta }\right)}{{\rm{\Gamma }}\left(\tfrac{\alpha +4}{\beta }\right)}\displaystyle \frac{{\rm{\Delta }}{t}_{* ,\pm }}{{T}_{\mathrm{surv}}}\right],\end{eqnarray} \tag{ 24 }$$

where ${\rm{\Delta }}{t}_{* ,\pm }=\tfrac{32{\pi }^{2}}{{H}_{0}}\tilde{r}({z}_{s}){\left(\tfrac{{\sigma }_{* }}{c}\right)}^{4}{y}_{\pm ,\max }$ and T_surv is the mission duration of the telescope. Here we obtain the event rates of the lensed GW of PBHs and ABHs by considering T_surv = 4 yr as planned by DECIGO and B-DECIGO. The results are shown in Table 4. The differential lensing rates are also illustrated in Figures 4 and 5. The high-redshift proportion of PBHs detected by both DECIGO and B-DECIGO will increase as shown in Figures 2, 4 and 3, 5. So, the distributions of ABH and PBH events are more distinguishable. These differences are reflected in the detection numbers as well. The lensed signals have a higher ratio between PBH and ABH than no-lensed signals as shown in Tables 4 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper, we yield predictions for the detection rates of the GW signals from coalescences of BBHs using DECIGO and B-DECIGO. These will provide us a complementary tool to distinguish the ABHs and PBHs, especially using the lensed signals.

The detection capability of DECIGO is significantly improved, and it is likely to detect 10⁵–10⁶ BBHs every year. This will enable us to distinguish the differences in event rates and redshift distributions of GW from ABHs and PBHs. It suggests that the predicted event number of GW signals from the coalescence of PBHs is ∼4 times larger than the one from the coalescence of ABHs. This difference implies that we can probably explore the origins of BBHs, i.e., PBHs or ABHs, with the detection of GW signals from the coalescence of them. Though the redshift distributions of both ABHs and PBHs peak at z ∼ 2, the distribution of ABHs declines rapidly and ends at z ∼ 6. For the distribution of PBHs, it shows a more flattened shape and can reach z > 10.

With these numerous and high-redshift detections, we further calculated the results of gravitationally lensed signals detected by DECIGO. The differences between ABH and PBH scenarios are more visible. DECIGO will report ∼5 and ${ \mathcal O }(100)$ lensed GW signals from ABHs and PBHs, respectively. That is, the event rates of lensed GW signals from PBHs might become about 20 times larger than that of ABHs. Meanwhile, the redshift distributions of lensed GW signals from PBHs will become more flattened because of the magnification effect. Moreover, the difference between ABH and PBH scenarios in the high-redshift range will thus be more significant.

B-DECIGO will also be able to provide differences both in the event rates and the distributions, though the differences are not as significant as DECIGO due to its lower sensitivity. The results also become more distinguishable in the lensed events. With the help of the gravitational lensing effect, we expect to explore more distant coalescences of BBHs by (B-)DECIGO and other future detectors.

We would like to thank Shaoqi Hou and Marek Biesiada for their helpful discussions. This work was supported by the National Key Research and Development Program of China grant No. 2021YFC2203001, the National Natural Science Foundation of China under grant Nos. 12275021, 11920101003, 11722324, 11603003, 11633001, and U1831122, the science research grants from the China Manned Space Project with No. CMS-CSST-2021-B11, the Strategic Priority Research Program of the Chinese Academy of Sciences, grant No. XDB23040100, and the Interdiscipline Research Funds of Beijing Normal University.