# Massive binary black holes from Population II and III stars

Guglielmo Costa1,2,3,4 , Michela Mapelli1,2,3, Giuliano Iorio1,2,3, Filippo
Santoliquido1,2, Gastón J. Escobar1,2, Ralf S. Klessen5, and Alessandro
Bressan6,3
1Dipartimento di Fisica e Astronomia Galileo Galilei, Università di Padova,
Vicolo dell’Osservatorio 3, I–35122 Padova, Italy
2INFN-Padova, Via Marzolo 8, I–35131 Padova, Italy
3INAF-Padova, Vicolo dell’Osservatorio 5, I–35122 Padova, Italy
4Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230 Saint-Genis-Laval, France
5Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Str. 2, D–69120 Heidelberg, Germany
6SISSA, via Bonomea 365, I–34136 Trieste, Italy E-mail:
guglielmo.costa.astro<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Population III stars, born from the primordial gas in the Universe, lose a
negligible fraction of their mass via stellar winds and possibly follow a top-
heavy mass function. Hence, they have often been regarded as the ideal
progenitors of massive black holes (BHs), even above the pair instability mass
gap. Here, we evolve a large set of Population III binary stars (metallicity
$Z=10^{-11}$) with our population-synthesis code sevn, and compare them with
Population II binary stars ($Z=10^{-4}$). In our models, the lower edge of the
pair-instability mass gap corresponds to a BH mass of $\approx{86}$
($\approx{91}$) M⊙ for single Population III (II) stars. Overall, we find only
mild differences between the properties of binary BHs (BBHs) born from
Population III and II stars, especially if we adopt the same initial mass
function and initial orbital properties. Most BBH mergers born from Population
III and II stars have primary BH mass below the pair-instability gap, and the
maximum secondary BH mass is $<50$ M⊙. Only up to $\approx{3.3}$%
($\approx{0.09}$%) BBH mergers from Population III (II) progenitors have
primary mass above the gap. Unlike metal-rich binary stars, the main formation
channel of BBH mergers from Population III and II stars involves only stable
mass transfer episodes in our fiducial model.

###### keywords:

black hole physics – stars: Population II – stars: Population III –
gravitational waves – methods: numerical

††pubyear: 2023††pagerange: Massive binary black holes from Population II and
III stars–Massive binary black holes from Population II and III stars

##  1 Introduction

Population III (hereafter, Pop. III) stars formed from metal-free primordial
gas in the early Universe, and have eluded any attempt to observe them to date
(e.g., Bromm & Larson, 2004; Klessen & Glover, 2023, for a review). Their
initial mass function (IMF) is commonly considered to be more top-heavy than
that of metal-rich stars, mostly because molecular hydrogen is an inefficient
coolant with respect to dust (e.g., Bromm & Larson, 2004; Schneider et al.,
2006; Stacy & Bromm, 2013; Susa et al., 2014; Hirano et al., 2014, 2015;
Wollenberg et al., 2020; Chon et al., 2021; Tanikawa et al., 2021b; Jaura et
al., 2022; Prole et al., 2022). Massive Pop. III stars lose a negligible
fraction of their mass during their life, because stellar winds are highly
inefficient for a nearly metal-free chemical composition (e.g., Woosley et
al., 2002; Volpato et al., 2023). If Pop. III stars avoid pair instability
(Woosley, 2017), they might thus end their life with a direct collapse,
leading to the formation of massive black holes (BHs, e.g., Woosley et al.,
2002). For this reason, Pop. III stars have been extensively studied (e.g.,
Hartwig et al., 2016; Belczynski et al., 2017; Kinugawa et al., 2020; Tanikawa
et al., 2021b, 2022a) as possible progenitors of the most massive BHs observed
by the LIGO–Virgo–KAGRA (LVK) collaboration (Abbott et al., 2020a, b, 2021a,
2021b).

Population II (hereafter, Pop. II) stars formed from material that was already
enriched in metals by Pop. III stars. With a metallicity ranging from
$Z\sim{10^{-6}}$ to a few $\times{}10^{-4}$ (here and in the following, we
define $Z$ as the mass fraction of elements heavier than helium, in absolute
values) Pop. II stars are way more common in the Universe than Pop. III stars
(e.g., Smith et al., 2015): we observe them in metal-poor globular clusters,
as well as in the halo of the Milky Way and in some metal-poor dwarf galaxies
(e.g., Frebel et al., 2007; Frebel & Norris, 2015). It is still unclear
whether Pop. II stars follow the same IMF as metal-rich stars (e.g., Schneider
et al., 2012; Chiaki et al., 2018; Chon et al., 2021; Sharda & Krumholz,
2022). Their metal content is still sufficiently low that stellar winds are
heavily quenched in Pop. II stars, too (e.g., Chen et al., 2015). Thus,
massive Pop. II stars might also collapse leaving massive compact remnants at
the end of their life, but their contribution to the population of BHs and
intermediate-mass BHs has been less investigated than that of Pop. III stars
(Spera & Mapelli, 2017; Renzo et al., 2020).

Both massive Pop. III and Pop. II stars are expected to undergo pair
instability or pulsational pair instability if their central temperature and
density lead to an efficient production of electron and positron pairs (Fowler
& Hoyle, 1964; Barkat et al., 1967; Rakavy & Shaviv, 1967; Woosley et al.,
2007). If the helium core mass grows to $\sim{64-135}$ M⊙ at the end of carbon
burning, the star is expected to be completely disrupted by a pair instability
supernova, leaving no compact remnant, while higher He-core masses enable the
direct collapse of the star to a BH (Woosley et al., 2002). For He-core masses
in the range $\sim{32-64}$ M⊙ (Woosley, 2017, 2019; Marchant et al., 2019;
Farmer et al., 2019), pair instability triggers pulsations of the star, which
enhance mass loss and, in the end, allow the star to find a stable
configuration. While the boundaries of pair instability and the final compact
remnant masses are still highly uncertain (e.g., Leung et al., 2019; Farmer et
al., 2019, 2020; Marchant et al., 2019; Stevenson et al., 2019; Renzo et al.,
2020; Marchant et al., 2021; Mapelli et al., 2020; Costa et al., 2021; Woosley
& Heger, 2021; Vink et al., 2021), this process has a key impact on the final
population of binary BHs (BBHs) born from metal-free and metal-poor stars.

Here, we model a population of BHs and BBHs born from Pop. III and Pop. II
stars. Our Pop. III (II) star models assume a metallicity $Z=10^{-11}$
($Z=10^{-4}$). We probe a large range of initial configurations for the IMF
and orbital parameters of Pop. III and II binary stars. We show that the
differences between the two BH populations are subtle. Both Pop. II and III
stars can give birth to very massive BHs above the pair instability mass gap.
However, most BBH mergers born via isolated binary evolution host BHs below
the pair-instability mass gap. When the initial semi-major axis distribution
is skewed toward small values ($<10^{3}$ R⊙), the vast majority of BBH mergers
originate from Pop. III and Pop. II binary stars that evolve only via stable
mass transfer, without common envelope. In a companion paper (Santoliquido et
al., 2023), we explore the impact of these models on the cosmic merger rate of
BBHs.

This paper is structured as follows. In Section 2 we describe our stellar
tracks and population synthesis simulations. Section 3 summarizes our main
results, that we discuss in Section 4 by considering the main formation
channels of the simulated BBH mergers. We draw our main conclusions in Section
5.

##  2 Methods

###  2.1 Binary population synthesis code (sevn)

In this work, we use the sevn code version 2 (Iorio et al., 2022). sevn
integrates the evolution of stellar properties (e.g., total mass, photospheric
radius, He and carbon-oxygen, CO, core mass and radius, luminosity) by
interpolating a set of stellar tracks (Spera & Mapelli, 2017), and models the
main binary evolution processes (mass transfer via stellar winds, Roche lobe
overflow, common envelope evolution, tides, and gravitational-wave decay)
according to the semi-analytic formalism presented in Hurley et al. (2002),
with several updates described in Iorio et al. (2022). In the following, we
adopt the same set up as the fiducial model presented in Iorio et al. (2022).
In particular, Roche-lobe overflow mass transfer is always stable for main
sequence (MS) and Hertzsprung gap donor stars, while we follow the
prescriptions by Hurley et al. (2002) in all the other cases. We model the
mass accretion rate during Roche-lobe overflow as

$\dot{M}_{\mathrm{a}}=\left\\{\begin{array}[]{ll}\min{(\dot{M}_{\rm
Edd},\,{}-f_{\mathrm{MT}}\,{}\dot{M}_{\mathrm{d}})}&\mathrm{if}\,{}{\rm
the}\,{}\mathrm{accretor}\,{}\mathrm{is}\,{}\mathrm{a}\,{}\mathrm{compact}\,{}\mathrm{object},\\\
-f_{\mathrm{MT}}\,{}\dot{M}_{\mathrm{d}}&\mathrm{otherwise},\end{array}\right.$
(1)

where $\dot{M}_{\rm Edd}$ is the Eddington rate (Eq. 67 of Hurley et al.
2002), $\dot{M}_{\mathrm{d}}$ is the mass-loss rate of the donor star, and
$f_{\rm MT}\in[0,1]$ is the mass accretion efficiency; here, we use $f_{\rm
MT}=0.5$. Furthermore, we assume that the mass not accreted during the Roche-
lobe overflow is lost from the vicinity of the accretor as an isotropic wind
(isotropic re-emission). At the onset of Roche-lobe overflow, sevn
circularises the orbit at periastron.

To model a common-envelope phase, we assume an efficiency parameter
$\alpha_{\rm CE}=1$ and estimate the envelope binding energy using the same
$\lambda_{\mathrm{CE}}$ formalism as in Claeys et al. (2014).

We model the final fate of intermediate-mass and high-mass stars as described
in Iorio et al. (2022). In particular, we use the rapid formalism by Fryer et
al. (2012) for core-collapse supernovae, we model electron-capture supernovae
as in Giacobbo & Mapelli (2019), and (pulsational) pair-instability supernovae
according to Mapelli et al. (2020). Compact objects receive a natal kick at
their birth. In our models, we randomly draw the natal kick magnitude as
(Giacobbo & Mapelli, 2020):

$V_{\rm kick}=f_{\rm H05}\,{}\frac{\langle{}M_{\rm NS}\rangle{}}{M_{\rm
rem}}\,{}\frac{M_{\rm ej}}{\langle{}M_{\rm ej}\rangle},$ (2)

where $\langle{}M_{\rm NS}\rangle{}$ and $\langle{}M_{\rm ej}\rangle$ are the
average neutron star mass and ejecta mass from single stellar evolution,
respectively, while $M_{\rm rem}$ and $M_{\rm ej}$ are the compact object mass
and the ejecta mass (Giacobbo & Mapelli, 2020). The term $f_{\rm H05}$ is a
random number drawn from a Maxwellian distribution with one-dimensional root
mean square $\sigma_{\mathrm{kick}}=265\ \mathrm{km}\,{}\mathrm{s}^{-1}$,
coming from a fit to the proper motions of 73 young pulsars ($<3$ Myr) in the
Milky Way (Hobbs et al., 2005). In this formalism, stripped and ultra-stripped
supernovae result in lower kicks with respect to the other explosions, owing
to the lower amount of ejected mass $M_{\rm ej}$ (Bray & Eldridge, 2016,
2018). BHs originating from a direct collapse receive zero natal kicks from
this mechanism. We report the sevn input parameter list in Zenodo at
https://doi.org/10.5281/zenodo.7736310 (Costa et al., 2023).

###  2.2 Tracks and single star evolution

Figure 1: Hertzsprung–Russell (HR) diagram of Pop. III (left) and Pop. II
(right) tracks. The black thick lines show the evolution of selected tracks
with $M_{\mathrm{ZAMS{}{}}}$= 2, 5, 10, 14, 24, 40, 100, 200, and 600 M⊙. All
the other tracks are shown with solid grey lines. The dashed black line
indicates the ZAMS. The red dashed line indicates the end of the MS. The
orange stars (circles) mark the beginning (end) of core He burning. The blue
circles indicate the final position of the star in the diagram before the
supernova. Diagonal grey dashed lines correspond to constant stellar radii in
R⊙. Figure 2: Radius versus age of some selected massive stars
($\mbox{$M_{\mathrm{ZAMS{}{}}}$}\geq 10~{}\mbox{M${}_{\odot}$}$). The solid
black (dashed red) lines indicate Pop. III (Pop. II) stars. The orange stars
(circles) mark the beginning (end) of core He burning. The blue circles
indicate the final radius of the star before the supernova.

We generated sets of Pop. III and Pop. II stellar tracks with the parsec code
(Bressan et al., 2012; Costa et al., 2021, 2022; Nguyen et al., 2022). Pop.
III stars have typical behaviours, which Pop. II stars do not show in their
evolution (Cassisi & Castellani, 1993; Marigo et al., 2001; Murphy et al.,
2021). For instance, during the early-MS, Pop. III stars cannot ignite the CNO
tri-cycle because of the initial lack of carbon, nitrogen and oxygen. Pop. III
stars need very high central temperatures to reach pressure support just with
the energy provided by the proton-proton (pp) chain. Depending on the stellar
mass, the central temperature becomes so high that some carbon could be
synthesized from the triple-$\alpha$ reaction (i.e. helium burning), even in
the MS. This leads the CNO tri-cycle to ignite and suddenly replace the pp
chain as the main source of energy of the star (Marigo et al., 2001; Murphy et
al., 2021). Moreover, due to the high central temperatures reached at the end
of the MS, Pop. III stars have a smoother transition to the core helium
burning (CHeB) phase with respect to more metal-rich stars. These
characteristic features of Pop. III stars evolution arise at metallicity
$Z\lesssim 10^{-10}$ (Cassisi & Castellani, 1993). Hence, here we assume a
metallicity $Z=10^{-11}$ for Pop. III stars (see also Tanikawa et al., 2021b).
For Pop. II stars, we take a metallicity $Z=10^{-4}$.

We adopt the Caffau et al. (2011) solar partition of chemical elements. Each
set has an initial mass at the zero-age main sequence (ZAMS),
$M_{\mathrm{ZAMS{}{}}}$, which ranges from 2 to 600 M⊙. All tracks evolve
until advanced evolutionary phases. Stars with $2\leq M_{\rm
ZAMS}/\mbox{M${}_{\odot}$}{}<10$ reach the early- asymptotic giant branch
(AGB) phase (post-core helium-burning phase characterized by the burning of
helium in a shell above the CO core). Stars with $M_{\rm ZAMS}\geq
10\,{}\mbox{M${}_{\odot}$}$ evolve until advanced phases of the core O burning
or the beginning of the pair-instability regime. The Pop. III set of tracks
extends the collection of parsec tracks, already used in sevn (Iorio et al.,
2022), and will soon be made available
online111http://stev.oapd.inaf.it/PARSEC.

We also computed new tracks of pure-He stars to extend the database to lower
metallicities. We use tracks of pure-He stars to model naked-He stars formed
via stripping during mass transfer at low metallicity (Kruckow et al., 2018;
Spera et al., 2019; Mapelli et al., 2020; Iorio et al., 2022; Agrawal et al.,
2023). The metallicity adopted is $Z=10^{-6}$, and the masses range from
$M_{\rm ZAMS,He}=0.36$ to $350\,{}\mbox{M${}_{\odot}$}$. This metallicity is
similar to the metal content we find in He cores of Pop. III stars after the
MS phase. For instance, depending on the initial mass, the carbon mass
fraction in the He cores goes from $7\times 10^{-7}$ to $2.6\times 10^{-6}$
for 2 M⊙ and 600 M⊙, respectively. Therefore, we do not expect to have
completely metal-free pure-He stars.

All the new tracks are computed with the same physical set-up as described in
Costa et al. (2021, 2022) for stellar winds, nuclear reaction network,
opacities and equation of state. In Pop. III tracks, stellar winds are
naturally quenched due to the lack of metals (particularly iron) which we
accounted for (see Chen et al., 2015). Concerning the stellar convection, we
adopt the Schwarzschild criterion (Schwarzschild, 1958) for defining the
unstable region and the mixing-length theory (Böhm-Vitense, 1958) with a
solar-calibrated value for the parameter $\alpha_{\mathrm{MLT}}=1.74$ (Bressan
et al., 2012). Above the convective core, we adopt a penetrative overshooting
with a characteristic parameter of $\lambda_{\mathrm{ov}}=0.5$ in units of
pressure scale height, computed with the ballistic approach (Bressan et al.,
1981). In this framework, $\lambda_{\mathrm{ov}}$ is the mean free path of the
unstable element across the border of the convective region, and its value
corresponds to about $f_{\mathrm{ov}}=0.025$ in the exponential decay
overshooting formalism (in the diffusive scheme, Herwig et al., 1997). We also
included undershooting at the bottom of the convective envelope, with a
characteristic distance of $\Lambda_{\mathrm{env}}=0.7$ in pressure scale
heights. The undershooting could play a role in the ending fate of massive
stars, triggering dredge-up episodes which may stabilize the star against the
pair-instability (Costa et al., 2021). In the interpolation methods used in
sevn (described in detail in Iorio et al., 2022), we cut the evolution just
before the early-AGB or the ignition of the core C burning.

Figure 1 shows the two sets of tracks used in this work, Pop. III and II
stars, in the Hertzsprung–Russell (HR) diagram. Pop. III stars begin their
life as metal-free objects and, in the ZAMS, are more compact and hotter than
their Pop. II counterparts.

Both Pop. III and Pop. II stars evolve towards the red part of the diagram
during the MS. Figure 1 also shows a clear trend of the star position at the
end of the MS as a function of the initial stellar mass. Pop. III stars with
an $\mbox{$M_{\mathrm{ZAMS{}{}}}$}\leq 200$ M⊙ end the MS as blue super-giant
(BSG) stars, while, stars with $\mbox{$M_{\mathrm{ZAMS{}{}}}$}>200$ M⊙
complete the MS as yellow super-giant (YSG) or red super-giant (RSG) stars.
This trend is similar in Pop. II stars, but with a lower transition mass
(about 150 M⊙). Such a trend for Pop. III stars has also been found by other
authors (e.g., Tanikawa et al., 2021b), but for higher initial masses
($\mbox{$M_{\mathrm{ZAMS{}{}}}$}>600$ M⊙). Stars that become RSG during the MS
develop large convective envelopes, differently from stars that remain BSG,
which have mostly radiative envelopes. The transition mass that separates the
two evolutionary pathways depends on the convection treatment adopted. This
peculiar evolution of the most massive Pop. III and Pop. II stars can
dramatically affect the evolution of a binary system since they become giant
stars with very large radii during the MS, in which there is still no well-
defined transition between the core and the envelope. Therefore, binary
interactions in such cases may lead to early mergers.

The post-MS evolution of Pop. III stars shows several features in the HR
diagram, which depend on $M_{\mathrm{ZAMS{}{.}}}$ Due to the high central
temperatures during the MS, all the tracks ignite helium as BSG stars shortly
after the end of the MS. Stars with a mass $\mbox{$M_{\mathrm{ZAMS{}{}}}$}\leq
40~{}\mbox{M${}_{\odot}$}$ end the CHeB phase in the blue side of the HR
diagram. After the CHeB phase, stars with an initial mass between
$2\leq\mbox{$M_{\mathrm{ZAMS{}{}}}$}/\mbox{M${}_{\odot}$}\leq 10$ evolve to
the asymptotic-giant branch. Stars above 10 M⊙ evolve through all the advanced
phases up to the oxygen burning, but they die with different final
configurations. Stars in the mass range
$10<\mbox{$M_{\mathrm{ZAMS{}{}}}$}/\mbox{M${}_{\odot}$}<14$ move to the red
part of the HR diagram and explode as RSG. Stars in the mass range
$14\leq\mbox{$M_{\mathrm{ZAMS{}{}}}$}/\mbox{M${}_{\odot}$}\leq 40$ die as BSG.
Similar behaviour in this mass range has been found by other authors (e.g.
Marigo et al., 2001; Tanikawa et al., 2021b). Stars in the mass range
$40<\mbox{$M_{\mathrm{ZAMS{}{}}}$}/\mbox{M${}_{\odot}$}\leq 100$ deplete
helium in the red part of the diagram and finish their evolution as RSGs.
Finally, stars with $\mbox{$M_{\mathrm{ZAMS{}{}}}$}>100$ M⊙ ignite helium as
RSGs, become BSGs during the CHeB, and remain BSGs until their final fate.

Concerning Pop. II stars, intermediate-mass stars
($2\leq\mbox{$M_{\mathrm{ZAMS{}{}}}$}/\mbox{M${}_{\odot}$}<10$) do a blue loop
during the CHeB, and then, after helium depletion, move to the AGB phase.
Stars in the mass range
$10<\mbox{$M_{\mathrm{ZAMS{}{}}}$}/\mbox{M${}_{\odot}$}\leq 100$ burn helium
as BSGs before moving to the red part of the diagram. While stars with an
initial mass $\mbox{$M_{\mathrm{ZAMS{}{}}}$}\geq 100~{}\mbox{M${}_{\odot}$}$
ignite helium as RSGs, burn it as BSGs (blue loop again), and later move back
to the red. All massive stars ($\mbox{$M_{\mathrm{ZAMS{}{}}}$}\geq
10~{}\mbox{M${}_{\odot}$}$) burn all the elements up to oxygen and finally
explode as RSGs.

Figure 2 shows the comparison between the radius evolution of Pop. III and
Pop. II stars. Pop. III stars with $\mbox{$M_{\mathrm{ZAMS{}{}}}$}\leq
100~{}\mbox{M${}_{\odot}$}$ evolve and reach the RSG stage before Pop. II
stars. Pop. III and Pop. II stars with a mass of $\geq{}200$ M⊙ become RSGs
with comparable timescales. Figure 2 shows that Pop. III stars with mass
$\mbox{$M_{\mathrm{ZAMS{}{}}}$}>100~{}\mbox{M${}_{\odot}$}$ reach the pre-
supernova stage as compact BSGs stars, whereas Pop. II stars explode as very
large RSG stars. Pop. III stars evolve to the final pre-supernova stage faster
than Pop. II stars.

###  2.3 Binary initial conditions

In this Section, we describe the initial conditions used in this work for
computing binary-population catalogues.

Table 1: Initial conditions.

Model | $M_{\rm ZAMS,1}$ | $M_{\rm ZAMS}$ | Mass ratio $q$ | Period $P$ | Eccentricity $e$ |  $N$ [$\times{}10^{7}$] | Total mass [$\times{}10^{9}$ M⊙]
---|---|---|---|---|---|---|---
LOG1 | Flat in log | – | S12 | S12 | S12 | 1.45 | 2.59
LOG2 | Flat in log | – | S12 | SB13 | Thermal | 1.45 | 2.58
LOG3 | – | Flat in log | Sorted | S12 | S12 | 1.38 | 3.19
LOG4 | Flat in log | – | SB13 | S12 | Thermal | 1.53 | 2.60
LOG5 | Flat in log | – | SB13 | SB13 | Thermal | 1.53 | 2.60
KRO1 | K01 | – | S12 | S12 | S12 | 5.23 (2.00$\dagger$) | 1.35 (0.89$\dagger$)
KRO5 | K01 | – | SB13 | SB13 | Thermal | 6.11 (2.00$\dagger$) | 1.52 (0.93$\dagger$)
LAR1 | L98 | – | S12 | S12 | S12 | 2.00 | 1.20
LAR5 | L98 | – | SB13 | SB13 | Thermal | 2.27 (2.00$\dagger$) | 1.30 (1.24$\dagger$)
TOP1 | Top heavy | – | S12 | S12 | S12 | 1.05 | 4.16
TOP5 | Top heavy | – | SB13 | SB13 | Thermal | 1.07 | 4.03

Column 1 reports the model name. Columns 2 describes how we generate the ZAMS
mass of the primary star (i.e., the most massive of the two members of the
binary system). Column 3 describes how we generate the ZAMS mass of the
overall stellar population (without differentiating between primary and
secondary stars). We follow this procedure only for model LOG3 (see the text
for details). Columns 3, 4, and 5 specify the distributions we used to
generate the mass ratios, the orbital periods and the orbital eccentricity.
See Section 2.3 for a detailed description of such distributions. The last two
columns report the total number and the total mass of the of simulated
binaries, respectively. $\dagger$The ICs for such models are under-sampled,
the actual number of simulated systems and their total mass is reported in
parenthesis. See main text for additional details.

####  2.3.1 Initial mass function (IMF)

There is still no consensus about the initial mass function (IMF) of Pop. III
stars, although several papers suggest that it should be rather top-heavy with
respect to that of local stars (e.g., Abel et al., 2002; Bromm & Larson, 2004;
Yoshida et al., 2006; Bromm, 2013; Glover, 2013). The transition between a
top-heavy and a bottom-heavy mass function likely happened in the metallicity
range of Pop. II stars (Chon et al., 2021; Sharda & Krumholz, 2022). Here,
given the uncertainties, we consider the same set of possible IMFs for both
Pop. III and Pop. II stars, as follows.

Figure 3: Initial conditions for the models LOG1–5. From upper left to bottom
right: ZAMS mass of the primary star $M_{\rm ZAMS,1}$, mass ratio $q=M_{\rm
ZAMS,2}/M_{\rm ZAMS,1}$, initial orbital period $P$, and initial orbital
eccentricity $e$. In the upper-left panel, we do not show models LOG2, LOG4
and LOG5 for simplicity, because they follow the same distribution as LOG1
(i.e., $M_{\rm ZAMS,1}$ is sampled from a flat-in-log distribution), while in
model LOG3 we sample the entire stellar population (i.e., both primary and
secondary stars) from a flat-in-log distribution. In the upper-righ panel
($q$), model LOG2 follows the same distribution as LOG1 (S12), while model
LOG5 follows the same distribution as LOG4 (SB13). In the lower-left panel
($P$), models LOG3 and LOG4 follow the same distribution as LOG1 (S12), while
models LOG2 and LOG5 are sampled from SB13. Finally, in the lower-right panel
($e$), models LOG1 and LOG3 follow S12, while the other models adopt a thermal
distribution (we show only LOG5 for simplicity). Table 1 describes the models
in detail. Figure 4: Distribution of the primary ZAMS mass (left) and mass
ratio (right) of models LOG1, KRO1, LAR1, and TOP1. These models follow the
same initial distribution for the other orbital parameters ($P$) and ($e$).
Table 1 describes the models in detail.

  * •

A flat-in-log distribution (Stacy & Bromm, 2013; Susa et al., 2014; Hirano et
al., 2014, 2015; Wollenberg et al., 2020; Chon et al., 2021; Tanikawa et al.,
2021b; Jaura et al., 2022; Prole et al., 2022):

$\xi(M_{\rm ZAMS})\propto M_{\rm ZAMS}^{-1}.$ (3)

This IMF will be our fiducial model for Pop. III stars.

  * •

A Kroupa (2001) IMF (hereafter, K01):

$\xi(M_{\rm ZAMS})\propto M_{\rm ZAMS}^{-2.3}.$ (4)

This mass function is often adopted for stars in the low-redshift Universe and
will be our fiducial model for Pop. II stars. With respect to the original
K01, which has a flatter slope for $M_{\rm ZAMS}<0.5$ M⊙, here we assume a
single slope because we do not generate ZAMS masses $<5$ M⊙ from this
distribution.

  * •

A Larson (1998) distribution (hereafter, L98):

$\xi(M_{\rm ZAMS})\propto M_{\rm ZAMS}^{-2.35}e^{-M_{\mathrm{cut1}}/M_{\rm
ZAMS}},$ (5)

where $M_{\rm cut1}=20$ M⊙ (Valiante et al., 2016).

  * •

A top-heavy distribution (Stacy & Bromm, 2013; Jaacks et al., 2019; Liu &
Bromm, 2020a):

$\xi(M_{\rm ZAMS})\propto M_{\rm ZAMS}^{-0.17}e^{-M_{\mathrm{cut2}}^{2}/M_{\rm
ZAMS}^{2}},$ (6)

where $M_{\mathrm{cut2}}=20$ M⊙.

We label our input catalogues by taking the IMF name and adding a number that
indicates the mass ratios, periods and eccentricity combinations. Therefore
all the LOG catalogues adopt the flat-in-log IMF, the KRO ones use the K01,
the LAR ones use the L98, and finally, the TOP catalogues adopt the top-heavy
IMF. Table 1 lists the properties of our input catalogues. In all our models
but LOG3 (Table 1), we randomly sample the ZAMS mass of the primary star
$M_{\rm ZAMS,1}$ (i.e., the ZAMS mass of the most massive member) of the
binary system in the range $[5,\,{}550]$ M⊙, according to one of the four
aforementioned distributions. We then randomly sample the ZAMS mass of the
secondary star ($M_{\rm ZAMS,2}$) based on the mass ratio distributions
described in the following section (Section 2.3.2). We assume that the
secondary mass can be as low as $2.2$ M⊙.

In model LOG3, we instead randomly sample the entire IMF in the range $M_{\rm
ZAMS}\in{[5,\,{}550]}$ M⊙, according to the LOG distribution. We then randomly
pair the generated stellar masses. The primary (secondary) star is thus the
component with the higher (lower) initial mass (see model sorted in Section
2.3.2).

####  2.3.2 Mass ratio and secondary mass

We draw the ZAMS mass of the secondary star ($M_{\rm ZAMS,\,{}2}$) according
to three different distributions.

  * •

We use the distribution of the mass ratio ($q=M_{\rm ZAMS,2}/M_{\rm ZAMS,1}$)
from Sana et al. (2012, hereafter S12):

$\xi(q)\propto q^{-0.1}\;{\rm with}\;q\in[0.1,1]\;{\rm and}\;M_{\rm
ZAMS,2}\geq 2.2~{}{\rm M}_{\odot}.$ (7)

This distribution is a fit to the mass ratio of O- and B-type binary stars in
the local Universe (Sana et al., 2012).

  * •

In the sorted distribution, we draw the ZAMS mass of the entire star
population from the same IMF, and then we randomly pair two stars from this
distribution, imposing that $M_{\rm ZAMS,2}\leq{}M_{\rm ZAMS,1}$. In this
model, the minimum mass of the secondary is equal to that of the primary (5
M⊙) by construction.

  * •

The mass ratio distribution by Stacy & Bromm (2013, hereafter SB13):

$\xi(q)\propto{}q^{-0.55}\;{\rm with}\;q\in[0.1,1]\;{\rm and}\;M_{\rm
ZAMS,2}\geq 2.2~{}{\rm M}_{\odot}.$ (8)

This distribution was obtained from a fit to Pop. III stars formed in
cosmological simulations (SB13).

####  2.3.3 Orbital period

We consider two different distributions for the initial orbital period ($P$):

  * •

The distribution derived by S12 for O- and B- stars in the local Universe:

$\xi(\pi)\propto{}\pi^{-0.55}\quad{\rm with}\quad\pi=\log(P/{\rm
day})\in[0.15,5.5].$ (9)

  * •

The period distribution found by SB13:

$\xi(\pi)\propto{}\exp{\left[-(\pi-\mu)^{2}/(2\,{}\sigma{}^{2})\right]}.$ (10)

This is a Gaussian distribution with $\mu{}=5.5$, and $\sigma{}=0.85$,
favouring long periods with respect to the S12 distribution. While this
distribution is likely affected by numerical resolution, which reduces the
number of systems with short orbital periods, we decide to consider it as a
robust upper limit to the orbital period of Pop. III and II binary stars.

####  2.3.4 Eccentricity

We compare two distributions for the orbital eccentricity:

  * •

the distribution obtained by S12 and based on a sample of O- and B-type stars
in the local Universe:

$\xi(e)\propto e^{-0.42}\quad{\rm with}\;e\in[0,1).$ (11)

  * •

The thermal distribution, adopted for Pop. III binaries by, e.g., Kinugawa et
al. (2014); Hartwig et al. (2016); Tanikawa et al. (2021b):

$\xi(e)=2\,{}e\quad{\rm with}\;e\in[0,1).$ (12)

####  2.3.5 Input catalogues

We build 11 different input catalogues by varying the aforementioned
distributions of the IMF, $q$, $P$, and $e$. For each of these catalogues, we
consider the two metallicities for Pop. III and Pop. II, i.e. $Z=10^{-11}$ and
$10^{-4}$, respectively.

We set the total number of generated binaries to obtain $10^{7}$ binaries in
the high-mass regime ($M_{\rm ZAMS,2}\geq 10~{}{\rm M}_{\odot}$, and $M_{\rm
ZAMS,1}\geq 10~{}{\rm M}_{\odot}$ by construction). For the models that draw
the primary mass from K01 and L98 (Section 2.3.1), we limit the total number
of generated binaries to 20M: 10M in the high-mass range, and 10M in the low-
mass range ($M_{\rm ZAMS,2}\leq 10~{}{\rm M}_{\odot}$). As a consequence, the
low-mass regime is under-sampled by a factor of $\approx$ 4-5 for K01 and
$\lessapprox$1.2 for L98. We take into account the incomplete sampling of the
initial conditions by performing an a posteriori over-sampling of the
simulated binaries with $M_{\rm ZAMS,2}\leq 10~{}{\rm M}_{\odot}$.

This ensures a good sampling of the high-mass regime and reduces the
stochastic fluctuations (e.g., Iorio et al., 2022). Table 1 lists the
properties of our input catalogues. Hereafter, we assume the models LOG1 and
KRO1 as our fiducial case for Pop. III and Pop. II stars, respectively.

Figure 3 shows the distributions of IMF, $q$, initial orbital period $P$, and
initial orbital eccentricity $e$ for the models with a flat-in-log IMF. Figure
4 shows the distribution of $M_{\rm ZAMS,1}$ and $q$ for models LOG1, KRO1,
LAR1, and TOP1. These models have the same initial distributions for mass
ratios, orbital periods and eccentricity. They differ only by the ZAMS mass
function of the primary star (Table 1).

##  3 Results

###  3.1 Black holes from single star evolution

Figure 5: The upper and lower panel show the mass of the compact remnant, and
the mass of the He core at the onset of core collapse as a function of the
initial mass, $M_{\mathrm{ZAMS{}{}}}$. Black solid and red dashed lines refer
to Pop. III and Pop. II stars, respectively.

Figure 5 shows the mass of the compact remnant ($M_{\mathrm{rem{}{}}}$) as a
function of the ZAMS mass ($M_{\mathrm{ZAMS{}{}}}$) for Pop. II and Pop. III
stars evolved via single stellar evolution. Pop. II and III stars evolving via
single star evolution produce similar He core masses and, thus, similar
compact remnant masses. The only differences are (i) in the range between
pulsational pair instability and pair instability ($M_{\rm
ZAMS}\in[60,\,{}170]$ M⊙), where envelope overshooting can cause dredge-up
episodes, and (ii) at extremely high BH masses ($M_{\rm rem}>400$ M⊙), where
Pop. II stars suffer from slightly higher mass loss rates.

In the region between pulsational pair instability and pair instability
($M_{\rm ZAMS}\in[60,\,{}170]$ M⊙), the He-core mass does not grow
monotonically, especially in the case of Pop. II stars. The core decrease in
some mass ranges is caused by dredge-up episodes triggered by envelope
undershooting (Costa et al., 2021). Different choices for the convection
parameters, such as the core overshooting ($\lambda{}_{\rm ov}=0.5$ in our
models), can change the behaviour and the occurrence of dredge-up episodes.
For instance, Pop. II stars (Z = 0.0001) with $\lambda{}_{\rm ov}=0.4$ show a
monotonic trend of the He core mass (see discussion in Iorio et al., 2022). In
the models presented in this work, for ZAMS mass $M_{\rm ZAMS}\in{}[145-160]$
M⊙ we expect an "island" of massive BH formation for Pop. II stars inside a
region of pair instability. This happens because a dredge-up episode reduces
the mass of the He and CO core below the threshold for pair instability in
this range for Pop. II stars, but not for Pop. III stars.

The maximum mass of a Pop. III BH below the pair-instability mass gap is 86 M⊙
for the adopted pair-instability model. Similarly, the maximum mass of a Pop.
II BH below the mass gap is 91 M⊙. In both cases, this mass is reached for a
ZAMS mass $\approx{100-105}$ M⊙. Below the mass gap, our models predict
several sharp features in the BH mass spectrum because of dredge-up episodes
that affect the He core mass in this range. The mass spectrum in this region
is maximally sensitive to several details of the input physics that are highly
uncertain (e.g., core overshooting, nuclear reaction rates), as already
discussed in previous papers (e.g., Leung et al., 2019; Farmer et al., 2019,
2020; Mapelli et al., 2020; Costa et al., 2021; Woosley & Heger, 2021; Vink et
al., 2021).

In our models, the upper edge of the mass gap is at $M_{\rm ZAMS}\approx{242}$
M⊙ and $\approx{236}$ M⊙ for Pop. III and II stars, respectively. Above the
mass gap, both Pop. III and II stars produce intermediate-mass BHs from direct
collapse. The mass of a BH born from a Pop. III star in this regime is very
similar to that of a BH formed by a Pop. II star with the same ZAMS mass,
because stellar winds are already extremely quenched at $Z=10^{-4}$.

The maximum BH mass in our models is $M_{\rm rem}\approx 545$ M⊙
($\approx{510}$ M⊙) for Pop. III (Pop. II) stars, corresponding to a ZAMS mass
$M_{\rm ZAMS}=550$ M⊙. We obtain these masses with the optimistic assumption
that the residual H-rich envelope of the progenitor star collapses to a BH
entirely when the star collapses. A fraction of the H-rich envelope mass might
be lost even in the case of a failed explosion, because of shocks induced by
the emission of neutrinos (e.g., Fernández et al., 2018; Renzo et al., 2020).

###  3.2 Binary evolution

Figures 6 and 7 show the secondary BH mass (i.e. the mass of the least massive
BH) versus the primary BH mass (i.e. the mass of the most massive BH) for all
our simulated BBH mergers. The masses of BBHs born from Pop. III stars are
qualitatively similar to those of BBHs born from Pop. II stars, for all the
considered models.

Mergers with at least one component above the pair-instability mass gap are
not common, and mergers inside the gap are even rarer. In our binary
simulations, it is even difficult to identify sharp edges for the pair
instability mass gap, because of dredge-up episodes and mass transfer (Figure
8). Assuming that the pair-instability mass gap spans from 85 to 230 M⊙, we
find that BBH mergers with primary BH masses above the gap are up to $3.3$%
(LOG3) and up to $0.09$% (TOP5) for Pop. III and Pop. II stars, respectively.
With the same definition, mergers with primary BH mass inside the gap are up
to $1.9$% (LOG3) and up to $2.4$% (TOP5) for Pop. III and Pop. II stars,
respectively. In general, Pop. II stars seem to produce more BBH mergers with
primary BH mass above the gap with respect to Pop. III stars, with the
exception of model LOG3 (Fig. 8). Furthermore, no secondary BHs in our BBH
mergers have mass inside or above the gap.

The most common primary BH masses are around 30–40 M⊙ (Fig. 8). There is a
dearth of low-mass primary BHs ($8-10$ M⊙) with respect to LVK mergers (Abbott
et al., 2019, 2021d, 2021c) in all of our runs, even KRO1. This is a
consequence of the negligible mass loss rate and relatively compact stellar
radii at such low metallicity.

Figure 6: Distribution of secondary ($M_{\mathrm{s}}$) versus primary
($M_{\mathrm{p}}$) mass of all BBH mergers in our simulations LOG1–5. Left
(Right): Pop. III (II) stars. The colour bar indicates the number of BBHs in
each cell. Figure 7:  Same as Fig. 6 but for models KRO1, KRO5, LAR1, LAR5,
TOP1, and TOP5.

Our models show a preference for equal mass systems but also a non-negligible
contribution from unequal mass mergers (Figs. 6 and 7). The secondary BH mass
is always $M_{\rm s}\leq{}45$ M⊙ for both Pop. II and III stars. Figure 9
highlights some differences between Pop. II and III BBHs. For example, the
most common secondary BH mass for Pop. III stars is $\sim{20}$ M⊙, while for
Pop. II stars it is either $\sim{8-10}$ M⊙ (KRO1, LAR1), or $\sim{35-38}$ M⊙
(LOG2, LOG5, TOP5), depending on the model.

We find another interesting difference between Pop. III and II BBHs if we look
at the mass ratio (Fig. 10). In the case of Pop. II stars, equal-mass BBHs are
the most common systems regardless of the model, even if models LOG5, KRO1,
LAR5, and TOP5 show a mild secondary peak for $q\sim{0.4-0.5}$. In contrast,
for Pop. III stars, the most common BBH mass ratio is $\sim{0.8-0.9}$ for the
models LOG1, LOG3, LOG4, KRO1, LAR1, and TOP1, i.e. for all the models
adopting the S12 initial period distribution. This is a consequence of the
dominant evolutionary channels in such models (see Section 4.1).

Finally, the distribution of delay times $t_{\rm del}$ (i.e., the time elapsed
between the formation of the binary system and the BBH merger) shows another
difference between Pop. III and II BBHs (Fig. 11). All Pop. III models seem to
match the trend $t_{\rm del}\propto{}t^{-1}$ between 3 and $10^{4}$ Myr. In
contrast, some of the Pop. II models (LOG2, LOG5, KRO5, LAR5, and TOP5) show
an excess of short delay times ($3-10$ Myr). The models showing this excess
share the SB13 orbital period distribution. This feature is another signature
of the formation channel, as we discuss in Section 4.1.

Figure 8:  Distribution of the primary BH mass (i.e., the most massive BH in
each binary system) in the simulated BBH mergers. Upper (lower) panel: Pop.
III (Pop. II) stars.  Figure 9: Same as Fig. 8, but for the secondary BH mass
(i.e., the least massive BH of each binary system). Figure 10: Same as Fig. 8,
but for the mass ratio $q_{\rm BH}=M_{\rm s}/M_{\rm p}$ between the secondary
and primary BH. Figure 11:  Distribution of delay times, i.e. the time elapsed
from the formation of the binary star to the merger of the two BHs. The gray
dashed line shows the $\propto{}t^{-1}$ predicted trend (e.g., Dominik et al.,
2012).

##  4 Discussion

###  4.1 Formation channels of BBH mergers

The features of BBH mergers we described in the previous section (Section 3.2)
can be interpreted by looking at the formation channels of our BBHs. Figure 12
summarizes the main formation channels of BBH mergers from Pop. III and Pop.
II stars. Tables 2 and 3 report the percentages in detail.

We distinguish five main channels, following the definition by Iorio et al.
(2022). BBH mergers that go through Channel 0 do not undergo any mass transfer
episodes during the evolution of their progenitors. Systems belonging to this
channel are always very rare ($\ll{1}\%$). Since this channel is so uncommon,
we do not show it in Figure 12 and in the following Figures.

Channel I is often referred to as the "traditional" formation channel of BBH
mergers: the two progenitor stars undergo stable mass transfer before the
formation of the first BH. Then, after the formation of the first BH, the
system evolves through at least one common envelope.

In channel II, the system evolves only via stable mass transfer episodes.
Finally, in both channel III and IV, the system undergoes at least one common
envelope before the formation of the first BH. The only difference between
channel III and IV is that in the former the companion star preserves a
residual of the original H-rich envelope at the time of the formation of the
first BH, while in the latter the companion has already been stripped of its
envelope when the first BH forms.

Figure 12:  Percentage of BBH mergers that evolve via each of the four
channels considered in this work. The x-axis refers to the simulation set,
while the y-axis specifies the channel. The upper (lower) panel refers to BBHs
that form from Pop. III (II) stars.  Table 2: Percentage of BBH mergers from
Pop. III stars.

Model | BBHm | Ch. 0 | Ch. I | Ch. II | Ch. III | Ch. IV
---|---|---|---|---|---|---
| [%] | [%] | [%] | [%] | [%] | [%]
LOG1 | 11.25 | 0.1 | 3.65 | 74.81 | 7.62 | 13.71
LOG2 | 0.75 | 0.22 | 3.3 | 1.27 | 58.77 | 36.11
LOG3 | 9.33 | 0.09 | 3.35 | 70.41 | 17.31 | 8.76
LOG4 | 11.57 | 0.12 | 3.26 | 67.57 | 18.13 | 10.69
LOG5 | 0.68 | 0.23 | 3.3 | 1.52 | 64.6 | 30.07
KRO1 | 14.66 | 0.19 | 1.93 | 85.82 | 4.08 | 7.7
KRO5 | 0.85 | 0.34 | 2.04 | 2.18 | 65.11 | 29.35
LAR1 | 14.34 | 0.16 | 2.13 | 83.57 | 4.89 | 9.03
LAR5 | 0.91 | 0.38 | 2.26 | 1.98 | 64.97 | 29.74
TOP1 | 6.47 | 0.06 | 5.65 | 66.17 | 10.39 | 17.67
TOP5 | 0.36 | 0.19 | 4.47 | 1.26 | 64.81 | 29.14

Column 1: Model; column2: percentage of BBH mergers with respect to all
simulated BBHs; columns 3, 4, 5, 6, and 7: BBH mergers formed via channel 0,
I, II, III, and IV, respectively.

Table 3: Percentage of BBH mergers from Pop. II stars.

Model | BBHm | Ch. 0 | Ch. I | Ch. II | Ch. III | Ch. IV
---|---|---|---|---|---|---
| [%] | [%] | [%] | [%] | [%] | [%]
LOG1 | 13.53 | 0.02 | 5.94 | 70.48 | 7.59 | 15.44
LOG2 | 0.97 | 0.12 | 1.3 | 1.44 | 48.11 | 45.88
LOG3 | 10.88 | 0.03 | 5.98 | 70.55 | 12.34 | 10.71
LOG4 | 14.46 | 0.02 | 4.03 | 67.06 | 13.37 | 14.79
LOG5 | 0.86 | 0.11 | 1.48 | 1.53 | 54.96 | 39.05
KRO1 | 16.23 | 0.03 | 6.16 | 79.64 | 2.85 | 10.01
KRO5 | 1.15 | 0.23 | 0.79 | 1.66 | 41.21 | 45.93
LAR1 | 16.15 | 0.04 | 5.28 | 78.93 | 3.37 | 11.35
LAR5 | 1.16 | 0.22 | 0.91 | 1.6 | 43.11 | 46.57
TOP1 | 8.52 | 0.02 | 7.42 | 59.61 | 14.73 | 17.99
TOP5 | 0.48 | 0.1 | 1.9 | 1.44 | 64.8 | 30.76

Column 1: Model; column2: percentage of BBH mergers with respect to all
simulated BBHs; columns 3, 4, 5, 6, and 7: BBH mergers formed via channel 0,
I, II, III, and IV, respectively.

Figure 13:  Main properties of BBH mergers and their progenitors in model LOG1
for Pop. III stars. From top to bottom and from left to right: ZAMS mass of
the progenitor of the primary BH $M_{\rm ZAMS}(M_{\rm p})$, mass ratio of the
progenitors $q_{\rm ZAMS}=M_{\rm ZAMS}(M_{\rm s})/M_{\rm ZAMS}(M_{\rm p})$
(i.e. the ratio between the ZAMS mass of the progenitor of the secondary BH
and the ZAMS mass of the progenitor of the primary BH), mass of the primary BH
($M_{\rm p}$), mass ratio of the two BHs ($q_{\rm BH}=M_{\rm s}/M_{\rm p}$),
initial semi-major axis ($a$), initial eccentricity ($e$). Light-blue line:
channel I, blue line: channel II, red line: channel III, black line: channel
IV. Figure 14:  Same as Fig. 13 but for model LOG3 for Pop. III stars.  Figure
15:  Same as Figure 13 but for model LOG5 for Pop. III stars.  Figure 16:
Same as Figure 13 but here we show the KRO1 model and compare BBHs from Pop.
III (solid lines) and Pop. II stars (dashed lines). We show only channel II
and IV because they are the two most important channels (especially channel
II, Fig. 12) and to make the plot more readable.  Figure 17:  Delay time
distribution for Pop. II and Pop. III BBHs (dashed thick and solid thin lines,
respectively). From left to right and from top to bottom: models LOG1, KRO1,
LOG5 and KRO5. In each panel, we show only the two most important channels for
each considered model: channels II and IV in the upper panels, and channels
III and IV in the lower panels.

Channel II (i.e. stable mass transfer) is the dominant channel for models
LOG1, LOG3, LOG4, KRO1, LAR1, and TOP1, i.e. for all the systems that adopt
the initial S12 distribution of the orbital periods. Figures 13 and 14 show
that most channel II systems are associated with short initial semi-major axis
$a=10-10^{3}$ R⊙. These short initial semi-major axes are heavily suppressed
with the orbital period distribution by SB13 (Fig. 15).

Figures 13 and 14 also show that channel II is associated with relatively low
mass ratios in the ZAMS ($q_{\rm ZAMS}\sim{0.5-0.9}$) and relatively low mass
ratios between the two final BHs ($q_{\rm BH}=0.75-0.9$) for Pop. III stars.
Hence, the predominance of channel II in models LOG1, LOG3, LOG4, KRO1, LAR1,
and TOP1 (adopting the S12 orbital period distribution) explains why these
models have a preference for BBH mass ratios $q_{\rm BH}=0.75-0.9$ in the case
of Pop. III stars, as discussed in the previous Section (Fig. 10). In these
systems, the mass difference between primary and secondary stars is
sufficiently large that the system undergoes the first mass transfer while the
secondary is still on the MS.

In contrast, channels III and IV are the dominant channels for models LOG2,
LOG5, KRO5, LAR5, and TOP5, i.e. models that adopt the SB13 orbital period
distribution. As shown by, e.g., Fig. 15, the large initial semi-major axes of
distribution SB13 suppress systems with initial orbital separation $a<10^{3}$
R⊙, hence suppressing channel II. Channel III and IV preferentially arise when
$a\sim{10^{3}-10^{5}}$ R⊙. In this case, mass transfer takes place only when
the radii of the two stars become very large, i.e. in the late evolutionary
stages. Channel IV is the preferred channel of equal-mass stars, that evolve
nearly at the same time and strip off each other’s envelopes. It mainly leads
to the formation of equal-mass BBHs, explaining the preference of these models
for equal-mass mergers (Fig. 10). In contrast, channel III has a preference
for markedly unequal-mass systems, explaining the population of BBHs with
$q\leq{}0.6$ in models LOG5, KRO5, LAR5, and TOP5 (especially for Pop. II
stars, Fig. 10).

Figure 16 compares the properties of BBH mergers from Pop. III and Pop. II
stars in the case of model KRO1 (the fiducial model for Pop. II stars). Since
channels I and III are less important than channels II and IV in this model,
we show only the latter channels for simplicity. This Figure shows that Pop.
III and Pop. II stars have a very similar behaviour in the case of channel IV.
As for channel II, we see three main differences: Pop. II stars have a
preference for i) higher ZAMS mass $M_{\rm ZAMS}(M_{\rm p})$, ii) larger mass
ratios in the ZAMS $q_{\rm ZAMS}$, and iii) and larger BBH mass ratios $q_{\rm
BH}$ with respect to Pop. III stars. The correlation between these three
properties explains why Pop. II stars tend to produce equal-mass BBHs, while
Pop. III stars produce BBHs with a mass ratio peaking at $q_{\rm
BH}\approx{0.8}$ (Fig. 10).

Figure 17 compares the delay time distribution of Pop. III and Pop. II binary
systems if we consider models KRO1 (fiducial model for Pop. II), LOG1
(fiducial model for Pop. III), LOG5 and KRO5. Channel IV is skewed toward the
shortest delay times, for both Pop. III and II binary stars, because it is
associated with the most efficient orbital shrinking during common envelope.
The shrinking is more efficient for Pop. II stars (especially in models LOG5
and KRO5), because they reach even larger radii in their late evolutionary
stages (see Figures 1 and 2). This explains why the overall delay time
distribution of Pop. II systems (Figure 11) has an excess at very low values
($t_{\rm del}\leq{10}$ Myr) in models LOG2, LOG5, KRO5, LAR5, and TOP5. In
contrast, channel II is associated with long delay times, because stable mass
transfer is not as efficient as common envelope in reducing the orbital
separation.

Tables 2 and 3, and Figure 12 show that both Pop. III and Pop. II stars behave
in a very different way from more metal rich binary systems. In fact, only
$\lesssim{5.7}\%$ ($\lesssim{7.5}\%$) of all BBH mergers evolve via channel I
in the case of Pop. III (II) binary stars. For comparison, Iorio et al. (2022)
show that between 50 and 80% of all BBH mergers evolve via channel I at
metallicity between $Z=2\times{}10^{-3}$ and $Z=10^{-2}$ (see Fig. 14 of Iorio
et al. 2022 for $\alpha=1$).

The abundance of channel II systems and the dearth of channel I systems for
Pop. III and II binary stars with respect to metal-rich binary systems
($Z\in{}[2\times{}10^{-3}-10^{-2}]$) are a consequence of the large BH masses
at such low metallicity. In fact, both channel I and II systems go through a
stable mass transfer before the formation of the first compact object, and
then undergo a second mass transfer after the formation of the first compact
object. The only difference between the two channels is that the mass transfer
episode between the companion star and the first-born BH becomes unstable in
channel I and remains stable in channel II. In sevn (as in most binary
population-synthesis codes) the stability of mass transfer is evaluated
through a critical mass ratio $q_{\rm crit}$ between the donor and the
accretor: the systems with mass ratio $q\geq{}q_{\rm crit}$ undergo a
dynamically unstable mass transfer (i.e., a common-envelope episode), while
mass transfer remains stable in the other cases (Hurley et al., 2002). Since
most our Pop. II and III BHs in tight binary systems are relatively massive
($\gtrsim{}20$ M⊙, Fig. 5), we have that $q<q_{\rm crit}$ in most binary
systems (where $q$ is the mass ratio between the donor star and the first born
BH), ensuring the stability of most mass-transfer episodes.

Furthermore, the abundance of channel II versus channel I systems depends on
the assumed value of common-envelope efficiency $\alpha_{\rm CE}$. We have re-
run the fiducial case LOG1 with $\alpha_{\rm CE}=3$. We find that for Pop. III
stars, the percentage of channel II systems decreases from 75% for
$\alpha_{\rm CE}=1$ down to 51% for $\alpha_{\rm CE}=3$. This happens because,
when $\alpha_{\rm CE}$ is large, the common envelope process is less efficient
in shrinking the orbital separation. In contrast, the relative abundance
between channel I and II systems is not significantly affected by our
assumption that mass transfer is always stable for MS and Hertzsprung-gap
donors. In fact, both channels I and II undergo a stable mass transfer when
the primary star is still a MS or an Hertzsprung-gap star. Relaxing the
aforementioned assumption has a more sizable impact on the evolution of
channels III and IV systems, which evolve via common envelope before the
formation of the first-born BH.

###  4.2 BBH mergers with primary above the mass gap

Channel III is the key to interpret the abundance of BBHs with primary mass
above the gap in model LOG3. These systems are associated with a population of
binary stars with very low $q_{\rm ZAMS}$ (Figure 14), mostly following the
relation $q\approx 23\,{}M_{\odot}/M_{\mathrm{ZAMS,1}}$ with
$M_{\mathrm{ZAMS,1}}\gtrsim 250\ \mbox{M${}_{\odot}$}$ and large initial
separations $a\gtrsim 10^{3}\ \mbox{M${}_{\odot}$}$. Such systems evolve
through channel III triggering a Roche-lobe overflow episode that becomes
unstable when the primary star enters the red super-giant phase. Systems with
smaller initial separation merge due to a double Roche-lobe overflow. Only
systems with initial $q_{\rm ZAMS}<0.15$ can evolve through this channel. This
is the reason why most of the high-mass primaries come from the model LOG3,
followed by LOG4 and LOG1, LOG2 (Fig. 3). The fact we do not see this feature
in the model LOG3 of Pop. II stars is related to the evolution of radius.
High-mass Pop II stars reach the radius to start the interaction earlier than
Pop. III stars (see Fig. 2).

###  4.3 Comparison with previous work

Several authors have explored the formation of BBHs from Pop. III stars (e.g.,
Kinugawa et al., 2014, 2016, 2020; Hartwig et al., 2016; Belczynski et al.,
2017; Liu & Bromm, 2020b; Tanikawa et al., 2021b, 2022a, 2022b; Wang et al.,
2022). The ZAMS mass range we consider in this work is comparable to the one
explored by Tanikawa et al. (2021b). In our models, we do not have any mergers
with both BHs above the mass gap, while these are very common in their
fiducial model. This discrepancy springs from the intrinsic differences
between our single star evolution models. In fact, the very massive stars
($M_{\rm ZAMS}>200$ M⊙) considered by Tanikawa et al. (2021b) end the MS as
compact blue super-giant stars, while our very massive stars expand during the
MS and become red super-giant stars already at the end of the MS. As a
consequence, the very massive binary systems by Tanikawa et al. (2021b)
undergo stable mass transfer and leave BHs above the mass gap, while our very
massive binary systems start an unstable common envelope phase as soon as they
leave the MS and merge prematurely, before becoming BHs.

The same line of reasoning explains why our delay times (Figs. 11 and Figs.
17) are generally much shorter than the one presented, e.g., in Figure 3 from
Tanikawa et al. (2021a). Almost all BBH mergers from Pop. III stars evolve via
stable mass transfer in the models by Tanikawa et al. (2021a) and thus have
long delay times, while our channel IV mergers (which do not form in Tanikawa
et al. 2021a) have very short delay times. These results confirm the key role
of single star evolution (including uncertainties about core overshooting,
convection and rotation) for the formation of merging BBHs.

##  5 Summary and conclusions

We have presented a new set of Pop. III stars ($Z=10^{-11}$) obtained with the
stellar evolution code parsec (Bressan et al., 2012; Costa et al., 2019; Costa
et al., 2021). Our Pop. III stars range from 2 to 600 M⊙. With respect to Pop.
II stars ($Z=10^{-4}$), Pop. III stars with initial mass $M_{\rm
ZAMS}\in{[14,40]}$ M⊙ evolve with much more compact radii ($R\ll{}10^{2}$ R⊙).
Furthermore, the most massive Pop. III stars
($\mbox{$M_{\mathrm{ZAMS{}{}}}$}>{100}$ M⊙) end their life as blue super-giant
stars, whereas Pop. II stars in the same mass range die as yellow and red
super-giant stars.

We use these tracks as input tables for our fast binary population synthesis
code sevn (Spera et al., 2019; Mapelli et al., 2020; Iorio et al., 2022), in
order to study the population of BHs and BBHs born from Pop. II and Pop. III
stars. We explore a variety of initial conditions for our binary stars,
including a flat-in-log, a Kroupa, a Larson, and a top-heavy IMF.

We estimate similar BH masses from the evolution of single massive Pop. II and
Pop. III stars. In our fiducial model, the maximum BH mass below the pair-
instability gap is 91 and 86 M⊙ for Pop. II and III stars, respectively. Above
the gap, both Pop. II and Pop. III stars produce BHs more massive than
$\approx{230}$ M⊙ if they can achieve a ZAMS mass of $\gtrsim{240}$ M⊙ (Figure
5). Assumptions on core overshooting, envelope undershooting, and stellar
rotation can significantly affect this result, because they influence the mass
of the He and CO cores, hence the central temperature and density.
Furthermore, these results spring from the assumption that the residual H-rich
envelope is not ejected during the failed supernova (e.g., Costa et al.,
2022).

Most BBH mergers from both Pop. II and Pop. III stars have primary BH mass
below the mass gap. In order to populate the region above the gap, we need
very compact stellar radii, that can be achieved either with fast rotation
(chemically homogeneous evolution, de Mink & Belczynski 2015) or by
suppressing core overshooting (Tanikawa et al., 2022b). With our evolutionary
models, we find no mergers with secondary BH mass above the gap. We expect
that only dynamics of dense stellar systems can pair up BHs with both primary
and secondary mass above the gap, and populate the gap as well (e.g., Mapelli
et al., 2022; Wang et al., 2022).

The mass ratios are one of the main signatures of Pop. III versus Pop. II
BBHs. In most of our models, Pop. II BBHs are predominantly equal-mass
systems, whereas Pop. III BBHs have a peak at mass ratio $q_{\rm BH}=0.8-0.9$.
This difference is too subtle for current detectors, even at a population
level, but it can be investigated with next-generation ground-based detectors.

A distinctive signature of Pop. III and II BBHs with respect to BBHs born from
metal-rich stars are the evolutionary channels. Assuming the orbital period
distribution from Sana et al. (2012), the vast majority ($60-80$%) of Pop. III
and II progenitor stars of BBH mergers evolve via channel II, i.e. just stable
mass transfer, with no common envelope. In contrast, at higher metallicity
($Z\sim{2\times{}10^{-3}-10^{-2}}$) and with the same set-up for binary
evolution, the dominant evolutionary pathway ($50-80$% BBH mergers) becomes
channel I, characterized by common envelope between the first-born BH and its
companion star (Iorio et al., 2022).

Overall, our models show that Pop. III and II stars produce a similar BBH
population, especially if we adopt the same IMF and initial orbital
properties. The actual IMF and maximum mass of metal-poor and metal-free stars
are two of the main sources of uncertainty. Since the peak of the star
formation rate density of Pop. II stars is about one order of magnitude higher
than that of Pop. III stars and occurs at a lower redshift ($z\sim{4-5}$
versus $z\sim{12-16}$, Hartwig et al. 2022, e.g.,), we expect that the
contribution of Pop. II binary stars to the population of massive BBHs (with
total mass $\gtrsim{50}$ M⊙) detectable by LVK is dominating with respect to
the contribution of Pop. III binary stars.

## Acknowledgements

We are grateful to Simon Glover, Tilman Hartwig, Tomoya Kinugawa, Mario Spera,
and Ataru Tanikawa for their enlightening comments. GC, GI, MM, and FS
acknowledge financial support from the European Research Council for the ERC
Consolidator grant DEMOBLACK, under contract no. 770017\. This research made
use of NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), IPython
(Perez & Granger, 2007). For the plots we used Matplotlib (Hunter, 2007).

## Data Availability

The data underlying this article and the sevn configurations files are
available on Zenodo at https://doi.org/10.5281/zenodo.7736310 (Costa et al.,
2023). sevn is publicly available at https://gitlab.com/sevncodes/sevn.git:
the version used in this work is the commit 0f9ae3bf in the branch
Costa23popIII (https://gitlab.com/sevncodes/sevn/-/tree/Costa23popIII).
Further data will be shared on reasonable request to the corresponding
authors.

## References

  * Abbott et al. (2019) Abbott B. P., et al., 2019, ApJ, 882, L24
  * Abbott et al. (2020a) Abbott R., et al., 2020a, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2020b) Abbott R., et al., 2020b, ApJ, 900, L13
  * Abbott et al. (2021a) Abbott R., et al., 2021a, arXiv e-prints, p. arXiv:2108.01045
  * Abbott et al. (2021b) Abbott R., et al., 2021b, arXiv e-prints, p. arXiv:2111.03606
  * Abbott et al. (2021c) Abbott R., et al., 2021c, arXiv e-prints, p. arXiv:2111.03634
  * Abbott et al. (2021d) Abbott R., et al., 2021d, ApJ, 913, L7
  * Abel et al. (2002) Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93
  * Agrawal et al. (2023) Agrawal P., Hurley J., Stevenson S., Rodriguez C. L., Szecsi D., Kemp A., 2023, arXiv e-prints, p. arXiv:2303.10187
  * Barkat et al. (1967) Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379
  * Belczynski et al. (2017) Belczynski K., Ryu T., Perna R., Berti E., Tanaka T. L., Bulik T., 2017, MNRAS, 471, 4702
  * Böhm-Vitense (1958) Böhm-Vitense E., 1958, Z. Astrophys., 46, 108
  * Bray & Eldridge (2016) Bray J. C., Eldridge J. J., 2016, MNRAS, 461, 3747
  * Bray & Eldridge (2018) Bray J. C., Eldridge J. J., 2018, MNRAS, 480, 5657
  * Bressan et al. (1981) Bressan A. G., Chiosi C., Bertelli G., 1981, A&A, 102, 25
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Bromm (2013) Bromm V., 2013, Reports on Progress in Physics, 76, 112901
  * Bromm & Larson (2004) Bromm V., Larson R. B., 2004, ARA&A, 42, 79
  * Caffau et al. (2011) Caffau E., Ludwig H. G., Steffen M., Freytag B., Bonifacio P., 2011, Sol. Phys., 268, 255
  * Cassisi & Castellani (1993) Cassisi S., Castellani V., 1993, ApJS, 88, 509
  * Chen et al. (2015) Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068
  * Chiaki et al. (2018) Chiaki G., Susa H., Hirano S., 2018, MNRAS, 475, 4378
  * Chon et al. (2021) Chon S., Omukai K., Schneider R., 2021, MNRAS, 508, 4175
  * 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
  * Costa et al. (2019) Costa G., Girardi L., Bressan A., Marigo P., Rodrigues T. S., Chen Y., Lanza A., Goudfrooij P., 2019, MNRAS, 485, 4641
  * Costa et al. (2021) Costa G., Bressan A., Mapelli M., Marigo P., Iorio G., Spera M., 2021, MNRAS, 501, 4514
  * Costa et al. (2022) Costa G., Ballone A., Mapelli M., Bressan A., 2022, MNRAS, 516, 1072
  * Costa et al. (2023) Costa G., Mapelli M., Iorio G., Santoliquido F., Escobar G. J., Bressan A., 2023, Data from "Massive binary black holes from Population II and III stars", doi:10.5281/zenodo.7736310, https://doi.org/10.5281/zenodo.7736310
  * Dominik et al. (2012) Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2012, ApJ, 759, 52
  * Farmer et al. (2019) Farmer R., Renzo M., de Mink S. E., Marchant P., Justham S., 2019, ApJ, 887, 53
  * Farmer et al. (2020) Farmer R., Renzo M., de Mink S. E., Fishbach M., Justham S., 2020, ApJ, 902, L36
  * Fernández et al. (2018) Fernández R., Quataert E., Kashiyama K., Coughlin E. R., 2018, MNRAS, 476, 2366
  * Fowler & Hoyle (1964) Fowler W. A., Hoyle F., 1964, ApJS, 9, 201
  * Frebel & Norris (2015) Frebel A., Norris J. E., 2015, ARA&A, 53, 631
  * Frebel et al. (2007) Frebel A., Johnson J. L., Bromm V., 2007, MNRAS, 380, L40
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Giacobbo & Mapelli (2019) Giacobbo N., Mapelli M., 2019, MNRAS, 482, 2234
  * Giacobbo & Mapelli (2020) Giacobbo N., Mapelli M., 2020, ApJ, 891, 141
  * Glover (2013) Glover S., 2013, in Wiklind T., Mobasher B., Bromm V., eds, Astrophysics and Space Science Library Vol. 396, The First Galaxies. p. 103 (arXiv:1209.2509), doi:10.1007/978-3-642-32362-1_3
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hartwig et al. (2016) Hartwig T., Volonteri M., Bromm V., Klessen R. S., Barausse E., Magg M., Stacy A., 2016, MNRAS, 460, L74
  * Hartwig et al. (2022) Hartwig T., et al., 2022, ApJ, 936, 45
  * Herwig et al. (1997) Herwig F., Bloecker T., Schoenberner D., El Eid M., 1997, A&A, 324, L81
  * Hirano et al. (2014) Hirano S., Hosokawa T., Yoshida N., Umeda H., Omukai K., Chiaki G., Yorke H. W., 2014, ApJ, 781, 60
  * Hirano et al. (2015) Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W., 2015, MNRAS, 448, 568
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Iorio et al. (2022) Iorio G., et al., 2022, arXiv e-prints, p. arXiv:2211.11774
  * Jaacks et al. (2019) Jaacks J., Finkelstein S. L., Bromm V., 2019, MNRAS, 488, 2202
  * Jaura et al. (2022) Jaura O., Glover S. C. O., Wollenberg K. M. J., Klessen R. S., Geen S., Haemmerlé L., 2022, MNRAS, 512, 116
  * Kinugawa et al. (2014) Kinugawa T., Inayoshi K., Hotokezaka K., Nakauchi D., Nakamura T., 2014, MNRAS, 442, 2963
  * Kinugawa et al. (2016) Kinugawa T., Miyamoto A., Kanda N., Nakamura T., 2016, MNRAS, 456, 1093
  * Kinugawa et al. (2020) Kinugawa T., Nakamura T., Nakano H., 2020, MNRAS, 498, 3946
  * Klessen & Glover (2023) Klessen R. S., Glover S. C. O., 2023, arXiv e-prints, p. arXiv:2303.12500
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kruckow et al. (2018) Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G., 2018, MNRAS, 481, 1908
  * Larson (1998) Larson R. B., 1998, MNRAS, 301, 569
  * Leung et al. (2019) Leung S.-C., Nomoto K., Blinnikov S., 2019, ApJ, 887, 72
  * Liu & Bromm (2020a) Liu B., Bromm V., 2020a, MNRAS, 495, 2475
  * Liu & Bromm (2020b) Liu B., Bromm V., 2020b, MNRAS, 495, 2475
  * Mapelli et al. (2020) Mapelli M., Spera M., Montanari E., Limongi M., Chieffi A., Giacobbo N., Bressan A., Bouffanais Y., 2020, ApJ, 888, 76
  * Mapelli et al. (2022) Mapelli M., Bouffanais Y., Santoliquido F., Arca Sedda M., Artale M. C., 2022, MNRAS, 511, 5797
  * Marchant et al. (2019) Marchant P., Renzo M., Farmer R., Pappas K. M. W., Taam R. E., de Mink S. E., Kalogera V., 2019, ApJ, 882, 36
  * Marchant et al. (2021) Marchant P., Pappas K. M. W., Gallegos-Garcia M., Berry C. P. L., Taam R. E., Kalogera V., Podsiadlowski P., 2021, A&A, 650, A107
  * Marigo et al. (2001) Marigo P., Girardi L., Chiosi C., Wood P. R., 2001, A&A, 371, 152
  * Murphy et al. (2021) Murphy L. J., et al., 2021, MNRAS, 501, 2745
  * Nguyen et al. (2022) Nguyen C. T., et al., 2022, A&A, 665, A126
  * Perez & Granger (2007) Perez F., Granger B. E., 2007, Computing in Science Engineering, 9, 21
  * Prole et al. (2022) Prole L. R., Clark P. C., Klessen R. S., Glover S. C. O., 2022, MNRAS, 510, 4019
  * Rakavy & Shaviv (1967) Rakavy G., Shaviv G., 1967, ApJ, 148, 803
  * Renzo et al. (2020) Renzo M., Farmer R., Justham S., Götberg Y., de Mink S. E., Zapartas E., Marchant P., Smith N., 2020, A&A, 640, A56
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Santoliquido et al. (2023) Santoliquido F., Mapelli M., Iorio G., Costa G., Glover S. C. O., Hartwig T., Klessen R. S., Merli L., 2023, arXiv e-prints, p. arXiv:2303.15515
  * Schneider et al. (2006) Schneider R., Omukai K., Inoue A. K., Ferrara A., 2006, MNRAS, 369, 1437
  * Schneider et al. (2012) Schneider R., Omukai K., Bianchi S., Valiante R., 2012, MNRAS, 419, 1566
  * Schwarzschild (1958) Schwarzschild M., 1958, Structure and evolution of the stars.. Princeton, Princeton University Press, 1958.
  * Sharda & Krumholz (2022) Sharda P., Krumholz M. R., 2022, MNRAS, 509, 1959
  * Smith et al. (2015) Smith B. D., Wise J. H., O’Shea B. W., Norman M. L., Khochfar S., 2015, MNRAS, 452, 2822
  * Spera & Mapelli (2017) Spera M., Mapelli M., 2017, MNRAS, 470, 4739
  * Spera et al. (2019) Spera M., Mapelli M., Giacobbo N., Trani A. A., Bressan A., Costa G., 2019, MNRAS, 485, 889
  * Stacy & Bromm (2013) Stacy A., Bromm V., 2013, MNRAS, 433, 1094
  * Stevenson et al. (2019) Stevenson S., Sampson M., Powell J., Vigna-Gómez A., Neijssel C. J., Szécsi D., Mandel I., 2019, ApJ, 882, 121
  * Susa et al. (2014) Susa H., Hasegawa K., Tominaga N., 2014, ApJ, 792, 32
  * Tanikawa et al. (2021a) Tanikawa A., Kinugawa T., Yoshida T., Hijikawa K., Umeda H., 2021a, MNRAS, 505, 2170
  * Tanikawa et al. (2021b) Tanikawa A., Susa H., Yoshida T., Trani A. A., Kinugawa T., 2021b, ApJ, 910, 30
  * Tanikawa et al. (2022a) Tanikawa A., Chiaki G., Kinugawa T., Suwa Y., Tominaga N., 2022a, PASJ, 74, 521
  * Tanikawa et al. (2022b) Tanikawa A., Yoshida T., Kinugawa T., Trani A. A., Hosokawa T., Susa H., Omukai K., 2022b, ApJ, 926, 83
  * Valiante et al. (2016) Valiante R., Schneider R., Volonteri M., Omukai K., 2016, MNRAS, 457, 3356
  * Vink et al. (2021) Vink J. S., Higgins E. R., Sander A. A. C., Sabhahit G. N., 2021, MNRAS, 504, 146
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Volpato et al. (2023) Volpato G., Marigo P., Costa G., Bressan A., Trabucchi M., Girardi L., 2023, ApJ, 944, 40
  * Wang et al. (2022) Wang L., Tanikawa A., Fujii M., 2022, MNRAS, 515, 5106
  * Wollenberg et al. (2020) Wollenberg K. M. J., Glover S. C. O., Clark P. C., Klessen R. S., 2020, MNRAS, 494, 1871
  * Woosley (2017) Woosley S. E., 2017, ApJ, 836, 244
  * Woosley (2019) Woosley S. E., 2019, ApJ, 878, 49
  * Woosley & Heger (2021) Woosley S. E., Heger A., 2021, ApJ, 912, L31
  * Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Yoshida et al. (2006) Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ, 652, 6
  * de Mink & Belczynski (2015) de Mink S. E., Belczynski K., 2015, ApJ, 814, 58

  *[MS]: main sequence
  *[ZAMS]: zero-age main sequence
  *[CHeB]: core helium burning
  *[AGB]: asymptotic giant branch
  *[BSG]: blue super-giant
  *[YSG]: yellow super-giant
  *[RSG]: red super-giant
  *[RSGs]: red super-giant
  *[BSGs]: blue super-giant
  *[IMF]: initial mass function