# The merger rate of primordial black hole binaries as a probe of Hubble
parameter

Qianhang Ding<EMAIL_ADDRESS>Cosmology, Gravity and Astroparticle Physics
Group, Center for Theoretical Physics of the Universe,
Institute for Basic Science, Daejeon 34126, Republic of Korea

###### Abstract

We propose that the merger rate of primordial black hole (PBH) binaries can be
a probe of Hubble parameter by constraining PBH mass function in the
redshifted mass distribution of PBH binaries. In next-generation gravitational
wave (GW) detectors, the GWs from PBH binaries would be detected at high
redshifts, which gives their redshifted mass and luminosity distances. From a
number of detected events, the redshifted mass distribution of PBH binaries
can be statistically obtained, and it depends on PBH mass function and
redshift distribution of detected PBH binaries. The PBH mass function can be
inversely solved by applying the gradient descent method in the relation
between redshifted mass distribution and redshift distribution. However, the
construction of redshift distribution requires an assumed Hubble parameter in
a background cosmology to extract redshift from luminosity distances, which
causes solved PBH mass function also depends on assumed Hubble parameter. To
determine the Hubble parameter, the merger rate of PBH binaries constrains on
this Hubble parameter-dependent PBH mass function by comparing calculated
merger rate distribution with observed one, and the best-fit result produces
an approximate mass distribution of the physical PBH mass function and pins
down the Hubble parameter.

##  I Introduction

The value of Hubble parameter is of central importance in modern cosmology.
However, cosmological observations report an inconsistency in the measurement
value of Hubble parameter. The local measurement of Hubble parameter gives a
local value of $H_{0}=73.04\pm
1.04\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$ [1] (also see [2] for a
review), which is $5\sigma$ larger than the measurement value of Hubble
parameter $H_{0}=67.4\pm 0.5\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$
from Planck 2018 result under Lambda-cold dark matter ($\Lambda$CDM) model
[3]. This is so-called Hubble tension and its $5\sigma$ inconsistency has
become a cosmological crisis in $\Lambda$CDM cosmology [4, 5, 6].

To reveal Hubble tension and toward a concordance cosmological model, the
dynamical evolution of Hubble parameter from the early universe to the late
universe should be well measured. Currently, for early universe probes, cosmic
microwave background (CMB) provides a Hubble parameter measurement at redshift
$z\simeq 1100$ by fitting CMB power spectrum [3]. For late universe probes,
Type Ia supernovae (SNe) provide a Hubble parameter measurement at redshift
$z<\mathcal{O}(1)$ by constraining redshift-luminosity distance relation [7,
8], which works as standard candles. And baryon acoustic oscillation (BAO) can
measure Hubble parameter at redshift $z\sim\mathcal{O}(1)$ by calibrating
redshift-angular diameter distance relation [9, 10, 11], which is so-called
standard rulers. In the intermediate epoch between early universe and late
universe, gravitational waves (GWs) from compact binaries and their
electromagnetic counterparts can be a potential Hubble parameter probe up to
redshift $z\sim\mathcal{O}(10)$ as a standard siren in next-generation GW
detectors [12, 13, 14, 15]. However, Hubble parameter measurement between
redshift $z\sim\mathcal{O}(10)$ and $z\sim\mathcal{O}(1000)$ is still not
promising, due to lack of enough cosmological signals in the dark age [16],
which causes difficulty in constructing redshift-distance relation to extract
Hubble parameter.

The possible existence of primordial black hole (PBH) binaries changes the
story. PBHs are hypothetical astrophysical objects, which were born from
primordial perturbations in highly overdense regions by gravitational collapse
[17, 18, 19]. The formation of PBH binaries occurs in early universe epoch,
like radiation-domination era or matter-radiation equality, depends on PBH
mass and abundance [20, 21]. The GWs emitted from PBH binaries could propagate
cosmological distance from the early epoch to the present, where luminosity
distance is encoded in their GW waveform. Although redshift information is
still missing, due to no electromagnetic counterpart of PBH binaries in the
dark age [16], and causes difficulty in breaking the mass-redshift degeneracy
in GW waveform [12]. The unique statistic properties of PBH binaries open
another window in cosmic probes, the statistic distribution of PBH binaries,
such as PBH mass function and probability distribution on orbital parameters,
can provide an additional constraint on degenerate parameters to help extract
their intrinsic values, for instance, studying the evolution of probability
distribution on semi-major axis and eccentricity of PBH binaries from the
initial state to a later state can construct a redshift-time relation, which
constrains the Hubble parameter, this is one kind of standard timers [22, 23].
With possible existence of PBH binaries in cosmic history, their unique
statistic properties can help PBH binaries become a potential probe in
tracking cosmic evolution.

In this paper, we focus on using the merger rate of PBH binaries to put a
constraint on PBH mass function, which leads to a measurement of Hubble
parameter. From a number of detected PBH binaries at different redshifts, the
redshifted mass of PBHs and luminosity distance can be extracted from the GW
waveform of each PBH binary. We can statistically obtain a redshifted mass
distribution of PBH binaries from extracted redshifted PBH mass. Meanwhile,
cosmic redshifts can be solved from obtained luminosity distances under an
assumed Hubble parameter in a given cosmology such as $\Lambda$CDM model. Then
a redshift distribution of detected PBH binaries is statistically constructed
from solved redshifts. Combine the redshifted PBH mass distribution and
redshift distribution, a Hubble parameter-dependent PBH mass function can be
solved. To determine the Hubble parameter, the merger rate of PBH binaries put
the constraint on PBH mass function [24], we can calculate a merger rate
distribution of PBH binaries based on derived PBH mass function and compared
it with observed one. The fitting result of merger rate distribution of PBH
binaries depends on this Hubble parameter-dependent PBH mass function, when
solved PBH mass function is approaching its physical mass distribution, this
comparison can approach the best-fit result, which breaks the degeneracy
between PBH mass function and redshift distribution, and pins down the correct
Hubble parameter.

Although, we have not observed any PBH binaries so far, there still existing
some GW events which may hint the existence of PBH binaries, e.g., GW190425
and GW190814 show a companion mass of compact binary is smaller than
$3\,M_{\odot}$ [25, 26], which could be candidates of PBHs [27], GW190521
reports a compact binary event with mass in astrophysical BH mass gap [28],
which can be explained in PBH scenarios [27, 29]. In the next-generation GW
detectors, such as Laser Interferometer Space Antenna (LISA) [30] and Big Bang
Observer (BBO) [31], can detect GWs at redshift $z>20$, where astrophysical
BHs have not formed, therefore the detection of high redshift compact binaries
could be a smoking gun for the existence of PBH binaries [32, 33].

This paper is organized as follows. In Sec. II, we construct the connection
among redshifted mass distribution of PBH binaries, PBH mass function and
redshift distribution of detected PBH binaries. In Sec. III, we apply the
gradient descent method in redshifted mass distribution of PBH binaries to
inversely solve PBH mass function. In Sec. IV, we calculate the dependence of
redshift distribution of detected PBH binaries on Hubble parameter, which
leads to a Hubble parameter-dependent PBH mass function in the inverse
process. In Sec. V, we use the merger rate distribution of PBH binaries to
constrain the Hubble parameter-dependent PBH mass function, which works as a
Hubble parameter probe. In Sec. VI, the conclusion and discussions on this
Hubble parameter probe are given.

##  II Redshifted mass distribution of PBH binaries

In next-generation GW detectors, the detectable redshift range of BH binaries
would be larger than $20$, which provides an observational window for PBH
binaries [33]. From GW waveform of PBH binaries, key properties like
redshifted chirp mass $\mathcal{M}_{z}$, mass ratio $q$, luminosity distance
$d_{L}$, etc., can be extracted. Then an observed redshifted mass distribution
of PBH binaries $P(m_{1}^{z},m_{2}^{z})$ can be statistically obtained, where
$m_{1}^{z}$ and $m_{2}^{z}$ are the redshifted mass of PBHs in binaries, that
can be expressed by observable quantities as follows,

$\displaystyle
m_{1}^{z}=\mathcal{M}_{z}q^{2/5}(1+q)^{1/5}~{},~{}~{}~{}m_{2}^{z}=\mathcal{M}_{z}\frac{(1+q)^{1/5}}{q^{3/5}}~{}.$
(1)

The observed redshifted mass distribution of PBH binaries in GW detectors
depends on PBH mass function, sensitivity of GW detectors, and redshift
distribution of detected PBH binaries. It can be calculated from the
cumulative distribution of redshifted PBH mass $C(m_{1}^{z},m_{2}^{z})$, which
can be expressed as follows,

$\displaystyle C(m_{1}^{z},m_{2}^{z})$
$\displaystyle=\int_{0}^{m_{1}^{z}}\int_{0}^{m_{2}^{z}}P(m_{a},m_{b})\,dm_{a}\,dm_{b}$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{\frac{m_{1}^{z}}{1+z}}\int_{0}^{\frac{m_{2}^{z}}{1+z}}n(m_{1})n(m_{2})W(m_{1},m_{2};z)p(z)\,dm_{1}\,dm_{2}\,dz~{}.$
(2)

Here, $n(m)$ is the PBH mass function, $W(m_{1},m_{2};z)$ is the detectable
window function of PBH binaries with mass $m_{1}$ and $m_{2}$ at redshift $z$,
and $p(z)$ is the redshift distribution of detected PBH binaries. From Eq.
(II), we can obtain the redshifted mass distribution of PBH binaries by
derivating $C(m_{1}^{z},m_{2}^{z})$ with respect to redshifted mass
$m_{1}^{z}$ and $m_{2}^{z}$ as

$\displaystyle P(m_{1}^{z},m_{2}^{z})$
$\displaystyle=\frac{\partial^{2}C(m_{1}^{z},m_{2}^{z})}{\partial
m_{1}^{z}\,\partial m_{2}^{z}}$
$\displaystyle=\int_{0}^{\infty}n(\frac{m_{1}^{z}}{1+z})n(\frac{m_{2}^{z}}{1+z})W(\frac{m_{1}^{z}}{1+z},\frac{m_{2}^{z}}{1+z};z)\frac{p(z)}{(1+z)^{2}}\,dz~{}.$
(3)

In order to connect redshifted mass distribution of PBH binaries
$P(m_{1}^{z},m_{2}^{z})$ with PBH mass function $n(m)$, we need to figure out
the detectable window function $W(m_{1},m_{2};z)$ and redshift distribution
$p(z)$. Then the relation between $P(m_{1}^{z},m_{2}^{z})$ and $n(m)$ can be
obtained and further can be applied in an inverse process to solve $n(m)$ from
$P(m_{1}^{z},m_{2}^{z})$.

Window function $W(m_{1},m_{2};z)$ denotes the detection probability of the
observed PBH binaries among all existing PBH binaries with mass $m_{1}$ and
$m_{2}$ at redshift $z$, it comes from the selection effect due to limited
sensitivity and frequency band of GW detectors, which causes only a part of
PBH binaries with suitable orbital parameters, semi-major axis $a$ and
eccentricity $e$, can be observed by given GW detectors [34, 35, 36]. The
window function can be expressed as follows,

$\displaystyle
W(m_{1},m_{2};z)=\frac{N_{\mathrm{obs}}(m_{1},m_{2};z)}{N_{\mathrm{tot}}(m_{1},m_{2};z)}~{},$
(4)

where $N_{\mathrm{obs}}(m_{1},m_{2};z)$ and $N_{\mathrm{tot}}(m_{1},m_{2};z)$
are the number of observed and total PBH binaries with mass $m_{1}$ and
$m_{2}$ at redshift $z$, respectively. The detection of a PBH binary requires
a signal-to-noise ratio ($\mathrm{SNR}$) should exceed a threshold for a high
confidence level, where we set the threshold value of $8$ [37], which makes
sure the detection probability is $>95\%$, and a false alarm probability
$<0.1\%$. The SNR can be calculated as shown in [37, 38],

$\displaystyle\mathrm{SNR}=\sqrt{4\int_{f_{\mathrm{min}}}^{f_{\mathrm{max}}}\frac{|\tilde{h}(f)|^{2}}{S_{n}(f)}df}~{},$
(5)

where $[f_{\mathrm{min}},f_{\mathrm{max}}]$ is the observational frequency
range. $S_{n}(f)$ is the noise strain of the detector, and $S_{n}(f)$ for each
detector can be found from Fig. A2 in [39]. In following parts, we adopt the
$S_{n}(f)$ of BBO in calculations. $\tilde{h}(f)$ is the Fourier
transformation of GW amplitude $h(t)$, which can be expressed as following
form under the stationary approximation [40],

$\displaystyle\tilde{h}(f)=\sqrt{\frac{5}{24}}\frac{(G\mathcal{M}_{c}(1+z))^{5/6}}{\pi^{2/3}c^{3/2}d_{L}(z)}f^{-7/6}~{},$
(6)

where $G$ is the Newton’s constant, $c$ is the speed of light,
$\mathcal{M}_{c}=\mathcal{M}_{z}/1+z$ is the intrinsic chirp mass of PBH
binaries, and $d_{L}$ is the luminosity distance to the PBH binary.

In order to evaluate Eq. (4), we need integrate the PBH binaries probability
distribution on semi-major axis and eccentricity $P(a,e;z)$ over detectable
parameter region as follows,

$\displaystyle
W(m_{1},m_{2};z)=\int_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}\int_{e_{\mathrm{min}}}^{e_{\mathrm{max}}}P(a,e;z)dade=\int_{a_{\mathrm{min}}^{\mathrm{ini}}}^{a_{\mathrm{max}}^{\mathrm{ini}}}\int_{e_{\mathrm{min}}^{\mathrm{ini}}}^{e_{\mathrm{max}}^{\mathrm{ini}}}P(a,e;z_{\mathrm{ini}})dade~{},$
(7)

where $a_{\mathrm{min}}$ ($e_{\mathrm{min}}$) and $a_{\mathrm{max}}$
($e_{\mathrm{max}}$) are minimal and maximal detectable orbit parameters,
respectively. Their values can be numerically solved from the detection
requirement $\mathrm{SNR}>8$ in Eqs. (5) and (6). We mainly focus on the
detection of small eccentricity orbits ($e\sim 0$) due to lack of suitable
high-eccentricity GW waveform template [41]. The superscript $\mathrm{ini}$ on
$a$ and $e$ denotes the initial value of semi-major axis and eccentricity when
PBH binaries form. The initial value of $a$ and $e$ can be solved by the
formula in [42] as follows,

$\displaystyle\frac{da}{dt}=-\frac{64}{5}\frac{G^{3}m_{1}m_{2}(m_{1}+m_{2})}{c^{5}a^{3}(1-e^{2})^{7/2}}\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right)~{},~{}~{}~{}\frac{de}{dt}=-\frac{304}{15}\frac{G^{3}m_{1}m_{2}(m_{1}+m_{2})}{c^{5}a^{4}(1-e^{2})^{5/2}}e\left(1+\frac{121}{304}e^{2}\right)~{}.$
(8)

Then we solve Eq. (8) backward in time from the observable parameters $(a,e)$
at redshift $z$ to the initial parameters
$(a^{\mathrm{ini}},e^{\mathrm{ini}})$ at PBH binary formation redshift
$z_{\mathrm{ini}}$ in Eq. (25). The probability distribution on semi-major
axis and eccentricity at initial redshift $z_{\mathrm{ini}}$ for PBH binaries
with mass $m_{1}$ and $m_{2}$ can be approximated as [20, 21]

$\displaystyle
P(a,e,z_{\mathrm{ini}})=\frac{3}{4}f_{b}^{3/2}\frac{a^{1/2}}{\bar{x}^{3/2}}\frac{e}{(1-e^{2})^{3/2}}~{},$
(9)

where $f_{b}$ is the energy density fraction of PBH binaries with mass $m_{1}$
and $m_{2}$ in the dark matter and $\bar{x}$ is the physical mean separation
of PBH binaries with mass $m_{1}$ and $m_{2}$ at matter-radiation equality,
see Eq. (27) and (28). Combine Eqs. (4 – 9), window function
$W(m_{1},m_{2};z)$ for PBH binaries with mass $m_{1}$ and $m_{2}$ at redshift
$z$ can be numerically calculated.

Redshift distribution $p(z)$ describes the distribution of detected PBH
binaries at different redshifts, it depends on detectable event rate of PBH
binaries at different redshifts, comoving volume as a function of redshift,
and selection effect due to the redshift detection efficiency in GW detectors.
In order to obtain this redshift distribution $p(z)$, we can statistically
count detected PBH binary events at different redshifts in observations.
Theoretically, we use normalized redshift distribution of PBH merger rate
population as an approximation as shown in [43],

$\displaystyle p(z)\propto\frac{\dot{n}(z)}{1+z}\frac{dV_{c}}{dz}~{},$ (10)

where we set a detectable redshift range $20\leq z\leq 100$ as an
approximation of the selection effect on redshift. $\dot{n}(z)$ is the
volumetric merger rate density of PBHs, which is power-law dependent on the
age of the universe $t(z)$ as $\dot{n}(z)\propto(t(z)/t_{0})^{-34/37}$ [44],
where $t_{0}$ is the current age of the universe. $1+z$ term comes from the
cosmic time dilation. $dV_{c}/dz$ is the differential comoving volume, which
can be calculated as $dV_{c}/dz=4\pi D_{H}D_{M}^{2}/E_{z}(z)$ in flat
universe, where $D_{H}=c/H_{0}$ is Hubble distance, $D_{M}$ is transverse
comoving distance, and
$E_{z}(z)=\sqrt{\Omega_{\mathrm{M}}(1+z)^{3}+\Omega_{\Lambda}}$ in flat
$\Lambda$CDM model [45].

Taking account Eqs. (7) and (10) in Eq. (II), and we consider two typical PBH
mass functions, log-normal distribution [46] and power-law distribution [47]

$\displaystyle n_{\mathrm{LN}}(m)=\frac{1}{\sqrt{2\pi}\sigma
m}\exp\left(-\frac{\log^{2}(m/m_{c})}{2\sigma^{2}}\right)~{},~{}~{}~{}n_{\mathrm{PL}}(m)=\frac{\alpha-1}{M}\left(\frac{m}{M}\right)^{-\alpha}~{},$
(11)

where $\sigma$ denotes the mass width of log-normal distribution, $m_{c}$ is
the critical mass in log-normal distribution, and $\alpha$ is the power index
in power-law distribution and $M$ is the minimal mass in power-law
distribution. Then we can obtain the redshifted mass distribution of PBH
binaries as shown in Fig. 1.

Figure 1:  Left panel: The normalized redshifted mass distribution of PBH
binaries for a log-normal PBH mass function with $m_{c}=30\,M_{\odot}$ and
$\sigma=1$. Right panel: The normalized redshifted mass distribution of PBH
binaries for a power-law PBH mass function with $\alpha=1.5$ and
$M=2M_{\odot}$. The assumed cosmology for both panels is a flat $\Lambda$CDM
model with Planck 2018 parameters $\Omega_{\mathrm{M}}=0.315$,
$\Omega_{\Lambda}=0.685$ and
$H_{0}=67.4\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$ [3], and the
redshifted mass $m_{1}^{z}$ and $m_{2}^{z}$ is set to follow $m_{1}^{z}\leq
m_{2}^{z}$.

It shows the redshifted mass distribution of PBH binaries strongly depends on
PBH mass functions, and we can find the redshifted mass distribution with a
log-normal PBH mass function spreads broadly, while a power-law PBH mass
function produces a narrow redshifted mass distribution that is centered
around a small value of redshifted mass. This indicates the redshifted mass
distribution of PBH binaries could be a potential tool in distinguishing
different PBH mass functions.

##  III PBH mass function from inverse redshifted mass distribution

Theoretically, we can calculate the redshifted mass distribution of PBH
binaries based on PBH mass function $n(m)$, detector-dependent window function
$W(m_{1},m_{2};z)$ and the redshift distribution of detected PBH binaries
$p(z)$. Observationally, we can only detect the redshifted mass distribution
of PBH binaries. The relation between redshifted mass distribution of PBH
binaries and PBH mass function follows Eq. (II) as

$\displaystyle
P_{O}(m_{1}^{z},m_{2}^{z})=\int_{0}^{\infty}n_{p}(\frac{m_{1}^{z}}{1+z})n_{p}(\frac{m_{2}^{z}}{1+z})W(\frac{m_{1}^{z}}{1+z},\frac{m_{2}^{z}}{1+z};z)\frac{p(z)}{(1+z)^{2}}\,dz~{},$
(12)

where $P_{O}(m_{1}^{z},m_{2}^{z})$ with subscript $O$ denotes the observed
redshifted mass distribution and $n_{p}(m)$ describes the physical mass
distribution function of PBHs. Apart from $n_{p}(m)$,
$P_{O}(m_{1}^{z},m_{2}^{z})$ only depends on detector-dependent window
function $W(m_{1},m_{2};z)$ and redshift distribution of detected PBH binaries
$p(z)$. For a given GW detectors, the selection effect in window function is
fixed. Also, redshift distribution can be extracted from observed PBH
binaries, where GW waveform of PBH binaries provides the luminosity distances
between PBH binaries and the earth, if we know the physical cosmological model
with a fixed Hubble parameter, redshift can be determined from obtained
luminosity distance, then redshift distribution of detected PBH binaries can
be statistically obtained. After knowing the contribution of
$W(m_{1},m_{2};z)$ and $p(z)$ in Eq. (12), the physical PBH mass function
$n_{p}(m)$ can be inversely solved from observed redshifted mass distribution
$P_{O}(m_{1}^{z},m_{2}^{z})$.

To obtain PBH mass function from inverse redshifted mass distribution, we
apply the gradient descent method (see [48] for detailed introductions) in Eq.
(12). There are several steps in the gradient descent method as follows,

  * 1.

Give an initial mass distribution as an initial condition of PBH mass
function.

  * 2.

Calculate theoretical redshifted mass distribution and compare it with
observed one to obtain the error function.

  * 3.

Update PBH mass function by calculating the gradient of error function with
respect to the variation of PBH mass function.

  * 4.

Iterate PBH mass function until finding the minimal value of error function,
which is a best-fit approximation of physical PBH mass function.

In order to decrease time complexity of this algorithm, we start from
discretizing PBH mass function $n(m)$ to a PBH mass function vector
$\mathbf{n}$, where the component of vector follows
$\mathbf{n}_{\mathrm{i}}=n(m_{i})$ and a set of $\\{m_{i}|1\leq i\leq N\\}$ is
$N$ mass points taken from the possible PBH mass range. Then we provide an
initial value of this PBH mass function vector $\mathbf{n}$. This initial
value could be an arbitrary mass function, however, choosing a more realistic
and reasonable initial mass function could effectively decrease the time
complexity of the gradient descent method. After giving the initial value of
$\mathbf{n}$, an initial PBH mass function can be approximately constructed by
interpolating the initial PBH mass function vector.

For a given assumed PBH mass function, we can combine it with window function
$W(m_{1},m_{2};z)$ and redshift distribution $p(z)$ in Eq. (12) to calculate a
theoretical redshifted mass distribution. Then we can define an error function
$E(\mathbf{n})$ by comparing theoretical value with the observed value of
redshifted mass distribution at a selected set of PBH binary mass pairs as
follows (see [49] for more details),

$\displaystyle E(\mathbf{n})\equiv\sqrt{\frac{\sum_{1\leq i\leq
j,\,j=1}^{j=N}[P_{T}(m_{i}^{z},m_{j}^{z})-P_{O}(m_{i}^{z},m_{j}^{z})]^{2}}{(1+N)N/2}}~{},$
(13)

where $P_{T}(m_{1}^{z},m_{2}^{z})$ and $P_{O}(m_{1}^{z},m_{2}^{z})$ are
theoretical and observed redshifted mass distribution, respectively, and $N$
is the number of selected redshifted mass points for $m_{1}^{z}$ and
$m_{2}^{z}$.

The error function $E(\mathbf{n})$ only depends on the assumed value of PBH
mass function vector $\mathbf{n}$, so any further derivation of $\mathbf{n}$
from the physical PBH mass function $n_{p}(m)$ would enhance the difference
between $P_{T}(m_{1}^{z},m_{2}^{z})$ and $P_{O}(m_{1}^{z},m_{2}^{z})$ and
increase the value of error function. It indicates that, in order to minimize
the difference between assumed PBH mass function vector and the physical PBH
mass function, minimizing the error function $E(\mathbf{n})$ could be an
indicator to make $\mathbf{n}$ approaching the physical PBH mass function. To
minimize the error function $E(\mathbf{n})$, we use the following iteration
equation to update $\mathbf{n}$,

$\displaystyle
n_{\mathrm{k+1}}(m_{i})=n_{\mathrm{k}}(m_{i})-\gamma\frac{\partial
E(\mathbf{n_{\mathrm{k}}})}{\partial n_{\mathrm{k}}(m_{i})}~{},$ (14)

where $\mathbf{n}_{\mathrm{k}}$ is the $k\,$th iterated PBH mass vector and
$n_{\mathrm{k}}(m_{i})$ is its $i\,$th component. $\gamma$ is so-called
learning rate in the gradient descent method, which describes the length of
updated step in each iteration. With the update of $\mathbf{n}$, error
function $E(\mathbf{n})$ is approaching its minimal value during iterations,
the gradient of error function $\nabla E(\mathbf{n})$ is approaching zero
vector, and the updated PBH mass function vector
$\mathbf{n}_{\mathrm{k+1}}\approx\mathbf{n}_{\mathrm{k}}$. In this case, the
iterated $\mathbf{n}$ can be viewed as a discretized approximation of the
physical PBH mass function $n_{p}(m)$.

Figure 2:  Upper panel: The data points are solved PBH mass function vector
via the gradient descent method. The solid curve is the assumed log-normal PBH
mass function with $m_{c}=30\,M_{\odot}$ and $\sigma=1$. Lower panel: The
relative error between solved PBH mass function vector and assumed PBH mass
function. The background cosmology is a flat $\Lambda$CDM model with Planck
2018 parameters [3].

In Fig. 2, we use the gradient descent method to inversely solve PBH mass
function from a redshifted mass distribution of PBH binaries. This redshifted
mass distribution is calculated under an assumed log-normal PBH mass function
with $m_{c}=30\,M_{\odot}$ and $\sigma=1$, which is shown as the solid curve,
and the background cosmology is assumed as the flat $\Lambda$CDM model with
Planck 2018 parameters [3]. The data points are solved discrete PBH mass
function vector, where we discretize PBH mass function into $50$ data points.
In this iteration process, we obtain the minimal value in error function after
$\mathcal{O}(50)$ times iterations, and we can find the solved PBH mass
function vector can fit assumed PBH mass function to precision
$\mathcal{O}(10^{-2})$. However, the relative error at small mass and large
mass is relatively large compared with the relative error at medium mass
range. This is because the small mass tail and large mass tail only occupy a
small part in PBH mass function and further provide a small contribution in
derived error function. This causes the constraints from derivating the error
function on small mass and large mass range is relatively weak in the gradient
descent method.

##  IV Hubble parameter–dependent PBH mass function

In extracting PBH mass function via inverse process from observed redshifted
PBH mass function, the redshift distribution of detected PBH binaries $p(z)$
plays an important role in determining obtained PBH mass function. However,
the observables we can extract from PBH binaries in the GW channel are their
redshifted mass $m_{z}$ and luminosity distance $d_{L}$, which cannot produce
a redshift distribution of detected PBH binaries due to lack of redshift
information of PBH binaries to break mass-redshift degeneracy [12, 50].

Although the redshift is degenerate with intrinsic mass in the observables of
PBH binaries, we can still break this degeneracy by fixing some terms in the
observables. As we can find that the luminosity distance $d_{L}(z)$ depends on
the redshift of PBH binary and the cosmological model, if we fix an assumed
cosmological model and its cosmological parameters, the redshift can be
extracted under such an assumption. Then, a model-dependent redshift
distribution $p(z)$ can be constructed from the observables. The redshift and
cosmological model dependence of luminosity distance in a flat cosmology [45]
is

$\displaystyle
d_{L}(z)=(1+z)\int_{0}^{z}\frac{c}{H(z^{\prime})}dz^{\prime}=\frac{1+z}{H_{0}}\int_{0}^{z}\frac{c}{E_{z}(z^{\prime})}dz^{\prime}~{},$
(15)

where $H(z)$ is the Hubble parameter at different redshifts, $H_{0}$ is the
present Hubble parameter, $E_{z}(z)=H(z)/H_{0}$ describes the evolution of
Hubble parameter with respect to redshift. For any given cosmology with
corresponding parameters, $E_{z}(z)$ would be fixed, then luminosity distance
become a function of redshift and Hubble parameter $d_{L}\equiv
d_{L}(z,H_{0})$, we can also generalize this Hubble parameter dependence as
$d_{L}\equiv d_{L}(z,H)$, where $H$ is the Hubble parameter at redshift in the
range of $[0,z]$. In following discussion, we mainly focus on the dependence
of present Hubble parameter $H_{0}$.

Then redshift becomes a function of the luminosity distance and Hubble
parameter $z(d_{L},H_{0})$, and can be extracted from the luminosity distance
by fixing a Hubble parameter. From a number of detected PBH binaries, we have
a set of luminosity distance $\\{d_{L}(z_{i},H_{0})|1\leq i\leq N\\}$, where
$N$ is the number of detected PBH binaries. For any given Hubble parameter
$H_{0}$, a set of redshift can be obtained as
$\\{z_{i}=z(d_{L}^{\,i},H_{0})|1\leq i\leq N\\}$, which constructs a Hubble
parameter-dependent redshift distribution $p(z;H_{0})$. The dependence of
redshift distribution on Hubble parameter behaves as follows,

$\displaystyle
p(z;\widetilde{H}_{0})=\frac{dN}{d\tilde{z}}=\frac{dN}{dz}\frac{dz}{d\tilde{z}}=p(z;H_{0})\frac{dz}{d\tilde{z}}~{},$
(16)

where $\widetilde{H}_{0}$ and $\tilde{z}$ are assumed Hubble parameter and its
corresponding derived redshift, respectively. The differential relation of
redshift $dz/d\tilde{z}$ can be found by derivating the relation
$d_{L}(z,H_{0})=d_{L}(\tilde{z},\widetilde{H}_{0})$, which gives,

$\displaystyle\frac{dz}{d\tilde{z}}=\frac{\widetilde{D}_{c}+\widetilde{D}_{H}(1+\tilde{z})/E_{z}(\tilde{z})}{D_{c}+D_{H}(1+z)/E_{z}(z)}~{}.$
(17)

Here, $D_{c}=D_{H}\int_{0}^{z}dz^{\prime}/E_{z}(z)$ is the line-of-sight
comoving distance [45] and $\widetilde{D}_{c}$ is the line-of-sight comoving
distance with an assumed Hubble parameter $\widetilde{H}_{0}$. Then, we can
obtain the relation between redshift distribution with assumed Hubble
parameter $p(z;\widetilde{H}_{0})$ and the physical redshift distribution
$p(z;H_{0})$ from Eq. (16).

In GW observation of PBH binaries, the intrinsic PBH mass cannot be determined
from their observable, redshifted mass $m_{z}$ and luminosity distance
$d_{L}$, however, if we can know the physical cosmological model and Hubble
parameter, the redshift can be determined from the luminosity distance, and
the intrinsic PBH mass can be pinned down from mass-redshift degeneracy as
follows,

$\displaystyle m=\frac{m_{z}}{1+z(d_{L};H_{0})}~{}.$ (18)

It shows that PBH mass inferred from derived redshift in luminosity distance
also depends on Hubble parameter, which indicates the PBH mass function should
also be Hubble parameter-dependent. This can be found as follows,

$\displaystyle
P_{O}(m_{1}^{z},m_{2}^{z})=\int_{0}^{\infty}\tilde{n}(\frac{m_{1}^{z}}{1+z})\tilde{n}(\frac{m_{2}^{z}}{1+z})W(\frac{m_{1}^{z}}{1+z},\frac{m_{2}^{z}}{1+z};z)\frac{p(z;\widetilde{H}_{0})}{(1+z)^{2}}\,dz~{},$
(19)

where $p(z;\widetilde{H}_{0})$ is a Hubble parameter-dependent redshift
distribution under the assumed Hubble parameter and $\tilde{n}(m)$ is the
unknown PBH mass function under the assumed Hubble parameter, then we apply
the gradient descent method in Eq. (19) to inversely solve PBH mass function
from observed redshifted mass distribution of PBH binaries as discussed in
Sec. III, and obtained PBH mass function $\tilde{n}(m)$ is also a Hubble
parameter-dependent PBH mass function $n(m;\widetilde{H}_{0})$.

Figure 3:  The Hubble parameter-dependent PBH mass functions and assumed
physical PBH mass function. The solid curve is assumed physical PBH mass
function with Hubble parameter
$H_{0}=67.4\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$, this physical
PBH mass function is a log-normal type with $m_{c}=30\,M_{\odot}$ and
$\sigma=1$. The dashed curves are inversely solved normalized PBH mass
functions with different assumed Hubble parameters
$h=H_{0}/(100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1})$ as
annotations. The background cosmology is a flat $\Lambda$CDM model with Planck
2018 parameters [3].

In Fig. 3, we can find the solved PBH mass function is Hubble parameter-
dependent as we have discussed above. A larger assumed Hubble parameter would
cause a lighter mass distribution function, this can be intuitively understood
as a larger value of Hubble parameter gives a larger value of solved redshift
in Eq. (15), and this larger redshift would cause a smaller value of intrinsic
mass in Eq. (18). This figure shows that if we can have a constraint on the
PBH mass function (we will discuss on it in Sec. V), we can obtain a
measurement on the Hubble parameter. In such a measurement, measurement
uncertainty is the key issue, which can be estimated based on the measurement
precision of redshifted mass, intrinsic mass and luminosity distance in Eqs.
(15) and (18) as $\delta H_{0}\sim\delta m_{z}-\delta m-\delta d_{L}$, where
$\delta$ denotes the relative error. Currently, the measurement uncertainty of
$m_{z}$ and $d_{L}$ is around $\mathcal{O}(10^{-1})$ in LIGO and Virgo (see
TABLE VI of [51]), if the measurement precision of intrinsic mass can achieve
the same order as redshifted mass, this would causes the $\delta
H_{0}\sim\mathcal{O}(10^{-1})$. To further improve the measurement precision
of the Hubble parameter, so that we can reveal the Hubble tension, the
measurement uncertainty of Hubble parameter should achieve $\delta H_{0}<0.1$,
this requires a large number of high SNR GW events to put a strong constraint
on PBH mass function. With the development of GW detectors, especially space-
based low frequency GW detectors, the measurement precision of extracted
parameters can be further improved, which provides a higher measurement
precision on Hubble parameter in this method and helps us understand the
Hubble tension.

##  V Merger rate of PBH binaries as a probe of Hubble parameter

The measurement of Hubble parameter requires a connection between different
observables, where the scale of the universe and cosmic dynamics should be
encoded [52], such as redshift-luminosity distance relation in standard
candles [53] and standard sirens [12, 13], redshift-cosmic time relation in
cosmic chronometers [54] and standard timers [23, 22]. In Eq. (19), such a
relation is also connected, the redshift is encoded in observed
$P_{O}(m_{1}^{z},m_{2}^{z})$, while cosmic dynamics is hidden in inferred
$p(z;H_{0})$, which requires an assumed cosmology and Hubble parameter to
solve redshifts from the observed luminosity distances. In connecting the
relation between $P_{O}(m_{1}^{z},m_{2}^{z})$ and $p(z;H_{0})$, we need to
know the physical PBH mass function. However, the PBH mass function is
currently unknown, this causes the degeneracy between PBH mass function and
redshift distribution, that is for an observed $P_{O}(m_{1}^{z},m_{2}^{z})$,
different assumed Hubble parameters would produce various PBH mass functions
and redshift distributions. Theoretically, we can calculate PBH mass function
based on inflation models and PBH formation mechanisms, a lot of work has been
done to derive PBH mass function [47, 55, 56, 57, 58]. Observationally, we
have not proved the existence of the PBH, even though some GW events could be
potential clues for PBHs [25, 26, 28]. To determine the Hubble parameter from
above relation, we need to construct a relation bwtween an observable and PBH
mass function and use this observable to constrain the PBH mass function. This
relation must be independent on relation in Eq. (19), then it can be used to
break the degeneracy between PBH mass function and redshift distribution and
pin down the Hubble parameter.

The introduction of merger rate of PBH binaries $R$ changes the story, which
provides another connection between an observable PBH merger rate and unknown
PBH mass function. Combine this observable PBH merger rate with the relation
in Eq. (19), the PBH mass function can be constrained and the Hubble parameter
can be pinned down as a probe of Hubble parameter. It works as follows,

  * 1.

Give an assumed Hubble parameter $\widetilde{H}_{0}$ in a given cosmology.

  * 2.

Use obtained $d_{L}$ from detected PBH binaries with $\widetilde{H}_{0}$ to
construct $p(z;\widetilde{H}_{0})$.

  * 3.

Combine $p(z;\widetilde{H}_{0})$ with observed $P_{O}(m_{1}^{z},m_{2}^{z})$ to
solve $n(m;\widetilde{H}_{0})$ via the gradient descent method.

  * 4.

Calculate $R$ based on $n(m;\widetilde{H}_{0})$ and compare it with observed
PBH merger rate.

  * 5.

Vary $\widetilde{H}_{0}$ until obtained $R$ matches the observed result, which
pins down the correct Hubble parameter.

From the discussion in Sec. (II–IV), we can use the gradient descent method to
solve a Hubble parameter-dependent PBH mass function $n(m;\widetilde{H}_{0})$
from an observed redshifted mass distribution $P_{O}(m_{1}^{z},m_{2}^{z})$.
Then we can calculate PBH binary merger rate based on $n(m;\widetilde{H}_{0})$
and compare it with observed PBH merger rate $R_{O}$ to find the best-fit
Hubble parameter. To calculate the merger rate of PBH binaries, a number of
works have been done, see [20, 21, 59, 24, 60, 44] for details. We follow [24]
and use the merger rate of PBH binaries $R$ with an extended mass function as,

$\displaystyle
R_{ij}=\rho_{\mathrm{PBH}}\min\left(\frac{n(m_{i})}{m_{i}},\frac{n(m_{j})}{m_{j}}\right)\Delta_{m}\frac{dP}{dt}~{},$
(20)

where, $\rho_{\mathrm{PBH}}$ is the present PBH density, expressed as
$\rho_{\mathrm{PBH}}=3H_{0}^{2}\Omega_{\mathrm{DM}}f_{\mathrm{PBH}}/8\pi G$
and $f_{\mathrm{PBH}}$ is the energy density fraction of PBHs in dark matter.
$n(m_{i})$ and $\Delta_{m}$ are the discretized PBH mass function and mass
interval, respectively, which follow the normalization of probability as
$\sum_{m_{\mathrm{min}}}^{m_{\mathrm{max}}}n(m_{i})\Delta_{m}=1$. $dP/dt$ is
the probability distribution of merger time, see Appendix. A for detailed
derivation.

In Eq. (20), the merger rate of PBH binaries depends on PBH mass function
differently compared with the PBH mass function dependence in redshifted mass
distribution, which indicates the merger rate of PBH binaries can provide an
independent constraint on PBH mass function. As we have discussed above, the
solved PBH mass function from inversing the redshifted mass distribution
requires an assumed Hubble parameter, then we can use the solved Hubble
parameter-dependent PBH mass function $n(m;\widetilde{H}_{0})$ to calculate
the merger rate distribution of PBH binaries theoretically, which should also
be Hubble parameter-dependent $R(m_{i},m_{j};\widetilde{H}_{0})$. Then we use
calculated observed PBH mass function at a given redshift $z$, where the
observed merger rate is a redshift-diluted intrinsic merger rate, it follows
$R_{O}=R/(1+z)$, to compare with the physical observed merger rate
distribution. Only if we assumed the correct Hubble parameter, the theoretical
result can match the observed result, and this comparison can help pin down
the correct value of Hubble parameter, which is shown in the left panel of
Fig. 4. Also, we can use the merger rate population of PBH binaries to
constrain Hubble parameter-dependent PBH mass function, which can be
calculated as follows,

$\displaystyle\dot{N}=\int_{z_{\mathrm{min}}}^{z_{\mathrm{max}}}\frac{R}{1+z}\frac{dV_{c}}{dz}dz~{},$
(21)

where $\dot{N}$ denotes the merger rate population of PBH binaries between
redshift $z_{\mathrm{min}}$ and $z_{\mathrm{max}}$, $dV_{c}/dz$ is the
differential comoving volume as defined in Eq. (10). The comparison of merger
rate population of PBH binaries between assumed physical population
distribution and calculated population distribution is shown in the right
panel of Fig. 4.

Figure 4:  Left panel: The comparison of observed merger rate distribution of
PBH binaries at redshift $z=20$ between assumed physical merger rate and
calculated Hubble parameter-dependent merger rate. Right panel: The comparison
of merger rate population distribution of PBH binaries within redshift $z=20$
and $z=100$ between assumed physical merger rate population and calculated
Hubble parameter-dependent merger rate population. The label of horizontal
axis $m_{\mathrm{H}}$ denotes the mass of heavier PBH companion in PBH
binaries. The solid curves are assumed physical merger rate distribution and
merger rate population distribution, respectively. The assumed physical merger
rate and population are calculated based on a log-normal PBH mass function
with $m_{c}=30\,M_{\odot}$ and $\sigma=1$. The dashed curves are calculated
merger rate distribution and merger rate population distribution,
respectively. The Hubble parameters denoted in dashed curves are
$h=H_{0}/(100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1})$. The
background cosmology for both panels is a flat $\Lambda$CDM model with Planck
2018 parameters [3] and the energy density fraction of PBHs in dark matter
$f_{\mathrm{PBH}}$ is set as $f_{\mathrm{PBH}}=0.001$.

In Fig. 4, we calculate the distribution of merger rate and population of PBH
binaries with respect to the mass of heavier PBH companion, which follow as
$R(m_{\mathrm{H}})=\sum_{i=1}^{i=\mathrm{H}}R(m_{i},m_{\mathrm{H}})$ and
$\dot{N}(m_{\mathrm{H}})=\sum_{i=1}^{i=\mathrm{H}}\dot{N}(m_{i},m_{\mathrm{H}})$.
We can find that the calculated merger rate distribution and population
distribution depend on the assumed value of Hubble parameter, where a larger
value of Hubble parameter produces more merger rates of PBH binaries . This
can be intuitively understood as a large value of Hubble parameter causes the
intrinsic PBH mass tend to a smaller value as we have discussed in Sec. IV,
then a smaller intrinsic mass in PBH mass function would produce a larger
merger rate in Eq. (20). The difference in assumed value of Hubble parameter
would produce distinct merger rate distributions, which can be used to pin
down the correct Hubble parameter. As we can find in the both panels of Fig.
4, under an assumed flat $\Lambda$CDM model with
$H_{0}=67.4\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$, we apply
different values of Hubble parameter to inversely solve PBH mass function from
the redshifted mass distribution of PBH binaries and calculate their merger
rate distributions, only when Hubble parameter is approaching
$67.4\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$, we can obtain the
best-fit result. However, we still need to notice that this is a model-
dependent method to measure the Hubble parameter, various background
cosmologies would pin down different Hubble parameters. To achieve a better
performance in probing Hubble parameter, a combination of this probe with
other cosmic probes would pin down more precise Hubble parameter.

##  VI Conclusion and Discussions

To summarize, we propose that the merger rate of PBH binaries can constrain
the PBH mass function in the redshifted mass distribution of PBH binaries, and
such a constraint on PBH mass function can help break the degeneracy between
PBH mass function and redshift distribution and pin down the Hubble parameter.
In GW observations on PBH binaries, a redshifted mass distribution of PBH
binaries can be statistically obtained from redshifted mass of PBH binaries,
and this redshifted mass distribution depends on the PBH mass function and the
redshift distribution of detected PBH binaries, where the redshift
distribution is not an observable in GW observations. To construct redshift
distribution of detected PBH binaries, redshifts need to be extracted from
luminosity distances by assuming a value of the Hubble parameter in a given
cosmology, which produces a Hubble parameter-dependent redshift distribution.
In the connection among redshifted mass distribution, PBH mass function and
the Hubble parameter-dependent redshift distribution, Hubble parameter cannot
be extracted, due to the unknown PBH mass function. To achieve it, we use the
merger rate distribution of PBH binaries to provide an independent constraint
on PBH mass function, which breaks this degeneracy between PBH mass function
and redshift distribution and pins down the Hubble parameter.

In calculation, we focus on the redshifted mass distribution of PBH binaries
at redshifts within $20\leq z\leq 100$, where astrophysical BHs have not
formed and contributed, such a high redshift detection can be achieved in the
next-generation GW detectors [30, 31]. Then we can construct redshift
distribution from a series of observed luminosity distances under an assumed
Hubble parameter in a given cosmology (we use $\Lambda$CDM model in this
work). To obtain PBH mass function from redshifted mass distribution, we apply
the gradient descent method in the redshifted mass distribution to inversely
solve the PBH mass function. Also, this solved PBH mass function depends on
assumed Hubble parameter, due to the dependence of redshift distribution on
Hubble parameter. Based on solved Hubble parameter-dependent PBH mass
function, we calculate the merger rate distribution of PBH binaries
theoretically and compare it with observed result, where the compared
observational merger rate distribution is calculated based on a log-normal PBH
mass function with $m_{c}=30\,M_{\odot}$ and $\sigma=1$ and $\Lambda$CDM model
with Planck 2018 parameters [3]. To find the physical value of Hubble
parameter, we vary the value of assumed Hubble parameter until the calculated
merger rate distribution can match the observed result, where the best-fit
Hubble parameter can be an approximation of its physical value.

In above discussion, we mainly focus on the contribution of intrinsic PBH mass
function in the redshifted PBH mass function without considering the dynamical
evolution of PBH mass function due to accretion effect [61, 62, 63] and the
multiple mergers of PBH binaries [64]. Although some dynamical effect can be
ignored due to their negligible contribution [65], a realistic Hubble
parameter probe requires adding an extra contribution term in the physical PBH
mass function in Eq. (12). Also, the measurement precision of Hubble parameter
is the key issue in this method, its measurement uncertainty depends on the
measurement uncertainty of redshifted mass, luminosity distance and PBH merger
rate, that causes the precision of this method is not comparable with current
methods, like standard candles and CMB power spectrum, as we have discussed in
Sec. IV. For an applicable Hubble parameter measurement in this probe, it
requires a measurement uncertainty for extracted parameters from GWs of PBH
binaries, such as redshifted mass and luminosity distance, should achieve
around $\mathcal{O}(10^{-2})$ in next-generation GW detectors.

To measure Hubble parameter, a number of methods could be applied in PBH
binaries, such as gravitational-wave lensing [66, 67, 68], black hole
population [69, 70] and statistic mass distribution in compact binaries [71],
with these new methods and development of GW observations, PBH binaries would
bring us a high redshift measurement on Hubble parameter to reveal Hubble
tension. In addition, PBH binaries can also provide a window in exploring the
fundamental physics, e.g., the nature of dark matter [72, 73, 74], Hawking
radiation [75, 76], cosmic evolution [77], etc. Therefore, pursuing these
potential studies in PBH binaries could shed light on the nature of the
universe and shows the great value of PBH binaries as a cosmic probe.

## Acknowledgements

I thank Masahide Yamaguchi, Otto A. Hannuksela, Xiao-Dong Li, Chao Chen, Leo
W.H. Fung for useful discussions. I would like to thank Dr. Xingwei Tang for
nice discussions in the gradient descent method. This work is supported by the
Institute for Basic Science under the project code, IBS-R018-D3.

##  Appendix A The formation and merger rate of PBH binaries

In this section, we give a brief review on the formation and the merger rate
of PBH binaries, see [21, 24] for more details.

The PBH binary forms when two neighboring BHs are close enough and decouple
from the Hubble flow. Given the equation of motion of two-point masses $M$ at
rest with an initial separation $x$, the proper separation $r$ along the axis
of motion evolves as

$\displaystyle\ddot{r}-(\dot{H}+H^{2})r+\frac{2M}{r^{2}}\frac{r}{|r|}=0~{},$
(22)

where dots represent the differentiation with respect to proper time. In order
to describe the early-universe evolution, we use $s\equiv a/a_{\mathrm{eq}}$
the scale factor normalized to unity at matter-radiation equality to express
the Hubble parameter as $H(s)=(8\pi\rho_{\mathrm{eq}}/3)^{1/2}h(s)$, where $h$
is defined as $h(s)\equiv\sqrt{s^{-3}+s^{-4}}$ and $\rho_{\mathrm{eq}}$ is the
matter density at equality. Then, we can rewrite Eq. (22) by introducing
$\chi\equiv r/x$ as

$\displaystyle\chi^{\prime\prime}+\frac{sh^{\prime}+h}{s^{2}h}(s\chi^{\prime}-\chi)+\frac{1}{\lambda}\frac{1}{(sh)^{2}}\frac{1}{\chi^{2}}\frac{\chi}{|\chi|}=0~{},$
(23)

where primes denote differentiation with respect to $s$. Here, the
dimensionless parameter $\lambda$ is defined as $\lambda\equiv
4\pi\rho_{\mathrm{eq}}x^{3}/3M$. The initial condition is given in the
condition that the two neighboring BHs follow the Hubble flow $\chi(s)=s$, the
initial conditions are

$\displaystyle\chi(0)=0~{},~{}~{}~{}\chi^{\prime}(0)=1~{}.$ (24)

Then, the numerical solution in [21] shows that the binary effectively
decouples from the Hubble flow at $s\approx\lambda/3$. The corresponding
redshift is

$\displaystyle z=\frac{3(1+z_{\mathrm{eq}})}{\lambda}-1~{}.$ (25)

In order to get the merger rate of PBH binaries with mass distribution $P(m)$
at cosmic time $t$, we follow [24] and start from discretizing the PBH mass
function $n(m)$ to define a binned mass distribution $n(m_{i})$ and mass
interval $\Delta_{m}$, which follows

$\displaystyle\sum_{m_{\mathrm{min}}}^{m_{\mathrm{max}}}n(m_{i})\Delta_{m}=1~{}.$
(26)

Then, the average distance between two nearby PBHs is

$\displaystyle\langle
x_{ij}\rangle=(\bar{x}_{i}^{-3}+\bar{x}_{j}^{-3})^{-1/3}=\mu_{ij}^{1/3}\bar{x}_{ij}~{},$
(27)

where, $\mu_{ij}$ and $\bar{x}_{ij}$ are defined as

$\displaystyle\mu_{ij}=\frac{2m_{i}m_{j}f_{b}}{m_{b}f(n(m_{j})m_{i}+n(m_{i})m_{j})}~{},~{}~{}~{}\bar{x}_{ij}^{3}=\frac{3}{8\pi}\frac{m_{b}}{\rho_{\mathrm{eq}}f_{b}\Delta_{m}}~{},$
(28)

where, $f$ is the fraction of PBH energy density in matter,
$f_{b}=f(n(m_{i})+n(m_{j}))$ and $m_{b}=m_{i}+m_{j}$.
$\rho_{\mathrm{eq}}=\Omega_{\mathrm{M}}\rho_{\mathrm{crit}}(1+z_{\mathrm{eq}})^{3}$
is the total energy density of matter at matter-radiation equality. After the
formation of PBH binaries, the orbit of PBH binary shrinks due to the emission
of GWs, which leads the merger of PBH binary, the coalescence time of the PBH
binary is given in [42] by,

$\displaystyle
t=\frac{3}{85}\frac{a^{4}}{G^{3}m_{1}m_{2}(m_{1}+m_{2})}j^{7}~{},$ (29)

From Eq. (29), the merger rate of PBH binaries can be obtained from the
initial semi-major axis $a$ distribution and the initial dimensionless angular
momentum $j$ distribution (the dimensionless angular momentum can be
represented in term of eccentricity of PBH binary $e$ as $j=\sqrt{1-e^{2}}$).
The initial semi-major axis $a$ after the binary formation is numerically
given in [21],

$\displaystyle a\approx 0.1\lambda
x=\frac{0.1\bar{x}_{ij}}{f_{b}\Delta_{m}}X^{4/3},$ (30)

where, $X\equiv x^{3}/\bar{x}_{ij}^{3}$. Therefore, the initial semi-major
axis distribution is determined by the separation $x$ distribution, we follow
[20, 21, 24] that assuming PBHs follow a random distribution, the probability
distribution of the separation $x$ is

$\displaystyle\frac{dP}{dX}=\mu_{ij}^{-1}e^{-X\frac{4\pi}{3}\bar{x}_{ij}^{3}n_{T}}~{},$
(31)

where $n_{T}\equiv
f_{\mathrm{PBH}}\rho_{\mathrm{DM}}(1+z_{\mathrm{eq}})^{3}\int_{0}^{\infty}\frac{n(m)}{m}dm$.
Considering fixed $X$, the dimensionless angular momentum $j$ can be given by
Eq. (29) and Eq. (30),

$\displaystyle
j(t;X)=\left(\frac{3}{85}\frac{G^{3}m_{1}m_{2}(m_{1}+m_{2})(f_{b}\Delta_{m})^{4}}{(0.1\bar{x}_{ij})^{4}X^{16/3}}t\right)^{1/7}~{}.$
(32)

The differential probability distribution of $(X,t)$ is given by

$\displaystyle\frac{d^{2}P}{dXdt}=\frac{dP}{dX}\left(\frac{\partial
j}{\partial
t}\frac{dP}{dj}\bigg{|}_{X}\right)_{j(t;X)}=\frac{\mu_{ij}^{-1}}{7t}e^{-X\frac{4\pi}{3}\bar{x}_{ij}^{3}n_{T}}\mathcal{P}(j/j_{X})~{},$
(33)

where $j_{X}=0.5fX/f_{b}\Delta_{m}$,
$\mathcal{P}(j/j_{X})=(j/j_{X})^{2}/(1+j^{2}/j_{X}^{2})^{3/2}$. Integrating
Eq. (33) gives the merger time probability distribution

$\displaystyle\frac{dP}{dt}=\frac{\mu_{ij}^{-1}}{7t}\int
dXe^{-X\frac{4\pi}{3}\bar{x}_{ij}^{3}n_{T}}\mathcal{P}(j/j_{X})~{}.$ (34)

Then, the comoving merger rate $R_{ij}$ for PBH binaries at time $t$ is

$\displaystyle
R_{ij}(t)=\rho_{\mathrm{PBH}}\mathrm{min}\left(\frac{n(m_{i})}{m_{i}},\frac{n(m_{j})}{m_{j}}\right)\Delta_{m}\frac{dP}{dt}~{}.$
(35)

## References

  * [1] A. G. Riess et al., “A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km/s/Mpc Uncertainty from the Hubble Space Telescope and the SH0ES Team,” Astrophys. J. Lett. 934 no. 1, (2022) L7, arXiv:2112.04510 [astro-ph.CO].
  * [2] A. G. Riess and L. Breuval, “The Local Value of H0,”  8, 2023.  arXiv:2308.10954 [astro-ph.CO].
  * [3] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [4] L. Perivolaropoulos and F. Skara, “Challenges for $\Lambda$CDM: An update,” New Astron. Rev. 95 (2022) 101659, arXiv:2105.05208 [astro-ph.CO].
  * [5] E. Abdalla et al., “Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies,” JHEAp 34 (2022) 49–211, arXiv:2203.06142 [astro-ph.CO].
  * [6] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, “In the realm of the Hubble tension—a review of solutions,” Class. Quant. Grav. 38 no. 15, (2021) 153001, arXiv:2103.01183 [astro-ph.CO].
  * [7] J. D. Fernie, “The period-luminosity relation: A historical review,” Publications of the Astronomical Society of the Pacific 81 (Dec, 1969) 707.
  * [8] Pan-STARRS1 Collaboration, D. M. Scolnic et al., “The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample,” Astrophys. J. 859 no. 2, (2018) 101, arXiv:1710.00845 [astro-ph.CO].
  * [9] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, “The 6df galaxy survey: baryon acoustic oscillations and the local hubble constant,” Monthly Notices of the Royal Astronomical Society 416 no. 4, (2011) 3017–3032.
  * [10] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, “The clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at $z=0.15$,” Mon. Not. Roy. Astron. Soc. 449 no. 1, (2015) 835–847, arXiv:1409.3242 [astro-ph.CO].
  * [11] BOSS Collaboration, T. Delubac et al., “Baryon acoustic oscillations in the Ly$\alpha$ forest of BOSS DR11 quasars,” Astron. Astrophys. 574 (2015) A59, arXiv:1404.1801 [astro-ph.CO].
  * [12] B. F. Schutz, “Determining the Hubble Constant from Gravitational Wave Observations,” Nature 323 (1986) 310–311.
  * [13] D. E. Holz and S. A. Hughes, “Using gravitational-wave standard sirens,” Astrophys. J. 629 (2005) 15–22, arXiv:astro-ph/0504616.
  * [14] LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119 no. 16, (2017) 161101, arXiv:1710.05832 [gr-qc].
  * [15] LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER Collaboration, B. P. Abbott et al., “A gravitational-wave standard siren measurement of the Hubble constant,” Nature 551 no. 7678, (2017) 85–88, arXiv:1710.05835 [astro-ph.CO].
  * [16] J. Miralda-Escude, “The dark age of the universe,” Science 300 (2003) 1904–1909, arXiv:astro-ph/0307396.
  * [17] S. Hawking, “Gravitationally collapsed objects of very low mass,” Mon. Not. Roy. Astron. Soc. 152 (1971) 75.
  * [18] Y. B. Zel’dovich and I. D. Novikov, “The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model,” Soviet Astron. AJ (Engl. Transl. ), 10 (1967) 602.
  * [19] B. J. Carr and S. W. Hawking, “Black holes in the early Universe,” Mon. Not. Roy. Astron. Soc. 168 (1974) 399–415.
  * [20] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, “Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914,” Phys. Rev. Lett. 117 no. 6, (2016) 061101, arXiv:1603.08338 [astro-ph.CO]. [Erratum: Phys.Rev.Lett. 121, 059901 (2018)].
  * [21] Y. Ali-Haïmoud, E. D. Kovetz, and M. Kamionkowski, “Merger rate of primordial black-hole binaries,” Phys. Rev. D 96 no. 12, (2017) 123523, arXiv:1709.06576 [astro-ph.CO].
  * [22] Q. Ding, “Toward cosmological standard timers in primordial black hole binaries,” Phys. Rev. D 108 no. 2, (2023) 023514, arXiv:2206.03142 [astro-ph.CO].
  * [23] Y.-F. Cai, C. Chen, Q. Ding, and Y. Wang, “Cosmological standard timers from unstable primordial relics,” Eur. Phys. J. C 83 no. 10, (2023) 913, arXiv:2112.10422 [astro-ph.CO].
  * [24] Z.-C. Chen and Q.-G. Huang, “Merger Rate Distribution of Primordial-Black-Hole Binaries,” Astrophys. J. 864 no. 1, (2018) 61, arXiv:1801.10327 [astro-ph.CO].
  * [25] LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$,” Astrophys. J. Lett. 892 no. 1, (2020) L3, arXiv:2001.01761 [astro-ph.HE].
  * [26] LIGO Scientific, Virgo Collaboration, R. Abbott et al., “GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object,” Astrophys. J. Lett. 896 no. 2, (2020) L44, arXiv:2006.12611 [astro-ph.HE].
  * [27] S. Clesse and J. Garcia-Bellido, “GW190425, GW190521 and GW190814: Three candidate mergers of primordial black holes from the QCD epoch,” Phys. Dark Univ. 38 (2022) 101111, arXiv:2007.06481 [astro-ph.CO].
  * [28] LIGO Scientific, Virgo Collaboration, R. Abbott et al., “GW190521: A Binary Black Hole Merger with a Total Mass of $150M_{\odot}$,” Phys. Rev. Lett. 125 no. 10, (2020) 101102, arXiv:2009.01075 [gr-qc].
  * [29] V. De Luca, V. Desjacques, G. Franciolini, P. Pani, and A. Riotto, “GW190521 Mass Gap Event and the Primordial Black Hole Scenario,” Phys. Rev. Lett. 126 no. 5, (2021) 051101, arXiv:2009.01728 [astro-ph.CO].
  * [30] P. Bender, A. Brillet, I. Ciufolini, A. Cruise, C. Cutler, K. Danzmann, F. Fidecaro, W. Folkner, J. Hough, P. McNamara, et al., “Lisa pre-phase a report,” Max-Planck Institut für Quantenoptik (1998) .
  * [31] G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner, and E. S. Phinney, “Laser interferometry for the big bang observer,” Class. Quant. Grav. 23 (2006) 4887–4894. [Erratum: Class.Quant.Grav. 23, 7361 (2006)].
  * [32] T. Nakamura et al., “Pre-DECIGO can get the smoking gun to decide the astrophysical or cosmological origin of GW150914-like binary black holes,” PTEP 2016 no. 9, (2016) 093E01, arXiv:1607.00897 [astro-ph.HE].
  * [33] Q. Ding, “Detectability of primordial black hole binaries at high redshift,” Phys. Rev. D 104 no. 4, (2021) 043527, arXiv:2011.13643 [astro-ph.CO].
  * [34] I. Mandel, W. M. Farr, and J. R. Gair, “Extracting distribution parameters from multiple uncertain observations with selection biases,” Mon. Not. Roy. Astron. Soc. 486 no. 1, (2019) 1086–1093, arXiv:1809.02063 [physics.data-an].
  * [35] S. Vitale, D. Gerosa, W. M. Farr, and S. R. Taylor, “Inferring the properties of a population of compact binaries in presence of selection effects,” arXiv:2007.05579 [astro-ph.IM].
  * [36] D. Gerosa, G. Pratten, and A. Vecchio, “Gravitational-wave selection effects using neural-network classifiers,” Phys. Rev. D 102 no. 10, (2020) 103020, arXiv:2007.06585 [astro-ph.HE].
  * [37] E. E. Flanagan and S. A. Hughes, “Measuring gravitational waves from binary black hole coalescences: 1. Signal-to-noise for inspiral, merger, and ringdown,” Phys. Rev. D 57 (1998) 4535–4565, arXiv:gr-qc/9701039.
  * [38] P. A. Rosado, P. D. Lasky, E. Thrane, X. Zhu, I. Mandel, and A. Sesana, “Detectability of Gravitational Waves from High-Redshift Binaries,” Phys. Rev. Lett. 116 no. 10, (2016) 101102, arXiv:1512.04950 [astro-ph.CO].
  * [39] C. J. Moore, R. H. Cole, and C. P. L. Berry, “Gravitational-wave sensitivity curves,” Class. Quant. Grav. 32 no. 1, (2015) 015014, arXiv:1408.0740 [gr-qc].
  * [40] S. Droz, D. J. Knapp, E. Poisson, and B. J. Owen, “Gravitational waves from inspiraling compact binaries: Validity of the stationary phase approximation to the Fourier transform,” Phys. Rev. D 59 (1999) 124016, arXiv:gr-qc/9901076.
  * [41] V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczepanczyk, I. Bartos, M. Campanelli, S. Klimenko, C. O. Lousto, and R. O’Shaughnessy, “Eccentricity estimate for black hole mergers with numerical relativity simulations,” Nature Astron. 6 no. 3, (2022) 344–349, arXiv:2009.05461 [astro-ph.HE].
  * [42] P. C. Peters, “Gravitational Radiation and the Motion of Two Point Masses,” Phys. Rev. 136 (1964) B1224–B1232.
  * [43] K. K. Y. Ng, G. Franciolini, E. Berti, P. Pani, A. Riotto, and S. Vitale, “Constraining High-redshift Stellar-mass Primordial Black Holes with Next-generation Ground-based Gravitational-wave Detectors,” Astrophys. J. Lett. 933 no. 2, (2022) L41, arXiv:2204.11864 [astro-ph.CO].
  * [44] M. Raidal, C. Spethmann, V. Vaskonen, and H. Veermäe, “Formation and Evolution of Primordial Black Hole Binaries in the Early Universe,” JCAP 02 (2019) 018, arXiv:1812.01930 [astro-ph.CO].
  * [45] D. W. Hogg, “Distance measures in cosmology,” arXiv:astro-ph/9905116.
  * [46] A. Dolgov and J. Silk, “Baryon isocurvature fluctuations at small scales and baryonic dark matter,” Phys. Rev. D 47 (1993) 4244–4255.
  * [47] B. J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201 (1975) 1–19.
  * [48] S. Ruder, “An overview of gradient descent optimization algorithms,” arXiv preprint arXiv:1609.04747 (2016) .
  * [49] L. Ciampiconi, A. Elwood, M. Leonardi, A. Mohamed, and A. Rozza, “A survey and taxonomy of loss functions in machine learning,” arXiv preprint arXiv:2301.05579 (2023) .
  * [50] X. Chen, S. Li, and Z. Cao, “Mass–redshift degeneracy for the gravitational-wave sources in the vicinity of supermassive black holes,” Mon. Not. Roy. Astron. Soc. 485 no. 1, (2019) L141–L145, arXiv:1703.10543 [astro-ph.HE].
  * [51] LIGO Scientific, Virgo Collaboration, R. Abbott et al., “GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run,” Phys. Rev. X 11 (2021) 021053, arXiv:2010.14527 [gr-qc].
  * [52] M. Moresco et al., “Unveiling the Universe with emerging cosmological probes,” Living Rev. Rel. 25 no. 1, (2022) 6, arXiv:2201.07241 [astro-ph.CO].
  * [53] D. Scolnic, A. G. Riess, J. Wu, S. Li, G. S. Anand, R. Beaton, S. Casertano, R. I. Anderson, S. Dhawan, and X. Ke, “CATS: The Hubble Constant from Standardized TRGB and Type Ia Supernova Measurements,” Astrophys. J. Lett. 954 no. 1, (2023) L31, arXiv:2304.06693 [astro-ph.CO].
  * [54] D. Stern, R. Jimenez, L. Verde, M. Kamionkowski, and S. A. Stanford, “Cosmic Chronometers: Constraining the Equation of State of Dark Energy. I: H(z) Measurements,” JCAP 02 (2010) 008, arXiv:0907.3149 [astro-ph.CO].
  * [55] B. V. Lehmann, S. Profumo, and J. Yant, “The Maximal-Density Mass Function for Primordial Black Hole Dark Matter,” JCAP 04 (2018) 007, arXiv:1801.00808 [astro-ph.CO].
  * [56] T. Suyama and S. Yokoyama, “A novel formulation of the primordial black hole mass function,” PTEP 2020 no. 2, (2020) 023E03, arXiv:1912.04687 [astro-ph.CO].
  * [57] V. De Luca, G. Franciolini, and A. Riotto, “On the primordial black hole mass function for broad spectra,” Phys. Lett. B 807 (2020) 135550, arXiv:2001.04371 [astro-ph.CO].
  * [58] X.-M. Bi, L. Chen, and K. Wang, “Primordial Black Hole Mass Function with Mass Gap,” arXiv:2210.16479 [astro-ph.CO].
  * [59] M. Raidal, V. Vaskonen, and H. Veermäe, “Gravitational Waves from Primordial Black Hole Mergers,” JCAP 09 (2017) 037, arXiv:1707.01480 [astro-ph.CO].
  * [60] V. Vaskonen and H. Veermäe, “Lower bound on the primordial black hole merger rate,” Phys. Rev. D 101 no. 4, (2020) 043015, arXiv:1908.09752 [astro-ph.CO].
  * [61] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, “The evolution of primordial black holes and their final observable spins,” JCAP 04 (2020) 052, arXiv:2003.02778 [astro-ph.CO].
  * [62] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, “Primordial Black Holes Confront LIGO/Virgo data: Current situation,” JCAP 06 (2020) 044, arXiv:2005.05641 [astro-ph.CO].
  * [63] V. De Luca, G. Franciolini, A. Kehagias, P. Pani, and A. Riotto, “Primordial black holes in matter-dominated eras: The role of accretion,” Phys. Lett. B 832 (2022) 137265, arXiv:2112.02534 [astro-ph.CO].
  * [64] L. Liu, Z.-K. Guo, and R.-G. Cai, “Effects of the merger history on the merger rate density of primordial black hole binaries,” Eur. Phys. J. C 79 no. 8, (2019) 717, arXiv:1901.07672 [astro-ph.CO].
  * [65] Y. Wu, “Merger history of primordial black-hole binaries,” Phys. Rev. D 101 no. 8, (2020) 083008, arXiv:2001.03833 [astro-ph.CO].
  * [66] O. A. Hannuksela, T. E. Collett, M. Çalışkan, and T. G. F. Li, “Localizing merging black holes with sub-arcsecond precision using gravitational-wave lensing,” Mon. Not. Roy. Astron. Soc. 498 no. 3, (2020) 3395–3402, arXiv:2004.13811 [astro-ph.HE].
  * [67] Z.-F. Wu, L. W. L. Chan, M. Hendry, and O. A. Hannuksela, “Reducing the impact of weak-lensing errors on gravitational-wave standard sirens,” Mon. Not. Roy. Astron. Soc. 522 no. 3, (2023) 4059–4077, arXiv:2211.15160 [astro-ph.CO].
  * [68] S.-J. Huang, Y.-M. Hu, X. Chen, J.-d. Zhang, E.-K. Li, Z. Gao, and X.-Y. Lin, “Measuring the Hubble constant using strongly lensed gravitational wave signals,” JCAP 08 (2023) 003, arXiv:2304.10435 [astro-ph.CO].
  * [69] Z.-C. Chen, S.-S. Du, Q.-G. Huang, and Z.-Q. You, “Constraints on primordial-black-hole population and cosmic expansion history from GWTC-3,” JCAP 03 (2023) 024, arXiv:2205.11278 [astro-ph.CO].
  * [70] LIGO Scientific, Virgo, KAGRA Collaboration, R. Abbott et al., “Constraints on the Cosmic Expansion History from GWTC–3,” Astrophys. J. 949 no. 2, (2023) 76, arXiv:2111.03604 [astro-ph.CO].
  * [71] L. W. H. Fung, T. Broadhurst, and G. F. Smoot, “Pure Gravitational Wave Estimation of Hubble’s Constant using Neutron Star-Black Hole Mergers,” arXiv:2308.02440 [astro-ph.CO].
  * [72] N. Speeney, A. Antonelli, V. Baibhav, and E. Berti, “Impact of relativistic corrections on the detectability of dark-matter spikes with gravitational waves,” Phys. Rev. D 106 no. 4, (2022) 044027, arXiv:2204.12508 [gr-qc].
  * [73] S. Pilipenko, M. Tkachev, and P. Ivanov, “Evolution of a primordial binary black hole due to interaction with cold dark matter and the formation rate of gravitational wave events,” Phys. Rev. D 105 no. 12, (2022) 123504, arXiv:2205.10792 [astro-ph.CO].
  * [74] A. Akil and Q. Ding, “A dark matter probe in accreting pulsar-black hole binaries,” JCAP 09 (2023) 011, arXiv:2304.08824 [astro-ph.HE].
  * [75] Y.-F. Cai, C. Chen, Q. Ding, and Y. Wang, “Ultrahigh-energy gamma rays and gravitational waves from primordial exotic stellar bubbles,” Eur. Phys. J. C 82 no. 5, (2022) 464, arXiv:2105.11481 [astro-ph.CO].
  * [76] H. Li and J. Wang, “Towards the merger of Hawking radiating black holes,” Int. J. Mod. Phys. D 30 no. 08, (2021) 2150060, arXiv:2002.08048 [gr-qc].
  * [77] T. Nakama, “Gravitational waves from binary black holes as probes of the structure formation history,” Phys. Dark Univ. 28 (2020) 100476, arXiv:2001.01407 [astro-ph.CO].