# On the population III binary black hole mergers with intermediate mass black
holes: dependence on common envelope parameter

Kotaro Hijikawa,1 Tomoya Kinugawa,2 Ataru Tanikawa,3 Takashi Yoshida4,1 and
Hideyuki Umeda1
1Department of Astronomy, Graduate School of Science, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, 113-0033 Tokyo, Japan
2Institute for Cosmic Ray Research, The University of Tokyo, Kashiwa, Chiba,
Japan
3Department of Earth Science and Astronomy, College of Arts and Sciences, The
University of Tokyo, Meguro-ku, Tokyo, Japan
4Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa
Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The current gravitational wave (GW) detectors have successfully observed many
binary compact objects, and the third generation ground-based GW detectors
such as Einstein telescope and space-borne detectors such as LISA will start
their GW observation in a decade. Ahead of the arrival of this new era, we
perform a binary population synthesis calculation for very massive ($\sim$
100–1000 $M_{\odot}$) population (Pop.) III stars, derive the various property
of binary black hole (BBH) mergers with intermediate mass black holes (IMBHs)
and investigate the dependence on common envelope parameter $\alpha\lambda$
which is still not a well understood parameter. We find that the maximum mass
of primary BH mass is larger for smaller value of common envelope parameter.
In this study, we adopt double power law initial mass function (IMF) for Pop.
III stars, and put some constraints on Pop. III IMF by comparing our obtained
merger rate density at the local Universe with that derived from gravitational
wave (GW) observation. We compute the detection rate and show that the third
generation ground-based GW detector, Einstein telescope, have a potential to
detect $\sim$ 10–1000 BBHs with IMBHs per year. We also find that we may be
able to obtain the insight into $\alpha\lambda$ if a BBH with total mass
$\gtrsim 500M_{\odot}$ are detected by advanced LIGO (O4) or LISA.

###### keywords:

stars: Population III, binaries: general, black hole mergers, gravitational
waves

††pubyear: 2022††pagerange: On the population III binary black hole mergers
with intermediate mass black holes: dependence on common envelope parameter–On
the population III binary black hole mergers with intermediate mass black
holes: dependence on common envelope parameter

##  1 Introduction

Many binary black hole (BBH) mergers have been discovered by the gravitational
wave (GW) detectors, advanced LIGO and Virgo. The catalog, GWTC-3, was
released last year by the LIGO Scientific Collaboration, the Virgo
Collaboration and the KAGRA Collaboration (The LIGO Scientific Collaboration
et al., 2021b), and the number of discovered BBH mergers is approaching 100.
Detected BBHs have relatively high mass ($\sim$ 20–50 $M_{\odot}$) compared to
Galactic BHs, and have low spin ($\sim$ 0). Origin of such discovered BBH
mergers have been suggested by many authors. For example, field population
(Pop.) I or II isolated binary stars (e.g., Bethe & Brown, 1998; Marchant et
al., 2016; Mandel & de Mink, 2016; Stevenson et al., 2017; Eldridge et al.,
2017; Eldridge et al., 2019; Mapelli et al., 2017, 2019; Mapelli et al., 2021;
Kruckow et al., 2018; Belczynski et al., 2020; Bavera et al., 2021; Riley et
al., 2021), Pop. III isolated binary stars (e.g., Kinugawa et al., 2014, 2016,
2020, 2021b; Hartwig et al., 2016; Belczynski et al., 2017; Tanikawa et al.,
2021b, 2022), a dynamical formation in stellar clusters (e.g. Portegies Zwart
& McMillan, 2000; O’Leary et al., 2006; Sadowski et al., 2008; Banerjee et
al., 2010; Banerjee, 2017, 2018a, 2018b, 2021a, 2021b; Downing et al., 2010,
2011; Tanikawa, 2013; Bae et al., 2014; Ziosi et al., 2014; Mapelli, 2016;
Rodriguez et al., 2016, 2018; Rodriguez et al., 2019; Fujii et al., 2017; Park
et al., 2017; Hong et al., 2018, 2020; Samsing, 2018; Samsing et al., 2018; Di
Carlo et al., 2019; Di Carlo et al., 2020; Kumamoto et al., 2019, 2020;
Kumamoto et al., 2021; Anagnostou et al., 2020; Kremer et al., 2020; Fragione
& Banerjee, 2021; Trani et al., 2021; Weatherford et al., 2021), a Pop. III
dynamical formation (e.g. Liu & Bromm, 2020, 2021; Liu et al., 2021; Wang et
al., 2022), a formation in the galactic nuclei (e.g. O’Leary et al., 2009;
Antonini & Perets, 2012; Antonini & Rasio, 2016; Stone et al., 2017; Leigh et
al., 2018; McKernan et al., 2018; Rasskazov & Kocsis, 2019; Yang et al., 2019;
Arca Sedda, 2020; Tagawa et al., 2020a; Tagawa et al., 2020b, 2021a, 2021b),
or triple or multiple star systems (e.g. Antonini et al., 2017; Silsbee &
Tremaine, 2017; Rodriguez & Antonini, 2018; Fragione & Kocsis, 2019; Fragione
et al., 2020; Hamers et al., 2021; Trani et al., 2022; Vynatheya & Hamers,
2022).

In addition to the typical population, a more massive BBH, GW190521, with
total mass $\sim 150M_{\odot}$ is also discovered (Abbott et al., 2020). Since
its primary BH mass is about 90 $M_{\odot}$ and it lies in the pair-
instability mass gap 111Note that GW190521 may be a pair of two BHs above and
below the pair-instability mass gap (Fishbach & Holz, 2020; Nitz & Capano,
2021)., it is considered to be difficult to form such BBHs through a general
isolated binary formation channel, and the hierarchical merger process in the
dynamical formation scenario can explain the existence of $\sim 90M_{\odot}$
BH. However, there are some isolated binary scenarios to explain the formation
of GW190521. For example, by reducing
${}^{12}\mathrm{C}(\alpha,\gamma)^{16}\mathrm{O}$ reaction rate within its
uncertainty, the pair-instability mass gap shifts toward the heavier side, and
GW190521-like BBHs can be formed through isolated binary formation channel
(Farmer et al., 2019, 2020; Belczynski, 2020). Furthermore, a massive
extremely metal-poor star or a Pop. III star may be able to become a BH with
65–100 $M_{\odot}$, and thus the pair-instability mass gap may be filled, and
GW190521-like BBHs can be formed through isolated binary formation channel
(Umeda et al., 2020; Kinugawa et al., 2021a; Tanikawa et al., 2021a, b;
Farrell et al., 2021; Vink et al., 2021). From the discovery of GW190521 and
studies on this observation, the equivalent objects or more massive BBH
mergers, i.e., the mergers of an intermediate mass black hole (IMBH; $\gtrsim
100M_{\odot}$) with a stellar mass BH or an IMBH are expected to be detected
by the future GW observations.

So far, a number of stellar mass BHs ($\lesssim 100M_{\odot}$) have been found
by GW observations or X-ray binary observations. On the other hand, few IMBHs
have yet been discovered (one can find more information about the searchs for
IMBHs in the review, Greene et al., 2020). The most reliable one is the merger
remnant of the BBH merger, GW190521 (Abbott et al., 2020). The X-ray sources
HLX-1 is also considered to be an IMBH (Farrell et al., 2009). The presence of
the supermassive BHs (SMBHs; $\gtrsim 10^{6}M_{\odot}$) at the center of
galaxies has been confirmed, and recently the Event Horizon Telescope
successfully got images of SMBHs at the center of M87 (Event Horizon Telescope
Collaboration et al., 2019a, b, c, d, e, f) and the Milky Way (Event Horizon
Telescope Collaboration et al., 2022a, b, c, d, e, f). The formation and
evolution of SMBHs are still not well understood. According to some models,
the seed of SMBHs has repeatedly merged with other BHs multiple times and
increases its mass. Since there are few IMBHs that have been found, it is
difficult to constrain the formation theory and to know where the SMBHs came
from. Therefore, it is important to know what mass of IMBHs existed at each
redshift in order to impose restrictions on the formation model of the SMBH.
As well as the probe for SMBH, IMBHs may be useful for other problems. For
example, the findings of IMBHs just above the pair-instability mass gap can be
used as the “standard silen” for the cosmology (Farr et al., 2019; Farmer et
al., 2019, 2020). From the GW observation, the luminosity distance and the
redshifted mass can be obtained. If we have a knowledge of the upper end of
the pair-instability mass gap, then we can get the redshift of the sources
from the true and redshifted mass. Furthermore, the pair of stellar mass BH
and IMBH, or the intermediate mass ratio inspiral (also know as IMRI) will
allow us to test the gravitational theory (e.g., Sopuerta & Yunes, 2009;
Amaro-Seoane, 2018).

In the future, the third generation ground-based GW detector, Einstein
telescope (Punturo et al., 2010; Sathyaprakash et al., 2012) or Cosmic
Explorer (Reitze et al., 2019), and space-borne GW detector, LISA (Amaro-
Seoane et al., 2017), TianQin (Luo et al., 2016; Wang et al., 2019), TianGO
(Kuns et al., 2020), or DECIGO (Seto et al., 2001) will start the construction
in ten years. Since their sensitivity will improve and their horizon redshift
will be $z\sim$ 10–1000 for BBH mergers with total mass $M_{\mathrm{t}}\sim$
10–10 000 $M_{\odot}$, the inspiral and merger of more massive BBHs, or binary
IMBHs, located in the earlier Universe are expected to be discovered.
Previously, the GWs from IMBH-IMBH mergers during the runaway merging process
were investigated (e.g., Matsubayashi et al., 2004; Shinkai et al., 2017). As
to the isolated binary, the BBH mergers above the pair-instability mass gap
originated from isolated Pop. I/II binaries were investigated by Mangiagli et
al. (2019). Pop. III star formation simulations have been performed (e.g.,
Hirano et al., 2014, 2015; Susa et al., 2014), and it is claimed that the
typical Pop. III stellar mass is 10–100 $M_{\odot}$. Furthermore, very massive
($\sim$ 100–1000 $M_{\odot}$) Pop. III stars might be formed (e.g., Hirano et
al., 2015). A Pop. III star with 100–200 $M_{\odot}$ finally causes the pair-
instability supernova, and leaves no remnant. On the other hand, more massive
Pop. III star than this directly collapses and finally becomes an IMBH.
Therefore, in our previous paper (Hijikawa et al., 2021), we focused on the
very massive ($\sim$ 100–1000 $M_{\odot}$) Pop. III stars, performed binary
population synthesis calculation, derived the merger rate density of BBHs with
IMBHs assuming single power law initial mass function (IMF), and compared this
with the observation. Thus far, no BBH with IMBHs have been discovered, and
hence the upper limit on the merger rate density is obtained. We compared our
results with this upper limit, and found that the power of single power law
IMF needs to be lower than 2.8. However, this is inconsistent with the results
of Pop. III star formation simulation claiming a relatively top heavy IMF
where $\lesssim 100M_{\odot}$.

In this study, we adopt the double power law IMF with extensive parameter
range, perform a binary population synthesis calculation, and investigate what
types of BBH mergers with IMBHs will be observed by current and future GW
detectors when using a reasonable IMF which meets the Pop. III star formation
simulation and current GW observation. In addition, we investigate the
dependence of the properties of BBHs with IMBHs on the common envelope
parameter $\alpha\lambda$.

This paper is constructed as follows. In Section 2, we describe our method of
binary population synthesis calculations. In Section 3, we show the formation
channel, primary BH mass distribution, mass ratio distribution, effective spin
distribution, delay time distribution and the redshift evolution of merger
rate density for all sub-populations of BBH mergers with IMBHs. We also show
in detail the dependence of such distribution on the common envelope parameter
$\alpha\lambda$. We compare our obtained merger rate density with the upper
limit and impose restrictions on IMF and $\alpha\lambda$ models, and compute
the detection rate and detected BBH mass distribution in Section 4. Finally,
we summarize our paper in Section 5.

##  2 Method

In order to calculate the merger rate density of Pop. III BBH mergers for
various Pop. III IMFs, we perform binary population synthesis calculations.
Our binary population synthesis code is identical with that of Hijikawa et al.
(2021). In section 2.1, we overview the stellar and binary evolution
implemented in our code. In section 2.2, we describe IMFs and other initial
conditions we adopt. How we compute the merger rate density of Pop. III BBH
mergers is described in Section 2.3. Finally, in Section 2.4, we summarize our
parameter sets.

###  2.1 Stellar and Binary Evolution Physics

Our binary population synthesis code is based on BSE (Hurley et al., 2000,
2002) 111http://astronomy.swin.edu.au/j̃hurley/. Some additions and changes
have been made to this, and these are described below.

In our code, the fitting formulae for $Z=10^{-8}Z_{\odot}$ stars were
implemented (Tanikawa et al., 2020, 2021a, 2021b, 2022), published as
BSEEMP222https://github.com/atrtnkw/bseemp. $Z$ and $Z_{\odot}=0.02$ are the
stellar metallicity and sollar metallicity, respectively. In their studies,
the evolution of $Z=10^{-8}Z_{\odot}$ stars up to 1280 $M_{\odot}$. are
calculated by the 1D stellar evolution code, HOSHI (Takahashi et al., 2016,
2019, 2018; Yoshida et al., 2019). In our study, we treat $Z=10^{-8}Z_{\odot}$
stars as Pop. III stars ($Z=0$). There is little difference in the stellar
evolution between $Z=10^{-8}Z_{\odot}$ stars and $Z=0$ stars for the following
reason. As the abundance of heavy elements in a $Z=10^{-8}Z_{\odot}$ star is
extremely low, the CNO cycle does not operate initially. The abundance of
carbon, nitrogen and oxygen are increased by the triple alpha reaction, and
then the CNO cycle commences. These evolutionary properties are same as those
of $Z=0$ stars. Therefore, the difference in metallicity, $Z=10^{-8}Z_{\odot}$
and $Z=0$, has little effect on the stellar evolution. Our code has two
stellar models in which the convective overshooting is more effective and less
so, L and M models (Yoshida et al., 2019; Tanikawa et al., 2020, 2021a, 2021b,
2022). In this paper, we use only L model.

We adopt the ‘rapid’ model in Fryer et al. (2012) with the modification for
the pulsational pair-instability and pair-instability supernova (see equations
5–7 in Tanikawa et al., 2021b). We assume that the mass of a BH formed through
pulsational pair-instability is 45 $M_{\odot}$, and the minimum helium core
mass of direct collapse is 135 $M_{\odot}$ in the same way as Belczynski et
al. (2016). Therefore, the mass range of pair-instability mass gap is 45–135
$M_{\odot}$.

In a binary system, orbital angular momenta and spin angular momenta are
exchanged by tidal interaction. In our code, the tidal coefficient factor $E$
for the dynamical tide is implemented as

$E=\left\\{\begin{aligned}
&10^{-0.42}\left(\frac{R_{\mathrm{con}}}{R}\right)^{7.5}&&\text{for H-rich
stars},\\\ &10^{-0.93}\left(\frac{R_{\mathrm{con}}}{R}\right)^{6.7}&&\text{for
He-rich stars},\end{aligned}\right.$ (1)

(Qin et al., 2018), where $R$ is the stellar radius and $R_{\mathrm{con}}$ is
the convective core radius, respectively. We use the values
$R_{\mathrm{con}}=1R_{\odot}$ and $0.5R_{\odot}$ for main-sequence stars and
naked helium stars, respectively (Tanikawa et al., 2020; Kinugawa et al.,
2020).

When one of the star in the binary system fills its Roche lobe, the stellar
mass overflows. If this mass transfer is unstable, the binary system enters a
common envelope phase; if not, the stable mass transfer occurs. In our code,
the mass transfer rate of stable mass transfer is calculated as

$\displaystyle\dot{M}_{\mathrm{don}}=-\frac{f(\mu)M_{\mathrm{don}}}{\sqrt{R_{\mathrm{don}}^{3}/GM_{\mathrm{don}}}}\left(\frac{\Delta
R_{\mathrm{don}}}{R_{\mathrm{don}}}\right)^{n+3/2}d_{n},$ (2)

where

$\displaystyle
f(\mu)=\frac{4\mu\sqrt{\mu}\sqrt{1-\mu}}{(\sqrt{\mu}+\sqrt{1-\mu})^{4}}\left(\frac{a}{R_{\mathrm{don}}}\right)^{3},$
(3)

$\mu=M_{\mathrm{don}}/M_{\mathrm{tot}}$, $\Delta
R_{\mathrm{don}}=R_{\mathrm{don}}-R_{\mathrm{L,\mathrm{don}}}$, $G$,
$M_{\mathrm{don}}$, $R_{\mathrm{don}}$, $n$, $d_{n}$, $a$, $M_{\mathrm{tot}}$,
and $R_{\mathrm{L,\mathrm{don}}}$ are the gravitational constant, the donor
mass, the donor radius, the polytropic index of the donor star, the
normalization factor depending on $n$, orbital separation, the total mass, and
the donor Roche lobe radius, respectively (e.g., Paczyński & Sienkiewicz,
1972; Savonije, 1978; Edwards & Pringle, 1987). In our study, we assume that
the mass transfer is conservative, unless the accretor is a compact object.

If a mass transfer becomes unstable
($\zeta_{\mathrm{L}}>\zeta_{\mathrm{ad}}$), the binary system enters a common
envelope phase. $\zeta_{\mathrm{L}}=\mathrm{d}\log
R_{\mathrm{L}}/\mathrm{d}\log M$ is the response of the Roche lobe radius
$R_{\mathrm{L}}$ to the change in stellar mass $M$ and
$\zeta_{\mathrm{ad}}=(\mathrm{d}\log R/\mathrm{d}\log M)_{\mathrm{ad}}$ is the
response of stellar radius $R$ to the change in stellar mass within the
dynamical timescale. We use the same $\zeta_{\mathrm{ad}}$ value as Kinugawa
et al. (2014), and this criteria is roughly consistent with Pavlovskii et al.
(2017). In this study, we use the $\alpha\lambda$ formalism for common
envelope evolution (Webbink, 1984), and the orbital separation just after the
common envelope phase $a_{\mathrm{f}}$ is calculated by the following energy
budget if the accretor star does not have a clear core-envelope structure.

$\displaystyle\alpha\left(\frac{GM_{\mathrm{c,don}}M_{\mathrm{acc}}}{2a_{\mathrm{f}}}-\frac{GM_{\mathrm{don}}M_{\mathrm{acc}}}{2a_{\mathrm{i}}}\right)=\frac{GM_{\mathrm{don}}M_{\mathrm{env,don}}}{\lambda
R_{\mathrm{don}}}.$ (4)

Here, $M_{\mathrm{c,don}}$ and $M_{\mathrm{env,don}}$ are the stellar core and
envelope mass of the donor star, $M_{\mathrm{acc}}$ is the accretor stellar
mass and $a_{\mathrm{i}}$ is the orbital separation just before the common
envelope phase. $\alpha$ is the efficiency parameter showing how much orbital
energy is used to strip the stellar envelope. $\lambda$ is the binding energy
parameter. The value of the common envelope parameters, $\alpha$ and
$\lambda$, are still not well known. Thus, in this study, we use various
values for the product of the common envelope parameters, $\alpha\lambda$:
$\alpha\lambda=$ 0.01, 0.1, 0.5, 1.0, 2.0 and 5.0. We call the common envelope
phase during which the accretor star does not have a core-envelope structure
“single common envelope”.

When the accretor star which does not fill its Roche lobe is also a giant
star, the orbital energy is used to strip the envelope of both stars. We call
this common envelope “double common envelope”, and the orbital separation just
after the common envelope phase $a_{\mathrm{f}}$ is calculated as

$\displaystyle\alpha\left(\frac{GM_{\mathrm{c,don}}M_{\mathrm{c,acc}}}{2a_{\mathrm{f}}}-\frac{GM_{\mathrm{don}}M_{\mathrm{acc}}}{2a_{\mathrm{i}}}\right)=$
$\displaystyle\frac{GM_{\mathrm{don}}M_{\mathrm{env,don}}}{\lambda
R_{\mathrm{don}}}+$
$\displaystyle\frac{GM_{\mathrm{acc}}M_{\mathrm{env,acc}}}{\lambda
R_{\mathrm{acc}}}.$ (5)

$M_{\mathrm{c,acc}}$, $M_{\mathrm{env,acc}}$ and $R_{\mathrm{acc}}$ are the
core and envelope mass and stellar radius of accretor star. In our
calculation, when both stars fill their Roche lobes, the binary system enters
a common envelope phase even though $\zeta_{\mathrm{L}}>\zeta_{\mathrm{ad}}$
is not met. We call this common envelope “contact common envelope”, and the
orbital separation just after the common envelope phase $a_{\mathrm{f}}$ is
calculated as equation (2.1).

###  2.2 Initial Conditions

In this study, we adopt various IMFs and calculate the merger rate density of
BBH mergers. In Section 2.2.1, we briefly review the recent Pop. III star
formation simulation, and describe our IMFs for Pop. III stars. Other initial
conditions such as initial mass ratio function, initial orbital separation
function and initial eccentricity function are shown in Section 2.2.2.

####  2.2.1 Initial Mass Functions

Pop. III stars are formed from metal free gas. Star forming clouds are formed
in the mini-halo and at the center the gas cloud collapses due to the
$\mathrm{H}_{2}$ and HD cooling. The cooling efficiency is weaker than in the
Pop. I/II case, and hence the Jeans mass of the gas cloud is larger. According
to Pop. III star formation simulation (Hirano et al., 2014, 2015; Susa et al.,
2014), the ZAMS mass distribution of Pop. III stars is roughly flat for
$\lesssim 100M_{\odot}$ and there is break around at $\sim 100M_{\odot}$. In
this study, in order to reproduce this structure, we adopt a double power law
IMF for the Pop. III single stars and primary stars of Pop. III binaries. Our
double power law IMFs $f(m)$ have three parameters as shown in the following
equations.

$f(m)\propto\left\\{\begin{aligned} &m^{-\gamma_{1}}&&(m_{\mathrm{min}}\leq
m\leq m_{\mathrm{crit}}),\\\ &m_{\mathrm{crit}}^{\gamma_{2}-\gamma_{1}}\times
m^{-\gamma_{2}}&&(m_{\mathrm{crit}}<m\leq
m_{\mathrm{max}}),\end{aligned}\right.$ (6)

where $m$ is the stellar mass and $m_{\mathrm{min}}=$ 10 $M_{\odot}$ and
$m_{\mathrm{max}}=$ 1500 $M_{\odot}$ are the minimum mass and maximum mass in
our calculation, respectively. At the critical mass $m_{\mathrm{crit}}$, the
slope of the IMF changes from $\gamma_{1}$ to $\gamma_{2}$. In order to
connect the IMF continuously at $m=m_{\mathrm{crit}}$, the coefficient
$m_{\mathrm{crit}}^{\gamma_{2}-\gamma_{1}}$ is multiplied. According to Pop.
III star formation simulations (Susa et al., 2014; Hirano et al., 2014, 2015),
we set the critical mass $m_{\mathrm{crit}}$ to 100, 200 and 300 $M_{\odot}$.
The slope of lower mass region ($m<m_{\mathrm{crit}}$), $\gamma_{1}$, is set
to 0 and 1 to reproduce a top heavy function in the lower mass region. The
slope of higher mass region ($m>m_{\mathrm{crit}}$), $\gamma_{2}$, is set to
1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5 and 5.0.

####  2.2.2 Other Initial Conditions

In order to determine the initial condition, we set the initial mass ratio,
initial orbital separation, and initial eccentricity distribution in addition
to the IMF. The initial mass ratio distribution is flat
($\propto\mathrm{const}.$), and the mass ratio range is $q_{\mathrm{min}}$–1,
where $q_{\mathrm{min}}=10/M_{1}$ and $M_{1}$ is the mass of initially heavier
star. The initial orbital separation distribution is logarithmically flat
distribution, and the orbital separation range is $a_{\mathrm{min}}$–
$10^{5}R_{\odot}$, where $a_{\mathrm{min}}$ is determined so that the ZAMS
radius does not exceed the Roche lobe radius. In this study, we assume that
the initial eccentricity is zero (circular orbit).

###  2.3 Rate Calculation

In this study, we sample $N_{\mathrm{sam}}=10^{6}$ binaries consisting of ZAMS
stars with the metallicity $10^{-8}Z_{\odot}$ and primary mass
$\in[10,1500]M_{\odot}$ according to various initial functions and follow the
evolution of $N_{\mathrm{sam}}$ binaries for each parameter sets. The detail
of parameter sets are described in Section 2.4. The merger rate
$\mathcal{R}(z)$ at the redshift $z$ is calculated as follows.

$\displaystyle\mathcal{R}(z(t))=\int_{0}^{t}f_{\mathrm{b}}\frac{\mathrm{SFR}(\tau)}{\langle
m\rangle}\frac{n(t-\tau)}{N_{\mathrm{sam}}}~{}\mathrm{d}\tau.$ (7)

$t$ is the age of the Universe at the redshift $z$ and is calculated assuming
the $\Lambda$-CDM cosmology. The Hubble constant is
$H_{0}=67.4~{}\mathrm{km}~{}\mathrm{s}^{-1}~{}\mathrm{Mpc}^{-1}$ and the
matter density parameter is $\Omega_{\mathrm{m0}}=0.315$ (Planck Collaboration
et al., 2020). $\mathrm{SFR}(\tau)$ is the star formation rate, and we adopt
the star formation history in de Souza et al. (2011) for Pop. III stars, but
reduce it by a factor of three (Inayoshi et al., 2016; Kinugawa et al., 2020).
We assume that the binary fraction $f_{\mathrm{b}}$ is 0.5.

$\displaystyle\langle
m\rangle=\int_{I}mf(m)~{}\mathrm{d}m+f_{\mathrm{b}}\int_{I}\mathrm{d}m~{}f(m)\int_{J(m)}mqg(q;m)~{}\mathrm{d}q,$
(8)

is the mean mass of Pop. III stellar objects (single stars and binaries).
$q\in(0,1]$ is the mass ratio,
$I=[m_{\mathrm{min}},m_{\mathrm{max}}]=[10,1500]M_{\odot}$ is the primary mass
range, $g(q;m)=1/(1-m_{\mathrm{min}}/m)$ is the initial mass ratio function
when the primary mass is $m$ and $J(m)=[m_{\mathrm{min}}/m,1]$ is the initial
mass ratio range.

$\displaystyle n(t)=\frac{\text{the number of BBHs whose delay time is in
}[t,t+\delta t]}{\delta t},$ (9)

is the delay time distribution function and we set $\delta t=0.1$ Gyr. As long
as we use such a small value for $\delta t$, the results do not change.

For all parameter sets, we additionally sample the $N_{\mathrm{sam},J}=10^{6}$
binaries in the primary initial mass range $J=[150,600]M_{\odot}\subset I$ in
order to increase the number of samples efficiently and to obtain the more
precise results. In such a case, the merger rate density $\mathcal{R}(z)$ is
calculated as follows.

$\displaystyle\mathcal{R}(z(t))=\int_{0}^{t}$ $\displaystyle
f_{\mathrm{b}}\frac{\mathrm{SFR}(\tau)}{\langle
m\rangle}\left(A_{J}\frac{n_{I}^{(J)}(t-\tau)+n_{J}(t-\tau)}{N_{\mathrm{sam},I}^{(J)}+N_{\mathrm{sam},J}}+\right.$
$\displaystyle(1-A_{J})\left.\frac{n_{I}(t-\tau)-n_{I}^{(J)}(t-\tau)}{N_{\mathrm{sam},I}-N_{\mathrm{sam},I}^{(J)}}\right)~{}\mathrm{d}\tau,$
(10)

$N_{\mathrm{sam},I}=N_{\mathrm{sam}}$ is the number of samples for the
sampling in the primary initial mass range $I$.
$N_{\mathrm{sam},I}^{(J)}<N_{\mathrm{sam},I}$ is the number of samples whose
primary initial mass is in the range $J\subset I$ for the sampling in the
primary initial mass range $I$. $n_{I}(t)=n(t)$ and $n_{J}(t)$ are the delay
time distribution for the sampling in the primary initial mass range $I$ and
$J$, respectively. $n^{(J)}_{I}(t)<n_{I}(t)$ is the delay time distribution of
BBHs whose primary initial mass is in the range $J$ for the sampling in the
primary initial mass range $I$. Finally,

$\displaystyle A_{J}=\int_{J}f(m)~{}\mathrm{d}m<\int_{I}f(m)~{}\mathrm{d}m=1,$
(11)

is calculated from an IMF.

###  2.4 Parameter Sets

As written in Section 2.1 and Section 2.2.1, we use $\alpha\lambda=$ 0.01,
0.1, 0.5, 1.0, 2.0 and 5.0 for common envelope evolution, and $\gamma_{1}=$
0.0, 1.0, $\gamma_{2}=$ 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5 and 5.0, and
$m_{\mathrm{crit}}/M_{\odot}=$ 100, 200 and 300 for the Pop. III IMF.
Therefore, the number of IMFs in our study is $2\times 8\times 3=48$, and thus
the number of parameter sets in our study is $6\times 48=288$. For these sets,
we perform the binary population synthesis calculation and commpute the merger
rate density of Pop. III BBH mergers. Finally, we summarize our parameter sets
in Table 1.

Table 1: Parameter sets Parameter | Value
---|---
common envelope parameter, $\alpha\lambda$ | 0.01, 0.1, 0.5, 1.0, 2.0, 5.0
the slope of the lower mass region, $\gamma_{1}$ | 0.0, 1.0
the slope of the higher mass region, $\gamma_{2}$ | 1.5, 2.0, 2.5, $\ldots$, 4.5, 5.0
critical mass, $m_{\mathrm{crit}}$ | 100, 200, 300 $M_{\odot}$

##  3 Results

The BHs in the mass range 45–135 $M_{\odot}$ cannot be formed due to the pair-
instability effect (see Section 2.1). Therefore, the BBHs with IMBHs ($>100$
$M_{\odot}$) are classified into two types. In this paper, we call a BH whose
mass is lower than the lower edge of pair-instability mass gap (45
$M_{\odot}$) “low mass BH” and call a BH whose mass is higher than the upper
edge of pair-instability mass gap (135 $M_{\odot}$) “high mass BH”. We
describe the pair of “low mass BH” and “high mass BH” as “low mass $+$ high
mass” and describe the pair of “high mass BH”s as “high mass $+$ high mass”.
In Section 3.1 and 3.2, we show the primary BH mass, mass ratio, effective
spin , delay time and merger rate density of “high mass $+$ high mass” and
“low mass $+$ high mass”. We also discuss the dependence of these distribution
on common envelope parameter $\alpha\lambda$.

###  3.1 “high mass $+$ high mass”

####  3.1.1 formation channel

Figure 1: An example of “double common envelope channel” when
$\alpha\lambda=1.0$. “CHeB” and “nHe” stand for the core helium burning phase
and naked helium star, respectively. The dashed circles in phase (B) indicate
the Roche lobes of both stars. Figure 2: An example of “contact common
envelope channel” when $\alpha\lambda=1.0$. Figure 3: An example of “single
common envelope channel (2)” when $\alpha\lambda=0.01$. “ShHeB” stands for the
shell helium burning phase.

First, we introduce the dominant formation channel of “high mass $+$ high
mass” for each common envelope parameter $\alpha\lambda$. In our calculation,
when $\alpha\lambda=$ 0.5, 1.0, 2,0 and 5.0, “high mass $+$ high mass” BBHs
are formed through “double common envelope channel” and “contact common
envelope channel”. The definitions of these evolutionary channels are that a
binary system pass through one “double common envelope” and “contact common
envelope”, respectively (see Section 2.1). An example of “double common
envelope channel” and “contact common envelope channel” are shown in Figure 1
and 2 (for $\alpha\lambda=1.0$). When $\alpha\lambda=$ 0.1, “single common
envelope channel (2)” emerges and becomes dominant in addition to “double
common envelope channel” and “contact common envelope channel”. “Single common
envelope channel (2)” is the evolutionary channel that passes through one
“single common envelope” whose donor star is the secondary star at the ZAMS
time, i.e., the initially lighter star. When $\alpha\lambda=0.01$, only
“single common envelope channel (2)” can form a merging BBH. An example
channel is shown in Figure 3 (for $\alpha\lambda=0.01$). Almost all BBH formed
through “single common envelope channel (2)” are mass inverted.

In general, it is easier to shrink through “double common envelope channel”
and “contact common envelope channel” than through “single common envelope
channel (2)” because more orbital energy are used to strip the stellar
envelope through the former common envelope than the latter common envelope.
Therefore, if $\alpha\lambda$ is too small such as 0.01, the binary system
entering “double common envelope channel” and “contact common envelope
channel” coalesces after the common envelope phase, and thus they cannot form
a merging BBH. On the other hand, when $\alpha\lambda=$ 0.5, 1.0, 2.0, it is
impossible to shrink the orbital separation sufficiently to merge within a
Hubble time through “single common envelope channel (2)”. When
$\alpha\lambda=0.1$, both formation channels can occur.

####  3.1.2 primary BH mass distribution

Figure 4: The primary BH mass distributions of “high mass $+$ high mass” BBH
mergers which merge within a Hubble time. From the upper left to upper right
panel, $\alpha\lambda=$ 0.01, 0.1 and 0.5, and from the lower left to lower
right panel, $\alpha\lambda=$ 1.0, 2.0 and 5.0. Each line corresponds to each
IMF model. The normal lines and relatively transparent lines correspond to
$\gamma_{1}=0.0$ and 1.0, respectively. The red, orange, yellow, green, blue,
navy, purple and pink lines correspond to $\gamma_{2}=1.5$, 2.0, 2.5, 3.0,
3.5, 4.0, 4.5 and 5.0, respectively. The solid, dashed and dotted lines
correspond to $m_{\mathrm{crit}}=100M_{\odot}$, $200M_{\odot}$ and
$300M_{\odot}$, respectively.

The primary BH mass distributions of “high mass $+$ high mass” BBH mergers
which merge within a Hubble time are shown in Figure 4. For simplicity, we
only show models with $\gamma_{1}=0.0$, but models with $\gamma_{1}=1.0$ also
have a similar distribution. The most distinctive feature is that the primary
BH mass distributions have an upper limit. This upper limit decreases as the
common envelope parameter $\alpha\lambda$ increases. The reason for this is as
follows. If the ZAMS mass of Pop. III stars is larger than $\sim
600M_{\odot}$, they expand extremely, become a red super giant, and have a
convective envelope when they still in their MS phase. Therefore, if such a
massive Pop. III MS star fills its Roche lobe, the mass transfer may be
unstable, and the binary system enters a common envelope phase. After that,
the binary system always coalesces because the core–envelope structure of MS
star is not clear, and it will no longer evolve to a BBH. In order to avoid
this common envelope episode, the initial orbital separation needs to be large
enough.

Hereafter, in order to understand why there is an upper limit on the primary
ZAMS mass of merging “high mass $+$ high mass” BBHs, we consider the
dependence of the merging timescale of the binary system on the primary ZAMS
mass. The merging timescale $t$ is proportional to
$a^{4}m_{\mathrm{BH,p}}^{-1}m_{\mathrm{BH,s}}^{-1}m_{\mathrm{tot}}^{-1}$,
where $a$ is the orbital separation when BBH system is formed,
$m_{\mathrm{BH,p}}$ is the primary BH mass, $m_{\mathrm{BH,s}}$ is the
secondary BH mass and $m_{\mathrm{tot}}$ is the total mass of the BBH. Here,
if we fix the secondary BH mass, $m_{\mathrm{BH,s}}$, then the merging
timescale $t$ is roughly proportional to $a^{4}m_{\mathrm{BH,p}}^{-2}$. The
orbital separation, $a$, is proportional to the radius $r$ of a Pop. III giant
star heavier than $\simeq 600M_{\odot}$, and the primary BH mass,
$m_{\mathrm{BH,p}}$, is proportional to the ZAMS mass, $m$. Since there is the
relation $r\propto m^{0.6}$, $t\propto m^{0.4}$. Therefore, as the primary
ZAMS mass increases, the merging timescale also increases, and exceeds a
Hubble time at a critical ZAMS mass. That is why there is an upper limit on
the primary ZAMS mass of merging BBHs.

If the value of $\alpha\lambda$ is smaller, the efficiency of envelope
stripping through common envelope gets lower, and thus the binary system can
shrink more during the common envelope phase. Therefore, the more massive Pop.
III stars can merge within a Hubble time with the smaller $\alpha\lambda$, and
the upper limit of primary BH mass increases when the value of $\alpha\lambda$
decreases.

####  3.1.3 mass ratio distribution

Figure 5: The mass ratio distributions of “high mass $+$ high mass” BBH
mergers which merge within a Hubble time. From the upper left to upper right
panel, $\alpha\lambda=$ 0.01, 0.1 and 0.5, and from the lower left to lower
right panel, $\alpha\lambda=$1.0, 2.0 and 5.0. The legends are same as that of
Figure 4.

The mass ratio distributions of “high mass $+$ high mass” BBH mergers which
merge within a Hubble time are shown in Figure 5. When $\alpha\lambda\gtrsim$
0.5, “high mass $+$ high mass” BBH mergers with the mass ratio $>0.8$ and
$<0.8$ are formed through “contact common envelope channel” and “double common
envelope channel”, respectively. If the binary is initially symmetric (initial
mass ratio $\gtrsim 0.8$), the evolutionary timescale is similar to each other
and initially lighter star fills its Roche lobe while the initially heavier
star fills its Roche lobe and stable mass transfer proceeds. After that, a
“contact common envelope” occurs, only helium cores remain and the binary
system becomes a binary naked helium star. Since the helium core mass ratio is
roughly equal to initial mass ratio, the mass ratio of BBH formed through
“contact common envelope channel” is $\gtrsim 0.8$.

If the initial mass ratio is $\lesssim 0.8$, the evolutionary timescale of
both stars are somewhat different. Therefore, while the initially heavier star
fills its Roche lobe and stable mass transfer proceeds, the initially lighter
star does not yet expand significantly. After the initially lighter star
becomes the post MS star but before it fills its Roche lobe, the initially
heavier star evolves and its envelope becomes a convective, and then the
binary system enters a common envelope phase. Like “contact common envelope
channel”, the binary system becomes a binary naked helium star after the
common envelope phase, and thus the mass ratio of BBHs formed through “double
common envelope channel” is roughly same as initial mass ratio. If the initial
mass ratio gets smaller, the wider orbital separation is needed for “double
common envelope channel”. The reason is as follows. In order to cause a
“double common envelope”, the initially heavier star fills its Roche lobe
after the initially lighter star becomes a giant. If the initial mass ratio
gets smaller, the difference between each evolutionary timescale gets larger,
and thus the initially lighter star takes a relatively longer period to become
a giant star. In order for the initially heavier star not to fill its Roche
lobe during this relatively longer period, the wider initial orbital
separation is needed. If the initial orbital separation is too wide, the delay
time exceeds a Hubble time. In this way the minimum mass ratio of BBHs formed
through “double common envelope channel” is determined. If the $\alpha\lambda$
gets larger, the orbital separation just after the common envelope gets wider
(see Equation 2.1), and thus the minimum mass ratio gets larger. When
$\alpha\lambda=0.1$, unlike the $\alpha\lambda\gtrsim 0.5$ case, there are
“high mass $+$ high mass” BBHs formed through “single common envelope channel
(2)”, but the mechanism determining the minimum value of mass ratio $q$ is the
same. When $\alpha\lambda=0.01$, all “high mass $+$ high mass” BBHs which
merge within a Hubble time are formed through “single common envelope channel
(2)” (see Figure 3). Almost all “high mass $+$ high mass” BBHs are mass
inverted for $\alpha\lambda=0.01$, and have the mass ratio $q$ roughly equal
to the ratio of the helium core mass of initially heavier star to the helium
core mass of initially lighter star which gains significant mass from the
initially heavier star (phase B in Figure 3). Therefore, as the initial mass
ratio increases, the mass ratio $q$ decreases. Since the initial mass ratio
has an upper limit ($=1$), the mass ratio $q$ also has lower limit.

####  3.1.4 effective spin distribution

Figure 6: The effective spin distribution of “high mass $+$ high mass” BBH
mergers which merge within a Hubble time. From the upper left to upper right
panel, $\alpha\lambda=$ 0.01, 0.1 and 0.5, and from the lower left to lower
right panel, $\alpha\lambda=$ 1.0, 2.0 and 5.0. The legends are same as that
of Figure 4.

The effective spin distributions of “high mass $+$ high mass” BBH mergers
which merge within a Hubble time are shown in Figure 6. The effective spin,

$\displaystyle\chi_{\mathrm{eff}}=\frac{m_{\mathrm{BH,p}}\vec{\chi}_{\mathrm{BH,p}}+m_{\mathrm{BH,s}}\vec{\chi}_{\mathrm{BH,s}}}{m_{\mathrm{BH,p}}+m_{\mathrm{BH,s}}}\cdot\frac{\vec{L}_{\mathrm{orb}}}{|\vec{L}_{\mathrm{orb}}|},$
(12)

is the mass weighted average of the dimensionless spin components
perpendicular to the orbital plane. Here, $m_{\mathrm{BH,p}}$ and
$m_{\mathrm{BH,s}}$ are the primary and secondary BH mass,
$\vec{\chi}_{\mathrm{BH,p}}$ and $\vec{\chi}_{\mathrm{BH,s}}$ are the primary
and secondary BH spin vector and $\vec{L}_{\mathrm{orb}}$ is the orbital
angular momentum vector. Since we do not assume that the kick is imparted
during the BH formation in this study, BHs always aligned with the orbital
angular momentum vector. Therefore, the effective spin can be written as

$\displaystyle\chi_{\mathrm{eff}}$
$\displaystyle=\frac{m_{\mathrm{BH,p}}\chi_{\mathrm{BH,p}}+m_{\mathrm{BH,s}}\chi_{\mathrm{BH,s}}}{m_{\mathrm{BH,p}}+m_{\mathrm{BH,s}}}$
$\displaystyle=\frac{\chi_{\mathrm{BH,p}}+q\chi_{\mathrm{BH,s}}}{1+q}.$ (13)

When $\alpha\lambda\gtrsim$ 0.5, the “high mass $+$ high mass” BBH mergers are
formed through “contact common envelope channel” or “double common envelope
channel”. Through these channels, both stars lose their hydrogen envelope
during the common envelope phase and lose spin angular momenta at the same
time. Therefore, both BHs have very low spin ($\sim 0$) and the effective spin
distribution sharply peaks at 0. If the orbital separation after the common
envelope phase is too short, the naked helium star can be spun up by tidal
interaction. In such a case, BHs have a non-zero spin and the effective spin
of BBH mergers have a non-zero value. If $\alpha\lambda$ gets larger, the
orbital separation after the common envelope gets relatively wider, and the
degree of spin-up by tidal interaction gets smaller. Therefore, the maximum
effective spin gets smaller when $\alpha\lambda$ gets larger. When
$\alpha\lambda<$ 0.1, BBH mergers formed through “single common envelope (2)”
emerges. Through this channel, the progenitor of one of the BH lose its
hydrogen envelope by the common envelope, and thus this BH have a very low
spin ($\sim 0$). On the other hand, the progenitor of the other BH did not
undergo a significant mass loss process, and thus the spin of this BH is high
($\sim 1$). Therefore, the effective spin is written as

$\displaystyle\chi_{\mathrm{eff}}\sim\frac{q}{1+q}\in[0,0.5].$ (14)

Therefore, the effective spin distribution of “high mass $+$ high mass” BBH
mergers assemble around non-zero value, 0.35–0.5 for $\alpha\lambda=$0.01 and
0.25–0.5 for $\alpha\lambda=$0.1.

####  3.1.5 delay time distribution

Figure 7: The delay time distributions of “high mass $+$ high mass” BBH
mergers which merge within a Hubble time. From the upper left to upper right
panel, $\alpha\lambda=$ 0.01, 0.1 and 0.5, and from the lower left to lower
right panel, $\alpha\lambda=$ 1.0, 2.0 and 5.0. The legends are same as that
of Figure 4.

The delay time distributions of “high mass $+$ high mass” BBH mergers which
merge within a Hubble time are shown in Figure 7. The delay time
$t_{\mathrm{d}}$ is roughly proportional to
$a^{4}m_{\mathrm{BH,p}}^{-1}m_{\mathrm{BH,s}}^{-1}m_{\mathrm{tot}}^{-1}$. If
the $\alpha\lambda$ gets larger, the orbital separation after the common
envelope gets wider, and thus the delay time gets longer. Therefore, the
minimum delay time increases when $\alpha\lambda$ gets larger, and exceeds 10
and 100 Myr when $\alpha\lambda=$2.0 and 5.0. When the value of
$\alpha\lambda$ is too small such as 0.01, the orbital separation after the
common envelope is too short, and thus the maximum delay time is lower than
100 Myr.

####  3.1.6 merger rate density

Figure 8: The redshift evolution of merger rate density of “high mass $+$ high
mass” BBH mergers. From the upper left to upper right panel, $\alpha\lambda=$
0.01, 0.1 and 0.5, and from the lower left to lower right panel,
$\alpha\lambda=$ 1.0, 2.0 and 5.0. The legends are same as that of Figure 4.
The lines at the bottom of each panel are the redshift evolution of BBH merger
rate density. The relatively thin lines at the top of each panel indicate the
binary formation rate [$\mathrm{Gpc}^{-3}$ $\mathrm{yr}^{-1}$]. The binary
formation rate is the number density of Pop. III binaries born at each
redshift and is calculated from star formation rate.

The redshift evolution of “high mass $+$ high mass” BBH merger rate density is
shown in Figure 8. The merger rate densities roughly trace the star formation
rate history (de Souza et al., 2011), and the peaks of merger rate densities
are located near the redshift $z\sim 10$. However, when $\alpha\lambda=$ 2.0,
the minimum delay time exceeds 10 Myr, hence the peak of merger rate density
is shifted toward low redshift compared to that of $\alpha\lambda\lesssim$
1.0. When $\alpha\lambda=$ 5.0, the minimum delay time exceeds 100 Myr, hence
the peak of merger rate density is significantly shifted toward low redshift
compared to other cases. Furthermore, the merger rate density becomes
significantly smaller at the redshift $z\gtrsim$ 15–20. When
$\alpha\lambda=0.01$, the maximum delay time is so short that the merger rate
density rapidly decreases after the end of Pop. III star formation. Therefore,
the merger rate densities are strictly 0 when $z\lesssim 3$.

###  3.2 low mass $+$ high mass

####  3.2.1 formation channel

Figure 9: The example channel of “low mass $+$ high mass” when
$\alpha\lambda=1.0$.

The almost all “low mass $+$ high mass” BBH mergers are formed through “single
common envelope channel (2)”, and example formation channel is shown in Figure
9. Both stars would cause the pair-instability supernova if they were not in
the binary system, but thanks to the binary interaction such as stable mass
transfer one of the star gains the mass and becomes a “high mass BH”. On the
other hand, the other star loses the mass and becomes a “low mass BH”. In
order to significantly transfer the stellar mass to avoid the pair-instability
supernova, the mass transfer rate needs to be high and the initial orbital
separation needs to be short. If the initial orbital separation is short and
$\alpha\lambda$ is smaller, the binary system coalesces after the common
envelope phase, and thus “low mass $+$ high mass” BBH mergers cannot be formed
when $\alpha\lambda\lesssim$ 0.5.

####  3.2.2 primary BH mass distribution

Figure 10: The primary BH mass distribution of “low mass $+$ high mass” BBH
mergers which merge within a Hubble time. From the left to right panel,
$\alpha\lambda=$ 1.0, 2.0 and 5.0. The legends are same as that of Figure 4.

The primary BH mass distributions of “low mass $+$ high mass” BBH mergers
which merge within a Hubble time are shown in Figure 10. The primary BH mass
distributions are almost same for the different $\alpha\lambda$ values. Since
both stars are initially 100–200 $M_{\odot}$ and the progenitor of primary BH
can only increase its mass up to $\sim 300M_{\odot}$, the primary BH mass is
from 135 $M_{\odot}$ to $\sim 200M_{\odot}$. 135 $M_{\odot}$ is the upper end
of pair-instability mass gap in this study.

####  3.2.3 mass ratio distribution

Figure 11: The mass ratio distribution of “low mass $+$ high mass” BBH mergers
which merge within a Hubble time. From the left to right panel,
$\alpha\lambda=$ 1.0, 2.0 and 5.0. The legends are same as that of Figure 4.

The mass ratio distributions of “low mass $+$ high mass” BBH mergers which
merge within a Hubble time are shown in Figure 11. The mass ratio
distributions are roughly same for the different $\alpha\lambda$ values. The
maximum mass ratio is equal to the ratio of lower end of pair-instability mass
gap ($=45M_{\odot}$) to the upper end ($=135M_{\odot}$), so that 0.33. When
the $\alpha\lambda=$ 2.0 and 5.0, the orbital separation does not shrink much
through common envelope phase compared to when $\alpha\lambda=1.0$. Therefore,
the initially shorter orbital separation are allowed for $\alpha\lambda=$ 2.0
and 5.0, and in these cases the stellar mass is transferred more significantly
(phase B and D in Figure 9). That is why minimum mass ratio gets lower when
$\alpha\lambda=$ 2.0 and 5.0.

####  3.2.4 effective spin distribution

Figure 12: The effective spin distribution of “low mass $+$ high mass” BBH
mergers which merge within a Hubble time. From the left to right panel,
$\alpha\lambda=$ 1.0, 2.0 and 5.0. The legends are same as that of Figure 4.

The effective spin distributions of “low mass $+$ high mass” BBH mergers which
merge within a Hubble time are shown in Figure 12. The progenitor of secondary
BH loses its spin angular momentum through significant mass transfer process
(phase B and D in Figure 9), and thus the BH also has very low spin ($\sim
0$). The progenitor of the primary BH also loses its spin angular momentum
during the common envelope phase (phase G in Figure 9). However, the orbital
separation after the common envelope is so short that the progenitor of the
primary BH is highly spun up by tidal interaction. Thus, the spin of primary
BH is high. As the value of $\alpha\lambda$ increases, the degree of spin-up
decreases, and thus the effective spin distribution shifted toward lower
value.

####  3.2.5 delay time distribution

Figure 13: The delay time distribution of “low mass $+$ high mass” BBH mergers
which merge within a Hubble time. From the left to right panel,
$\alpha\lambda=$ 1.0, 2.0 and 5.0. The legends are same as that of Figure 4.

The delay time distribution of “low mass $+$ high mass” BBH mergers which
merge within a Hubble time are shown in Figure 13. As the value of
$\alpha\lambda$ increases, the orbital separation after the common envelope
widens, and thus the delay time distribution shifted toward the higher value.
Since the initial orbital separation is short in order to form “low mass $+$
high mass” BBHs, the delay time also becomes short, and does not exceeds
$10^{3}$ Myr for $\alpha\lambda=$ 1.0, 2.0 and 5.0.

####  3.2.6 merger rate density

Figure 14: The redshift evolution of merger rate density of “low mass $+$ high
mass” BBH. From the left to right, for $\alpha\lambda=$ 1.0, 2.0 and 5.0. The
legends are same as that of Figure 4. The lines at the bottom of each panel
are the redshift evolution of BBH merger rate density. The relatively thin
lines at the top of each panel indicate the binary formation rate
[$\mathrm{Gpc}^{-3}$ $\mathrm{yr}^{-1}$].

The redshift evolution of “low mass $+$ high mass” BBH merger rate density is
shown in Figure 14. The merger rate densities roughly trace the star formation
rate history (de Souza et al., 2011), and the peaks of merger rate densities
are located near the redshift $z\sim 10$. However, when $\alpha\lambda=$ 5.0,
the minimum delay time exceeds 10 Myr, hence the peak of merger rate density
is shifted toward low redshift compared to that of $\alpha\lambda\lesssim$
2.0. The maximum delay time is so short that the merger rate density rapidly
decreases after the end of Pop. III star formation. Therefore, the merger rate
densities are strictly 0 when $z\lesssim 3$.

##  4 Discussion

So far, the BBH mergers with IMBHs have not been discovered by GW observation.
Thus, the upper limit on the merger rate density of the BBH mergers with IMBHs
is obtained, and the value is $\sim 0.01$ $\mathrm{Gpc}^{-3}$
$\mathrm{yr}^{-1}$ (The LIGO Scientific Collaboration et al., 2021a; Ezquiaga
& Holz, 2021). In Section 4.1, we compare the merger rate density of BBHs with
IMBHs with this upper limit and investigate what parameter set meets the
criterion obtained from the non-detection. In Section 4.2, using these
parameter sets, we infer the detection rate by future GW detectors.
Furthermore, we show the likelihood ratio for each model assuming a BBH with
IMBHs is observed and explain how we can limit the value of $\alpha\lambda$
and IMF models.

###  4.1 comparison with the current GW observation

The merger rate density of BBHs with IMBHs at the redshift $z=0$ are shown in
Table 2–6 for $\alpha\lambda=$ 0.1–5.0. The value lower than the upper limit,
0.01 $\mathrm{Gpc}^{-3}$ $\mathrm{yr}^{-1}$, is written in bold. The BBHs with
IMBHs are classified as “high mass $+$ high mass” BBHs and “low mass $+$ high
mass” BBHs. The merger rate density of “low mass $+$ high mass” BBHs is 0 at
the redshift $z=0$ (see Section 3.2.6), hence the merger rate density shown in
the Table 2–6 are that of “high mass $+$ high mass” BBH mergers. When
$\alpha\lambda=0.01$, the merger rate density of “high mass $+$ high mass”
BBHs is also 0 at the redshift $z=0$ (see Section 3.1.6), hence the merger
rate density of BBHs with IMBHs is strictly 0.

Since the delay time is so short that the Pop. III BBHs with IMBHs cannot
merge at the local Universe and the merger rate is 0 when
$\alpha\lambda=0.01$, any IMF is allowed. Therefore it is impossible to impose
any restrictions on Pop. III IMF. When $\alpha\lambda=0.1$, the contribution
from the Pop. III star heavier than $\sim 400M_{\odot}$ is significant.
Therefore, an IMF with the sufficiently large $\gamma_{2}$ is allowed for
$\gamma_{1}=$ 0 and 1 and $m_{\mathrm{crit}}=$ 100–300 $M_{\odot}$ (see Table
2). When $\alpha\lambda\gtrsim 0.5$, the contribution from the Pop. III star
$\sim 200M_{\odot}$ is significant. In this study, the Pop. III IMF is
relatively top heavy below the critical mass $m_{\mathrm{crit}}$. Therefore,
when $m_{\mathrm{crit}}=$ 200 and 300 $M_{\odot}$, it is impossible to
decrease the merger rate density by steeping the IMF slope (i.e., increasing
the value of $\gamma_{2}$) above the critical mass, and thus the Pop. III IMF
with $m_{\mathrm{crit}}=$ 200 and 300 $M_{\odot}$ are not allowed (see Table
3–6).

Table 2: The merger rate density [$\mathrm{Gpc}^{-3}$ $\mathrm{yr}^{-1}$] of BBHs with IMBHs at the redshift $z=0$ for $\alpha\lambda=0.1$. The value lower than the upper limit, 0.01 $\mathrm{Gpc}^{-3}$ $\mathrm{yr}^{-1}$, is written in bold. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | 2.44e-02 | 1.09e-02 | 9.83e-03 | 1.23e-02 | 3.12e-03 | 4.77e-03 | 2.40e-04 | 3.79e-04
0.0 | 200 | 1.82e-02 | 2.33e-02 | 1.42e-02 | 2.29e-02 | 1.34e-02 | 6.04e-03 | 8.70e-03 | 2.38e-03
0.0 | 300 | 2.29e-02 | 2.74e-02 | 2.91e-02 | 2.60e-02 | 2.30e-02 | 1.50e-02 | 9.06e-03 | 9.98e-03
1.0 | 100 | 2.69e-02 | 1.36e-02 | 1.43e-02 | 6.08e-03 | 8.34e-03 | 1.50e-04 | 8.01e-04 | 5.68e-05
1.0 | 200 | 2.42e-02 | 2.08e-02 | 1.86e-02 | 1.60e-03 | 6.40e-03 | 2.74e-03 | 2.41e-03 | 1.79e-03
1.0 | 300 | 2.18e-02 | 1.49e-02 | 1.91e-02 | 1.29e-02 | 1.08e-02 | 8.18e-03 | 9.28e-03 | 7.72e-03
Table 3: Same as Table 2 but for $\alpha\lambda=0.5$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | 3.97e-02 | 3.77e-02 | 3.79e-02 | 3.38e-02 | 1.38e-02 | 8.73e-03 | 6.69e-03 | 3.16e-03
0.0 | 200 | 4.37e-02 | 5.06e-02 | 5.33e-02 | 4.76e-02 | 3.99e-02 | 4.95e-02 | 3.81e-02 | 3.04e-02
0.0 | 300 | 5.43e-02 | 6.29e-02 | 7.86e-02 | 7.33e-02 | 8.46e-02 | 7.89e-02 | 7.03e-02 | 9.14e-02
1.0 | 100 | 3.80e-02 | 3.95e-02 | 2.67e-02 | 1.93e-02 | 9.67e-03 | 7.76e-03 | 4.48e-03 | 2.09e-03
1.0 | 200 | 4.28e-02 | 3.89e-02 | 3.57e-02 | 4.34e-02 | 3.19e-02 | 2.36e-02 | 2.10e-02 | 1.31e-02
1.0 | 300 | 3.95e-02 | 4.72e-02 | 4.85e-02 | 6.53e-02 | 6.30e-02 | 4.98e-02 | 5.99e-02 | 4.95e-02
Table 4: Same as Table 2 but for $\alpha\lambda=1.0$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | 5.32e-02 | 5.30e-02 | 4.16e-02 | 3.03e-02 | 2.51e-02 | 1.32e-02 | 1.03e-02 | 4.85e-03
0.0 | 200 | 5.68e-02 | 7.11e-02 | 6.59e-02 | 6.69e-02 | 6.61e-02 | 6.22e-02 | 5.57e-02 | 3.76e-02
0.0 | 300 | 5.92e-02 | 7.65e-02 | 9.86e-02 | 1.08e-01 | 1.16e-01 | 1.21e-01 | 1.17e-01 | 1.19e-01
1.0 | 100 | 4.19e-02 | 4.32e-02 | 3.90e-02 | 2.22e-02 | 1.59e-02 | 1.02e-02 | 6.22e-03 | 3.90e-03
1.0 | 200 | 5.04e-02 | 5.36e-02 | 3.86e-02 | 4.56e-02 | 4.56e-02 | 3.35e-02 | 3.69e-02 | 3.32e-02
1.0 | 300 | 4.85e-02 | 6.47e-02 | 7.73e-02 | 7.43e-02 | 7.41e-02 | 6.96e-02 | 7.63e-02 | 7.73e-02
Table 5: Same as Table 2 but for $\alpha\lambda=2.0$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | 4.28e-02 | 4.95e-02 | 4.76e-02 | 4.19e-02 | 2.05e-02 | 1.50e-02 | 1.00e-02 | 5.87e-03
0.0 | 200 | 4.92e-02 | 6.40e-02 | 7.69e-02 | 6.86e-02 | 6.32e-02 | 6.37e-02 | 5.83e-02 | 4.89e-02
0.0 | 300 | 5.19e-02 | 7.84e-02 | 8.36e-02 | 1.06e-01 | 1.13e-01 | 1.20e-01 | 1.13e-01 | 1.33e-01
1.0 | 100 | 4.27e-02 | 4.21e-02 | 3.67e-02 | 2.33e-02 | 1.79e-02 | 1.10e-02 | 6.06e-03 | 3.30e-03
1.0 | 200 | 5.44e-02 | 6.12e-02 | 5.09e-02 | 5.03e-02 | 4.99e-02 | 4.07e-02 | 4.03e-02 | 2.92e-02
1.0 | 300 | 5.48e-02 | 6.41e-02 | 6.45e-02 | 7.90e-02 | 7.25e-02 | 7.56e-02 | 8.26e-02 | 8.28e-02
Table 6: Same as Table 2 but for $\alpha\lambda=5.0$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | 3.32e-02 | 4.47e-02 | 3.36e-02 | 3.16e-02 | 1.90e-02 | 9.64e-03 | 8.18e-03 | 6.29e-03
0.0 | 200 | 3.59e-02 | 4.99e-02 | 6.03e-02 | 5.30e-02 | 5.99e-02 | 5.56e-02 | 5.61e-02 | 4.78e-02
0.0 | 300 | 4.39e-02 | 5.44e-02 | 7.18e-02 | 8.86e-02 | 1.03e-01 | 1.07e-01 | 1.01e-01 | 1.08e-01
1.0 | 100 | 3.73e-02 | 3.90e-02 | 3.33e-02 | 2.28e-02 | 1.39e-02 | 1.06e-02 | 6.21e-03 | 2.76e-03
1.0 | 200 | 3.12e-02 | 3.67e-02 | 4.56e-02 | 4.11e-02 | 4.06e-02 | 3.71e-02 | 2.95e-02 | 2.88e-02
1.0 | 300 | 3.96e-02 | 4.88e-02 | 4.80e-02 | 6.68e-02 | 7.08e-02 | 8.16e-02 | 7.09e-02 | 9.28e-02

###  4.2 future observation

####  4.2.1 detection rate

Figure 15: The horizon redshift of future GW detectors. The blue, red, green, orange and yellow lines correspond to advanced LIGO (O4), Einstein telescope, LISA, Cosmic Explorer and TianQin respectively. Table 7: The detection rate [$\mathrm{yr}^{-1}$] of BBHs with IMBHs for $\alpha\lambda=0.01$. The detection rate by advanced LIGO (O4), Einstein telescope and LISA are written in roman, italic and bold font, respectively. If an IMF model does not meet the upper limit from the non-detection, only “-” is written in the cell. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 |  | 0.0
---
23.6
213.8
| 0.0
---
24.0
140.7
| 0.0
---
22.1
79.3
| 0.0
---
22.1
79.3
| 0.0
---
11.5
21.8
| 0.0
---
6.0
5.4
| 0.0
---
4.6
1.4
| 0.0
---
3.4
1.0
0.0 | 200 |  | 0.0
---
25.3
240.5
| 0.0
---
30.6
192.4
| 0.0
---
34.0
115.2
| 0.0
---
35.1
87.0
| 0.0
---
36.7
54.1
| 0.0
---
36.9
27.5
| 0.0
---
30.6
15.4
| 0.0
---
29.6
7.8
0.0 | 300 |  | 0.0
---
25.1
263.3
| 0.0
---
30.2
231.9
| 0.0
---
33.9
170.3
| 0.0
---
42.9
125.6
| 0.0
---
49.6
93.5
| 0.0
---
58.3
64.6
| 0.0
---
58.8
38.5
| 0.0
---
57.4
32.4
1.0 | 100 |  | 0.0
---
18.9
197.3
| 0.0
---
22.0
138.5
| 0.0
---
16.8
71.2
| 0.0
---
12.4
32.6
| 0.0
---
11.2
17.1
| 0.0
---
6.7
3.3
| 0.0
---
3.4
3.4
| 0.0
---
1.6
0.0
1.0 | 200 |  | 0.0
---
25.2
228.0
| 0.0
---
27.0
160.4
| 0.0
---
25.5
104.8
| 0.0
---
28.9
46.6
| 0.0
---
23.5
37.6
| 0.0
---
29.5
23.2
| 0.0
---
21.1
9.1
| 0.0
---
19.1
6.4
1.0 | 300 |  | 0.0
---
20.8
240.3
| 0.0
---
26.6
174.9
| 0.0
---
27.7
137.9
| 0.0
---
32.8
85.1
| 0.0
---
38.1
68.0
| 0.0
---
39.9
39.5
| 0.0
---
42.4
30.0
| 0.0
---
43.8
15.9
Table 8: Same as Table 7 but for $\alpha\lambda=0.1$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | - | - |  | 0.32
---
604
76
- |  | 0.10
---
308
15
| 0.052
---
195
5.4
| 0.033
---
118
1.5
| 0.010
---
71
0.68
0.0 | 200 | - | - | - | - | - |  | 0.20
---
785
21
| 0.14
---
693
13
| 0.082
---
587
8.3
0.0 | 300 | - | - | - | - | - | - |  | 0.34
---
1600
31
| 0.50
---
1646
21
1.0 | 100 | - | - | - |  | 0.13
---
350
25
| 0.057
---
222
6.5
| 0.017
---
137
0.77
| 0.012
---
84
5.1
| 0.005
---
47
0.014
1.0 | 200 | - | - | - |  | 0.30
---
713
69
| 0.22
---
623
28
| 0.20
---
543
22
| 0.13
---
456
7.5
| 0.069
---
378
6.9
1.0 | 300 | - | - | - | - | - |  | 0.32
---
1077
39
| 0.26
---
1083
23
| 0.28
---
1085
17
Table 9: Same as Table 7 but for $\alpha\lambda=0.5$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | - | - | - | - | - |  | 1.6
---
212
1.6
| 0.8
---
134
0.7
| 0.5
---
83
0.4
0.0 | 200 | - | - | - | - | - | - | - | -
0.0 | 300 | - | - | - | - | - | - | - | -
1.0 | 100 | - | - | - | - |  | 1.5
---
240
1.9
| 1.0
---
153
0.9
| 0.6
---
90
0.5
| 0.4
---
56
0.3
1.0 | 200 | - | - | - | - | - | - | - | -
1.0 | 300 | - | - | - | - | - | - | - | -
Table 10: Same as Table 7 but for $\alpha\lambda=1.0$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | - | - | - | - | - | - | - |  | 0.6
---
143
0.2
0.0 | 200 | - | - | - | - | - | - | - | -
0.0 | 300 | - | - | - | - | - | - | - | -
1.0 | 100 | - | - | - | - | - | - |  | 0.7
---
136
0.2
| 0.5
---
93
0.1
1.0 | 200 | - | - | - | - | - | - | - | -
1.0 | 300 | - | - | - | - | - | - | - | -
Table 11: Same as Table 7 but for $\alpha\lambda=2.0$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | - | - | - | - | - | - | - |  | 0.65
---
164
0.12
0.0 | 200 | - | - | - | - | - | - | - | -
0.0 | 300 | - | - | - | - | - | - | - | -
1.0 | 100 | - | - | - | - | - | - |  | 0.71
---
147
0.12
| 0.41
---
104
0.02
1.0 | 200 | - | - | - | - | - | - | - | -
1.0 | 300 | - | - | - | - | - | - | - | -
Table 12: Same as Table 7 but for $\alpha\lambda=5.0$. | $\gamma_{2}$
---|---
$\gamma_{1}$ |  $m_{\mathrm{crit}}$ [$M_{\odot}$] | 1.5 | 2.0 | 2.5 | 3.0 | 3.5 | 4.0 | 4.5 | 5.0
0.0 | 100 | - | - | - | - | - |  | 1.28
---
216
0.08
| 0.84
---
145
0.05
| 0.53
---
133
0.03
0.0 | 200 | - | - | - | - | - | - | - | -
0.0 | 300 | - | - | - | - | - | - | - | -
1.0 | 100 | - | - | - | - | - | - |  | 0.56
---
114
0.03
| 0.37
---
85.3
0.02
1.0 | 200 | - | - | - | - | - | - | - | -
1.0 | 300 | - | - | - | - | - | - | - | -

In this section, we compute the detection rate by the second generation GW
detector, advanced LIGO (O4), third generation detector, Einstein telescope
and the space-borne detector, LISA for the models which meet the upper limit
from the non-detection. Figure 15 shows the horizon of each detector. The
detection rate of BBHs with IMBHs by advanced LIGO (O4), Einstein telescope
and LISA are shown in Table 7–12 for $\alpha\lambda=$ 0.01–5.0. The detection
rate by advanced LIGO (O4), Einstein telescope and LISA are written in roman,
italic and bold font, respectively.

Next, we choose the IMF model with the highest detection rate for each
detector and for each $\alpha\lambda$, and show the detection rate as a
function of mass. If we choose different IMF model, the mass distribution does
not change significantly. Therefore, it is sufficient to show the mass
distribution of one IMF model for each detector and for each $\alpha\lambda$.

(a) $\alpha\lambda=0.01$

(b) $\alpha\lambda=0.1$

(c) $\alpha\lambda=0.5$

(d) $\alpha\lambda=1.0$

(e) $\alpha\lambda=2.0$

(f) $\alpha\lambda=5.0$

Figure 16: The detection rate by advanced LIGO for each $\alpha\lambda$. Note
that not all panels have the same color scale.

(a) $\alpha\lambda=0.01$

(b) $\alpha\lambda=0.1$

(c) $\alpha\lambda=0.5$

(d) $\alpha\lambda=1.0$

(e) $\alpha\lambda=2.0$

(f) $\alpha\lambda=5.0$

Figure 17: The detection rate by Einstein telescope for each $\alpha\lambda$.
Note that not all panels have the same color scale.

(a) $\alpha\lambda=0.01$

(b) $\alpha\lambda=0.1$

(c) $\alpha\lambda=0.5$

(d) $\alpha\lambda=1.0$

(e) $\alpha\lambda=2.0$

(f) $\alpha\lambda=5.0$

Figure 18: The detection rate by LISA for each $\alpha\lambda$. Note that not
all panels have the same color scale.

The detection rate by advanced LIGO (O4) are shown in Figure 16. The IMF
models with the highest detection rate by advanced LIGO (O4) are
$(\gamma_{1}=0.0,\gamma_{2}=1.5,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=0.0,\gamma_{2}=5.0,m_{\mathrm{crit}}=300M_{\odot})$,
$(\gamma_{1}=0.0,\gamma_{2}=4.0,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=1.0,\gamma_{2}=4.5,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=1.0,\gamma_{2}=4.5,m_{\mathrm{crit}}=100M_{\odot})$ and
$(\gamma_{1}=0.0,\gamma_{2}=4.0,m_{\mathrm{crit}}=100M_{\odot})$ for
$\alpha\lambda=$0.01, 0.1, 0.5, 1.0, 2.0 and 5.0. When the value of
$\alpha\lambda$ is small, the delay time is also short. Therefore, BBHs with
IMBHs are hard to be detected by the advanced LIGO with the horizon $z\sim 1$.
On the other hand, when $\alpha\lambda$ is larger than 0.5, BBH with IMBHs
merges within $z\sim 1$, and the advanced LIGO have a potential to detect the
BBHs with total mass $\sim 300M_{\odot}$ and mass ratio $\sim 1$.

The detection rate by Einstein telescope are shown in Figure 17. The IMF
models with the highest detection rate by Einstein telescope are
$(\gamma_{1}=0.0,\gamma_{2}=4.5,m_{\mathrm{crit}}=300M_{\odot})$,
$(\gamma_{1}=0.0,\gamma_{2}=5.0,m_{\mathrm{crit}}=300M_{\odot})$,
$(\gamma_{1}=1.0,\gamma_{2}=3.5,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=0.0,\gamma_{2}=5.0,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=0.0,\gamma_{2}=5.0,m_{\mathrm{crit}}=100M_{\odot})$ and
$(\gamma_{1}=0.0,\gamma_{2}=4.0,m_{\mathrm{crit}}=100M_{\odot})$ for
$\alpha\lambda=$0.01, 0.1, 0.5, 1.0, 2.0 and 5.0. Einstein telescope can
detect a BBH merger up to $z\sim$ 10–100 (see Figure 15). The total mass and
mass ratio of observable “high mass $+$ high mass” BBHs are $\sim
300M_{\odot}$ and $\sim 1$, respectively. When $\alpha\lambda$ is larger than
1, “low mass $+$ high mass” BBH mergers can occur (see Section refsubsec: lh),
and can be detected by Einstein telescope. Their total mass and mass ratio are
$\sim 200M_{\odot}$ and $\sim 0.2$, respectively. The detection rate by Cosmic
Explorer is roughly same as that by Einstein telescope.

The detection rate by LISA are shown in Figure 18. The IMF models with the
highest detection rate by LISA are
$(\gamma_{1}=0.0,\gamma_{2}=1.5,m_{\mathrm{crit}}=300M_{\odot})$,
$(\gamma_{1}=0.0,\gamma_{2}=2.5,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=1.0,\gamma_{2}=3.5,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=1.0,\gamma_{2}=4.5,m_{\mathrm{crit}}=100M_{\odot})$,
$(\gamma_{1}=1.0,\gamma_{2}=4.5,m_{\mathrm{crit}}=100M_{\odot})$ and
$(\gamma_{1}=0.0,\gamma_{2}=4.0,m_{\mathrm{crit}}=100M_{\odot})$ for
$\alpha\lambda=$0.01, 0.1, 0.5, 1.0, 2.0 and 5.0. LISA can detect BBHs with
total mass $\sim 1000M_{\odot}$ up to $z\sim 1000$, but is hard to detect BBHs
with total mass $\sim 100M_{\odot}$ (see Figure 15). Since the maximum mass
decreases when $\alpha\lambda$ increases (see Section 3.1.2), BBHs with IMBHs
are hard to be detected by LISA when $\alpha\lambda$ is larger than 0.5. On
the other hand, when $\alpha\lambda$ is smaller than 0.1, BBHs with total mass
$\sim 1000M_{\odot}$ can form and merge. LISA can detect such very massive BBH
mergers. The detection rate by Tianqin is roughly same as that by LISA.

####  4.2.2 constraint on the models

(a) advanced LIGO (O4)

(b) LISA

Figure 19: The likelihood ratio of $\alpha\lambda$ for each detector.

(a) advanced LIGO (O4)

(b) LISA

Figure 20: The likelihood ratio of $\gamma_{1}$ for each detector.

(a) advanced LIGO (O4)

(b) LISA

Figure 21: The likelihood ratio of $\gamma_{2}$ for each detector.

(a) advanced LIGO (O4)

(b) LISA

Figure 22: The likelihood ratio of $m_{\mathrm{crit}}$ for each detector.

From Table 7–12, Einstein telescope will be able to detect several tens to
thousands BBHs with IMBHs in the future. Therefore, we may be able to impose
restrictions on models by using mass distribution. On the other hand, advanced
LIGO (O4) and LISA may only see one to a few BBHs. Thus, in this section, in
order to impose restrictions on the model when only one BBH with IMBHs is
detected, we compute and show the likelihood ratio for each parameter
($\alpha\lambda$, $\gamma_{1}$, $\gamma_{2}$ and $m_{\mathrm{crit}}$) assuming
that one BBH with IMBHs with a certain total mass is observed by a certain
detector. We simply assume that IMF models occurs with equal probability and
$\alpha\lambda$ have a probability density,
$p(\log(\alpha\lambda))\propto\mathrm{const}$.

Figure 19 shows the likelihood ratio to the maximum likelihood for
$\alpha\lambda$. For example, from Figure 19 (a), if a BBH with total mass
$\sim 500M_{\odot}$ are detected by advanced LIGO (O4), we can find that most
likely model is $\alpha\lambda=0.5$ and about three times more likely than the
second most likely model, $\alpha\lambda=0.1$. As described in Section 3.1.2,
the smaller $\alpha\lambda$ value is, the larger the maximum primary BH mass.
Therefore, as the detected BBH total mass increases, $\alpha\lambda$ value
with the maximum likelihood decreases. We can impose restrictions on
$\alpha\lambda$ to some degree if a BBH with total mass $\gtrsim 500M_{\odot}$
are detected by advanced LIGO (O4) or LISA.

Figure 20–22 show the likelihood ratio of IMF parameters. It is hard to
distinguish between $\gamma_{1}=0.0$ and 1.0 (see Figure 20). If more massive
BBH is detected, smaller $\gamma_{2}$, i.e., shallower IMF, and larger
$m_{\mathrm{crit}}$ are preferred. However, the difference between maximum
likelihood and that of the other models is not so large that it is difficult
to distinguish models.

In Figure 19(a)–22(a), the likelihood in the range $M_{\mathrm{t}}=$600–800 is
not plotted. This is because, as can be seen from Figure 16, BBHs in this mass
range do not merge within the horizon of advanced LIGO. When
$\alpha\lambda=0.1$, the total mass of BBH mergers with short delay time is
$\sim$ 250–1500 $M_{\odot}$. On the other hand, BBHs with long delay time have
two sub-populations with slightly different formation channels, and total mass
is $\lesssim 600M_{\odot}$ and $\gtrsim 800M_{\odot}$. Since the horizon of
advanced LIGO is $z\sim 1$ (see Figure 15), only BBHs with long delay time can
be observed by advanced LIGO. Therefore, when $\alpha\lambda=0.1$, BBHs with
$M_{\mathrm{t}}=$ 600–800 $M_{\odot}$ cannot be detected by advanced LIGO. On
the other hand, since LISA can detect BBHs with short delay time, likelihood
ratios are plotted in the range $M_{\mathrm{t}}=$ 600–800 $M_{\odot}$.

##  5 Summary

In this paper, we perform a binary population synthesis calculation for very
massive ($\sim$ 100–1000 $M_{\odot}$) Pop. III stars and derive the various
property of BBH mergers with IMBHs. We adopt wide range of common envelope
parameter $\alpha\lambda$ from 0.01 to 5 and investigate the dependence of
$\alpha\lambda$ on the primary BH mass, mass ratio, effective spin and delay
time distribution of BBH mergers. We find that if $\alpha\lambda$ is smaller,
the maximum mass of primary BH mass is larger (Section 3).

We derive the redshift evolution of merger rate density and compare our
results with the observation. So far, no BBH with IMBHs have been discovered
by gravitational wave detector, and thus only upper limit on merger rate
density have been obtained. In this study, we adopt various IMF models and
impose restrictions on IMF models by using this upper limit obtained from non-
detection (Section 4.1).

We also compute the detection rate by advanced LIGO (O4), Einstein telescope
and LISA, and obtain the mass distribution of detected BBHs with IMBHs
(Section 4.2.1). We find that Einstein telescope is a promising detector for
the Pop.III BBH mergers with IMBHs, and the detection rate is from $\sim$ 10
$\mathrm{yr}^{-1}$ to 1000 $\mathrm{yr}^{-1}$.

Finally, we calculate the likelihood ratio for model parameter
($\alpha\lambda$, $\gamma_{1}$, $\gamma_{2}$ and $m_{\mathrm{crit}}$) assuming
that one BBH with IMBHs with a certain total mass is observed by certain
detector. We find that we can impose restrictions on $\alpha\lambda$ to some
degree if a BBH with total mass $\gtrsim 500M_{\odot}$ are detected by
advanced LIGO (O4) or LISA.

## Acknowledgements

This work was supported by JSPS KAKENHI grant No. 17H06360 and 19K03907 (A.
T.), 21K13915, 22K03630(T. K), 20H05249 (T. Y.), and 17H01130, 17K05380 and
21H01123 (H. U.), and by the University of Tokyo Young Excellent Researcher
program (T. K.). We use the Python packages, numpy (van der Walt et al., 2011)
and matplotlib (Hunter, 2007) for our analysis.

## Data Availability

Results will be shared on reasonable request to authors.

## References

  * Abbott et al. (2020) Abbott R., et al., 2020, Phys. Rev. Lett., 125, 101102
  * Amaro-Seoane (2018) Amaro-Seoane P., 2018, Phys. Rev. D, 98, 063018
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, arXiv e-prints, p. arXiv:1702.00786
  * Anagnostou et al. (2020) Anagnostou O., Trenti M., Melatos A., 2020, Publ. Astron. Soc. Australia, 37, e044
  * Antonini & Perets (2012) Antonini F., Perets H. B., 2012, ApJ, 757, 27
  * Antonini & Rasio (2016) Antonini F., Rasio F. A., 2016, ApJ, 831, 187
  * Antonini et al. (2017) Antonini F., Toonen S., Hamers A. S., 2017, ApJ, 841, 77
  * Arca Sedda (2020) Arca Sedda M., 2020, ApJ, 891, 47
  * Bae et al. (2014) Bae Y.-B., Kim C., Lee H. M., 2014, MNRAS, 440, 2714
  * Banerjee (2017) Banerjee S., 2017, MNRAS, 467, 524
  * Banerjee (2018a) Banerjee S., 2018a, MNRAS, 473, 909
  * Banerjee (2018b) Banerjee S., 2018b, MNRAS, 481, 5123
  * Banerjee (2021a) Banerjee S., 2021a, MNRAS, 500, 3002
  * Banerjee (2021b) Banerjee S., 2021b, MNRAS, 503, 3371
  * Banerjee et al. (2010) Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 371
  * Bavera et al. (2021) Bavera S. S., et al., 2021, A&A, 647, A153
  * Belczynski (2020) Belczynski K., 2020, ApJ, 905, L15
  * Belczynski et al. (2016) Belczynski K., et al., 2016, A&A, 594, A97
  * Belczynski et al. (2017) Belczynski K., Ryu T., Perna R., Berti E., Tanaka T. L., Bulik T., 2017, MNRAS, 471, 4702
  * Belczynski et al. (2020) Belczynski K., et al., 2020, A&A, 636, A104
  * Bethe & Brown (1998) Bethe H. A., Brown G. E., 1998, ApJ, 506, 780
  * Di Carlo et al. (2019) Di Carlo U. N., Giacobbo N., Mapelli M., Pasquato M., Spera M., Wang L., Haardt F., 2019, MNRAS, 487, 2947
  * Di Carlo et al. (2020) Di Carlo U. N., et al., 2020, MNRAS, 498, 495
  * Downing et al. (2010) Downing J. M. B., Benacquista M. J., Giersz M., Spurzem R., 2010, MNRAS, 407, 1946
  * Downing et al. (2011) Downing J. M. B., Benacquista M. J., Giersz M., Spurzem R., 2011, MNRAS, 416, 133
  * Edwards & Pringle (1987) Edwards D. A., Pringle J. E., 1987, MNRAS, 229, 383
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  * Eldridge et al. (2019) Eldridge J. J., Stanway E. R., Tang P. N., 2019, MNRAS, 482, 870
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration et al., 2019a, ApJ, 875, L1
  * Event Horizon Telescope Collaboration et al. (2019b) Event Horizon Telescope Collaboration et al., 2019b, ApJ, 875, L2
  * Event Horizon Telescope Collaboration et al. (2019c) Event Horizon Telescope Collaboration et al., 2019c, ApJ, 875, L3
  * Event Horizon Telescope Collaboration et al. (2019d) Event Horizon Telescope Collaboration et al., 2019d, ApJ, 875, L4
  * Event Horizon Telescope Collaboration et al. (2019e) Event Horizon Telescope Collaboration et al., 2019e, ApJ, 875, L5
  * Event Horizon Telescope Collaboration et al. (2019f) Event Horizon Telescope Collaboration et al., 2019f, ApJ, 875, L6
  * Event Horizon Telescope Collaboration et al. (2022a) Event Horizon Telescope Collaboration et al., 2022a, ApJ, 930, L12
  * Event Horizon Telescope Collaboration et al. (2022b) Event Horizon Telescope Collaboration et al., 2022b, ApJ, 930, L13
  * Event Horizon Telescope Collaboration et al. (2022c) Event Horizon Telescope Collaboration et al., 2022c, ApJ, 930, L14
  * Event Horizon Telescope Collaboration et al. (2022d) Event Horizon Telescope Collaboration et al., 2022d, ApJ, 930, L15
  * Event Horizon Telescope Collaboration et al. (2022e) Event Horizon Telescope Collaboration et al., 2022e, ApJ, 930, L16
  * Event Horizon Telescope Collaboration et al. (2022f) Event Horizon Telescope Collaboration et al., 2022f, ApJ, 930, L17
  * Ezquiaga & Holz (2021) Ezquiaga J. M., Holz D. E., 2021, ApJ, 909, L23
  * Farmer et al. (2019) Farmer R., Renzo M., de Mink S. E., Marchant P., Justham S., 2019, ApJ, 887, 53
  * Farmer et al. (2020) Farmer R., Renzo M., de Mink S. E., Fishbach M., Justham S., 2020, ApJ, 902, L36
  * Farr et al. (2019) Farr W. M., Fishbach M., Ye J., Holz D. E., 2019, ApJ, 883, L42
  * Farrell et al. (2009) Farrell S. A., Webb N. A., Barret D., Godet O., Rodrigues J. M., 2009, Nature, 460, 73
  * Farrell et al. (2021) Farrell E., Groh J. H., Hirschi R., Murphy L., Kaiser E., Ekström S., Georgy C., Meynet G., 2021, MNRAS, 502, L40
  * Fishbach & Holz (2020) Fishbach M., Holz D. E., 2020, ApJ, 904, L26
  * Fragione & Banerjee (2021) Fragione G., Banerjee S., 2021, ApJ, 913, L29
  * Fragione & Kocsis (2019) Fragione G., Kocsis B., 2019, MNRAS, 486, 4781
  * Fragione et al. (2020) Fragione G., Loeb A., Rasio F. A., 2020, ApJ, 895, L15
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Fujii et al. (2017) Fujii M. S., Tanikawa A., Makino J., 2017, PASJ, 69, 94
  * Greene et al. (2020) Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257
  * Hamers et al. (2021) Hamers A. S., Fragione G., Neunteufel P., Kocsis B., 2021, MNRAS, 506, 5345
  * Hartwig et al. (2016) Hartwig T., Volonteri M., Bromm V., Klessen R. S., Barausse E., Magg M., Stacy A., 2016, MNRAS, 460, L74
  * Hijikawa et al. (2021) Hijikawa K., Tanikawa A., Kinugawa T., Yoshida T., Umeda H., 2021, MNRAS, 505, L69
  * Hirano et al. (2014) Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., Yorke H. W., 2014, ApJ, 781, 60
  * Hirano et al. (2015) Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MNRAS, 448, 568
  * Hong et al. (2018) Hong J., Vesperini E., Askar A., Giersz M., Szkudlarek M., Bulik T., 2018, MNRAS, 480, 5645
  * Hong et al. (2020) Hong J., Askar A., Giersz M., Hypki A., Yoon S.-J., 2020, MNRAS, 498, 4287
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Inayoshi et al. (2016) Inayoshi K., Kashiyama K., Visbal E., Haiman Z., 2016, MNRAS, 461, 2722
  * Kinugawa et al. (2014) Kinugawa T., Inayoshi K., Hotokezaka K., Nakauchi D., Nakamura T., 2014, MNRAS, 442, 2963
  * Kinugawa et al. (2016) Kinugawa T., Miyamoto A., Kanda N., Nakamura T., 2016, MNRAS, 456, 1093
  * Kinugawa et al. (2020) Kinugawa T., Nakamura T., Nakano H., 2020, MNRAS, 498, 3946
  * Kinugawa et al. (2021a) Kinugawa T., Nakamura T., Nakano H., 2021a, MNRAS, 501, L49
  * Kinugawa et al. (2021b) Kinugawa T., Nakamura T., Nakano H., 2021b, MNRAS, 504, L28
  * Kremer et al. (2020) Kremer K., et al., 2020, ApJ, 903, 45
  * Kruckow et al. (2018) Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G., 2018, MNRAS, 481, 1908
  * Kumamoto et al. (2019) Kumamoto J., Fujii M. S., Tanikawa A., 2019, MNRAS, 486, 3942
  * Kumamoto et al. (2020) Kumamoto J., Fujii M. S., Tanikawa A., 2020, MNRAS, 495, 4268
  * Kumamoto et al. (2021) Kumamoto J., Fujii M. S., Trani A. A., Tanikawa A., 2021, arXiv e-prints, p. arXiv:2102.09323
  * Kuns et al. (2020) Kuns K. A., Yu H., Chen Y., Adhikari R. X., 2020, Phys. Rev. D, 102, 043001
  * Leigh et al. (2018) Leigh N. W. C., et al., 2018, MNRAS, 474, 5672
  * Liu & Bromm (2020) Liu B., Bromm V., 2020, MNRAS, 495, 2475
  * Liu & Bromm (2021) Liu B., Bromm V., 2021, MNRAS, 506, 5451
  * Liu et al. (2021) Liu B., Meynet G., Bromm V., 2021, MNRAS, 501, 643
  * Luo et al. (2016) Luo J., et al., 2016, Classical and Quantum Gravity, 33, 035010
  * Mandel & de Mink (2016) Mandel I., de Mink S. E., 2016, MNRAS, 458, 2634
  * Mangiagli et al. (2019) Mangiagli A., Bonetti M., Sesana A., Colpi M., 2019, ApJ, 883, L27
  * Mapelli (2016) Mapelli M., 2016, MNRAS, 459, 3432
  * Mapelli et al. (2017) Mapelli M., Giacobbo N., Ripamonti E., Spera M., 2017, MNRAS, 472, 2422
  * Mapelli et al. (2019) Mapelli M., Giacobbo N., Santoliquido F., Artale M. C., 2019, MNRAS, 487, 2
  * Mapelli et al. (2021) Mapelli M., et al., 2021, MNRAS, 505, 339
  * Marchant et al. (2016) Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, A&A, 588, A50
  * Matsubayashi et al. (2004) Matsubayashi T., Shinkai H.-a., Ebisuzaki T., 2004, ApJ, 614, 864
  * McKernan et al. (2018) McKernan B., et al., 2018, ApJ, 866, 66
  * Nitz & Capano (2021) Nitz A. H., Capano C. D., 2021, ApJ, 907, L9
  * O’Leary et al. (2006) O’Leary R. M., Rasio F. A., Fregeau J. M., Ivanova N., O’Shaughnessy R., 2006, ApJ, 637, 937
  * O’Leary et al. (2009) O’Leary R. M., Kocsis B., Loeb A., 2009, MNRAS, 395, 2127
  * Paczyński & Sienkiewicz (1972) Paczyński B., Sienkiewicz R., 1972, Acta Astron., 22, 73
  * Park et al. (2017) Park D., Kim C., Lee H. M., Bae Y.-B., Belczynski K., 2017, MNRAS, 469, 4665
  * Pavlovskii et al. (2017) Pavlovskii K., Ivanova N., Belczynski K., Van K. X., 2017, MNRAS, 465, 2092
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Punturo et al. (2010) Punturo M., et al., 2010, Classical and Quantum Gravity, 27, 194002
  * Qin et al. (2018) Qin Y., Fragos T., Meynet G., Andrews J., Sørensen M., Song H. F., 2018, A&A, 616, A28
  * Rasskazov & Kocsis (2019) Rasskazov A., Kocsis B., 2019, ApJ, 881, 20
  * Reitze et al. (2019) Reitze D., et al., 2019, in Bulletin of the American Astronomical Society. p. 35 (arXiv:1907.04833)
  * Riley et al. (2021) Riley J., Mandel I., Marchant P., Butler E., Nathaniel K., Neijssel C., Shortt S., Vigna-Gómez A., 2021, MNRAS, 505, 663
  * Rodriguez & Antonini (2018) Rodriguez C. L., Antonini F., 2018, ApJ, 863, 7
  * Rodriguez et al. (2016) Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D, 93, 084029
  * Rodriguez et al. (2018) Rodriguez C. L., Amaro-Seoane P., Chatterjee S., Rasio F. A., 2018, Phys. Rev. Lett., 120, 151101
  * Rodriguez et al. (2019) Rodriguez C. L., Zevin M., Amaro-Seoane P., Chatterjee S., Kremer K., Rasio F. A., Ye C. S., 2019, Phys. Rev. D, 100, 043027
  * Sadowski et al. (2008) Sadowski A., Belczynski K., Bulik T., Ivanova N., Rasio F. A., O’Shaughnessy R., 2008, ApJ, 676, 1162
  * Samsing (2018) Samsing J., 2018, Phys. Rev. D, 97, 103014
  * Samsing et al. (2018) Samsing J., Askar A., Giersz M., 2018, ApJ, 855, 124
  * Sathyaprakash et al. (2012) Sathyaprakash B., et al., 2012, Classical and Quantum Gravity, 29, 124013
  * Savonije (1978) Savonije G. J., 1978, A&A, 62, 317
  * Seto et al. (2001) Seto N., Kawamura S., Nakamura T., 2001, Phys. Rev. Lett., 87, 221103
  * Shinkai et al. (2017) Shinkai H.-a., Kanda N., Ebisuzaki T., 2017, ApJ, 835, 276
  * Silsbee & Tremaine (2017) Silsbee K., Tremaine S., 2017, ApJ, 836, 39
  * Sopuerta & Yunes (2009) Sopuerta C. F., Yunes N., 2009, Phys. Rev. D, 80, 064006
  * Stevenson et al. (2017) Stevenson S., Vigna-Gómez A., Mandel I., Barrett J. W., Neijssel C. J., Perkins D., de Mink S. E., 2017, Nature Communications, 8, 14906
  * Stone et al. (2017) Stone N. C., Metzger B. D., Haiman Z., 2017, MNRAS, 464, 946
  * Susa et al. (2014) Susa H., Hasegawa K., Tominaga N., 2014, ApJ, 792, 32
  * Tagawa et al. (2020a) Tagawa H., Haiman Z., Kocsis B., 2020a, ApJ, 898, 25
  * Tagawa et al. (2020b) Tagawa H., Haiman Z., Bartos I., Kocsis B., 2020b, ApJ, 899, 26
  * Tagawa et al. (2021a) Tagawa H., Kocsis B., Haiman Z., Bartos I., Omukai K., Samsing J., 2021a, ApJ, 907, L20
  * Tagawa et al. (2021b) Tagawa H., Kocsis B., Haiman Z., Bartos I., Omukai K., Samsing J., 2021b, ApJ, 908, 194
  * Takahashi et al. (2016) Takahashi K., Yoshida T., Umeda H., Sumiyoshi K., Yamada S., 2016, MNRAS, 456, 1320
  * Takahashi et al. (2018) Takahashi K., Yoshida T., Umeda H., 2018, ApJ, 857, 111
  * Takahashi et al. (2019) Takahashi K., Sumiyoshi K., Yamada S., Umeda H., Yoshida T., 2019, ApJ, 871, 153
  * Tanikawa (2013) Tanikawa A., 2013, MNRAS, 435, 1358
  * Tanikawa et al. (2020) Tanikawa A., Yoshida T., Kinugawa T., Takahashi K., Umeda H., 2020, MNRAS, 495, 4170
  * Tanikawa et al. (2021a) Tanikawa A., Kinugawa T., Yoshida T., Hijikawa K., Umeda H., 2021a, MNRAS, 505, 2170
  * Tanikawa et al. (2021b) Tanikawa A., Susa H., Yoshida T., Trani A. A., Kinugawa T., 2021b, ApJ, 910, 30
  * Tanikawa et al. (2022) Tanikawa A., Yoshida T., Kinugawa T., Trani A. A., Hosokawa T., Susa H., Omukai K., 2022, ApJ, 926, 83
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration the Virgo Collaboration the KAGRA Collaboration 2021a, arXiv e-prints, p. arXiv:2105.15120
  * The LIGO Scientific Collaboration et al. (2021b) The LIGO Scientific Collaboration et al., 2021b, arXiv e-prints, p. arXiv:2111.03606
  * Trani et al. (2021) Trani A. A., Tanikawa A., Fujii M. S., Leigh N. W. C., Kumamoto J., 2021, MNRAS, 504, 910
  * Trani et al. (2022) Trani A. A., Rastello S., Di Carlo U. N., Santoliquido F., Tanikawa A., Mapelli M., 2022, MNRAS, 511, 1362
  * Umeda et al. (2020) Umeda H., Yoshida T., Nagele C., Takahashi K., 2020, ApJ, 905, L21
  * Vink et al. (2021) Vink J. S., Higgins E. R., Sander A. A. C., Sabhahit G. N., 2021, MNRAS, 504, 146
  * Vynatheya & Hamers (2022) Vynatheya P., Hamers A. S., 2022, ApJ, 926, 195
  * Wang et al. (2019) Wang H.-T., et al., 2019, Phys. Rev. D, 100, 043003
  * Wang et al. (2022) Wang L., Tanikawa A., Fujii M., 2022, MNRAS, 515, 5106
  * Weatherford et al. (2021) Weatherford N. C., Fragione G., Kremer K., Chatterjee S., Ye C. S., Rodriguez C. L., Rasio F. A., 2021, ApJ, 907, L25
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Yang et al. (2019) Yang Y., Bartos I., Haiman Z., Kocsis B., Márka Z., Stone N. C., Márka S., 2019, ApJ, 876, 122
  * Yoshida et al. (2019) Yoshida T., Takiwaki T., Kotake K., Takahashi K., Nakamura K., Umeda H., 2019, ApJ, 881, 16
  * Ziosi et al. (2014) Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, MNRAS, 441, 3703
  * de Souza et al. (2011) de Souza R. S., Yoshida N., Ioka K., 2011, A&A, 533, A32
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science and Engineering, 13, 22