# Merger rate density of binary black holes through isolated Population I, II,
III and extremely metal-poor binary star evolution

Ataru Tanikawa  Department of Earth Science and Astronomy, College of Arts and
Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902,
Japan  Takashi Yoshida  Yukawa Institute for Theoretical Physics, Kyoto
University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan  Tomoya
Kinugawa  Institute for Cosmic Ray Research, The University of Tokyo, Kashiwa,
Chiba  Alessandro A. Trani  Department of Earth Science and Astronomy, College
of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo
153-8902, Japan  Okinawa Institute of Science and Technology, Quantum Gravity
unit, 1919-1 Tancha, Onna-son, Okinawa 904-0495, Japan  Takashi Hosokawa
Department of Physics, Graduate School of Science, Kyoto University, Sakyo,
Kyoto 606-8502, Japan  Hajime Susa  Department of Physics, Faculty of Science
and Engineering, Konan University, Kobe 658-8501, Japan  Kazuyuki Omukai
Astronomical Institute, Graduate School of Science, Tohoku University, Aoba,
Sendai 980-8578, Japan

###### Abstract

We investigate the formation of merging binary black holes (BHs) through
isolated binary evolution, performing binary population synthesis calculations
covering an unprecedentedly wide metallicity range of Population (Pop) I, II,
III, and extremely metal-poor (EMP) binary stars. We find that the predicted
merger rate density and primary BH mass ($m_{1}$) distribution are consistent
with the gravitational wave (GW) observations. Notably, Pop III and EMP
($<10^{-2}$ $Z_{\odot}$) binary stars yield most of the pair instability (PI)
mass gap events with $m_{1}=65$–$130$ $M_{\odot}$. Pop III binary stars
contribute more to the PI mass gap events with increasing redshift, and all
the PI mass gap events have the Pop III origin at redshifts $\gtrsim 8$. Our
result can be assessed by future GW observations in the following two points.
First, there are no binary BHs with $m_{1}=100$–$130$ $M_{\odot}$ in our
result, and thus the $m_{1}$ distribution should suddenly drop in the range of
$m_{1}=100$–$130$ $M_{\odot}$. Second, the PI mass gap event rate should
increase toward higher redshift up to $\sim 11$, since those events mainly
originate from the Pop III binary stars. We find that the following three
assumptions are needed to reproduce the current GW observations: a top-heavy
stellar initial mass function and the presence of close binary stars for Pop
III and EMP binary stars, and inefficient convective overshoot in the main-
sequence phase of stellar evolution. Without any of the above, the number of
PI mass gap events becomes too low to reproduce current GW observations.

stars: binaries: close – stars: black holes – stars: Population II – stars:
Population III – gravitational waves

††software: Open data from the first and second observing runs of Advanced
LIGO and Advanced Virgo (Abbott et al., 2021c); HOSHI (Takahashi et al., 2016,
2018, 2019; Yoshida et al., 2019); BSE (Hurley et al., 2000, 2002; Tanikawa et
al., 2020, 2021a); Matplotlib (Hunter, 2007); NumPy (van der Walt et al.,
2011).

##  1 Introduction

A growing number of binary black hole (BH) mergers have been discovered by
gravitational wave (GW) observatories, reaching nearly $50$ by the end of the
first half of the third observing run (Abbott et al., 2019, 2021a; The LIGO
Scientific Collaboration et al., 2021b). The number of binary BHs has grown to
about $90$ (The LIGO Scientific Collaboration et al., 2021c), while this paper
is under review. The origin of merging binary BHs is still an open question
(see review by Mandel & Farmer, 2018; Gerosa & Fishbach, 2021; Mapelli, 2021).
The progenitors can be isolated Population (Pop) I/II binary stars that go
through common envelope evolution (Bethe & Brown, 1998; Belczynski et al.,
2016a, 2020a; Stevenson et al., 2017; Eldridge et al., 2017, 2019; Mapelli et
al., 2017, 2019; Kruckow et al., 2018; Olejak et al., 2020), through
chemically homogeneous evolution (Mandel & de Mink, 2016; Marchant et al.,
2016; Riley et al., 2021), and through stable mass transfer (van den Heuvel et
al., 2017; Neijssel et al., 2019; Bavera et al., 2021). Isolated Pop III
binary stars are also suggested to form binary BHs mainly through stable mass
transfer (Kinugawa et al., 2014, 2016; Inayoshi et al., 2017), and partly
through common envelope evolution (Kinugawa et al., 2020, 2021a, 2021b;
Tanikawa et al., 2021a, b). Merging binary BHs can be formed through dynamical
capture in globular clusters (Portegies Zwart & McMillan, 2000; O’Leary et
al., 2006; Sadowski et al., 2008; Downing et al., 2010, 2011; Tanikawa, 2013;
Bae et al., 2014; Rodriguez et al., 2016a, b, 2018, 2019; Askar et al., 2017;
Fujii et al., 2017; Park et al., 2017; Samsing et al., 2018; Samsing, 2018;
Hong et al., 2018, 2020; Anagnostou et al., 2020a; Kremer et al., 2020;
Weatherford et al., 2021; Trani et al., 2021; Wang et al., 2021), in young
star clusters (Banerjee et al., 2010; Banerjee, 2017, 2018a, 2018b, 2021a,
2021b; Ziosi et al., 2014; Mapelli, 2016; Rastello et al., 2019, 2021;
Kumamoto et al., 2019, 2020, 2021; Di Carlo et al., 2019, 2020a, 2020b;
Santoliquido et al., 2020, 2021; Fragione et al., 2020a; Fragione & Banerjee,
2021; Arca Sedda et al., 2021), at galactic centers (O’Leary et al., 2009;
Antonini & Rasio, 2016; Stone et al., 2017; Leigh et al., 2018; McKernan et
al., 2018; Rasskazov & Kocsis, 2019; Trani et al., 2019; Yang et al., 2019;
Arca Sedda, 2020; Mapelli et al., 2021a, b; Tagawa et al., 2020a, b, 2021a,
2021b), and/or in Pop III star clusters (Liu & Bromm, 2020a, b, 2021; Liu et
al., 2021). BH mergers can be also triggered by secular evolution in
hierarchical multiple stellar systems (Antonini et al., 2014, 2017; Silsbee &
Tremaine, 2017; Rodriguez & Antonini, 2018; Fragione & Kocsis, 2019; Fragione
et al., 2020b; Hamers & Safarzadeh, 2020; Hamers et al., 2021; Vynatheya &
Hamers, 2021), and around massive BHs at galactic centers (Antonini & Perets,
2012; VanLandingham et al., 2016; Hamers et al., 2018; Hoang et al., 2018).
Merging binary BHs can be primordial origin (Franciolini et al., 2021).
Multiple origins have been suggested (Bouffanais et al., 2021; Zevin et al.,
2021; Ng et al., 2021; Hütsi et al., 2021).

The isolated binary scenario, however, has been challenged by the discovery of
GW190521 with a $\sim 90$ $M_{\odot}$ BH with the $90$ % credible interval
(Abbott et al., 2020a), which is prohibited to form by the pulsational pair
instability (PPI) and pair instability (PI) supernovae (Heger et al., 2003;
Belczynski et al., 2016b; Spera & Mapelli, 2017; Woosley, 2017, 2019; Giacobbo
et al., 2018), in contrast to the hierarchical merger scenario in star
clusters (Rodriguez et al., 2019; Di Carlo et al., 2020a; Liu & Bromm, 2020b;
Anagnostou et al., 2020b; Tagawa et al., 2021b; Kimball et al., 2021; Mapelli
et al., 2021b; Liu & Lai, 2021; Gayathri et al., 2021; Gerosa & Fishbach,
2021; Arca Sedda et al., 2021). There are, however, still some way outs to
accommodate the existence of GW190521 in the framework of the isolated binary
scenario, which can be classified into three categories. The first possibility
is that GW190521 actually straddles the PI mass gap; the primary and secondary
BHs have masses beyond and below the PI mass gap, respectively (Fishbach &
Holz, 2020; Nitz & Capano, 2021; Estellés et al., 2021). The second is that
the PI mass gap may shift upward because of microphysics in single star
evolution different from what has been thought: e.g., ${}^{12}{\rm
C}(\alpha,\gamma)^{16}{\rm O}$ reaction rate three times smaller than the
standard one (Farmer et al., 2019, 2020; Costa et al., 2021), or dark matter
effects (Croon et al., 2020; Ziegler & Freese, 2021). In particular,
Belczynski (2020) have shown that isolated binary stars can form the PI mass
gap events consistent with GW190521-like event rate, taking into account the
reduced ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction rate. The third
is that the lower PI mass gap, $\sim 65$–$100$ $M_{\odot}$, may disappear with
more detailed consideration of macrophysics in single star evolution, such as
metal-poor star evolution (Kinugawa et al., 2021a; Farrell et al., 2021; Vink
et al., 2021; Tanikawa et al., 2021a), strength of stellar wind mass loss
(Belczynski et al., 2020b), and stellar rotation effect (Mapelli et al.,
2020). In this paper, we focus on the third possibility, to be more precise
metal-poor star evolution, to accommodate the formation of GW190521-like
events in the isolated binary scenario.

Tanikawa et al. (2021a) have shown that Pop III stars can produce sufficient
PI mass gap events to be reconciled with the GW190521-like event rate inferred
by GW observations, if convective overshoot between their convective core and
radiative envelope is inefficient during the main-sequence (MS) phase. Pop III
stars, however, produce fewer binary BH mergers in the lower mass range of
$\lesssim 50$ $M_{\odot}$ than the expectation from GW observations. Thus, we
need to take into account not only Pop III stars but also Pop I/II stars to
explain all the observed BH mergers in the framework of the isolated binary
scenario.

For this purpose, we perform binary population synthesis (BPS) calculations
for isolated Pop I/II/III and extremely metal-poor (EMP) binary stars in an
unprecedentedly wide metallicity range in this paper. The reason for remarking
the “unprecedentedly wide metallicity range” is that we cover not only Pop
I/II stars but also EMP and Pop III stars, which are strictly defined at the
end of this section. Previous studies of binary BHs have focused on only Pop
I/II stars as most of BPS codes support only Pop I/II stars like SeBa
(Portegies Zwart & Verbunt, 1996; Toonen et al., 2012), BSE (Hurley et al.,
2002), binary_c (Izzard et al., 2009), MOBSE (Giacobbo et al., 2018), COSMIC
(Breivik et al., 2020), and COMPAS (Team COMPAS et al., 2021). Some studies
have examined Pop III stars (Kinugawa et al., 2014; Hartwig et al., 2016), and
StarTrack (Belczynski et al., 2002, 2017) can cover Pop I/II/III stars.
However, no studies have included EMP stars.

We investigate isolated Pop I/II binary evolution by our own calculations
rather than making use of the results of previous studies on isolated Pop I/II
binary stars as done by Kinugawa et al. (2021c). We do this for consistency.
We construct our own Pop III star evolution model as in Tanikawa et al.
(2021a), decreasing the metallicity of our own Pop I/II star evolution model,
rather than adopting the widely-used SSE Pop I/II star evolution model (Hurley
et al., 2000). It is unclear if our own Pop I/II star evolution models yield
similar binary BHs to the SSE Pop I/II star evolution models, or if our own
model Pop I/II star evolution models can produce binary BHs in the mass range
of $\lesssim 50$ $M_{\odot}$ consistently with the GW observations, and do not
overproduce the PI mass gap events compared to the GW observations. We expect
that the SSE and our own Pop I star evolution models would yield similar
binary BHs, since these models are similar because of similar calibration to
nearby Pop I stars including the Sun (Pols et al., 1998; Hurley et al., 2000;
Yoshida et al., 2019). On the other hand, this may not be the case for the SSE
and our own Pop II evolution models. Therefore, we use the SSE model for Pop I
stars, and our own model for Pop II stars by our own calculations. We could
have borrowed Pop I binary BHs from previous studies. However, we do not so.
This enables us to examine Pop I binary BHs as well as Pop II, III and EMP
binary BHs in more detail, since we obtain their properties by ourselves. For
example, we can investigate the redshift evolution of binary BH mergers formed
from Pop I/II/III and EMP binary stars.

In this paper, we examine the isolated binary scenario by the BPS
calculations, and demonstrate that from the isolated Pop I/II/III and EMP
binary stars the observed BH merger rate density and its primary BH mass
distribution over $\sim 5$–$100$ $M_{\odot}$ are successfully reproduced under
the assumption of inefficient convective overshoot. This means that isolated
Pop III (and EMP) binaries can form GW190521-like events, and that isolated
Pop I/II star evolution models can produce binary BHs in the mass range of
$\lesssim 50$ $M_{\odot}$ consistently with GW observations without
overproducing PI mass gap events regardless of the choice of Pop I/II star
evolution models. We also predict the redshift evolution of merging binary
BHs, which should help assessing the validity of our model by future GW
observations, e.g., by Cosmic Explorer (Reitze et al., 2019) and Einstein
Telescope (Punturo et al., 2010). We also point out the key assumptions
allowing the predicted BH properties to be consistent with the observations,
considering large uncertainties in the initial conditions and evolution of Pop
III stars due to its non detection so far (Frebel & Norris, 2015; Magg et al.,
2018, 2019).

The structure of this paper is as follows. We present our numerical method to
calculate the BH merger rate density at a given redshift in section 2. We
present and discuss the numerical results in section 3 and 4, respectively. We
summarize this paper in section 5. We adopt the cosmological parameters based
on the latest observation by Planck Collaboration et al. (2020). We define
metallicities of Pop I, Pop II, EMP, and Pop III stars as $Z/Z_{\odot}>0.16$,
$10^{-3}<Z/Z_{\odot}\leq 0.16$, $0<Z/Z_{\odot}\leq 10^{-3}$, and
$Z/Z_{\odot}=0$, respectively, where $Z$ is the metallicity, $Z_{\odot}$ is
the solar metallicity set to $Z_{\odot}=0.02$. We treat $Z/Z_{\odot}=10^{-8}$
stars as Pop III stars, $Z/Z_{\odot}=10^{-6}$ and $10^{-4}$ as EMP stars,
$Z/Z_{\odot}=10^{-2},0.025,0.05$ and $0.1$ as Pop II stars, and
$Z/Z_{\odot}=0.25,0.5$ and $1$ as Pop I stars.

##  2 Method

We describe initial conditions for our BPS calculations in section 2.1. We
overview our single and binary star models for our BPS calculations in section
2.2. As we take into account several different initial conditions and
single/binary star models, we describe our parameter sets for them in section
2.3. In section 2.4, we formulate how to obtain the BH merger rate density and
its derivatives, using our BPS calculation results.

###  2.1 Initial conditions

We account for the star formation rate (SFR) density with two components as

$\displaystyle\phi(z)=\phi_{\rm I/II}(z)+\phi_{\rm III}(z),$ (1)

where $z$ is the redshift, and $\phi_{\rm I/II}(z)$ and $\phi_{\rm III}(z)$
are the Pop I/II and Pop III SFR densities, respectively. We adopt the Pop
I/II SFR density of Madau & Fragos (2017), expressed as

$\displaystyle\phi_{\rm
I/II}(z)=0.01\frac{(1+z)^{2.6}}{1+\left[\left(1+z\right)/3.2\right]^{6.2}}\;M_{\odot}\;{\rm
yr}^{-1}\;{\rm Mpc}^{-3}.$ (2)

We assume the metallicity distribution of Pop I/II components at a given $z$
as a logarithmic normal distribution similarly to Belczynski et al. (2020a)
and Santoliquido et al. (2021), given by

$\displaystyle p_{\rm I/II}(z,Z)$ $\displaystyle=\frac{1}{(2\pi\sigma_{\rm
Z}^{2})^{1/2}}$
$\displaystyle\times\exp\left\\{-\frac{\left[\log(Z/Z_{\odot})-\log(\overline{Z}_{\rm
I/II}(z)/Z_{\odot})\right]^{2}}{2\sigma_{\rm Z}^{2}}\right\\},$ (3)

where $\overline{Z}_{\rm I/II}(z)$ is the average metallicity of the Pop I/II
component at a given redshift $z$, and $\sigma_{\rm Z}$ is the metallicity
dispersion set to $\sigma_{\rm Z}=0.35$. The average metallicity evolution in
the Pop I/II component $\overline{Z}_{\rm I/II}(z)$ is considered in the same
way as Madau & Fragos (2017):

$\displaystyle\log\left(\overline{Z}_{\rm
I/II}(z)/Z_{\odot}\right)=0.153-0.074z^{1.34}.$ (4)

We use the following Pop III SFR density, obtained by simplifying the
simulation result of Skinner & Wise (2020) as shown in their figure 4:

$\displaystyle\log(\phi_{\rm III}(t)/M_{\odot}\;{\rm yr}^{-1}\;{\rm
Mpc}^{-3})=$ $\displaystyle\left\\{\begin{array}[]{ll}0.0130(t/{\rm
Myr}-100)-5.00&(100\leq t/{\rm Myr}<200)\\\ -0.00261(t/{\rm
Myr}-200)-3.70&(200\leq t/{\rm Myr}<400)\\\ -0.0108(t/{\rm
Myr}-400)-4.22&(400\leq t/{\rm Myr}<500)\\\
0&(\mbox{otherwise})\end{array}\right.,$ (9)

where $t$ indicates the cosmic time since the Big Bang. The three formulae in
Eq. (9) are not continuous, since the numbers in the three formulae are
written to $3$ significant digits. We regard that $Z=10^{-8}Z_{\odot}$ stars,
which are the most metal-poor in our single star model, are in our simulations
equivalent to Pop III (i.e. zero-metal) stars. This is because their evolution
is almost identical. Thus, we suppose the metallicity distribution of the Pop
III component as

$\displaystyle p_{\rm III}(z,Z)$ $\displaystyle=\delta(Z-10^{-8}Z_{\odot})$
(10) $\displaystyle=\left\\{\begin{array}[]{ll}\infty&(Z/Z_{\odot}=10^{-8})\\\
0&({\rm otherwise})\end{array}\right.,$ (13)

where the $\delta$ function has the usual property as seen above. This means
that $Z=10^{-8}Z_{\odot}$ (i.e. Pop III or zero-metal) stars have the SFR
density described in Eq. (9) throughout all the times (or redshifts). Then, we
can get the SFR density differentiated by the metallicity:

$\displaystyle\frac{d\phi(z)}{dZ}=\phi_{\rm I/II}(z)p_{\rm
I/II}(z,Z)+\phi_{\rm III}(z)p_{\rm III}(z,Z).$ (14)

This equation is used in section 2.4.

We assume that the binary number fraction is $0.5$ for all metallicities,
where the fraction is close to an intrinsic binary fraction
$0.69^{+0.09}_{-0.09}$ found by Sana et al. (2012). The initial mass function
(IMF) of the single stars and primary stars of the binary stars is given by a
combination of the following two mass functions. The first is the Kroupa’s IMF
(Kroupa, 2001) written as

$\displaystyle f(m)dm\propto\left\\{\begin{array}[]{ll}m^{-1.3}&(0.08\leq
m/M_{\odot}<0.5)\\\ m^{-2.3}&(0.5\leq m/M_{\odot}\leq 150)\\\
\end{array}\right.,$ (17)

where $m$ is the stellar mass. The second is a function flat in the
logarithmic scale (hereafter, the logarithmically flat IMF), such that

$\displaystyle f(m)dm\propto m^{-1}$ $\displaystyle(10\leq m/M_{\odot}\leq
150).$ (18)

We add the above two mass functions with a weight such that more stars follow
the logarithmically flat IMF at lower metallicities. We adopt three different
methods to weight the two functions. The first is to adopt just the Kroupa’s
IMF for all the metallicities. We call this IMF model the “top-light IMF
model”. The second is to use the Kroupa’s IMF for $Z/Z_{\odot}\geq 10^{-5}$
and otherwise the logarithmically flat IMF. This is because the stellar IMF is
theoretically predicted to transition from top-light to top-heavy at
$Z/Z_{\odot}\sim 10^{-5}$ (Bromm & Larson, 2004; Omukai et al., 2005;
Schneider et al., 2006; Maio et al., 2010). We call this IMF model the
“stepwise IMF model”. In the third method, the IMF shifts gradually from
Kroupa’s IMF to the logarithmically flat one with decreasing metallicity. All
the stars follow the Kroupa’s IMF for $Z/Z_{\odot}>10^{-2}$. At
$Z/Z_{\odot}=10^{-2}$, $60$ % of them obey the logarithmically flat IMFs in
terms of mass, while this fraction increases to $90$ % at
$Z/Z_{\odot}=10^{-4}$. At metallicity as low as $Z/Z_{\odot}\leq 10^{-6}$, all
the stars follow the logarithmically flat distribution. This behavior is
inspired from the simulation results of Chon et al. (2021). We call this IMF
model the “transitional IMF model”.

For the initial distributions of binary parameters we refer to Sana et al.
(2012), who proposed the following expressions for the mass ratios ($q$),
periods ($P$), and orbital eccentricities ($e$) of binary stars:

$\displaystyle f(q)\propto q^{-0.1}\;(0.1\leq q\leq 1)$ (19) $\displaystyle
f(\log P)\propto(\log P)^{-0.55}\;(0.15<\log(P/{\rm day})<5.5)$ (20)
$\displaystyle f(e)\propto e^{-0.42}\;(0\leq e\leq 1),$ (21)

where $q$ is defined as the ratio of the secondary mass to the primary mass.
We call this binary star initial condition the “Sana’s model”. The period
range in Eq. (20) is the same as one adopted by de Mink & Belczynski (2015).
On the other hand, our binary number fraction ($0.5$) is different from that
of de Mink & Belczynski (2015) who have adopted $1$ for the binary number
fraction as default by careful consideration of the difference between their
and Sana’s period ranges.

In addition to this, we introduce another binary star initial condition, where
we exclude binary stars with pericenter distances of less than $200$
$R_{\odot}$ when $Z\leq 10^{-2}Z_{\odot}$. This condition is inspired by the
following consideration. $Z/Z_{\odot}\lesssim 10^{-2}$ stars tend to expand to
$\sim 100$ $R_{\odot}$ in radius during their protostellar evolution thanks to
the high mass accretion rate (Hosokawa & Omukai, 2009a, b; Omukai et al.,
2010). Even if a metal-poor star initially has a close companion star, they
may merge during their protostar phase, and may not make a close binary
system. We call this binary initial condition the “wide binary model”.
Nonetheless, we still adopt the Sana’s model as the fiducial one. This is
because the radial expansion depends on the geometry of protostellar
accretion, which is uncertain. If the accretion shifts from spherical to disk
accretion relatively early on, the maximum expansion radius would be less than
$100$ $R_{\odot}$ (Hosokawa et al., 2010). In addition, even with
$Z<10^{-2}Z_{\odot}$ close binaries can be formed through $N$-body interaction
after their maximal expansion phase (e.g. Liu et al., 2021).

###  2.2 Single and binary star models

We employ the SSE model for single stars with $5\times 10^{-3}\leq
Z/Z_{\odot}\leq 1.5$ (Hurley et al., 2000), and M and L models111The two
models are named after the Milky way and the Large Magellanic Cloud, whose
stars were used to calibrate the overshoot parameter of the M and L models,
respectively (Yoshida et al., 2019). for single stars with
$Z/Z_{\odot}=10^{-8},10^{-6},10^{-4}$, and $10^{-2}\leq Z/Z_{\odot}\leq 0.1$
(Tanikawa et al., 2020, 2021a). An L-model star experiences more efficient
convective overshoot than an M-model star does, for which the convective
overshoot happens between the convective core and radiative envelope in the MS
phase. The L-model star has a larger He core than the M-model star at the end
of their MS phase if these stars have the same initial mass. The L-model star
expands much more in radius than the M-model star at their post-MS phase. We
summarize the details of the M/L models, and compare them with the SSE models
in appendix A. We describe later in this section how to adopt the SSE, M, and
L models.

We take into account stellar wind mass-loss in the same way as Tanikawa et al.
(2020), based on the wind model of Belczynski et al. (2010) with modification
of metallicity-dependent luminous blue variable winds (Giacobbo et al., 2018).
For simplicity, we do not include rotational enhancement of stellar winds
(unlike as in Tanikawa et al., 2020).

Our supernova model is the rapid model (Fryer et al., 2012) with the PI mass
loss. The PI mass loss modifies a remnant mass given by the rapid model
($m_{\rm rapid}$) as follows:

$\displaystyle m_{\rm rem}=\left\\{\begin{array}[]{ll}m_{\rm rapid}&(m_{\rm
c,He}\leq m_{\rm c,He,PPI})\\\ m_{\rm c,He,PPI}&(m_{\rm c,He,PPI}<m_{\rm
c,He}\leq m_{\rm c,He,PISN})\\\ 0&(m_{\rm c,He,PISN}<m_{\rm c,He}\leq m_{\rm
c,He,DC})\\\ m_{\rm rapid}&(m_{\rm c,He}>m_{\rm c,He,DC})\\\
\end{array}\right.,$ (26)

where $m_{\rm rem}$ is the final remnant mass, $m_{\rm c,He}$ is the He core
mass, and $m_{\rm c,He,PPI}$, $m_{\rm c,He,PISN}$, and $m_{\rm c,He,DC}$ are
the lower He core mass limits of pulsational PI (PPI: Heger & Woosley, 2002;
Woosley et al., 2007; Yoshida et al., 2016; Leung et al., 2019), PI supernova
(PISN: Barkat et al., 1967; Fraley, 1968; Bond et al., 1984; El Eid & Langer,
1986; Fryer et al., 2001; Heger & Woosley, 2002; Umeda & Nomoto, 2002), and
direct collapse (DC), respectively. Note that, although $m_{\rm c,He,PPI}$ is
fixed to a constant value in the above equation, in reality it should depend
on $m_{\rm c,He}$ in some complex way (see figure 1 in Stevenson et al.,
2019). We consider two sets of the mass limits, since the mass limits strongly
depend on the uncertain ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction
rate (Farmer et al., 2020; Costa et al., 2021). The first and second sets are
$(m_{\rm c,He,PPI},m_{\rm c,He,PISN},m_{\rm
c,He,DC})=(45M_{\odot},65M_{\odot},135M_{\odot})$ and
$(90M_{\odot},90M_{\odot},180M_{\odot})$, which are similar to Belczynski et
al. (2016b) and Belczynski (2020), respectively. The first set is based on the
standard ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction rate with
respect to the STARLIB (Sallaska et al., 2013). Thus, we call the PI model
“the standard PI model”. On the other hand, the second set can be obtained
when the ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction rate is lower
than the standard one by $3\sigma$. Thus, we call the PI model “the $3\sigma$
PI model”.

Table 1: Summary of combinations of single star models and PI models.

Model | Pop III+EMP | Pop II | Pop I | PI
---|---|---|---|---
| ($10^{-8},10^{-6},10^{-4}Z_{\odot}$) | ($10^{-2},0.025,0.05,0.1Z_{\odot}$) | ($0.25,0.5,1Z_{\odot}$) |
M and SSE-Pop I with standard PI | M | M | SSE | Standard
L and SSE-Pop I with standard PI | L | L | SSE | Standard
L and SSE-Pop I with $3\sigma$ PI | L | L | SSE | $3\sigma$
M and SSE-Pop I/II with standard PI | M | SSE | SSE | Standard

Note. — Metallicities in parentheses indicate ones that we adopt for our BPS
calculations to realize Pop I, Pop II, EMP and Pop III stars. We regard
$10^{-8}Z_{\odot}$ as $0Z_{\odot}$.

We have four different combinations of single star models and PI models as
seen in Table 1. In all the combinations, we adopt the SSE model for Pop I
binary stars. In the first combination, we adopt the M model for Pop II and
Pop III+EMP stars, and the standard PI model. In the second and third
combinations, we adopt the L model for Pop II and Pop III+EMP stars, and the
standard and $3\sigma$ PI models, respectively. In the fourth combination, we
adopt the SSE and M models for Pop II and Pop III+EMP stars, respectively, and
the standard PI model. We do not consider the $3\sigma$ PI model when we adopt
the M model for Pop III+EMP stars. This is because we need either the
$3\sigma$ PI model or the M model in order to yield the PI mass gap event, but
not both.

We explain how to choose the SSE or M/L models for Pop I, Pop II, EMP, and Pop
III stars as seen in Table 1. We choose the SSE model for Pop I stars. The
practical reason for this choice is that $Z/Z_{\odot}>0.1$ star evolution
models of the M/L models are not implemented in our BPS code, although these
models are constructed in Yoshida et al. (2019). However, we expect that our
results do not depend on this choice. If we used the M/L models for Pop I
stars in our BPS calculations, these models would be similar to the SSE model.
This is because the M/L models for Pop I stars are based on the results of
Yoshida et al. (2019) which are calibrated to nearby Pop I stars including the
Sun, similarly to the SSE model (Pols et al., 1998; Hurley et al., 2000).
Rather, there is a possibility that our results are affected by the choice of
the SSE and M/L models for more metal-poor stars, i.e. Pop II stars. However,
we can see that there is little difference yielded by the choice of the SSE
and M/L models for Pop II stars as seen in appendix B. Despite that the SSE
model supports metallicity down to $5\times 10^{-3}Z_{\odot}$, we do not use
the SSE model with less than $10^{-2}Z_{\odot}$. The practical reason for this
is that we would like to adopt the same metallicities for Pop II stars among
the SSE and M/L models, and that the M/L models do not support $5\times
10^{-3}\leq Z/Z_{\odot}<10^{-2}$ star evolution. However, we can see that the
not using the SSE models for this range has no impact (see details in appendix
C). We adopt the M/L models for Pop III+EMP stars, since the SSE model does
not fully support Pop III+EMP star evolution. Especially, we regard the M/L
models with $Z/Z_{\odot}=10^{-8}$ as Pop III stars.

We model the BH spin angular momentum ($J_{\rm rem}$) as follows. We know spin
angular momenta of He core ($J_{\rm c,0}$) and H envelope ($J_{\rm e,0}$) for
each star at the pre-supernova phase from the stellar evolution model. We set

$\displaystyle J_{\rm rem}=J_{\rm c,0}(m_{\rm c,1}/m_{\rm c,0})+J_{\rm
e,0}(m_{\rm e,1}/m_{\rm e,0}),$ (27)

where $m_{\rm c,0}$ and $m_{\rm e,0}$ are the masses of He core and H envelope
at the pre-supernova phase, and $m_{\rm c,1}$ and $m_{\rm e,1}$ are the masses
of He core and H envelope contributing to the BH mass, i.e., the sum of
$m_{\rm c,1}$ and $m_{\rm e,1}$ is the BH mass. We obtain the above relation,
assuming that the He core and H envelope are in rigid rotation, and they each
have uniform densities. In this prescription, we may underestimate the BH spin
angular momentum, since stellar rotation and density become slower and smaller
outward, respectively. We can calculate the dimensionless BH spin as

$\displaystyle\chi_{\rm rem}=\min\left(\frac{J_{\rm rem}}{Gm_{\rm
rem}^{2}/c},1\right),$ (28)

where $G$ is the gravitational constant, and $c$ is the light speed.

We include natal kicks of neutron stars and BHs because of asymmetric
supernova explosions. The velocity distribution is given by a single
Maxwellian with $265$ ${\rm km~{}s}^{-1}$ (Hobbs et al., 2005) if the remnant
acquires no fallback mass. The natal kick velocity is reduced by $1-f_{\rm
b}$, where $f_{\rm b}$ is the fraction of the fallback mass (Fryer et al.,
2012). The direction of the natal kick is assumed to be isotropic.

Our binary star model is based on the BSE model (Hurley et al., 2002). We
overview the parameters of our choice, and the different points between the
BSE and our models. The BSE model includes wind accretion, tidal evolution,
stable mass transfer, common envelope evolution, magnetic braking, and orbital
decay through GW radiation. The wind accretion in our model is limited by the
Eddington limit expressed by Cameron & Mock (1967), while the wind accretion
in the BSE model is not. We adopt this limitation, since wind accretion may
exceed the Eddington limit when a star has a massive companion. The BSE and
our models have the same prescription for tidal evolution of stars with
convective envelopes. On the other hand, we adopt two different prescriptions
for tidal evolution of stars with radiative envelopes, one of which is the
same as the BSE model (described in detail later). The BSE and our models have
different criteria for whether a star has radiative or convective envelope. In
the BSE model, a massive star has radiative (convective) envelope at its core
(shell) He burning phases, respectively. On the other hand, in our model, a
massive star has radiative or convective envelope when its effective
temperature ($\log(T_{\rm eff}/\mbox{K})$) is above or below $3.65$,
respectively. The same criteria are applied also for the stability of mass
transfer. We need the different criteria, since the M- and L-model stars can
have radiative envelopes in their shell He burning phases, and convective
envelopes in their core He burning phases unlike the BSE model. In a stable
mass transfer phase, the maximum fraction of transferred mass is $0.5$ in our
model, while it is unity in the BSE model. Such non-conservative mass transfer
is also adopted by the latest BPS calculations (e.g. Belczynski et al., 2020a;
Kinugawa et al., 2020). Our model as well as the BSE model apply the $\alpha$
formalism for the common envelope evolution (Webbink, 1984). We set
$\alpha=1$, and adopt the formulae of Claeys et al. (2014) for $\lambda$. We
assume that when a star in the Hertzsprung gap phase fills its Roche lobe and
its mass transfer is unstable, the star merges with its companion star. This
is because a star in its Hertzsprung gap phase does not have steep density
gradient between the He core and H envelope (Ivanova & Taam, 2004). This
treatment is similar to other binary population synthesis models (Dominik et
al., 2012; Giacobbo et al., 2018). Note that such a star experiences common
envelope evolution in the BSE model.

We briefly describe the prescription for tidal evolution of stars with
radiative envelope, which is similar to Kinugawa et al. (2020). In this case,
the tidal evolution is subject to the radiative damping of the dynamical tide
(Zahn, 1975). The tidal coupling parameter $k/T$ is given by

$\displaystyle\frac{k}{T}$
$\displaystyle=\left(\frac{Gm}{R^{3}}\right)^{1/2}(1+q)^{5/6}E\left(\frac{R}{a}\right)^{5/2}$
(29) $\displaystyle=2.9\times 10^{-2}\;{\rm
yr^{-1}}\;\left(\frac{m}{M_{\odot}}\right)^{1/2}\left(\frac{R}{R_{\odot}}\right)$
$\displaystyle\times\left(\frac{a}{{\rm au}}\right)^{-5/2}(1+q)^{5/6}E$ (30)

(eq. (42) in Hurley et al., 2002, see also Zahn 1977; Hut 1981), where $R$ is
the stellar radius and $a$ is the binary semi-major axis. In the BSE model,
the variable $E$ depends on the stellar mass as

$\displaystyle E=1.1\times 10^{-6}\left(\frac{m}{10M_{\odot}}\right)^{2.84}.$
(31)

Recently, Yoon et al. (2010, see also ) proposed an alternative fitting
formula for the $E$ based on the ratio of the convective core radius ($R_{\rm
conv}$) to the stellar radius, such that

$\displaystyle E=\left\\{\begin{array}[]{ll}10^{-0.42}\left(R_{\rm
conv}/R\right)^{7.5}&\mbox{for H-rich stars}\\\ 10^{-0.93}\left(R_{\rm
conv}/R\right)^{6.7}&\mbox{for He-rich stars}\end{array}\right..$ (34)

We set the convective core radius as

$\displaystyle R_{\rm
conv}\sim\left\\{\begin{array}[]{ll}0.9R_{\odot}(Z/Z_{\odot})^{0.05}(m/10M_{\odot})^{0.8}&\mbox{for
H-rich stars}\\\ 0.5R_{\odot}&\mbox{for He-rich stars}\end{array}\right..$
(37)

We obtain the convective core radius for H-rich stars from our L-model star’s
data. We adopt the same convective core radius for He-rich stars as Kinugawa
et al. (2020), although this value is fitted to extremely metal-poor stars.
Hereafter, the prescriptions for tidal evolution expressed as Eqs. (31) and
(34) are called “original-tide” and “new-tide” models, respectively.

###  2.3 Parameter sets

We investigate eight parameter sets with different IMFs, binary star initial
conditions, and single and binary star models. We summarize the parameter sets
in Table 2.

We regard as the fiducial one the parameter set in which its IMF, binary star
initial condition, single star model, and binary star model are the
transitional IMF, the Sana’s model, the M and SSE-Pop I model with the
standard PI model, and new-tide model, respectively. We present the results of
the SSE-std set only in appendices B and C.

Table 2: Summary of parameter sets.

Parameter set | IMF | Binary IC | Single star | Binary star
---|---|---|---|---
Fiducial | Transitional | Sana | M and SSE-Pop I with standard PI | New tide
Top-light | Top-light | Sana | M and SSE-Pop I with standard PI | New tide
Stepwise | Stepwise | Sana | M and SSE-Pop I with standard PI | New tide
Wide-binary | Transitional | Wide binary | M and SSE-Pop I with standard PI | New tide
Original-tide | Transitional | Sana | M and SSE-Pop I with standard PI | Original tide
L-std | Transitional | Sana | L and SSE-Pop I with standard PI | New tide
L-$3\sigma$ | Transitional | Sana | L and SSE-Pop I with $3\sigma$ PI | New tide
SSE-std | Transitional | Sana | M and SSE-Pop I/II with standard PI | New tide

Note. — We describe IMF, initial condition (IC) of binary stars, and single
and binary star model for each model. For the parameter sets other than the
fiducial one, differences from the fiducial are accentuated by italic.

For parameter sets with the transitional IMF (e.g. the fiducial, wide-binary,
original-tide, L-std, L-$3\sigma$, and SSE-std parameter sets), we prepare
$12$ groups of $10^{6}$ binary stars, and perform BPS calculations for these
groups. The $12$ groups include $10$ different metallicities:
$Z/Z_{\odot}=10^{-8}$, $10^{-6}$, $10^{-4}$, $10^{-2}$, $0.25$, $0.05$, $0.1$,
$0.25$, $0.5$, and $1$. The two most metal-poor groups follow the
logarithmically flat IMF, and the six most metal-rich groups follow the
Kroupa’s IMF. The two groups with $Z/Z_{\odot}=10^{-4}$ and $10^{-2}$ have a
mixture of the Kroupa’s and logarithmically flat IMFs. For the top-light and
stepwise parameter sets, we prepare $10$ groups of $10^{6}$ binary stars with
the $10$ different metallicities. In the top-light parameter set, all the
groups follow the Kroupa’s IMF. For the stepwise parameter set, the two most
metal-poor groups follow the logarithmically flat IMF, and the other groups
follow the Kroupa’s IMF.

We generate $10^{6}$ binary stars without any restriction on the upper and
lower mass limits in the case of the logarithmically flat IMF. The total mass
of these binary stars is $8.2\times 10^{7}$ $M_{\odot}$. Since we assume the
binary number fraction to $0.5$ as described above, the $10^{6}$ binary star
formation should be accompanied by $10^{6}$ single star formation. The total
mass of the $10^{6}$ single stars and the $10^{6}$ binary stars is $1.3\times
10^{8}$ $M_{\odot}$. On the other hand, we generate $10^{6}$ binary stars
whose primary and secondary masses are more than $10$ $M_{\odot}$ for the
Kroupa’s IMF. This is because the number fraction of stars leaving behind BHs
is small for the Kroupa’s IMF. The total mass of these binary stars is
$4.2\times 10^{7}$ $M_{\odot}$. These binary stars should be accompanied by a
large number of single stars and binary stars whose primary and secondary
masses are less than $10$ $M_{\odot}$. The total mass of the single and binary
stars (including the $10^{6}$ binary stars) is $5.2\times 10^{8}$ $M_{\odot}$.

###  2.4 BH merger rate density

We use the results of BPS calculations, and obtain the BH merger rate density
(${\cal R}$) and its derivatives by binary BH properties (${\cal P}$), such as
primary mass, mass ratio, and their combination in the source frame. First, we
express the differential BH merger rate density at a given look-back time as

$\displaystyle\frac{d{\cal R}}{d{\cal P}}=\int dZ\int^{t_{\rm h}-t_{\rm
lb}}_{0}dt_{\rm d}\frac{d\phi(z(t_{\rm d}+t_{\rm lb}))}{dZ}\frac{d^{2}{\cal
N}(Z)}{dt_{\rm d}d{\cal P}},$ (38)

where $t_{\rm h}$ is the Hubble time, $t_{\rm lb}$ is the look-back time, and
$t_{\rm d}$ is the delay time of a BH merger from star formation. We can
regard ${\cal N}(Z)$ as the total number of BH mergers per stellar mass formed
at a given metallicity $Z$. We call ${\cal N}(Z)$ “production efficiency of BH
mergers”.222The production efficiency of BH mergers ${\cal N}(Z)$ is similar
to, but slightly different from the merger efficiency defined by Giacobbo &
Mapelli (2018) and the formation efficiency defined by Klencki et al. (2018).
This is because ${\cal N}(Z)$ includes BH mergers whose delay time is
infinity, or more than the Hubble time by definition. In section 3.1, we
define the production efficiency of BH mergers over the Hubble time $\eta(Z)$,
which is identical to the merger or formation efficiency.

We discretize the integral of $Z$ in Eq. (38), and obtain the following
equation:

$\displaystyle\frac{d{\cal R}}{d{\cal P}}\approx\sum_{i=1}^{N_{\rm
Z}}\int^{t_{\rm h}-t_{\rm lb}}_{0}dt_{\rm d}\Delta\phi_{i}(z(t_{\rm d}+t_{\rm
lb}))\frac{d^{2}{\cal N}(Z_{i})}{dt_{\rm d}d{\cal P}}.$ (39)

$N_{\rm Z}$ is the number of different metallicities (i.e. $10$), and
$Z_{i}/Z_{\odot}=10^{-8},10^{-6},\cdots,0.5,1$ for $i=1,2,\cdots,9,10$,
respectively. $\Delta\phi_{i}(z)$ is the SFR density at a given redshift $z$
in a metallicity range around metallicity $Z_{i}$. The lower and upper
boundaries of the metallicity range are the geometric means with its lower and
upper neighbor metallicities: $(Z_{i-1}Z_{i})^{1/2}$ and
$(Z_{i}Z_{i+1})^{1/2}$, respectively. We define $Z_{0}=0$ and $Z_{11}=\infty$.
Thus, we can calculate $\Delta\phi_{i}(z)$ as

$\displaystyle\Delta\phi_{i}(z)$
$\displaystyle=\int_{(Z_{i-1}Z_{i})^{1/2}}^{(Z_{i}Z_{i+1})^{1/2}}dZ\frac{d\phi(z)}{dZ},$
(40)

where $d\phi(z)/dZ$ can be expressed as Eq. (14). Note that we use the
$Z=Z_{\odot}$ result for this metallicity range, although
$(Z_{i}Z_{i+1})^{1/2}$ is infinite for $i=10$.

We can also discretize the integral in Eq. (39) as follows. The derivative in
Eq. (39) can be divided into two parts: binary stars subject to the Kroupa’s
and logarithmically flat IMFs, such that

$\displaystyle\frac{d^{2}{\cal N}(Z_{i})}{dt_{\rm d}d{\cal
P}}=\sum_{j=1}^{N_{\rm IMF}}f_{i,j}\frac{d^{2}{\cal N}_{i,j}}{dt_{\rm d}d{\cal
P}},$ (41)

where $N_{\rm IMF}$ is the number of IMFs (i.e. $N_{\rm IMF}=1$ or $2$),
$f_{i,j}$ is the mass fraction of stars formed in each IMF, and ${\cal
N}_{i,j}$ is production efficiency of BH mergers for each $Z_{i}$ and each
IMF. Note that $f_{i,1}+f_{i,2}=1$. Eq. (41) can be written as

$\displaystyle\frac{d^{2}{\cal N}_{i,j}}{dt_{\rm d}d{\cal P}}=\frac{1}{M_{{\rm
bps},i,j}}\sum_{k=1}^{N_{{\rm bps},i,j}}\delta(t_{\rm d}-t_{\rm
d}^{i,j,k})\delta({\cal P}-{\cal P}^{i,j,k}),$ (42)

where $N_{{\rm bps},i,j}$ is the number of simulated binary stars for each
$Z_{i}$ and each IMF (i.e. $N_{{\rm bps},i,j}=10^{6}$), $M_{{\rm bps},i,j}$ is
the total stellar mass required for forming the $N_{{\rm bps},i,j}$ binary
stars (i.e. $M_{{\rm bps},i,1}=5.2\times 10^{8}$ $M_{\odot}$ and $M_{{\rm
bps},i,2}=1.3\times 10^{8}$ $M_{\odot}$ as described in section 2.3), and
$t_{\rm d}^{i,j,k}$ and ${\cal P}^{i,j,k}$ are the delay time and binary BH
properties of each binary star. Using Eqs. (41) and (42), we can discretize
Eq. (39) as

$\displaystyle\frac{d{\cal R}}{d{\cal P}}\approx\sum_{i=1}^{N_{\rm
Z}}\sum_{j=1}^{N_{\rm IMF}}\sum_{k=1}^{N_{{\rm
bps},i,j}}\frac{f_{i,j}\Delta\phi_{i}(z(t_{\rm d}^{i,j,k}+t_{\rm
lb}))}{M_{{\rm bps},i,j}}\delta({\cal P}-{\cal P}^{i,j,k}).$ (43)

The BH merger rate density can also be obtained as

$\displaystyle{\cal R}\approx\sum_{i=1}^{N_{\rm Z}}\sum_{j=1}^{N_{\rm
IMF}}\sum_{k=1}^{N_{{\rm bps},i,j}}\frac{f_{i,j}\Delta\phi_{i}(z(t_{\rm
d}^{i,j,k}+t_{\rm lb}))}{M_{{\rm bps},i,j}}.$ (44)

##  3 Results

We present the results for the fiducial parameter set in section 3.1. We see
how the results change in the cases with other parameter sets in section 3.2.

###  3.1 Fiducial parameter set

Figure 1 shows the redshift evolution of the BH merger rate density for the
fiducial set. We find ${\cal R}\sim 20$ yr-1 Gpc-3 at $z=0$, comparable to the
observed rate, $19.3^{+15.1}_{-9.0}$ yr-1 Gpc-3 according to the LIGO-Virgo
Gravitational-Wave Transient Catalog 2 (GWTC-2) (Abbott et al., 2021a)333We
refer to the BH merger rate density, and BH and spin mass distribution based
on GWTC-2 (Abbott et al., 2021b), not the latest one on GWTC-3 (The LIGO
Scientific Collaboration et al., 2021a).. Note that our prediction is also
consistent with the GWTC-2.1 result (The LIGO Scientific Collaboration et al.,
2021b). The predicted merger rate density up to $z=1$ also falls within the
$90$ % credible interval inferred by GWTC-2.

Figure 1: Redshift evolution of the BH merger rate density in the case with
the fiducial parameter set. Black, red, green, blue, and cyan curves indicate
the BH merger rate densities for all $Z$, Pop I, Pop II, Pop III+EMP, and Pop
III, respectively. The gray shaded region shows the 90% credible interval
inferred by the LIGO-Virgo Gravitational-Wave Transient Catalog 2 (GWTC-2).
The dashed curve indicates the redshift evolution of the sum of the Pop I/II
and Pop III SFRs ($\phi_{\rm I/II}$ and $\phi_{\rm III}$, respectively) scaled
by an arbitrarily value $C$. The dashed-dotted curve indicates the redshift
evolution of the Pop III SFR scaled by $10C$.

The BH merger rate density monotonically increases with increasing redshift up
to $z\simeq 2.5$, and then decreases monotonically for higher redshift. The BH
merger rate density is dominated Pop I/II stars at redshifts of $z\lesssim 7$.
Pop III and EMP stars contribute only by $\sim 1$ ($10$) % to the BH merger
rate density at the redshift of $z=0$ ($z=5$, respectively). Their
contribution becomes dominant at redshifts higher than $z\sim 7$.

We find that the redshift at which the BH merger rate reaches the maximum
($z\sim 2.5$) is slightly higher than the redshift of the maximum SFR ($z\sim
2$). In naive expectations, the BH merger rate should achieve the maximum at
some epoch after the peak SFR as binary BHs merge after their progenitor
formation. However, this is not the case for the following reason. At $z<3$,
the SFR for $Z\sim 0.25$ $Z_{\odot}$ is lower than that for $Z\sim 0.5$ ($1$)
$Z_{\odot}$ by a factor of about $3$ ($10$). On the other hand, the production
efficiency of BH mergers over the Hubble time ($\eta(Z)=\int_{0}^{t_{\rm
h}}dt_{\rm d}d{\cal N}(Z)/dt_{\rm d}$)444As described in the previous
footnote, the production efficiency of BH mergers over the Hubble time
$\eta(Z)$ is the same as the merger efficiency called by Giacobbo & Mapelli
(2018) and the formation efficiency called by Klencki et al. (2018). for
$Z\sim 0.25$ $Z_{\odot}$ is higher than those for $Z\sim 0.5$ and $1$
$Z_{\odot}$ by a factor of about $30$ and $1000$, respectively, such that
$\eta(Z)=9.1\times 10^{-6}$ $M_{\odot}^{-1}$ for $Z=0.25$ $Z_{\odot}$,
$3.0\times 10^{-7}$ $M_{\odot}^{-1}$ for $Z=0.5$ $Z_{\odot}$, and $7.6\times
10^{-9}$ $M_{\odot}^{-1}$ for $Z=1$ $Z_{\odot}$. Such sharp dependence on
metallicity has also been pointed out by Giacobbo & Mapelli (2018) and Klencki
et al. (2018). Merging binary BHs are thus mainly yielded by stars of $Z\simeq
0.25$ $Z_{\odot}$ at $z<3$. Moreover, about half of them have much shorter
inspiral times than the Hubble time ($\lesssim 100$ Myr). The redshift
evolution of the BH merger rate density follows the SFR of $Z\simeq 0.25$
$Z_{\odot}$ stars, which achieves the maximum at $z\simeq 2.5$.

Figure 2: BH merger rate density differentiated by primary BH mass (top) and
mass ratio (bottom) at the redshift of $z=0$ in the case with the fiducial
set. The color codes are the same as Figure 1. The PI mass gap region is on
the right-hand side of the dashed line in the top panel. The dashed curve in
the bottom panel indicates $d{\cal R}/dq\propto q^{3}$.

We compare our obtaining BH masses with BH masses inferred by GWTC-2 (Abbott
et al., 2021b). Figure 2 shows the distributions of the primary BH masses
(top) and mass ratios (bottom) at $z=0$, which are defined as the BH merger
rate density differentiated by the primary BH mass and mass ratio,
respectively. The primary BH mass distribution is roughly consistent with the
observations in GWTC-2 in three respects. First, our primary BH mass
distribution has the global maximum at $m_{1}\sim 10$ $M_{\odot}$. Second, the
mass distribution suddenly decreases at $m_{1}\sim 45$ $M_{\odot}$ due to the
PPI effect. It should be noted that merging binary BHs with $m_{1}\sim 45$
$M_{\odot}$ are overproduced while those with $m_{1}\sim 35$ and $50$
$M_{\odot}$ are underproduced, relative to GWTC-2, because of the standard PI
modeling in which PPI produces BHs with a fixed mass ($45$ $M_{\odot}$). If we
had adopted a PPI model whose remnant mass has a broad distribution (e.g.
Yoshida et al., 2016; Leung et al., 2019; Umeda et al., 2020), this
discrepancy would have disappeared. Third, we find BHs within the PI mass gap,
$m_{1}\sim 65$–$90$ $M_{\odot}$. As shown in Tanikawa et al. (2021a), Pop III
stars of $65$–$90$ $M_{\odot}$ keep their small radii ($\lesssim 40$
$R_{\odot}$), and do not lose their H envelopes through binary interactions.
They have He cores with $\lesssim 45$ $M_{\odot}$, and thus leave behind
$65$–$90$ $M_{\odot}$ BHs without either PPI or PISN effect. As seen in
comparison between the blue and cyan curves in Figure 2, not only Pop III
binary stars but also EMP binary stars contribute to PI mass gap mergers. The
results in Figure 2 also eliminate our concerns described in section 1. In
other words, Pop I/II stars can produce binary BHs in the mass range of
$\lesssim 50$ $M_{\odot}$, and do not overproduce the PI mass gap events.

We can see that only Pop III and EMP stars can produce merging binary BHs
within the PI mass gap. In the mass range $65$–$90$ $M_{\odot}$, their radii
remain small, and thus they can keep their H envelope until the collapse to
BHs. They thus leave behind the PI mass gap BHs. On the other hand, Pop I/II
stars in the mass range $65$–$90$ $M_{\odot}$ expand greatly, and lose their H
envelope through binary interactions. They leave behind at most $\sim 45$
$M_{\odot}$ BHs. The radius evolution can be seen in the Appendix A. In
advance, we should remark that the formation of merging binary BHs within the
PI mass gap strongly depends on the choice of the single star model for Pop
III and EMP stars. If we choose the L model, we do not have such binary BHs as
seen in section 3.2.

The bottom panel of Figure 2 shows the mass ratio distribution of merging
binary BHs. The distribution roughly follows the power law $\propto q^{3}$
(dashed curve). According to Abbott et al. (2021b), the power-law index of the
GWTC-2 events is less than or equal to $3$. Our mass ratio distribution is
thus marginally consistent with it.

We find that Pop I/II binary stars preferentially result in equal mass ratio
merger events of $q\sim 1$. Because of the large radii in the post-MS phase,
they tend to experience binary interactions, such as common envelope and
stable mass transfer. After the binary interactions, which work to equalize
the binary component masses, they have similar BH masses. On the other hand,
Pop III and EMP stars hardly experience such binary interactions, and keep
their initial mass ratios and thus the initial flat $q$ distribution (see
section 2.1). In fact, their mass ratios follow nearly flat distribution as
seen in the lower panel in Figure 2 (blue).

Figure 3: Cumulative distribution of effective spin at the redshift of $z=0$
for the fiducial set. The color codes are the same as Figure 1. The magenta
solid and dashed curves indicate the cumulative distribution for
$m_{1}/M_{\odot}\geq 65$ at $z=0$ and $1$, respectively.

Figure 3 shows the cumulative distribution of effective spins for the fiducial
set of the parameters. The effective spin of a binary BH system is defined as

$\displaystyle\chi_{\rm
eff}=\frac{m_{1}\chi_{1}\cos\theta_{1}+m_{2}\chi_{2}\cos\theta_{2}}{m_{1}+m_{2}},$
(45)

where $\chi_{1}$ and $\chi_{2}$ are the dimensionless spin of the primary and
secondary BHs, respectively, and $\theta_{1}$ and $\theta_{2}$ are the
inclination angle between the primary and secondary spin vectors and the
orbital angular momentum vector of the binary BH, respectively. All of merging
binary BHs have zero or positive $\chi_{\rm eff}$. This means that they have
either no spin or spins aligned to their orbital angular momentum vectors,
since they have $\theta_{1}\sim 0$ and $\theta_{2}\sim 0$. Such spinning BHs
acquire their spin angular momenta in the progenitor phase through tidal
interactions with the companion stars.

About $10$ % of merging BHs have $\chi_{\rm eff}\sim 0.4$, which are formed
from Pop I binary stars in the following pathway. After the primary star in a
binary system collapses to a BH with little spin, the secondary star becomes a
giant, and fills its Roche lobe. The binary system then experiences the common
envelope evolution, and shrinks its separation down to $\lesssim 10$
$R_{\odot}$. The secondary star loses its H envelope and eventually becomes a
naked He star. Because of the small separation, it is efficiently spun up by
the tidal interaction (see the formulae in Eqs. (30), (34), and (37)), and
collapses to a BH conserving its spin angular momentum (see Eqs. (27) and
(28)). Since the resultant binary BHs have small separation, they merge after
a short delay time of $\sim 100$ Myr from their progenitor star formation.
This has been argued by Kushnir et al. (2016) and Hotokezaka & Piran (2017).

Pop II binary stars leave behind merging BHs with a wide range of positive
$\chi_{\rm eff}\sim 0$–$1$. These binary BHs especially dominate the
population of merging BHs with $\chi_{\rm eff}\gtrsim 0.4$. The reason for the
wide range of $\chi_{\rm eff}$ is because their progenitor stars are spun up
by tidal interactions, and lose their H envelope, i.e., large reservoir of
spin angular momenta, to some extent through stable mass transfer. Note that
their progenitors lose little mass through stellar winds because of their low
metallicity, and tend not to experience the common envelope evolution since
they become blue supergiant stars (see appendix A). The spin magnitudes of
these binary BHs depend on whether their progenitors keep the H envelope.

Most of Pop III and EMP binary stars also leave behind merging BHs with higher
$\chi_{\rm eff}$ than Pop II binary stars. Pop III and EMP stars with smaller
radii, tend to keep more massive H envelope than more metal-enriched stars
(see appendix A). Merging BHs resulting from Pop III and EMP binary stars
however, do not much affect the cumulative $\chi_{\rm eff}$ distribution of
all the merging binary BHs. This is because they account for only $1$ % of the
total BH merger rate density at $z=0$.

There are, however, discrepancies between our cumulative $\chi_{\rm
eff}$-distribution and observations in three respects. First, no binary BH
with negative $\chi_{\rm eff}$ is found in our sample. Second, our binary BHs
are fewer than observed in the range $\chi_{\rm eff}\sim 0$–$0.1$. Third,
there are more binaries BHs with $\chi_{\rm eff}\sim 0.4$ in our sample than
what is observed. We nevertheless expect that we can mitigate these deviations
once we account for observational uncertainties as in Bavera et al. (2020).

Figure 3 also shows the cumulative distribution of binary BHs with
$m_{1}/M_{\odot}\geq 65$, which are within the PI mass gap, at $z=0$ and $1$.
More than $90$ and $70$ % of binary BHs within the PI mass gap have $\chi_{\rm
eff}>0.5$ at the redshifts of $z=0$ and $1$, respectively. The events at
higher redshift have smaller $\chi_{\rm eff}$ since the progenitor binary
stars have smaller separations at higher redshift, and lose their H envelope
more promptly.

We have presented the cumulative distribution of binary BHs with PI mass gap
BHs both for $z=0$ and $1$, since some of the observed PI mass gap BHs are
located closer to $z=1$ rather than $z=0$. Here, we regard the observed binary
BHs with $m_{1}/M_{\odot}\geq 65$ as PI mass gap events even if they are not
conclusive. For example, GW190521 is conclusive, since its $m_{1}$ is in the
PI mass gap with the $90$ % credible intervals (Abbott et al., 2020a, b). On
the other hand, the followings are not: GW190519_153544, GW190602_175927,
GW190706_222641, and GW190929_012149 (Abbott et al., 2021a), and
GW190403_051519 and GW190426_190642 (The LIGO Scientific Collaboration et al.,
2021b). The conclusive event, GW190521, has small $\chi_{\rm eff}$
($0.03^{+0.32}_{-0.39}$), although the $90$ % credible interval is quite
large. On the other hand, the non-conclusive ones have a wide variety of
$\chi_{\rm eff}$. For example, GW190403_051519, has high $\chi_{\rm eff}$
($0.70^{+0.15}_{-0.27}$ with the $90$ % credible intervals). This suggests
that our model may still be consistent with the observation even considering
the high probability of PI mass gap BHs with high spins in our sample.

We do not try to compare our model with the observations in terms of the so-
called spin precession, i.e., the BH spin component perpendicular to the
binary orbital angular momentum vector. This is because the estimates of spin
precessions have larger uncertainties than those of effective spins, even
though meaningful constraints may be imposed on spin precessions of a few
events in GWTC-2 data (Abbott et al., 2021a), such as GW190412 (Abbott et al.,
2020c).

Figure 4: BH merger rate density as a function of primary BH mass at the
redshifts of $0$ (black), $2$ (blue), $4$ (green), and $8$ (red) for the
fiducial set. The gray, cyan, light-green, and pink curves indicate Pop III
contributions at the redshifts of $0$, $2$, $4$, and $8$, respectively. The
shaded gray regions show the 90% credible interval inferred by the LIGO-Virgo
Gravitational-Wave Transient Catalog 2 (GWTC-2), which is relevant for the
local redshift, $z<2$. The PI mass gap region is on the right-hand side of the
dashed line.

Figure 4 shows the redshift evolution of the primary BH mass distribution for
the fiducial set of the parameters. The global maximum at $m_{1}\sim 10$
$M_{\odot}$ increases from $z=0$ to $z=2-4$ by an order of magnitude, and then
again decreases from $z=4$ to $z=8$ by an order of magnitude. This redshift
evolution reflects that of the Pop I BH merger rate. Note that merging binary
BHs around the global maximum results from the Pop I progenitors. We find that
the BH merger rate density at $m_{1}\sim 10$ $M_{\odot}$ grows much more than
at $m_{1}\sim 15$ $M_{\odot}$ from $z=0$ to $z=2$ and $4$, despite that Pop I
binary stars can yield binary BHs with both $m_{1}\sim 10$ and $15$
$M_{\odot}$ at the redshift of $z=0$ (see Figure 2). This makes a spike in the
BH merger rate density at $m_{1}\sim 10$ $M_{\odot}$ at $z=2$ and $4$. This
spike comes from different formation pathways between Pop I binary BHs with
$m_{1}\sim 10$ and $15$ $M_{\odot}$. The former binary BHs are formed through
common envelope evolution, and typically have short delay time of $\sim 100$
Myr from their progenitor star formation. On the other hand, the latter binary
BHs are formed through stable mass transfer, and their typical merger time is
$\sim 10$ Gyr. Thus, the former merger rate density grow much more than the
latter from $z=0$ to $z=2$ and $4$. In advance, we note that this feature can
be seen in all the sets (see Figure 6).

Interestingly, the BH merger rate density in the mass range $m_{1}\sim
20$–$30$ $M_{\odot}$ does not change over the time interval $z=0$–$8$. This is
because they are originated mainly from Pop II binary stars (see Figure 2),
and their merger rate density does not change over those redshifts (see Figure
1).

The merger rate density of binary BHs within the PI mass gap ($\sim 65$–$90$
$M_{\odot}$) monotonically increases up to $z=8$. They are formed from the Pop
III and EMP progenitors (see Figure 2), whose merger rate increases with
redshift up to $z\sim 8$ (see Figure 1). This is because the Pop III SFR
increases toward higher redshifts up to $z\sim 20$.

We can see from Figure 4 that Pop III contributions in the PI mass gap become
larger with increasing redshift. At $z=0$, EMP binary stars contribute to the
PI mass gap equally to Pop III binary stars. However, Pop III binary stars
start dominating the PI mass gap from above. At $z=2$, $4$, and $8$, Pop III
binary stars have great roles on forming binary BHs with $m_{1}\gtrsim 90$,
$80$, and $70$ $M_{\odot}$ in the mass gap, respectively. We can expect that a
mass gap event has its origin in a Pop III binary star if it is detected at a
redshift $z\gtrsim 8$.

Figure 5: Cumulative distribution of effective spins for the fiducial set. The
color codes are the same as Figure 4.

Figure 5 shows the redshift evolution of the cumulative effective spin
distribution for the fiducial set. The contribution of binary BHs with
$\chi_{\rm eff}\sim 0.4$ increases from $z=0$ to $z=4$. As described above,
they are originated from Pop I binary stars, and have short delay time since
their progenitor formation. Thus, their merger events occur more frequently at
redshifts when star formation (especially of Pop I stars) is more active (i.e.
$z=2$ and $4$). At $z=8$, the fraction of binary BHs with $\chi_{\rm
eff}\simeq 0.4$ is low, while the fraction with higher $\chi_{\rm eff}$ is
high. This is because at this redshift binary BHs are originated dominantly
from the Pop II, Pop III and EMP progenitors, especially Pop III and EMP
progenitors.

###  3.2 Other parameter sets

Figure 6 shows the distributions of primary BH masses at the redshifts of
$z=0$, $2$, $4$, and $8$ for parameter sets other than the fiducial one. We
can see that in the case of the top-light set (top-left) the BH merger rate
for $m_{1}/M_{\odot}\gtrsim 50$ at $z=0$ is definitely smaller than observed,
although the rate for $m_{1}/M_{\odot}\lesssim 50$ at $z=0$ is comparable to
the observation. For the stepwise set (top-right), the BH merger rate of
$m_{1}/M_{\odot}\gtrsim 50$ at $z=0$ is consistent with the observation,
although the PI mass gap rate is smaller than that in the fiducial set. We
also find that there are few binary BHs within the PI mass gap in the case of
the wide-binary set (middle-left). These results demonstrate that the initial
conditions of metal-poor binary stars have great impacts on the formation of
binary BHs within the PI mass gap. Metal-poor binary stars should follow a
top-heavy IMF and have small initial pericenter distance in order to produce a
sufficient number of the PI mass gap events.

Figure 6: BH merger rate density differentiated by primary BH mass at the
redshifts of $0$ (black), $2$ (blue), $4$ (green), and $8$ (red) for parameter
sets other than the fiducial one. The gray, cyan, light-green, and pink curves
indicate Pop III contributions at the redshifts of $0$, $2$, $4$, and $8$,
respectively. The shaded gray regions show the 90% credible interval inferred
by the LIGO-Virgo Gravitational-Wave Transient Catalog 2 (GWTC-2).

The stepwise set (top-right) shows an interesting result as for Pop III binary
stars. Pop III binary stars dominate the PI mass gap events at any redshifts
from $z=0$ to $8$. In other words, EMP binary stars form little PI mass gap
event. This is because all the $Z/Z_{\odot}=10^{-4}$ stars have a top-light
IMF, and the star formation rate of $Z/Z_{\odot}=10^{-6}$ stars is small. If
IMFs depend on metallicity stepwisely like the stepwise set, the PI mass gap
events immediately indicate the presence of Pop III stars.

We next see the results of parameter sets with the L-model stars, where the
L-model stars experience efficient convective overshoot. The result for the
L-std set (bottom-left) shows that the primary BH mass distribution for
$m_{1}/M_{\odot}\lesssim 50$ is consistent with the observed one, while that
for $m_{1}/M_{\odot}\gtrsim 50$ is not. Even Pop III and EMP binary stars fail
to form binary BHs within the PI mass gap. This result is consistent with
Tanikawa et al. (2021b, a). In the L model, Pop III stars with $\sim 80$
$M_{\odot}$ expand in radius up to a few $10^{3}$ $R_{\odot}$. Their H
envelope is stripped through stable mass transfer or common envelope
evolution. Thus, they cannot form PI mass gap BHs through binary evolution.

In contrast to the L-std set, the L-$3\sigma$ set (bottom-right) can produce
binary BHs within the PI mass gap because of the different criteria for PPI
and PISN effects. The primary BH mass distribution at $z=0$ is consistent with
the observation, as is already pointed out by Belczynski (2020). The primary
BH mass distribution in the L-$3\sigma$ set has no sudden drop at
$m_{1}/M_{\odot}\sim 50$ unlike in the fiducial set. Recall that the sudden
drop is caused by PPI effects for the fiducial set. From this point of view,
the fiducial set might be slightly preferred over the L-$3\sigma$ set.

Figure 7: Redshift evolution of the BH merger rate density with
$m_{1}/M_{\odot}\geq 65$ in the case with the fiducial parameter set (top) and
the L-$3\sigma$ set (bottom). The color codes are the same as in Figure 1. Red
curves in both sets, and green curve in the fiducial set can not be seen,
since Pop I stars in both sets, and Pop II stars in the fiducial set do not
form binary BHs with $m_{1}/M_{\odot}\geq 65$.

The striking difference between the fiducial and L-$3\sigma$ sets is seen in
the redshift evolution of the BH merger rate density within the PI mass gap
(see Figure 4 and the bottom-left panel of Figure 6). The merger rate density
increases monotonically from $z=0$ to $z=8$ for the fiducial set, while it
stops increasing at $z=4$ for the L-$3\sigma$ set. This can be confirmed in
Figure 7. The merger rate density has a peak at a redshift of $\sim 11$ in the
fiducial set, and at a redshift of $\sim 6$ in the L-$3\sigma$ set. This is
because Pop III stars dominantly form the PI mass gap events in the fiducial
set, and do not in the L-$3\sigma$ set. Note that Pop III SFR has a peak at a
redshift of $\sim 20$. There are two reasons why the Pop III stars have
smaller contribution to PI mass gap events in the L-$3\sigma$ set than in the
fiducial set. First, Pop II stars can also form PI mass gap events. Second,
the number of Pop III stars forming PI mass gap BHs becomes smaller in the
L-$3\sigma$ set than in the fiducial set, since more massive Pop III stars
form PI mass gap BHs in the L-$3\sigma$ set. In the fiducial set, Pop III
stars with $65$–$90$ $M_{\odot}$ form the PI mass gap BHs, since they do not
lose their H envelopes through binary interactions, and collapse to BHs
without much mass loss. On the other hand, in the L-$3\sigma$ set, Pop III
stars with $\gtrsim 130$ $M_{\odot}$ form the PI mass gap BHs. They expand to
a few $10^{3}$ $R_{\odot}$, lose their H envelopes through binary
interactions, such as mass transfer or common envelope evolution, and become
totally naked He stars with $\gtrsim 65$ $M_{\odot}$. Finally, they leave
behind the PI mass gap BHs. Note that even stars with $65$–$90$ $M_{\odot}$ do
not cause PPI nor PISN in the L-$3\sigma$ set.

Surprisingly, we can find that Pop III binary stars have a significant impact
on the formation of the PI mass gap events at high redshift even in the
L-$3\sigma$ set. As seen in the bottom panel of Figure 6, 1 of 2–3 PI mass gap
events originates from Pop III binary stars at $z\gtrsim 3$. Even if stars
experience efficient convective overshoot, Pop III binary stars can produce
sufficient the PI mass gap events, depending on the PI model.

The primary BH mass distribution in the original-tide set is similar to that
in the fiducial one. Pop III contributions become more important with
increasing redshift, similarly to those in in the fiducial set. The
prescription for tidal interaction has no impact on the primary BH mass
distribution.

Figure 8: Cumulative distribution of effective spins for parameter sets other
than the fiducial one.. The color codes are the same as Figure 6, although Pop
III contributions are not shown.

Figure 8 shows the cumulative distribution of effective spins for parameter
sets other than the fiducial one. Except for the original-tide set, the
redshift evolution of the cumulative $\chi_{\rm eff}$ distributions is the
same as for the fiducial set. This is because merging binary BHs are dominated
by those with $m_{1}\lesssim 50$ $M_{\odot}$ for all the sets. These binary BH
progenitors acquire their spin angular momenta in the same process for all the
sets: Pop I binary stars acquire spins at the phase of naked He stars while
all but Pop I binary stars are tidally spun up keeping their H envelope.

The cumulative distribution in the original-tide set is different from that in
the fiducial set. In the original tidal prescription, naked He stars are less
spun up. Thus, binary BHs with $\chi_{\rm eff}\sim 0.4$ in the fiducial set
are shifted to $\chi_{\rm eff}\sim 0.2$–$0.4$ in the original-tide set. In the
original tidal prescription, blue supergiant stars are more strongly spun up.
Thus, binary BHs with $\chi_{\rm eff}\gtrsim 0.5$ in the original-tide set
contribute more than in the fiducial one. About $50$ % of binary BHs have
$\chi_{\rm eff}\sim 0$ in the fiducial set, while few binary BHs have
$\chi_{\rm eff}\sim 0$ in the original-tide set. This is because blue
supergiant stars are spun up more as described above.

In the original-tide set, the contribution of binary BHs with $\chi_{\rm
eff}\gtrsim 0.5$ is too large to be reconciled with the observation. The new
tidal prescription seems to have advantages in reproducing the GW
observations, which show the dominance of binary BHs with small $\chi_{\rm
eff}$ in GWTC-2 ($\sim 70$–$90$ %) (Abbott et al., 2021a, b; Galaudage et al.,
2021; Roulet et al., 2021).

##  4 Discussion

As described in section 3.1, the BH merger rate density and its derivatives
for the fiducial set are roughly consistent with the results of GWTC-2 and
GWTC-2.1. In this section, we compare the results of the fiducial set with
those of previous studies.

First, we compare our results with two other isolated binary scenarios,
Belczynski (2020) and Kinugawa et al. (2021c), who have claimed that isolated
binary stars can reproduce the BH mass distribution in the range of $5$–$100$
$M_{\odot}$ inferred by the GW observations. Belczynski (2020) have reported
that Pop I/II binary stars with metallicities $0.05$–$1.5$ $Z_{\odot}$ can
make merging binary BHs consistent with the GW observations with regard to the
BH mass and mass ratio distributions. Since we construct the $3\sigma$ PI
model with reference to the model of Belczynski (2020), the results of the
L-$3\sigma$ set are naturally similar to those of Belczynski (2020). Here, we
thus use the results of the L-$3\sigma$ set as a substitute for those of
Belczynski (2020) in discussing the differences between Belczynski (2020)’s
and our model predictions. As described in section 3.2, the BH merger rate
density with PI mass gap BHs evolve differently with redshift between the
models: it stops increasing at $z\sim 6$ in the L-$3\sigma$ set, while it is
still increasing up to $z\sim 11$ in the fiducial set. This is because Pop III
stars dominate the PI mass gap BHs in the fiducial set, and partially form
them in the L-$3\sigma$ set. With future GW observatories, such as Cosmic
Explorer (Reitze et al., 2019) and Einstein Telescope (Punturo et al., 2010),
which might be able to detect merging binary BHs in this mass range up to
$z\sim 10$, we can pick out the right scenario.

Kinugawa et al. (2021c) have argued that merging binary BHs with $m_{1}\gtrsim
30M_{\odot}$ are dominated by Pop III binary stars. On the other hand, not
only Pop III but also Pop II and EMP binaries dominantly yield merging binary
BHs with $m_{1}\gtrsim 50M_{\odot}$ in our model. The future GW observatories
can also distinguish these two models because of their capability to detect
binary BH mergers at high redshifts. In the result of Kinugawa et al. (2021c),
Pop III binary stars typically make binary BHs with $m_{1}\sim 30M_{\odot}$.
Thus, the BH merger rate density with $m_{1}\sim 30M_{\odot}$ should increase
up to at least $z\sim 8$. On the other hand, the BH merger rate density does
not change from $z=0$ to $8$ in our model (see Figure 4).

In section 3.2, we show that Pop III and EMP binary BHs strongly depend on
their IMFs, period distributions, and star evolution models. Here, we mention
that Pop I/II binary BHs are also sensitive to several uncertainties. For
example, the following two factors would make the local BH merger rate density
much lower than the observed one. First, the actual common envelope parameters
largely deviate from our choices, $\alpha=1$ and $\lambda\sim 0.1$ (Dominik et
al., 2012). Second, the metallicity distribution ($\sigma_{\rm Z}$) is much
smaller than our choice, $\sigma_{\rm Z}=0.35$ (Bouffanais et al., 2021). In
those cases, we need to consider dynamical capture scenarios described below.
In addition, the PI model (see Eq. (26)) can affect the BH mass distribution.
Even if the ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction rate is fixed
to the standard one, there are several possibilities of the PPI effects
(Stevenson et al., 2019). The PPI effects can change the BH mass distribution
around $m_{1}\sim 40$ $M_{\odot}$.

The observed primary BH mass distribution can also be explained by dynamical
capture scenarios in globular clusters (e.g. Rodriguez et al., 2019), open
clusters (e.g. Kumamoto et al., 2020; Di Carlo et al., 2020a; Santoliquido et
al., 2020), or active galactic nuclei (AGN) disks (e.g Tagawa et al., 2020b).
The difference between those models and ours is as follows. While binary BHs
with $m_{1}\sim 100-130M_{\odot}$ are in principle allowed to form in the
dynamical capture scenarios, they are not in our scenario. Note that binary
BHs with $m_{1}\gtrsim 130M_{\odot}$ can form if $\gtrsim 300$ $M_{\odot}$
stars are present (Mangiagli et al., 2019; Tanikawa et al., 2021b; Hijikawa et
al., 2021), since they overcome PPI and PISNe. Thus, the primary BH mass
distribution should smoothly extend beyond $m_{1}\sim 100M_{\odot}$ if the
dynamical capture is the dominant formation process of merging binary BHs. If
our scenario is correct, the BH merger rate density should suddenly decrease
at $m_{1}\sim 100-130M_{\odot}$. The fourth and fifth observing runs of the
current GW observatories possibly assess the validity of these two scenarios
if these observing runs detect binary BHs with $m_{1}\sim 100-130M_{\odot}$,
or put some constraints on the upper limit of a BH merger rate in this $m_{1}$
range.

##  5 Summary

We have derived the BH merger rate density and the distributions of merger
properties by means of BPS calculations for $Z/Z_{\odot}=0$–$1$ binary stars
equivalent to Pop I/II/III and EMP binary stars, which span an unprecedentedly
wide metallicity range. We adopted various initial conditions of binary stars
and single and binary star models. With the fiducial set of the parameters, we
successfully obtained the BH merger rate density, the primary BH mass in the
range $5$–$100$ $M_{\odot}$, and effective spin distributions consistent with
the current GW observations. In our results, Pop III and EMP binary stars are
responsible for the PI mass gap events like GW190521. Moreover, Pop III binary
stars become more responsible for the PI mass gap events with increasing
redshift, and contribute most of them at $z\gtrsim 8$. We have found that Pop
III and EMP binary stars should have a top-heavy IMF, small pericenter
distance at the initial time, and inefficient convective overshoot in order to
produce a sufficient number of the PI mass gap events to explain the PI mass
gap event rate inferred by GW observations.

We examine the dependence of binary BH properties on Pop III (and EMP) IMFs,
period distributions, and star evolution models with different convective
overshoot. We note that Pop III stars also contain other uncertain parameters,
such as Pop III SFR and binary fraction, since Pop III stars have not yet been
discovered. These uncertainties would directly affect our predicted BH merger
rate density, especially the PI mass gap event rate mainly yielded by Pop III
(and EMP) binary stars in our model. Apart from Pop III and EMP stars, Pop
I/II binary BHs can also depend on initial conditions (SFR and metallicity
distribution) and binary evolution parameters (common envelope parameters and
mass transfer parameters). If they largely deviate from our choice, the local
BH merger rate density can be much smaller (or larger) than the observed one.
However, examining the dependence is beyond of this paper.

Finally, we emphasize the Pop III contributions to the formation of the PI
mass gap events. In the parameter sets which reproduce the observed primary BH
mass distribution (the fiducial, stepwise, original-tide, and L-$3\sigma$
sets), Pop III binary stars have important roles in forming the PI mass gap
events. In the fiducial, stepwise, and original-tide sets, Pop III binary
stars dominate the PI mass gap at $z\gtrsim 8$. Even in the L-$3\sigma$ set, 1
of 2–3 PI mass gap events originates from Pop III binary stars at $z\gtrsim
3$. This means that future GW observations of the PI mass gap events will be
very useful for Pop III studies, if one of the above isolated binary scenarios
(Belczynski, 2020; Kinugawa et al., 2021c, and this paper) is correct.

The authors appreciate the anonymous referee for the thorough reading and many
fruitful suggestions. The authors wish to express their cordial gratitude to
Prof. Takahiro Tanaka, general PI of Innovative Area Grants-in-Aid for
Scientific Research ”Gravitational wave physics and astronomy: Genesis” for
his continuous interest and encouragement. This research could not been
accomplished without the support by Grants-in-Aid for Scientific Research
(17H01101, 17H01102, 17H01130, 17H02869, 17H06360, 17K05380, 19K03907,
19H01934, 20H00158, 20H05249, 21K13914) from the Japan Society for the
Promotion of Science. This research has made use of data, software and/or web
tools obtained from the Gravitational Wave Open Science Center
(https://www.gw-openscience.org/), a service of LIGO Laboratory, the LIGO
Scientific Collaboration and the Virgo Collaboration. LIGO Laboratory and
Advanced LIGO are funded by the United States National Science Foundation
(NSF) as well as the Science and Technology Facilities Council (STFC) of the
United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. Additional support for
Advanced LIGO was provided by the Australian Research Council. Virgo is
funded, through the European Gravitational Observatory (EGO), by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions
by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan,
Monaco, Poland, Portugal, Spain.

##  Appendix A Single star models

In this section, we show the evolutions of the M- and L-model stars with
various metallicities, and compare them with those of the SSE-model stars as a
reference. We remark again that the M- and L-model stars have inefficient and
efficient convective overshoot, respectively. Before showing their evolutions,
we overview how to construct evolution tracks of the M- and L-model stars,
which is described in detail in Tanikawa et al. (2020). We perform 1D
hydrodynamics simulations of stars with various metallicities and masses.
Based on the simulation results, we make fitting formulae of stellar evolution
tracks, and incorporate the fitting formulae into the BSE code. The main
components of the 1D hydrodynamics simulations are as follows. We use a 1D
stellar evolution code, HOSHI code (Takahashi et al., 2016, 2018, 2019;
Yoshida et al., 2019). We take into account convection, semiconvection, and
convective overshoot for chemical mixing, and do not account for rotation and
rotational mixing. We input stellar physics with nuclear reaction network with
49 species of nuclei (Takahashi et al., 2018), stellar equation of state
(Blinnikov et al., 1996; Vardya, 1960; Iben, 1963), OPAL, molecular, and
conductive opacities (Iglesias & Rogers, 1996; Ferguson et al., 2005; Cassisi
et al., 2007, respectively), and neutrino energy loss (Itoh et al., 1996).

The M- and L-model stars are modeled as stars with the less and more efficient
convective overshoot, respectively. In the HOSHI code, we treat convective
overshoot as a diffusive process above convective regions. The diffusion
coefficient of the convective overshoot exponentially decreases with the
distance from the convective boundary as

$\displaystyle D_{\rm cv}^{\rm ov}=D_{\rm cv,0}\exp\left(-2\frac{\Delta
r}{f_{\rm ov}H_{\rm P0}}\right),$ (A1)

where $D_{\rm cv,0}$ and $H_{\rm P0}$ are the diffusion coefficient and the
pressure scale height at the convective boundary, respectively, and $\Delta r$
is the distance from the convective boundary. The overshoot parameter $f_{\rm
ov}$ is set to be $0.01$ for the M-model stars, and $0.03$ for the L-model
stars, which are the same values of Sets MA and LA in Yoshida et al. (2019),
respectively. The overshoot parameter of the M-model stars is adjust to the MS
width of AB stars observed in open clusters in the Milky Way Galaxy (Ekström
et al., 2012). On the other hand, the overshoot parameter of the L-model stars
are determined based on the calculation to early-B type stars in the Large
Magellanic Cloud similarly to Stern model (Brott et al., 2011).

Figure 9: Hertzsprung-Russell diagram of the M model (our inefficient
convective overshoot model). Each metallicity is indicated in each panel. The
stellar masses are $20$, $40$, $80$, and $160$ $M_{\odot}$ from bottom to top
in each panel. Colors show stellar phases: main sequence (blue), Hertzsprung
gap (cyan), core helium burning (magenta), and shell helium burning (red)
phases. The gray dashed lines show stellar radii of $1$, $10$, $10^{2}$,
$10^{3}$, and $10^{4}$ $R_{\odot}$ from left to right. Stellar wind mass loss
is not considered for their evolution. Figure 10: The same as Figure 9, except
for the L model (our efficient convective overshoot model).

Figures 9 and 10 and show the Hertzsprung-Russell (HR) diagrams of the M-model
and L-model stars, respectively. We do not account for stellar wind mass loss
when making these stellar evolution tracks.

The most striking difference between the M- and L-model stars is the radius
evolution of Pop III and EMP stars with $m/M_{\odot}=80$–$160$. The L-model
stars expand to more than $10^{3}$ $R_{\odot}$ and evolve to red supergiant
stars, while the M-model stars expand up to $\sim 10^{2}$ $R_{\odot}$ and keep
blue supergiant stars. Since the L-model stars involve more efficient
convective overshoot, they transfer large amounts of H in their radiative
envelopes into their convective cores, and form large He cores at the end of
their MS phases. For example, the L-model stars with $m/M_{\odot}=80$ and
$Z/Z_{\odot}=10^{-8}$ (or Pop III) leave a $42$ $M_{\odot}$ He core. Such
large He cores create large luminosities in their post-MS phases, and expand
their H envelopes. On the other hand, the M-model stars form small He cores in
their MS phases. As a reference, the M-model stars with $m/M_{\odot}=80$ and
$Z/Z_{\odot}=10^{-8}$ (or Pop III) leave a $36$ $M_{\odot}$ He core. Since
such small He cores emit small luminosities, and their H envelopes have small
opacities because of their low metallicities, they do not expand much. Vink et
al. (2021) have also found similar results.

Both of the M- and L-model stars evolve to red supergiant stars if they are
Pop I/II stars with $m/M_{\odot}=40$–$160$. This is because their He cores are
large and the opacities of their H envelopes are large due to their high
metallicities. Both of the M- and L-model stars tend to keep blue supergiant
stars if they are Pop III and EMP stars with $m/M_{\odot}=20$–$40$. Their
luminosities are not large enough to expand their radii to red supergiant
stars.

Figure 11: The same as Figure 9, except for the SSE model. Note that the SSE
model does not support stars with $Z/Z_{\odot}=10^{-4},10^{-6}$, and
$10^{-8}$.

Figure 11 shows the HR diagrams of the SSE-model stars as a reference. Note
that the SSE model does not prepare stars with $Z/Z_{\odot}<5\times 10^{-3}$.
They expand to more than $10^{3}$ $R_{\odot}$ for all masses and $Z$. These
features are similar to the L-model stars.

##  Appendix B Comparison between the SSE- and M-model stars

In this section, we compare between the BH merger rate densities in the
fiducial and SSE-std sets as a reference. Additionally, we also make
comparison between the L-std and SSE-std sets. The difference between these
sets is single star models with Pop II metallicities. Thus, we focus only on
their results.

Figure 12: BH merger rate density differentiated by primary BH mass and
cumulative distribution of effective spins at the redshift of $z=0$ for the
fiducial and SSE-std sets (the top panels) and for the L-std and SSE-std sets
(the bottom panels). The color codes of the fiducial and L-std sets are the
same as Figures 2 and 3. The color codes of the SSE-std set are lighter than
the counterparts of the fiducial set.

The top panels of Figure 12 show the mass and spin distributions of merging
binary BHs at the redshift of $z=0$ in the fiducial and SSE-std sets. Roughly
speaking, the mass and spin distributions in the SSE-std set are in good
agreement with those in the fiducial set.

Especially, their spin distributions are quite similar. So, we focus only on
their mass distributions. The primary BH masses in the SSE-std set are cut off
at $m_{1}/M_{\odot}=45$ because of PPI and PISN effects, similarly to those in
the fiducial set. The SSE-model stars expand to more than $10^{3}$ $R_{\odot}$
and evolve to red supergiant stars. They lose their H envelopes through stable
mass transfer or common envelope evolution. Thus, they cannot form PI mass gap
BHs.

The BH merger rate density in the SSE-std set is a few times smaller than in
the fiducial set in the range of $m_{1}/M_{\odot}\sim 30$–$40$. The reason is
as follows. The SSE-model stars have larger radii in their MS phases than the
M-model stars. Thus, the SSE-model stars lose more masses in their MS phases
through stable mass transfer. Then, they leave smaller He cores. In their
post-MS phases, they further lose their masses through stable mass transfer,
and finally leave naked He stars. Note that they have companion stars close
enough to merge within the Hubble time. Thus, the resulting BHs in the SSE-std
set have smaller mass than in the fiducial set, and the BH merger rate density
with $m_{1}/M_{\odot}\sim 30$–$40$ is smaller.

In summary, the M-model and SSE-std stars generate consistent binary BHs,
although there are several small differences. Moreover, when it comes to Pop
II results, the M/L-model and SSE-std stars yield similar binary BHs.

##  Appendix C Effects of the metallicity choice in SSE-model stars

For the SSE-std set, we choose SSE-model stars with metallicities $0.25,0.5$,
and $1Z_{\odot}$ for Pop I stars, SSE-model stars with metallicities
$10^{-2},0.025,0.05$, and $0.1Z_{\odot}$ for Pop II stars, and M-model stars
with metallicities $10^{-8},10^{-6}$, and $10^{-4}Z_{\odot}$ for Pop III and
EMP stars as seen in Tables 1 and 2. We do not use SSE-model stars with
$5\times 10^{-3}Z_{\odot}$, despite that the SSE model supports metallicities
down to $5\times 10^{-3}Z_{\odot}$. In order to examine this, we prepare the
SSE-plus set. In the SSE-plus set, we add SSE-model stars with $5\times
10^{-3}Z_{\odot}$ for Pop II stars in addition to the metallicity list of the
SSE-std set.

Figure 13 compares the mass and spin distribution of binary BHs in the SSE-std
and SSE-plus sets. There is little difference between the two sets. Thus, we
can conclude that the not using SSE-model stars with $5\times
10^{-3}Z_{\odot}$ has little effect on our results.

Figure 13: (Left) The same as Figure 6 except that results of the SSE-std and
SSE-plus sets are shown in the top and bottom panels. (Right) The same as
Figure 8 except that results of the SSE-std and SSE-plus sets are shown in the
top and bottom panels.

## References

  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Physical Review X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, Phys. Rev. Lett., 125, 101102, doi: 10.1103/PhysRevLett.125.101102
  * Abbott et al. (2020b) —. 2020b, ApJ, 900, L13, doi: 10.3847/2041-8213/aba493
  * Abbott et al. (2020c) —. 2020c, Phys. Rev. D, 102, 043015, doi: 10.1103/PhysRevD.102.043015
  * Abbott et al. (2021a) —. 2021a, Physical Review X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2021b) —. 2021b, ApJ, 913, L7, doi: 10.3847/2041-8213/abe949
  * Abbott et al. (2021c) —. 2021c, SoftwareX, 13, 100658, doi: 10.1016/j.softx.2021.100658
  * Anagnostou et al. (2020a) Anagnostou, O., Trenti, M., & Melatos, A. 2020a, PASA, 37, e044, doi: 10.1017/pasa.2020.35
  * Anagnostou et al. (2020b) —. 2020b, arXiv e-prints, arXiv:2010.06161.  https://arxiv.org/abs/2010.06161
  * Antonini et al. (2014) Antonini, F., Murray, N., & Mikkola, S. 2014, ApJ, 781, 45, doi: 10.1088/0004-637X/781/1/45
  * Antonini & Perets (2012) Antonini, F., & Perets, H. B. 2012, ApJ, 757, 27, doi: 10.1088/0004-637X/757/1/27
  * Antonini & Rasio (2016) Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187, doi: 10.3847/0004-637X/831/2/187
  * Antonini et al. (2017) Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77, doi: 10.3847/1538-4357/aa6f5e
  * Arca Sedda (2020) Arca Sedda, M. 2020, ApJ, 891, 47, doi: 10.3847/1538-4357/ab723b
  * Arca Sedda et al. (2021) Arca Sedda, M., Rizzuto, F. P., Naab, T., et al. 2021, arXiv e-prints, arXiv:2105.07003.  https://arxiv.org/abs/2105.07003
  * Askar et al. (2017) Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36, doi: 10.1093/mnrasl/slw177
  * Bae et al. (2014) Bae, Y.-B., Kim, C., & Lee, H. M. 2014, MNRAS, 440, 2714, doi: 10.1093/mnras/stu381
  * Banerjee (2017) Banerjee, S. 2017, MNRAS, 467, 524, doi: 10.1093/mnras/stw3392
  * Banerjee (2018a) —. 2018a, MNRAS, 473, 909, doi: 10.1093/mnras/stx2347
  * Banerjee (2018b) —. 2018b, MNRAS, 481, 5123, doi: 10.1093/mnras/sty2608
  * Banerjee (2021a) —. 2021a, MNRAS, 500, 3002, doi: 10.1093/mnras/staa2392
  * Banerjee (2021b) —. 2021b, MNRAS, 503, 3371, doi: 10.1093/mnras/stab591
  * Banerjee et al. (2010) Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371, doi: 10.1111/j.1365-2966.2009.15880.x
  * Barkat et al. (1967) Barkat, Z., Rakavy, G., & Sack, N. 1967, Phys. Rev. Lett., 18, 379, doi: 10.1103/PhysRevLett.18.379
  * Bavera et al. (2020) Bavera, S. S., Fragos, T., Qin, Y., et al. 2020, A&A, 635, A97, doi: 10.1051/0004-6361/201936204
  * Bavera et al. (2021) Bavera, S. S., Fragos, T., Zevin, M., et al. 2021, A&A, 647, A153, doi: 10.1051/0004-6361/202039804
  * Belczynski (2020) Belczynski, K. 2020, ApJ, 905, L15, doi: 10.3847/2041-8213/abcbf1
  * Belczynski et al. (2010) Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217, doi: 10.1088/0004-637X/714/2/1217
  * Belczynski et al. (2016a) Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016a, Nature, 534, 512, doi: 10.1038/nature18322
  * Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407, doi: 10.1086/340304
  * Belczynski et al. (2017) Belczynski, K., Ryu, T., Perna, R., et al. 2017, MNRAS, 471, 4702, doi: 10.1093/mnras/stx1759
  * Belczynski et al. (2016b) Belczynski, K., Heger, A., Gladysz, W., et al. 2016b, A&A, 594, A97, doi: 10.1051/0004-6361/201628980
  * Belczynski et al. (2020a) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020a, A&A, 636, A104, doi: 10.1051/0004-6361/201936528
  * Belczynski et al. (2020b) Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020b, ApJ, 890, 113, doi: 10.3847/1538-4357/ab6d77
  * Bethe & Brown (1998) Bethe, H. A., & Brown, G. E. 1998, ApJ, 506, 780, doi: 10.1086/306265
  * Blinnikov et al. (1996) Blinnikov, S. I., Dunina-Barkovskaya, N. V., & Nadyozhin, D. K. 1996, ApJS, 106, 171, doi: 10.1086/192334
  * Bond et al. (1984) Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825, doi: 10.1086/162057
  * Bouffanais et al. (2021) Bouffanais, Y., Mapelli, M., Santoliquido, F., et al. 2021, MNRAS, 507, 5224, doi: 10.1093/mnras/stab2438
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71, doi: 10.3847/1538-4357/ab9d85
  * Bromm & Larson (2004) Bromm, V., & Larson, R. B. 2004, ARA&A, 42, 79, doi: 10.1146/annurev.astro.42.053102.134034
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115, doi: 10.1051/0004-6361/201016113
  * Cameron & Mock (1967) Cameron, A. G. W., & Mock, M. 1967, Nature, 215, 464, doi: 10.1038/215464a0
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819
  * Chon et al. (2021) Chon, S., Omukai, K., & Schneider, R. 2021, arXiv e-prints, arXiv:2103.04997.  https://arxiv.org/abs/2103.04997
  * Claeys et al. (2014) Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J., & Verbunt, F. W. M. 2014, A&A, 563, A83, doi: 10.1051/0004-6361/201322714
  * Costa et al. (2021) Costa, G., Bressan, A., Mapelli, M., et al. 2021, MNRAS, 501, 4514, doi: 10.1093/mnras/staa3916
  * Croon et al. (2020) Croon, D., McDermott, S. D., & Sakstein, J. 2020, arXiv e-prints, arXiv:2007.07889.  https://arxiv.org/abs/2007.07889
  * de Mink & Belczynski (2015) de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58, doi: 10.1088/0004-637X/814/1/58
  * Di Carlo et al. (2019) Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947, doi: 10.1093/mnras/stz1453
  * Di Carlo et al. (2020a) Di Carlo, U. N., Mapelli, M., Bouffanais, Y., et al. 2020a, MNRAS, 497, 1043, doi: 10.1093/mnras/staa1997
  * Di Carlo et al. (2020b) Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020b, MNRAS, 498, 495, doi: 10.1093/mnras/staa2286
  * Dominik et al. (2012) Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52, doi: 10.1088/0004-637X/759/1/52
  * Downing et al. (2010) Downing, J. M. B., Benacquista, M. J., Giersz, M., & Spurzem, R. 2010, MNRAS, 407, 1946, doi: 10.1111/j.1365-2966.2010.17040.x
  * Downing et al. (2011) —. 2011, MNRAS, 416, 133, doi: 10.1111/j.1365-2966.2011.19023.x
  * Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146, doi: 10.1051/0004-6361/201117751
  * El Eid & Langer (1986) El Eid, M. F., & Langer, N. 1986, A&A, 167, 274
  * Eldridge et al. (2019) Eldridge, J. J., Stanway, E. R., & Tang, P. N. 2019, MNRAS, 482, 870, doi: 10.1093/mnras/sty2714
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Estellés et al. (2021) Estellés, H., Husa, S., Colleoni, M., et al. 2021, arXiv e-prints, arXiv:2105.06360.  https://arxiv.org/abs/2105.06360
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36, doi: 10.3847/2041-8213/abbadd
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53, doi: 10.3847/1538-4357/ab518b
  * Farrell et al. (2021) Farrell, E., Groh, J. H., Hirschi, R., et al. 2021, MNRAS, 502, L40, doi: 10.1093/mnrasl/slaa196
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585, doi: 10.1086/428642
  * Fishbach & Holz (2020) Fishbach, M., & Holz, D. E. 2020, ApJ, 904, L26, doi: 10.3847/2041-8213/abc827
  * Fragione & Banerjee (2021) Fragione, G., & Banerjee, S. 2021, ApJ, 913, L29, doi: 10.3847/2041-8213/ac00a7
  * Fragione & Kocsis (2019) Fragione, G., & Kocsis, B. 2019, MNRAS, 486, 4781, doi: 10.1093/mnras/stz1175
  * Fragione et al. (2020a) Fragione, G., Loeb, A., & Rasio, F. A. 2020a, ApJ, 902, L26, doi: 10.3847/2041-8213/abbc0a
  * Fragione et al. (2020b) —. 2020b, ApJ, 895, L15, doi: 10.3847/2041-8213/ab9093
  * Fraley (1968) Fraley, G. S. 1968, Ap&SS, 2, 96, doi: 10.1007/BF00651498
  * Franciolini et al. (2021) Franciolini, G., Baibhav, V., De Luca, V., et al. 2021, arXiv e-prints, arXiv:2105.03349.  https://arxiv.org/abs/2105.03349
  * Frebel & Norris (2015) Frebel, A., & Norris, J. E. 2015, ARA&A, 53, 631, doi: 10.1146/annurev-astro-082214-122423
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91, doi: 10.1088/0004-637X/749/1/91
  * Fryer et al. (2001) Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372, doi: 10.1086/319719
  * Fujii et al. (2017) Fujii, M. S., Tanikawa, A., & Makino, J. 2017, PASJ, 69, 94, doi: 10.1093/pasj/psx108
  * Galaudage et al. (2021) Galaudage, S., Talbot, C., Nagar, T., et al. 2021, ApJ, 921, L15, doi: 10.3847/2041-8213/ac2f3c
  * Gayathri et al. (2021) Gayathri, V., Yang, Y., Tagawa, H., Haiman, Z., & Bartos, I. 2021, arXiv e-prints, arXiv:2104.10253.  https://arxiv.org/abs/2104.10253
  * Gerosa & Fishbach (2021) Gerosa, D., & Fishbach, M. 2021, Nature Astronomy, 5, 749, doi: 10.1038/s41550-021-01398-w
  * Giacobbo & Mapelli (2018) Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011, doi: 10.1093/mnras/sty1999
  * Giacobbo et al. (2018) Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959, doi: 10.1093/mnras/stx2933
  * Hamers et al. (2018) Hamers, A. S., Bar-Or, B., Petrovich, C., & Antonini, F. 2018, ApJ, 865, 2, doi: 10.3847/1538-4357/aadae2
  * Hamers et al. (2021) Hamers, A. S., Fragione, G., Neunteufel, P., & Kocsis, B. 2021, MNRAS, 506, 5345, doi: 10.1093/mnras/stab2136
  * Hamers & Safarzadeh (2020) Hamers, A. S., & Safarzadeh, M. 2020, ApJ, 898, 99, doi: 10.3847/1538-4357/ab9b27
  * Hartwig et al. (2016) Hartwig, T., Volonteri, M., Bromm, V., et al. 2016, MNRAS, 460, L74, doi: 10.1093/mnrasl/slw074
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288, doi: 10.1086/375341
  * Heger & Woosley (2002) Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532, doi: 10.1086/338487
  * Hijikawa et al. (2021) Hijikawa, K., Tanikawa, A., Kinugawa, T., Yoshida, T., & Umeda, H. 2021, MNRAS, 505, L69, doi: 10.1093/mnrasl/slab052
  * Hoang et al. (2018) Hoang, B.-M., Naoz, S., Kocsis, B., Rasio, F. A., & Dosopoulou, F. 2018, ApJ, 856, 140, doi: 10.3847/1538-4357/aaafce
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974, doi: 10.1111/j.1365-2966.2005.09087.x
  * Hong et al. (2020) Hong, J., Askar, A., Giersz, M., Hypki, A., & Yoon, S.-J. 2020, MNRAS, 498, 4287, doi: 10.1093/mnras/staa2677
  * Hong et al. (2018) Hong, J., Vesperini, E., Askar, A., et al. 2018, MNRAS, 480, 5645, doi: 10.1093/mnras/sty2211
  * Hosokawa & Omukai (2009a) Hosokawa, T., & Omukai, K. 2009a, ApJ, 691, 823, doi: 10.1088/0004-637X/691/1/823
  * Hosokawa & Omukai (2009b) —. 2009b, ApJ, 703, 1810, doi: 10.1088/0004-637X/703/2/1810
  * Hosokawa et al. (2010) Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478, doi: 10.1088/0004-637X/721/1/478
  * Hotokezaka & Piran (2017) Hotokezaka, K., & Piran, T. 2017, ApJ, 842, 111, doi: 10.3847/1538-4357/aa6f61
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543, doi: 10.1046/j.1365-8711.2000.03426.x
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897, doi: 10.1046/j.1365-8711.2002.05038.x
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Hütsi et al. (2021) Hütsi, G., Raidal, M., Vaskonen, V., & Veermäe, H. 2021, J. Cosmology Astropart. Phys, 2021, 068, doi: 10.1088/1475-7516/2021/03/068
  * Iben (1963) Iben, Jr., I. 1963, ApJ, 138, 452, doi: 10.1086/147658
  * Iglesias & Rogers (1996) Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Inayoshi et al. (2017) Inayoshi, K., Hirai, R., Kinugawa, T., & Hotokezaka, K. 2017, MNRAS, 468, 5020, doi: 10.1093/mnras/stx757
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411, doi: 10.1086/192264
  * Ivanova & Taam (2004) Ivanova, N., & Taam, R. E. 2004, ApJ, 601, 1058, doi: 10.1086/380561
  * Izzard et al. (2009) Izzard, R. G., Glebbeek, E., Stancliffe, R. J., & Pols, O. R. 2009, A&A, 508, 1359, doi: 10.1051/0004-6361/200912827
  * Kimball et al. (2021) Kimball, C., Talbot, C., Berry, C. P. L., et al. 2021, ApJ, 915, L35, doi: 10.3847/2041-8213/ac0aef
  * Kinugawa et al. (2014) Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & Nakamura, T. 2014, MNRAS, 442, 2963, doi: 10.1093/mnras/stu1022
  * Kinugawa et al. (2016) Kinugawa, T., Miyamoto, A., Kanda, N., & Nakamura, T. 2016, MNRAS, 456, 1093, doi: 10.1093/mnras/stv2624
  * Kinugawa et al. (2020) Kinugawa, T., Nakamura, T., & Nakano, H. 2020, MNRAS, 498, 3946, doi: 10.1093/mnras/staa2511
  * Kinugawa et al. (2021a) —. 2021a, MNRAS, 501, L49, doi: 10.1093/mnrasl/slaa191
  * Kinugawa et al. (2021b) —. 2021b, Progress of Theoretical and Experimental Physics, 2021, 021E01, doi: 10.1093/ptep/ptaa176
  * Kinugawa et al. (2021c) —. 2021c, MNRAS, 504, L28, doi: 10.1093/mnrasl/slab032
  * Klencki et al. (2018) Klencki, J., Moe, M., Gladysz, W., et al. 2018, A&A, 619, A77, doi: 10.1051/0004-6361/201833025
  * Kremer et al. (2020) Kremer, K., Spera, M., Becker, D., et al. 2020, ApJ, 903, 45, doi: 10.3847/1538-4357/abb945
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Kruckow et al. (2018) Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908, doi: 10.1093/mnras/sty2190
  * Kumamoto et al. (2019) Kumamoto, J., Fujii, M. S., & Tanikawa, A. 2019, MNRAS, 486, 3942, doi: 10.1093/mnras/stz1068
  * Kumamoto et al. (2020) —. 2020, MNRAS, 495, 4268, doi: 10.1093/mnras/staa1440
  * Kumamoto et al. (2021) Kumamoto, J., Fujii, M. S., Trani, A. A., & Tanikawa, A. 2021, arXiv e-prints, arXiv:2102.09323.  https://arxiv.org/abs/2102.09323
  * Kushnir et al. (2016) Kushnir, D., Zaldarriaga, M., Kollmeier, J. A., & Waldman, R. 2016, MNRAS, 462, 844, doi: 10.1093/mnras/stw1684
  * Leigh et al. (2018) Leigh, N. W. C., Geller, A. M., McKernan, B., et al. 2018, MNRAS, 474, 5672, doi: 10.1093/mnras/stx3134
  * Leung et al. (2019) Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, ApJ, 887, 72, doi: 10.3847/1538-4357/ab4fe5
  * Liu & Bromm (2020a) Liu, B., & Bromm, V. 2020a, MNRAS, 495, 2475, doi: 10.1093/mnras/staa1362
  * Liu & Bromm (2020b) —. 2020b, ApJ, 903, L40, doi: 10.3847/2041-8213/abc552
  * Liu & Bromm (2021) —. 2021, MNRAS, 506, 5451, doi: 10.1093/mnras/stab2028
  * Liu & Lai (2021) Liu, B., & Lai, D. 2021, MNRAS, 502, 2049, doi: 10.1093/mnras/stab178
  * Liu et al. (2021) Liu, B., Meynet, G., & Bromm, V. 2021, MNRAS, 501, 643, doi: 10.1093/mnras/staa3671
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, ApJ, 840, 39, doi: 10.3847/1538-4357/aa6af9
  * Magg et al. (2018) Magg, M., Hartwig, T., Agarwal, B., et al. 2018, MNRAS, 473, 5308, doi: 10.1093/mnras/stx2729
  * Magg et al. (2019) Magg, M., Klessen, R. S., Glover, S. C. O., & Li, H. 2019, MNRAS, 487, 486, doi: 10.1093/mnras/stz1210
  * Maio et al. (2010) Maio, U., Ciardi, B., Dolag, K., Tornatore, L., & Khochfar, S. 2010, MNRAS, 407, 1003, doi: 10.1111/j.1365-2966.2010.17003.x
  * Mandel & de Mink (2016) Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634, doi: 10.1093/mnras/stw379
  * Mandel & Farmer (2018) Mandel, I., & Farmer, A. 2018, arXiv e-prints, arXiv:1806.05820.  https://arxiv.org/abs/1806.05820
  * Mangiagli et al. (2019) Mangiagli, A., Bonetti, M., Sesana, A., & Colpi, M. 2019, ApJ, 883, L27, doi: 10.3847/2041-8213/ab3f33
  * Mapelli (2016) Mapelli, M. 2016, MNRAS, 459, 3432, doi: 10.1093/mnras/stw869
  * Mapelli (2021) —. 2021, Formation Channels of Single and Binary Stellar-Mass Black Holes, 4, doi: 10.1007/978-981-15-4702-7_16-1
  * Mapelli et al. (2017) Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422, doi: 10.1093/mnras/stx2123
  * Mapelli et al. (2019) Mapelli, M., Giacobbo, N., Santoliquido, F., & Artale, M. C. 2019, MNRAS, 487, 2, doi: 10.1093/mnras/stz1150
  * Mapelli et al. (2021a) Mapelli, M., Santoliquido, F., Bouffanais, Y., et al. 2021a, Symmetry, 13, 1678, doi: 10.3390/sym13091678
  * Mapelli et al. (2020) Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76, doi: 10.3847/1538-4357/ab584d
  * Mapelli et al. (2021b) Mapelli, M., Dall’Amico, M., Bouffanais, Y., et al. 2021b, MNRAS, 505, 339, doi: 10.1093/mnras/stab1334
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50, doi: 10.1051/0004-6361/201628133
  * McKernan et al. (2018) McKernan, B., Ford, K. E. S., Bellovary, J., et al. 2018, ApJ, 866, 66, doi: 10.3847/1538-4357/aadae5
  * Neijssel et al. (2019) Neijssel, C. J., Vigna-Gómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740, doi: 10.1093/mnras/stz2840
  * Ng et al. (2021) Ng, K. K. Y., Vitale, S., Farr, W. M., & Rodriguez, C. L. 2021, ApJ, 913, L5, doi: 10.3847/2041-8213/abf8be
  * Nitz & Capano (2021) Nitz, A. H., & Capano, C. D. 2021, ApJ, 907, L9, doi: 10.3847/2041-8213/abccc5
  * O’Leary et al. (2009) O’Leary, R. M., Kocsis, B., & Loeb, A. 2009, MNRAS, 395, 2127, doi: 10.1111/j.1365-2966.2009.14653.x
  * O’Leary et al. (2006) O’Leary, R. M., Rasio, F. A., Fregeau, J. M., Ivanova, N., & O’Shaughnessy, R. 2006, ApJ, 637, 937, doi: 10.1086/498446
  * Olejak et al. (2020) Olejak, A., Fishbach, M., Belczynski, K., et al. 2020, ApJ, 901, L39, doi: 10.3847/2041-8213/abb5b5
  * Omukai et al. (2010) Omukai, K., Hosokawa, T., & Yoshida, N. 2010, ApJ, 722, 1793, doi: 10.1088/0004-637X/722/2/1793
  * Omukai et al. (2005) Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627, doi: 10.1086/429955
  * Park et al. (2017) Park, D., Kim, C., Lee, H. M., Bae, Y.-B., & Belczynski, K. 2017, MNRAS, 469, 4665, doi: 10.1093/mnras/stx1015
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Pols et al. (1998) Pols, O. R., Schr”oder, K.-P., Hurley, J. R., Tout, C. A., & Eggleton, P. P. 1998, MNRAS, 298, 525, doi: 10.1046/j.1365-8711.1998.01658.x
  * Portegies Zwart & McMillan (2000) Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17, doi: 10.1086/312422
  * Portegies Zwart & Verbunt (1996) Portegies Zwart, S. F., & Verbunt, F. 1996, A&A, 309, 179
  * Punturo et al. (2010) Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Classical and Quantum Gravity, 27, 194002, doi: 10.1088/0264-9381/27/19/194002
  * Qin et al. (2018) Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28, doi: 10.1051/0004-6361/201832839
  * Rasskazov & Kocsis (2019) Rasskazov, A., & Kocsis, B. 2019, ApJ, 881, 20, doi: 10.3847/1538-4357/ab2c74
  * Rastello et al. (2019) Rastello, S., Amaro-Seoane, P., Arca-Sedda, M., et al. 2019, MNRAS, 483, 1233, doi: 10.1093/mnras/sty3193
  * Rastello et al. (2021) Rastello, S., Mapelli, M., Di Carlo, U. N., et al. 2021, MNRAS, 507, 3612, doi: 10.1093/mnras/stab2355
  * Reitze et al. (2019) Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, in BAAS, Vol. 51, 35.  https://arxiv.org/abs/1907.04833
  * Riley et al. (2021) Riley, J., Mandel, I., Marchant, P., et al. 2021, MNRAS, 505, 663, doi: 10.1093/mnras/stab1291
  * Rodriguez et al. (2018) Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., & Rasio, F. A. 2018, Phys. Rev. Lett., 120, 151101, doi: 10.1103/PhysRevLett.120.151101
  * Rodriguez & Antonini (2018) Rodriguez, C. L., & Antonini, F. 2018, ApJ, 863, 7, doi: 10.3847/1538-4357/aacea4
  * Rodriguez et al. (2016a) Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016a, Phys. Rev. D, 93, 084029, doi: 10.1103/PhysRevD.93.084029
  * Rodriguez et al. (2016b) Rodriguez, C. L., Haster, C.-J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016b, ApJ, 824, L8, doi: 10.3847/2041-8205/824/1/L8
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 100, 043027, doi: 10.1103/PhysRevD.100.043027
  * Roulet et al. (2021) Roulet, J., Chia, H. S., Olsen, S., et al. 2021, Phys. Rev. D, 104, 083010, doi: 10.1103/PhysRevD.104.083010
  * Sadowski et al. (2008) Sadowski, A., Belczynski, K., Bulik, T., et al. 2008, ApJ, 676, 1162, doi: 10.1086/528932
  * Sallaska et al. (2013) Sallaska, A. L., Iliadis, C., Champange, A. E., et al. 2013, ApJS, 207, 18, doi: 10.1088/0067-0049/207/1/18
  * Samsing (2018) Samsing, J. 2018, Phys. Rev. D, 97, 103014, doi: 10.1103/PhysRevD.97.103014
  * Samsing et al. (2018) Samsing, J., Askar, A., & Giersz, M. 2018, ApJ, 855, 124, doi: 10.3847/1538-4357/aaab52
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444, doi: 10.1126/science.1223344
  * Santoliquido et al. (2020) Santoliquido, F., Mapelli, M., Bouffanais, Y., et al. 2020, ApJ, 898, 152, doi: 10.3847/1538-4357/ab9b78
  * Santoliquido et al. (2021) Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021, MNRAS, 502, 4877, doi: 10.1093/mnras/stab280
  * Schneider et al. (2006) Schneider, R., Salvaterra, R., Ferrara, A., & Ciardi, B. 2006, MNRAS, 369, 825, doi: 10.1111/j.1365-2966.2006.10331.x
  * Silsbee & Tremaine (2017) Silsbee, K., & Tremaine, S. 2017, ApJ, 836, 39, doi: 10.3847/1538-4357/aa5729
  * Skinner & Wise (2020) Skinner, D., & Wise, J. H. 2020, MNRAS, 492, 4386, doi: 10.1093/mnras/staa139
  * Spera & Mapelli (2017) Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739, doi: 10.1093/mnras/stx1576
  * Stevenson et al. (2019) Stevenson, S., Sampson, M., Powell, J., et al. 2019, ApJ, 882, 121, doi: 10.3847/1538-4357/ab3981
  * Stevenson et al. (2017) Stevenson, S., Vigna-Gómez, A., Mandel, I., et al. 2017, Nature Communications, 8, 14906, doi: 10.1038/ncomms14906
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946, doi: 10.1093/mnras/stw2260
  * Tagawa et al. (2020a) Tagawa, H., Haiman, Z., Bartos, I., & Kocsis, B. 2020a, ApJ, 899, 26, doi: 10.3847/1538-4357/aba2cc
  * Tagawa et al. (2020b) Tagawa, H., Haiman, Z., & Kocsis, B. 2020b, ApJ, 898, 25, doi: 10.3847/1538-4357/ab9b8c
  * Tagawa et al. (2021a) Tagawa, H., Kocsis, B., Haiman, Z., et al. 2021a, ApJ, 907, L20, doi: 10.3847/2041-8213/abd4d3
  * Tagawa et al. (2021b) —. 2021b, ApJ, 908, 194, doi: 10.3847/1538-4357/abd555
  * Takahashi et al. (2019) Takahashi, K., Sumiyoshi, K., Yamada, S., Umeda, H., & Yoshida, T. 2019, ApJ, 871, 153, doi: 10.3847/1538-4357/aaf8a8
  * Takahashi et al. (2018) Takahashi, K., Yoshida, T., & Umeda, H. 2018, ApJ, 857, 111, doi: 10.3847/1538-4357/aab95f
  * Takahashi et al. (2016) Takahashi, K., Yoshida, T., Umeda, H., Sumiyoshi, K., & Yamada, S. 2016, MNRAS, 456, 1320, doi: 10.1093/mnras/stv2649
  * Tanikawa (2013) Tanikawa, A. 2013, MNRAS, 435, 1358, doi: 10.1093/mnras/stt1380
  * Tanikawa et al. (2021a) Tanikawa, A., Kinugawa, T., Yoshida, T., Hijikawa, K., & Umeda, H. 2021a, MNRAS, 505, 2170, doi: 10.1093/mnras/stab1421
  * Tanikawa et al. (2021b) Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., & Kinugawa, T. 2021b, ApJ, 910, 30, doi: 10.3847/1538-4357/abe40d
  * Tanikawa et al. (2020) Tanikawa, A., Yoshida, T., Kinugawa, T., Takahashi, K., & Umeda, H. 2020, MNRAS, 495, 4170, doi: 10.1093/mnras/staa1417
  * Team COMPAS et al. (2021) Team COMPAS, Riley, J., Agrawal, P., et al. 2021, arXiv e-prints, arXiv:2109.10352.  https://arxiv.org/abs/2109.10352
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration, The Virgo Collaboration, & The KAGRA Scientific Collaboration. 2021a, arXiv e-prints, arXiv:2111.03634.  https://arxiv.org/abs/2111.03634
  * The LIGO Scientific Collaboration et al. (2021b) The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2021b, arXiv e-prints, arXiv:2108.01045.  https://arxiv.org/abs/2108.01045
  * The LIGO Scientific Collaboration et al. (2021c) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021c, arXiv e-prints, arXiv:2111.03606.  https://arxiv.org/abs/2111.03606
  * Toonen et al. (2012) Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70, doi: 10.1051/0004-6361/201218966
  * Trani et al. (2019) Trani, A. A., Spera, M., Leigh, N. W. C., & Fujii, M. S. 2019, ApJ, 885, 135, doi: 10.3847/1538-4357/ab480a
  * Trani et al. (2021) Trani, A. A., Tanikawa, A., Fujii, M. S., Leigh, N. W. C., & Kumamoto, J. 2021, MNRAS, 504, 910, doi: 10.1093/mnras/stab967
  * Umeda & Nomoto (2002) Umeda, H., & Nomoto, K. 2002, ApJ, 565, 385, doi: 10.1086/323946
  * Umeda et al. (2020) Umeda, H., Yoshida, T., Nagele, C., & Takahashi, K. 2020, ApJ, 905, L21, doi: 10.3847/2041-8213/abcb96
  * van den Heuvel et al. (2017) van den Heuvel, E. P. J., Portegies Zwart, S. F., & de Mink, S. E. 2017, MNRAS, 471, 4256, doi: 10.1093/mnras/stx1430
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * VanLandingham et al. (2016) VanLandingham, J. H., Miller, M. C., Hamilton, D. P., & Richardson, D. C. 2016, ApJ, 828, 77, doi: 10.3847/0004-637X/828/2/77
  * Vardya (1960) Vardya, M. S. 1960, ApJS, 4, 281, doi: 10.1086/190051
  * Vink et al. (2021) Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146, doi: 10.1093/mnras/stab842
  * Vynatheya & Hamers (2021) Vynatheya, P., & Hamers, A. S. 2021, arXiv e-prints, arXiv:2110.14680.  https://arxiv.org/abs/2110.14680
  * Wang et al. (2021) Wang, L., Fujii, M. S., & Tanikawa, A. 2021, MNRAS, 504, 5778, doi: 10.1093/mnras/stab1157
  * Weatherford et al. (2021) Weatherford, N. C., Fragione, G., Kremer, K., et al. 2021, ApJ, 907, L25, doi: 10.3847/2041-8213/abd79c
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355, doi: 10.1086/161701
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244, doi: 10.3847/1538-4357/836/2/244
  * Woosley (2019) —. 2019, ApJ, 878, 49, doi: 10.3847/1538-4357/ab1b41
  * Woosley et al. (2007) Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390, doi: 10.1038/nature06333
  * Yang et al. (2019) Yang, Y., Bartos, I., Haiman, Z., et al. 2019, ApJ, 876, 122, doi: 10.3847/1538-4357/ab16e3
  * Yoon et al. (2010) Yoon, S. C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940, doi: 10.1088/0004-637X/725/1/940
  * Yoshida et al. (2019) Yoshida, T., Takiwaki, T., Kotake, K., et al. 2019, ApJ, 881, 16, doi: 10.3847/1538-4357/ab2b9d
  * Yoshida et al. (2016) Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS, 457, 351, doi: 10.1093/mnras/stv3002
  * Zahn (1975) Zahn, J. P. 1975, A&A, 41, 329
  * Zahn (1977) —. 1977, A&A, 500, 121
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152, doi: 10.3847/1538-4357/abe40e
  * Ziegler & Freese (2021) Ziegler, J., & Freese, K. 2021, Phys. Rev. D, 104, 043015, doi: 10.1103/PhysRevD.104.043015
  * Ziosi et al. (2014) Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703, doi: 10.1093/mnras/stu824