# Black hole—black hole total merger mass and the origin of LIGO/Virgo sources

Krzysztof Belczynski11affiliationmark:  , Zoheyr Doctor22affiliationmark:  ,
Michael Zevin33affiliationmark:  44affiliationmark:  , Aleksandra
Olejak11affiliationmark:  , Sambaran Banerjee44affiliationmark:
55affiliationmark:  , Debatri Chattopadhyay66affiliationmark:  1 Nicolaus
Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18,
00-716 Warsaw, Poland<EMAIL_ADDRESS><EMAIL_ADDRESS>
2 Center for Interdisciplinary Exploration and Research in Astrophysics
(CIERA) and Department of Physics & Astronomy, Northwestern University, 1800
Sherman Ave, Evanston, IL 60201, USA<EMAIL_ADDRESS>
3 Kavli Institute for Cosmological Physics, The University of Chicago, 5640
South Ellis Avenue, Chicago, Illinois 60637, USA
4 Enrico Fermi Institute, The University of Chicago, 933 East 56th Street,
Chicago, Illinois 60637, USA<EMAIL_ADDRESS>
5 Helmholtz-Instituts für Strahlen- und Kernphysik (HISKP), Nussallee 14-16,
D-53115 Bonn, Germany<EMAIL_ADDRESS>
6 Argelander-Institut für Astronomie (AIfA), Auf dem Hügel 71, D-53121 Bonn,
Germany
7 Gravity Exploration Institute, School of Physics and Astronomy, Cardiff
University, Cardiff, CF24 3AA, UK<EMAIL_ADDRESS>
###### Abstract

The LIGO-Virgo-KAGRA (LVK) Collaboration has reported nearly 100 BH-BH
mergers. The LVK provides estimates of rates, masses, effective spins, and
redshifts for these mergers. Yet, the formation channel(s) of the mergers
remains uncertain. One way to search for a formation site is to contrast
properties of detected BH-BH mergers with different models of BH-BH merger
formation. Our study is designed to investigate the usefulness of the total
BH-BH merger mass and its evolution with redshift in establishing the origin
of gravitational-wave sources. We find that the average intrinsic BH-BH total
merger mass shows exceptionally different behavior for the models that we
adopt for our analysis. In the local universe ($z=0$) the average merger mass
changes from $\overline{M}_{\rm tot,int}\sim 25{\rm~{}M}_{\odot}$ for CE
binary evolution and open clusters formation channels, to $\overline{M}_{\rm
tot,int}\sim 30{\rm~{}M}_{\odot}$ for the stable-RLOF binary channel, to
$\overline{M}_{\rm tot,int}\sim 45{\rm~{}M}_{\odot}$ for the globular cluster
channel. These differences are even more pronounced at larger redshifts.
However, these differences are diminished when considering LVK O3 detector
sensitivity. Comparison with LVK O3 data shows that none of our adopted models
can match the data despite large errors on BH-BH masses and redshifts. We
emphasize that our conclusions are derived from a small set of $6$ models that
are subject to numerous known uncertainties. We also note that BH-BH mergers
may originate from a mix of several channels and that other (than those
adopted here) BH-BH formation channels may exist.

###### Subject headings:

stars: black holes, neutron stars, x-ray binaries

##  1\. Introduction

The latest LIGO-Virgo-KAGRA (LVK) Collaboration catalog, GWTC-3 (The LIGO
Scientific Collaboration et al., 2021a), reports the detection of $90$ double
compact object mergers, $83$ of which are confidently classified as binary
black hole mergers. The LVK Collaboration also provided estimates of the
intrinsic population properties of these mergers such as their distributions
of rates, masses, and effective spins. The merger rate density of BH-BH
systems at a fiducial redshift of $z=0.2$ is
$17.3-45{\rm~{}Gpc}^{-3}{\rm~{}yr}^{-1}$ and the rate increases with redshift
(The LIGO Scientific Collaboration et al., 2021c). Masses of BH-BH mergers
show complex multi-feature behavior (Tiwari & Fairhurst, 2021; Edelman et al.,
2022; The LIGO Scientific Collaboration et al., 2021c). BH-BH merger effective
spins are mostly positive and small ($0<\chi_{\rm eff}<0.3$) with $\sim
10-30\%$ contribution of mergers with high positive ($\chi_{\rm eff}>0.3$)
effective spins (Roulet et al., 2021; Abbott et al., 2021, 2021).

Despite all the available information and rapidly increasing number of
detected BH-BH mergers, their origin and specific formation mechanism is not
known. Several formation channels are under close investigation, each
corresponding to a specific formation site of BH-BH mergers and specific
evolutionary processes involved in the evolution of BH-BH progenitors.
Evolution of binary systems in isolation (galactic fields) can provide ample
resources for BH-BH merger formation via several sub-channels; evolution that
involves a common envelope phase (CE, e.g., Belczynski et al. 2016a; Giacobbo
& Mapelli 2018; Kinugawa et al. 2020; García et al. 2021), evolution that only
requires binary interactions through stable Roche-lobe overflow (RLOF, e.g.,
van den Heuvel et al. 2017; Neijssel et al. 2019a; Bavera et al. 2021; Olejak
et al. 2021), and chemically homogeneous evolution that requires rapidly
rotating stars (e.g., Marchant et al. 2016; de Mink & Mandel 2016; Mandel & de
Mink 2016). Evolution of objects in dense stellar systems may also provide
satisfactory conditions for dynamical formation of BH-BH mergers; in young
open clusters (e.g., Banerjee (2018a); Di Carlo et al. (2020); Kremer et al.
(2020); Rastello et al. (2021); Chattopadhyay et al. (2022)), in dense
globular cluster (e.g., Rodriguez et al. 2016b; Askar et al. 2017; Antonini &
Gieles 2020), and in massive nuclear clusters (e.g., Arca-Sedda & Capuzzo-
Dolcetta 2017; Arca Sedda 2020; Fragione et al. 2021; Mapelli et al. 2021).
There are also formation channels that bridge the two above major classes of
isolated and dynamical BH-BH formation. Multiple systems, in particular
triples and quadruples, are also estimated to produce BH-BH mergers through
combination of typical stellar evolution assisted by dynamical interactions of
stars (e.g., Antonini et al. 2017; Safarzadeh et al. 2020). Alternatively, gas
assisted evolution in AGN disks may provide a way to form merging BH-BH
systems (e.g., Samsing et al. 2020; Ford & McKernan 2021).

Identifying the actual BH-BH merger formation channel or channels is hindered
by numerous modeling uncertainties (Arca Sedda & Benacquista, 2019; Belczynski
et al., 2022; Broekgaarden et al., 2021; Gallegos-Garcia et al., 2021). For
example, merger rates for double compact objects vary by orders of magnitude
for each formation channel (e.g., Mandel & Broekgaarden 2021), while predicted
mass distributions are extremely sensitive to various model assumptions (e.g.,
Olejak et al. 2021). Ideally, a given model should fit all physical properties
(rates, masses, spins) for all three double compact object types (BH-BH/BH-
NS/NS-NS mergers; NS: neutron star). Before this is accomplished, we need to
build understanding of how various model parameters change properties of
double compact objects within various formation channels. Numerous studies
were already dedicated to merger rates, NS/BH spins, or compact object masses
(e.g., Timmes et al. (1996); Lipunov et al. (1997); Belczynski et al. (2002);
Voss & Tauris (2003); Podsiadlowski et al. (2003); Zampieri & Roberts (2009);
O’Connor & Ott (2011); Dominik et al. (2012a); Mapelli (2016); Arca Sedda et
al. (2018); Stevenson et al. (2019); Spera et al. (2019); Olejak & Belczynski
(2021); van Son et al. (2021); Ford & McKernan (2021); Broekgaarden et al.
(2021)). Many of the above studies present BH-BH merger masses, but the
evolution of mass with redshift is usually not shown nor discussed in detail.
One notable exception is work by van Son et al. (2022) who present detailed
analysis of primary (more massive) BH mass and its evolution with redshift for
BH-BH mergers for two (RLOF and CE) sub-channels of isolated binary evolution.
None, however, have focused on average total BH-BH mass and its evolution with
redshift. Here we present and discuss the cosmic evolution of average total
mass of BH-BH mergers arising from various evolutionary channels. We note that
empirical studies have been performed to assess changes in the BH-BH mass
distribution shape with redshift, but thus far no definitive evidence of such
mass distribution evolution has been found (Fishbach et al., 2021; Tiwari,
2022; The LIGO Scientific Collaboration et al., 2021c).

##  2\. Models of BH-BH formation

###  2.1. Isolated Binary Evolution Models

We use three models of BH-BH formation from isolated binary stars that were
calculated with the StarTrack population synthesis code (Belczynski et al.,
2002, 2008). The standard evolution involves evolution of massive binary stars
with empirically derived initial conditions (Sana et al., 2012; de Mink &
Belczynski, 2015). In particular, for massive stars the initial mass function
(IMF) is adopted to be a power-law ($\propto M^{-2.3}$) and it extends to
initial star mass of $150{\rm~{}M}_{\odot}$. Input physics includes wind
losses for massive stars: O/B star and WR-type winds in addition to LBV winds
(Vink et al., 2001; Belczynski et al., 2010a). We treat the accretion onto
compact objects during RLOF and from stellar winds using approximations
presented by King et al. (2001); Mondal et al. (2020). CE evolution is treated
with the energy balance prescription of Webbink (1984) and we limit accretion
onto compact objects during a CE phase to $5\%$ of the Bondi rate (MacLeod et
al., 2017). We employ two models of core-collapse supernova (SN) engine in
NS/BH mass calculation (Fryer et al., 2012) that either allow (delayed SN
engine) or do not allow (rapid SN engine) for filling in the lower mass gap
(the dearth of compact objects with mass $\sim 2-5{\rm~{}M}_{\odot}$) with
compact objects between NSs and BHs (Belczynski et al., 2012; Zevin et al.,
2020). We assume the minimum BH mass of $2.5{\rm~{}M}_{\odot}$. In fact, we
use the updated variants of rapid and delayed SN engines from Olejak et al.
(2022). This study allows for gradual filling of the lower mass gap with
compact objects through a single parameter: $f_{\rm mix}$. Isolated binary
models employ either $f_{\rm mix}=0.5$ (close to original delayed SN engine)
and $f_{\rm mix}=4.0$ (close to the original rapid SN engine). For simplicity,
when discussing binary models we will refer to them as having delayed or rapid
SN engine physics. Pulsational pair-instability SNe (PPSN) and pair-
instability SNe (PSN) are assumed to be weak and allow for the formation of
BHs with mass up to $\sim 55{\rm~{}M}_{\odot}$ (Leung et al., 2019; Belczynski
et al., 2020). We use the fallback-decreased NS/BH natal kicks with
$\sigma=265{\rm~{}km}{\rm~{}s}^{-1}$ (Hobbs et al., 2005), and do not allow CE
survival for Hertzsprung gap donors (Belczynski et al., 2007). The results of
population synthesis calculations are post-processed with Madau & Fragos
(2017) star formation rate history in Universe and with evolving with redshift
metallicity distribution for Population I/II stars as presented in Belczynski
et al. (2020) to allow calculation of BH-BH merger population for redshift
range $z=0-15$.

The standard physics presented above produces a majority of BH-BH mergers
through the classical (CE) isolated binary evolution channel. In this channel,
binaries that undergo RLOF with a donor star that is approximately twice the
mass of the accretor or more ($q\gtrsim 2$) are subject to CE evolution that
leads to orbital decay, and may lead to BH-BH merger under favorable
conditions. We adopt model M380.B from Olejak et al. (2021) that produces
$99.5\%$ of BH-BH mergers through CE evolution. Delayed SN engine is employed
in this model.

In our second model, we adopt non-standard CE development criteria that allows
for stable RLOF to proceed even for systems with donors as massive as eight
times mass of the accretor ($q\lesssim 8$; Pavlovskii et al. 2017). Such
evolution, also under favorable conditions, may lead to orbital decay and
formation of BH-BH merger (van den Heuvel et al., 2017). We adopt model M480.B
from Olejak et al. (2021) that produces $94\%$ of BH-BH mergers through stable
RLOF evolution (and the rest through CE). Delayed SN engine is employed in
this model.

In our third model we adopt standard CE development criteria, but extend the
IMF to $200{\rm~{}M}_{\odot}$ and allow for less restrictive PPSN/PSN limits
on BH mass. It was argued that it is possible that the upper BH mass gap is
located in a different place on BH mass spectrum or may not exist at all
(Farmer et al., 2020; Costa et al., 2021; Farrell et al., 2021). We adopt a
model presented by Belczynski (2020) that avoids PPSN mass loss for helium
core masses up to $M_{\rm He}<90{\rm~{}M}_{\odot}$, with stars above this mass
threshold destroyed by PSNe. This allows for formation of BHs with mass up to
$90{\rm~{}M}_{\odot}$ from low-metallicity stars. Delayed SN engine is
employed in this model.

In our fourth model we adopt exactly the same input physics as in the third
model, but we employ the rapid SN engine instead.

This sequence of models with specific choices of input physics was design to
show the dependence of total BH-BH merger mass on some of the most important
and most uncertain parts of stellar and binary physics employed in our
isolated binary evolution models.

###  2.2. Globular Cluster Dynamical Models

We use globular cluster models from the BH-BH formation channel analysis of
Zevin et al. (2021), which have data relevant to this work publicly available
on Zenodo (Zevin, 2020). These models originate from a grid of $96$ $N$-body
models of collisionless star clusters simulated using the Hénon-style cluster
Monte Carlo code CMC (Hénon, 1971, 1971; Joshi et al., 2000; Pattabiraman et
al., 2013), as described in Rodriguez et al. (2019). Cluster birth times and
cluster initial metallicity follow cosmological models for globular cluster
formation (El-Badry et al., 2019), as outlined in (Rodriguez & Loeb, 2018) and
Rodriguez et al. (2018a). In particular, clusters with three metallicities are
used $Z=0.005,\ 0.001,\ 0.0002$. The initial binary fraction for stars born in
these cluster models is set to $10\%$ and is independent of stellar mass; the
true value for the binary fraction of stars born in globular clusters is
observationally unconstrained.

Stars and binaries that are born in GCs are evolved with rapid evolutionary
formulas (Hurley et al., 2000, 2002) with updated prescriptions for stellar
winds, compact object masses, and supernova natal kicks (see Rodriguez et al.,
2016a, 2018b). Initial masses of single stars and primaries in binary systems
are taken from a Kroupa IMF ($\propto M^{-1.3}$, in the range
$0.08-0.5{\rm~{}M}_{\odot}$ and $\propto M^{-2.3}$, in the range
$0.5-150{\rm~{}M}_{\odot}$), with secondary components in binaries drawn from
a uniform mass ratio distribution. Semimajor axes of initial binaries are
drawn flat in log ($P(a)da\propto 1/a$) between the point of stellar contact
and the hard-soft boundary, where the orbital velocity of the binary equals
the typical velocity of particles in the cluster (for a typical velocity of
$10{\rm~{}km}{\rm~{}s}^{-1}$ this translates to $\sim 1-5\times
10^{5}{\rm~{}R}_{\odot}$ for equal-mass binaries with component mass between
$30{\rm~{}M}_{\odot}$ and $150{\rm~{}M}_{\odot}$). Eccentricities are drawn
from a thermal distribution ($p(e)de=2ede$). Stellar winds for massive stars
are adopted from Vink et al. (2001) during the main sequence and from Vink &
de Koter (2005) for naked helium stars, with later evolutionary stages
following Belczynski et al. (2010a). Compact object masses are calculated
based on 1-D supernova models that allow for fallback and direct BH formation,
with a rapid supernova engine that produces a lower mass gap between NSs and
BHs (Fryer et al., 2012). Note that for isolated binary evolution BH-BH
formation we employed models with delayed SN engine. However, as we show in
Section 3.1 the choice of SN engine model does not impact average total BH-BH
merger mass and its evolution with redshift. Massive stars (with He cores
above $45{\rm~{}M}_{\odot}$ are subject to PPSN/PSN and stellar-origin natal
BH mass is limited to $\lesssim 45{\rm~{}M}_{\odot}$. BH natal kicks are taken
from a Maxwellian distribution with 1-D $\sigma=265{\rm~{}km}{\rm~{}s}^{-1}$,
and kicks are decreased by fallback (i.e., no natal kicks for massive BHs that
form through direct core collapse). This set of models assumes that BHs are
born with near-zero spin due to highly efficient angular momentum transport
(e.g., Fuller & Ma 2019; Qin et al. 2018; Belczynski et al. 2020). Spin-up of
massive stars (BH progenitors) by tidal interactions in close binaries is not
taken into account. Stars ($90\%$) and primordial binaries ($10\%$) in each
cluster form the first generation of BHs (with masses limited to $\lesssim
45{\rm~{}M}_{\odot}$).

Dynamical interactions and also mergers of BH-BH systems formed from
primordial binaries produce future generations of BHs. Obviously, these BHs
can be more massive than the first generation of BHs. The second generation of
BHs is typically formed with high spins ($a\sim 0.7$) due to the angular
momentum of the pre-merger double compact object orbit. Spins of the BH
remnant are calculated using the precession package (Gerosa & Kesden, 2016),
as described in Rodriguez et al. (2018a). Asymmetries in the BH components of
the merger (unequal masses, spins) will lead to a dispersion of spin
magnitudes around $a\sim 0.7$. Some second-generation BHs receive
gravitational-wave recoil kicks that may remove BHs from clusters, though
since the models we consider assume BHs are born with negligible spin, almost
$\sim 60\%$ of second-generation BHs are retained since recoil kicks depend
then entirely on mass ratio and are typically $\lesssim
100{\rm~{}km}{\rm~{}s}^{-1}$ (Rodriguez et al., 2019). Those that remain bound
continue to contribute to cluster evolution and to GC BH-BH merger
populations. Due to the relatively low escape velocities of globular clusters
($\sim 30-100{\rm~{}km}{\rm~{}s}^{-1}$), higher-generation merger products are
typically ejected from the cluster via gravitational-wave recoil kicks due to
their components having large spins as a result of their formation from
previous BH-BH mergers; $\lesssim 0.1\%$ mergers containing a second-
generation BH remain bound to the cluster (Rodriguez et al., 2019).

###  2.3. Open Cluster Dynamical Models

In this work, the long-term N-body evolutionary model set of star clusters as
described in Banerjee (2021c) is utilized. The computed models are described
in detail in Banerjee et al. (2020); Banerjee (2021b, c) and are summarized
below.

The model star clusters, initially, have mass within $2\times
10^{4}{\rm~{}M}_{\odot}\leq M_{\rm cl}\leq 10^{5}{\rm~{}M}_{\odot}$, size
(half-mass radius) $1{\rm~{}pc}\leq r_{\rm h}\leq 2{\rm~{}pc}$, and
metallicity $0.0001\leq Z\leq 0.02$, and they are subjected to a solar-
neighborhood-like external galactic field. The initial models are composed of
ZAMS stars of masses $0.08{\rm~{}M}_{\odot}\leq m_{\ast}\leq
150.0{\rm~{}M}_{\odot}$ that are distributed according to the standard IMF
(Kroupa, 2001). About half of the models have a primordial-binary population
(overall initial binary fraction $f_{\rm bin}\approx 5$% or 10%) where all
O-type stars (stars with ZAMS mass $\geq 16{\rm~{}M}_{\odot}$) are paired
among themselves (i.e, initial binary fraction among O-type stars is $f_{\rm
Obin}=100$%) following an observationally-deduced distribution of massive-star
binaries (Sana & Evans, 2011; Sana et al., 2013; Moe & Di Stefano, 2017). Such
cluster parameters and stellar membership are consistent with those observed
in ‘fully’-assembled, (near-)spherical, (near-)gas-free young massive clusters
(YMCs) which evolve into medium-mass open clusters (OCs) (Portegies Zwart et
al., 2010). They continue to form, evolve, and dissolve in the Milky Way and
other galaxies (as such, anywhere in the Universe) as a part of the galaxies’
ongoing star formation.

These model clusters evolve due to two-body relaxation (Spitzer, 1987), close
(relativistic) dynamical encounters (Heggie & Hut 2003, without applying any
gravitational softening), and stellar evolution. This is achieved using ${\tt
NBODY7}$, a state-of-the-art post-Newtonian (PN) direct N-body integration
code (Aarseth, 2012; Nitadori & Aarseth, 2012) that couples with the semi-
analytical stellar- and binary-evolutionary scheme ${\tt BSE}$ (Hurley et al.,
2000, 2002). As detailed in Banerjee et al. (2020), the integrated ${\tt BSE}$
is updated in regards to prescriptions of stellar wind mass loss and formation
of NSs and BHs. NSs and BHs form according to the ‘rapid’ or ‘delayed’ core-
collapse SN models of Fryer et al. (2012), and the PPSN and PSN models of
Belczynski et al. (2016b). The majority of the computed models of Banerjee
(2021c) employ the rapid-SN prescription although a few models employ the
delayed-SN prescription, for exploratory purposes. The dynamical evolution and
merger outcomes of the clusters are unlikely to be significantly affected by
this difference, as discussed in Banerjee et al. (2020); Banerjee (2020).

A newly formed NS or BH receives a natal kick that is modulated based on SN
matter fallback onto it, as in Belczynski et al. (2008). The adopted fallback
model is that of Fryer et al. (2012). For all core-collapse SN remnants, the
‘base’ natal kick is assigned from a Maxwellian distribution with one-
dimensional dispersion of $\sigma=265{\rm~{}km}{\rm~{}s}^{-1}$, i.e., the same
velocity dispersion as that of the Galactic single NS population (Hobbs et
al., 2005). However, due to momentum conservation, the material fallback slows
down the remnants, causing BHs of $\gtrsim 10{\rm~{}M}_{\odot}$ to be retained
in the clusters right after their birth. The material fallback also shapes the
mass distribution of NSs and BHs. NSs formed via electron-capture SN (ECS;
Podsiadlowski et al. 2004) also receive small natal kicks of a few
${\rm~{}km}{\rm~{}s}^{-1}$ (Gessner & Janka, 2018) and are retained in the
clusters at birth. The present remnant scheme forms first-generation of BHs in
the range $\approx 3-45{\rm~{}M}_{\odot}$ (the upper limit set by PPSN) and
exhibits a PSN mass gap between $\approx 45-120{\rm~{}M}_{\odot}$ (see Fig. 2
of Banerjee et al. 2020). See Banerjee et al. (2020) for further details and
the description of the implementations in ${\tt BSE}$ and ${\tt NBODY7}$.

In ${\tt NBODY7}$, the PN treatment is handled by the ${\tt ARCHAIN}$
algorithm (Mikkola & Tanikawa, 1999; Mikkola & Merritt, 2008). The treatment
allows for general relativistic (GR) evolution of the innermost NS- and/or BH-
containing binary of an in-cluster (i.e., gravitationally bound to the
cluster) triple or higher order compact subsystem, in tandem with the
Newtonian-dynamical evolution of the subsystem (Kozai-Lidov oscillation or
chaotic three-body interaction), leading to the binary’s (in-cluster) GR in-
spiral and potential merger. The PN treatment applies also to the GR evolution
of in-cluster NS/BH-containing binaries that are not a part of a higher-order
subsystem. As demonstrated in previous studies (Banerjee et al., 2010;
Banerjee, 2017, 2018b; Rastello et al., 2019; Banerjee, 2021b), the moderate
density ($10^{2}-10^{3}$ stars pc-3) and velocity dispersion ($\lesssim
10{\rm~{}km}{\rm~{}s}^{-1}$) in the model clusters make them efficient in
dynamically assembling PN subsystems, particularly, those composed of BHs. On
the other hand, the clusters’ low escape speed ($\sim
10{\rm~{}km}{\rm~{}s}^{-1}$) causes the dynamically-ejected BH-BH systems to
be typically wide (Portegies Zwart & McMillan, 2000), with a GR inspiral time
(Peters, 1964) often exceeding the Hubble time. These cause the majority of
the GR mergers from these computed clusters to be in-cluster BH-BH mergers.

In the computed models, star-star and star-remnant mergers happen due to
evolution of massive primordial binaries and three-body interactions. In the
present models, a star-star merger results in the loss of $\approx 10\%$ of
the total merging mass as suggested by hydrodynamical calculations (J. C.
Lombardi et al., 2002; Gaburov et al., 2008). In a BH-star merger (forming a
BH-Thorne-Zytkow Object, BH-TZO), $\approx 90\%$ of the star’s mass is taken
to be accreted onto the BH. The model grid used in this work comprises 74
long-term ($\sim 10$ Gyr) evolutionary cluster models (see Table A1 of
Banerjee 2021c).

##  3\. Results

Figure 1.—  Top: intrinsic average total BH-BH merger mass change as a
function of redshift for different BH-BH formation channels: globular cluster
and open cluster dynamical channels, and CE and RLOF isolated binary evolution
channels. For the RLOF channel we show three variations: first allows for
stars to form to $150{\rm~{}M}_{\odot}$ and PPSN/PSN limit BH masses to $\sim
55{\rm~{}M}_{\odot}$, second allows for star formation to
$200{\rm~{}M}_{\odot}$ and PPSN/PSN limit BH masses to $90{\rm~{}M}_{\odot}$,
while third is the same as the second variation but employs rapid SN engine
(all other isolated binary models employ delayed SN engine). We show the
approximate redshift range for design (future) LIGO/Virgo/KAGRA detection of
BH-BH mergers (the range for Einstein Telescope and Cosmic Explorer covers
entire plot). Bottom: detected average total BH-BH merger mass as a function
of redshift. The gravitational-wave instrumental and observational selection
effects are applied to the models from the top panel to show which BH-BH
mergers may be currently (O3 sensitivity of LVK) detectable. Note that rather
different intrinsic curves (top), become almost indistinguishable when
detection biases are applied (bottom).

###  3.1. Evolutionary Analysis

Our main result is shown in the top panel of Figure 1. It shows evolution of
the average intrinsic (not redshifted nor detection weighted) total BH-BH mass
($\overline{M}_{\rm tot,int}=M_{\rm bh1}+M_{\rm bh2}$) with redshift ($z$).
The evolution is shown for the six models (see Sec. 2) of BH-BH merger
formation: the globular cluster dynamical channel, the open cluster dynamical
channel, the classical CE isolated binary evolution channel, the RLOF isolated
binary evolution channel, the extended RLOF channel (with more massive stars
and BHs), and the extended RLOF channel with rapid SN engine.

Within the approximate LIGO-Virgo-KAGRA design redshift range for detection of
BH-BH mergers ($z<2$)111The maximum detectable redshift of a
50$M_{\odot}$-50$M_{\odot}$ BH merger with Advanced LIGO at design sensitivity
is $z_{h}=1.9$, as calculated with the gravitational wave distance calculator
described in Chen et al. (2021)., the CE channel produces the lowest mass BH-
BH mergers ($\overline{M}_{\rm tot,int}\sim 25{\rm~{}M}_{\odot}$, irrespective
of the redshift), the three RLOF channels produce more massive BH-BH mergers
that show a slow increase of the average mass with increasing redshift as the
metal-poor environments at higher redshifts produces more massive BHs
($\overline{M}_{\rm tot,int}\sim 30\rightarrow 35{\rm~{}M}_{\odot}$ for
$z=0\rightarrow 2$), and the globular cluster channel produces the most
massive BH-BH mergers with fast increase of the average mass with redshift
($\overline{M}_{\rm tot,int}\sim 45\rightarrow 60{\rm~{}M}_{\odot}$ for
$z=0\rightarrow 2$). The open cluster channel generates the most dramatic
evolution of the intrinsic average total BH-BH merger mass ($\overline{M}_{\rm
tot,int}\sim 25\rightarrow 50{\rm~{}M}_{\odot}$ for $z=0\rightarrow 2$).

In the entire presented redshift range ($z<15$) that will be accessible to
Einstein Telescope and Cosmic Explorer, the four isolated binary evolution
channels produce gradually more massive BH-BH mergers with increasing redshift
(at $z=15$ we find: $\overline{M}_{\rm tot,int}\sim 40{\rm~{}M}_{\odot}$ for
the CE channel, $\overline{M}_{\rm tot,int}\sim 55{\rm~{}M}_{\odot}$ for the
RLOF channel, and $\overline{M}_{\rm tot,int}\sim 70{\rm~{}M}_{\odot}$ for
both RLOF channels with extended IMF). For BH-BH mergers in the globular
cluster channel, we find that the average mass first increases with redshift
(reaching a maximum of $\overline{M}_{\rm tot,int}\sim 65{\rm~{}M}_{\odot}$ at
$z=3-5$) and then gradually decreases with increasing redshift. For the open
cluster channel, we find that the average BH-BH mass first increases with
redshift, reaching high values ($\overline{M}_{\rm tot,int}\sim
70{\rm~{}M}_{\odot}$ at $z=4-9$) and then sharply decreases for higher
redshifts.

The general increase of the intrinsic average BH-BH mass with redshift for all
four isolated binary evolution channels is easily understood in the framework
of cosmic chemical evolution. The chemical composition of stars (i.e.,
metallicity $Z$) is the primary factor setting the mass of stellar-origin BHs
(Belczynski et al., 2010a). At lower metallicities, stellar winds are weaker
(which removing mass from stars), and more massive BHs are able to form. Since
metallicity decreases with redshift, more massive stellar-origin BHs form at
higher redshifts. There is an exception from this general trend for the CE BH-
BH formation channel at low redshifts ($z\lesssim 2-3$), where the BH-BH
merger total mass remains constant and does not change with redshift. We
explain this behavior in Section 5.

The lowest average total BH-BH merger mass (at any redshift) is found for the
CE BH-BH binary formation channel. The actual value of the average mass is a
combination of the adopted stellar wind mass loss for massive stars and their
dependence on metallicity, the cosmic evolution of metallicity with redshift,
and the specific evolutionary sequences (set by stellar/binary input physics
assumptions) producing BH-BH mergers in the CE binary channel. The RLOF binary
channel produces on average more massive BH-BH mergers. Since for this channel
wind mass loss and cosmic metallicity evolution are the same as for the CE
channel, the entire difference in average total mass originates from different
stellar/binary evolutionary assumptions. We explain the (non-trivial) origin
of this difference in Section 6 Finally, the two RLOF models with high
(extended) IMF naturally produce the most massive BH-BH mergers among our
isolated binary models. In these two models, BHs may form up to
$90{\rm~{}M}_{\odot}$ as compared to maximum mass of a BH of $\sim
55{\rm~{}M}_{\odot}$ in the other two binary models. We note slightly larger
total BH-BH masses ($\lesssim$ few ${\rm~{}M}_{\odot}$) for model with the
rapid SN engine as the BHs that would form in the lower mass gap in the
delayed SN engine model are simply forming with masses just above the mass
gap. Note that the small differences between these two models with extended
IMF but with different SN engines are much smaller than differences associated
with the BH-BH formation mode (CE versus RLOF channels) or the change of IMF
and PPSN/PSN model.

The globular cluster channel produces on average significantly more massive
BH-BH mergers at low and intermediate redshifts ($z<5$) compared to any
presented isolated binary evolution model by about $\overline{M}_{\rm
tot,int}\sim 10-20{\rm~{}M}_{\odot}$ (see Fig. 1). This is due to the most
massive BHs in the cluster more readily synthesizing binaries through
dynamical interactions, as well as the contribution of mergers containing a
second-generation BH that was retained in the cluster following its first
merger. The average total merger mass drops at high redshifts ($z\gtrsim 7$)
when GC BH-BH mergers become more dominated by BH-BH mergers from primordial
binaries.

Unlike the isolated evolution channels considered, the globular cluster model
has a peak in the average mass of BH-BH mergers near the peak of globular
cluster formation, at redshifts of $z\sim 3-5$. The most massive BHs in the
cluster ($\gtrsim 20{\rm~{}M}_{\odot}$) are typically the first to be
dynamically processed, either merging in the cluster at early times or being
ejected and inspiraling over longer timescales. Less massive BHs born in the
clusters ($\sim 10-15{\rm~{}M}_{\odot}$) follow suit and typically merge at
later times than their more massive counterparts, leading to the decrease in
the average total mass lower redshifts ($z\lesssim 3$). At $z\gtrsim 5$,
primordial binaries in the cluster that evolve and lead to BH-BH mergers
formed in quasi-isolation (without much assistance from dynamical
interactions) have a larger contribution to the total population of BH-BH
mergers, driving the overall decline in the average mass at large redshifts.
These systems rapidly merge on timescales shorter than the mass segregation
timescale in their GC hosts of $\mathcal{O}(100\,\mathrm{Myr})$. However,
dynamically processed mergers still occur at these higher redshifts, which
combined with small number statistics (GCs are just beginning to form) leads
to the fluctuating behavior in the average total mass evolution at $z\gtrsim
5$. The relative contribution of BH mergers from quasi-isolated evolution
prior to the peak of globular cluster formation and their impact on the
average mass change with redshift depends on the metallicity evolution of
cluster formation, which is distinct from the metallicity evolution of
isolated binaries in galactic fields.

Interpretation along lines similar to the case of GCs applies to the evolution
of the average total mass with merger redshift for OCs (Figure 1; top panel).
However, a twist for OCs is that such systems, being 1-2 orders of magnitude
less populous than but of similar (initial) length scale as GCs, have a factor
of a few to tens shorter two-body relaxation time (Spitzer, 1987) compared to
GCs. In OCs, the most massive BHs/BH-progenitors segregate to the cluster
center in $\lesssim 10$ Myr (see Banerjee, 2021a). Also, before the formation
of BHs, a fraction of the BH progenitors already undergo complete mass
segregation due to the shorter relaxation/mass-segregation timescale. These
cause the dynamical processing and mergers among BHs to commence at earlier
ages (i.e. at larger redshifts) and closer to the quasi-isolated mergers in
OCs than in GCs. Also, due to the smaller member number and hence primordial-
binary population in OCs, the models produce a much smaller number of quasi-
isolated BH-BH mergers per cluster than the GC models (Fragione & Banerjee,
2020).

The effects of these differences in OCs can be seen in Figure 1. As opposed to
GCs, $\overline{M}_{\rm tot,int}$ for OCs rises sharply at $z\approx 9$ (fewer
quasi-isolated mergers), thereafter remaining nearly flat and peaking weakly
in between $5<z<6$ (pre/faster mass segregation), and drops off towards low
$z$ faster than GCs (faster dynamical processing due to shorter relaxation
times in OCs). It is also to be noted that the OCs are taken to follow a
different metallicity-dependent cosmic star formation history (Chruslinska &
Nelemans (2019) their ’moderate-Z’ relation; see Banerjee (2021c) for the
details) than those for both the isolated-binary and GC models presented here,
which would also cause differences in the redshift dependence of
$\overline{M}_{\rm tot,int}$. Although OCs are inefficient in producing
massive BHs via hierarchical mergers, the computed OC models adopt a small,
$\approx 10$%, mass loss in star-star mergers (based on hydrodynamic stellar
merger calculations; see Sec. 2.3), which is the main channel for forming the
most massive BHs lying within the PPSN/PSN mass gap in these models (Banerjee,
2021b, a). Furthermore, a $90\%$ BH-TZO accretion is adopted (Sec. 2.3).
These, along with the tendency of having lower metallicity at high redshifts
than globular cluster or binary models, results in a higher maximum
$\overline{M}_{\rm tot,int}$ for OCs (Figure 1).

###  3.2. Comparison to Observed Gravitational Waves

Although all models explored in this work show notable differences in the
average total mass of mergers and its evolution with redshift (the only
exception being the two RLOF models with extended IMF: one with rapid and one
with delayed SN engine), there is no guarantee these differences can be
measured with current gravitational-wave sensitivity. Firstly, gravitational-
wave sources are subject to strong selection effects, which not only limit the
distance to which BH-BH mergers can be detected, but also preferentially
enable detection of certain binary mass and spin configurations. As a result,
the properties of the detected population of BH-BH mergers will not match
those of the intrinsic population. Furthermore, gravitational-wave detectors
are subject to a variety of sources of noise, so the true masses, spins, and
distances of detected events can only be estimated, often with sizable
uncertainties. Finally, we are still in the era of small number statistics
when it comes to gravitational-wave detections, so even if the true properties
of a detected system were known, there still exists a counting uncertainty in
the number of similar events occurring in the Universe.

Given these combined effects as well as the poorly understood systematic
uncertainties in the six population synthesis models presented here, we take a
simple approach to assessing whether these models are distinguishable with
current detections, and whether any of them are a reasonable fit to the data,
in analog to the investigations in Fishbach et al. (2021). In Figure 2, we
show the average detected mass $\overline{M}_{\rm tot,det}$ as a function of
redshift under each of our six models by computing the average detected total
mass in 14 equally spaced redshift bins (i.e. 15 bin edges) from $z=0$ to
$z=1.5$. These curves are also reproduced in Figure 1 for comparison to the
intrinsic average masses. $\overline{M}_{\rm tot,det}$ is essentially
$\overline{M}_{\rm tot,int}$ filtered through the LVK detector and analysis
selection function. To simulate these effects of the detector and analysis
selection bias on the detected populations, we re-weight the detected software
injections provided in The LIGO Scientific Collaboration et al. (2021c) by the
intrinsic simulated populations. To perform the re-weighting, a histogram of
the population synthesis samples from each model is used to approximate the
intrinsic simulated population probability distribution over masses and
redshifts.

The model curves in Figure 2 show the average detected total mass increasing
as a function of redshift, which is expected given that less-massive sources
are detectable out to lower redshifts than the more massive sources. The rate
of change of the average detected total mass with redshift is approximately
the same for all our models, though the CE channel shows the strongest
increase with redshift. To compare these curves to those in Figure 1, consider
that the slopes of the $\overline{M}_{\rm tot,det}$ lines are strongly
influenced by the selection effects, which will increase the lowest detectable
mass with increasing redshift and therefore force models with low, constant
intrinsic average mass like the CE channel to appear to evolve with redshift
in the detected population. In contrast, for the open clusters, the average
intrinsic mass changes from $\sim 25$ to $\sim 50M_{\odot}$ between redshift 0
and 1 which matches what is seen in the detected average mass. This means that
the true average mass evolution with redshift in the open cluster model is
roughly following the O3 LVK selection effects. The similarity between slopes
of the $\overline{M}_{\rm tot,det}$ curves suggests that the changes in the
total mass distribution mean over redshift is not a good distinguisher of
these models given current sensitivity. However, the overall mass scale is
more telling: the GC model predicts a higher average detected mass than the
other models by $\sim 10-20M_{\odot}$, regardless of redshift.

To compare these simulated average detected masses to the masses of detected
BH-BHs, in gray we overlay the $90\%$ credible regions on detector-frame total
mass and redshift for the BH-BH systems in GWTC-3 from the third LVK observing
run (Abbott et al., 2021, 2021; The LIGO Scientific Collaboration et al.,
2021a). Only BH-BH mergers from O3 are selected for uniform sample. The
credible regions are generated via kernel density estimates (KDEs) of
posterior samples provided by the LVK, re-weighted to a prior that
approximates empirical measurements of the merger population distribution (The
LIGO Scientific Collaboration et al., 2021c; Roulet et al., 2021; The LIGO
Scientific Collaboration et al., 2021b). Specifically, we use (1) a broken
power law prior on the more massive BH in the binary’s mass, with
$p(m_{1})\propto m_{1}^{-1.7}$ for $m_{1}<40M_{\odot}$, and $p(m_{1})\propto
m_{1}^{-5}$ for $m_{1}\geq 40M_{\odot}$, (2) a $p(q)\propto q$ prior for the
mass ratio $q$, and (3) a uniform-in-comoving-time-volume prior on $z$, i.e.
$p(z)\propto\frac{dV}{dz}\frac{1}{1+z}$. Additionally, we estimate the mean
mass of these detections over redshift by computing a single total-mass mean
from the combined samples of all events in the same redshift bins used for the
six models’ curves. We show this estimate in gray in Figure 2.

At $z\sim 0.5$, all the models appreciably under-predict the estimated mean
detected mass by $\sim 25M_{\odot}$, except GC, which underpredicts by $\sim
10M_{\odot}$. This is consistent with the fact that the GC model is the only
model with a non-negligible fraction ($>1\%$) of systems with $M_{\rm
tot,int}>100M_{\odot}$ below $z=1$. At the highest redshifts, it is not clear
which model(s) may fit the data well given that there are few events detected
at these distances, and those events have large uncertainties and low
significance. At redshifts less than $\sim 0.2$, all the models over-predict
the average detected mass. This is simply because none of the models predict a
sufficiently high rate of low-mass mergers. Overall, these findings,
especially at low masses and redshifts, suggest that none of these models
alone (or even a mixture of them), are a good fit to GWTC-3 events. However,
the GC model offers the closest $\overline{M}_{\rm tot,det}$ prediction for
high redshift events, but would require either a different set of input
physics or another channel also with a different set of input physics than
explored in this work to fill in the low and high ends of the mass
distribution.

We emphasize that the explorations here are primarily illustrative. None of
the models herein have been tuned to fit the LVK data, nor have we performed a
simultaneous fit to all binary parameters in the LVK catalog accounting for
measurement uncertainties. However, comparisons such as the one shown in
Figure 2 will be invaluable for assessing the viability of BH-BH formation
channel models and enable a narrowing of the parameter space for easier
follow-up studies that perform a complete statistical analysis.

Figure 2.— The average detected total mass as a function of redshift for the
six models considered in this work. The models are overlaid on the gray $90\%$
credible interval parameter regions for BH-BH mergers from LVK O3, re-weighted
to an astrophysically realistic prior. The gray line shows an estimate of the
mean detected total mass of these events in different redshift bins for all
posterior samples under the aforementioned prior.

##  4\. Conclusions

We have employed several BH-BH merger formation models to show evolution of
the average intrinsic and detected total BH-BH merger mass with redshift. Our
main result is shown in Figures 1 and 2. In our analysis we have included the
classical isolated binary CE evolution formation channel, three stable RLOF
channels: one with a standard IMF (stars with masses $<150{\rm~{}M}_{\odot}$
and PPSN/PSN limiting BH masses to $<55{\rm~{}M}_{\odot}$), and two with an
extended IMF (stars with masses $<150{\rm~{}M}_{\odot}$ and PPSN/PSN limiting
BH masses to $<90{\rm~{}M}_{\odot}$) but with different SN engines (delayed
and rapid), and globular cluster and open cluster dynamical BH-BH formation
channels.

The model differences in average intrinsic BH-BH total mass are both
quantitative and qualitative. The lowest average masses are produced by the CE
binary channel ($\overline{M}_{\rm tot,int}\sim 25{\rm~{}M}_{\odot}$ at
$z=0-2$), with somewhat larger masses found for RLOF channels
($\overline{M}_{\rm tot,int}\sim 30{\rm~{}M}_{\odot}$ at $z=0$ and increasing
to $\sim 35{\rm~{}M}_{\odot}$ at $z=2$). The largest masses in
LIGO/Virgo/KAGRA range are produced by globular clusters ($\overline{M}_{\rm
tot,int}\sim 45{\rm~{}M}_{\odot}$ at $z=0$ with notable increase to $\sim
60{\rm~{}M}_{\odot}$ at $z=2$). Open cluster BH-BH mergers show the most
pronounced evolution of the average total mass in this redshift range
($\overline{M}_{\rm tot,int}\sim 25{\rm~{}M}_{\odot}$ at $z=0$ and $\sim
50{\rm~{}M}_{\odot}$ at $z=2$), bridging isolated binary and globular cluster
models.

These differences become more pronounced for higher redshifts. Isolated binary
channels show monotonic increase of the average intrinsic BH-BH merger mass
with redshift, with the exception that for the CE channel the average
intrinsic mass is constant for low redshifts ($z\lesssim 2-3$) and then
increases with increasing redshift. A different behavior is found for the
globular cluster model, for which the average intrinsic merger mass increases
with redshift to $z\sim 4$, and then decreases with redshift. The open cluster
model shows similar behavior to the globular cluster model, however the
largest intrinsic total BH-BH masses are produced at higher redshifts
($\overline{M}_{\rm tot,int}\gtrsim 70{\rm~{}M}_{\odot}$ for $z\sim 5-9$).
These trends are understood in the framework of input physics employed in the
models (see Sec. 3.1 for details).

We have shown that not only the formation channel (e.g., isolated binary
versus cluster BH-BH formation) significantly affects intrinsic average BH-BH
merger mass, but that some parts of highly uncertain input physics within a
given formation channel have big impact. For example, changing CE development
criteria for isolated binary evolution channel significantly impacts BH-BH
merger formation mode (CE versus RLOF) and the intrinsic total BH-BH merger
mass. Similar statement can be made about impact of treatment of PPSN/PSN and
the extend of IMF in binary evolution calculations. However, not all input
physics is that important for the intrinsic average total BH-BH merger mass
and its evolution with redshift. For example, SN engine model that regulates
the existence and numbers of compact objects with mass in the lower mass gap
($\sim 2-5{\rm~{}M}_{\odot}$) does not significantly impacts intrinsic average
total BH-BH merger mass. Since all these parts of physics of BH-BH merger
formation (CE, SN engine, PPSN/PSN) are highly uncertain we adopted models
with various assumptions on parameters regulating these processes.

The evolution of the averaged detected BH-BH merger total mass does not differ
as much between formation channels compared to the intrinsic mass (see Fig.
2). The four binary models and the open cluster model show a similar behavior
for redshifts in which mergers are currently detected ($z\lesssim 1$). Only
the globular cluster model results in moderately heavier detected BH-BH
mergers (by $\lesssim 10{\rm~{}M}_{\odot}$ for $z\lesssim 1$). This is caused
by LVK detection biases that affect different models in different ways (as
explained in Sec. 3.2). These differences in mass evolution with redshift are
too small to differentiate between models with the current LVK data due to
large detection errors. Additionally, none of the employed models matches the
average detected BH-BH merger masses over entire BH-BH detectable redshift
range. However, the adopted globular cluster model matches the current data
better than our isolated binary or open cluster models at moderate redshifts
$z\sim 0.5$.

The physics of BH-BH merger progenitors for isolated binary, globular cluster,
and open cluster formation channels is still subject to many uncertainties
(Belczynski et al., 2022; Santoliquido et al., 2021; Chatterjee et al., 2017;
Conroy & Gunn, 2010). Major overall uncertainties include the underlying
cosmic star formation rate of isolated binaries and stars/binaries in globular
and open clusters, and metallicity evolution with redshift (Sharda et al.,
2021) within each formation channel (Neijssel et al., 2019b). Basic stellar
evolution uncertainties include fusion reaction rates (Fields et al., 2018),
mixing in stellar interiors (Zhang, 2013), stellar wind mass loss rates
(Björklund et al., 2022), radial expansion, and core-collapse/supernova
physics (Schneider et al., 2019) associated with compact object formation
(Belczynski et al., 2022), as well as initial conditions such as the initial
mass function, mass ratio, orbital separation, and eccentricity (Klencki et
al., 2018). Binary evolution physics is highly uncertain in its treatment of
stellar interactions, whether these are stable or unstable (CE) RLOF events
(Ablimit & Maeda, 2018; Dominik et al., 2012b). In modeling of dynamical BH-BH
formation unknown initial conditions of globular and open clusters (size,
density, binary fractions) along with uncertain outcomes of stellar mergers
and BH (natal and gravitational-wave) kicks hinders their predictions.

Keeping these uncertainties in mind, we can deliver only a very limited
conclusion. For some choices of input physics in a given BH-BH formation
channel (six presented models) we predict that the average intrinsic total BH-
BH merger mass evolution with redshift changes from model to model, sometimes
in rather significant ways. However, due to BH-BH detection selection effects,
small statistics of detected sources, and significant uncertainties in
detected BH-BH masses and redshifts, the current LVK O3 data does not enable
us to identify any model as fitting the data better than another.
Additionally, no models (or their potential combinations) fit the entirety of
the data well due to their dearth of low-mass events relative to the detected
population. More data from LVK in O4 and O5 may change this conclusion. At
this moment, it is clear that the average total mass (and its evolution with
redshift) may be very different for various formation models, and this may
help in the future to eliminate some models and support others, and guide
theoretical work on the physics of BH-BH merger progenitors. The total BH-BH
merger mass and its evolution with redshift may provide, along other
observables, an extra leverage in search for the BH-BH formation site/channel.

KB and AO acknowledge support from the Polish National Science Center grant
Maestro (2018/30/A/ST9/00050). ZD is supported by the CIERA Board of Visitors
Research Professorship. Support for MZ is provided by NASA through the NASA
Hubble Fellowship grant HST-HF2-51474.001-A awarded by the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Incorporated, under NASA contract NAS5-26555. DC is
supported by the STFC grant ST/V005618/1. SB acknowledges support from the
Deutsche Forschungsgemeinschaft (DFG; German Research Foundation) through the
individual research grant “The dynamics of stellar-mass black holes in dense
stellar systems and their role in gravitational-wave generation” (BA 4281/6-1;
PI: S. Banerjee). SB acknowledges the generous support and efficient system
maintenance of the computing teams at the AIfA and HISKP. This work was
initiated and performed in part at Aspen Center for Physics, which is
supported by National Science Foundation grant PHY-1607611. This material is
based upon work supported by NSF’s LIGO Laboratory which is a major facility
fully funded by the National Science Foundation.

## References

  * Aarseth (2012) Aarseth, S. J. 2012, MNRAS, 422, 841
  * Abbott et al. (2021) Abbott, R., et al. 2021, Physical Review X, 11, 021053
  * Abbott et al. (2021) Abbott, R., et al. 2021, arXiv e-prints, arXiv:2108.01045
  * Ablimit & Maeda (2018) Ablimit, I., & Maeda, K. 2018, ApJ, 866, 151
  * Antonini & Gieles (2020) Antonini, F., & Gieles, M. 2020, Phys. Rev. D, 102, 123016
  * Antonini et al. (2017) Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77
  * Arca Sedda (2020) Arca Sedda, M. 2020, ApJ, 891, 47
  * Arca Sedda et al. (2018) Arca Sedda, M., Askar, A., & Giersz, M. 2018, ArXiv e-prints
  * Arca Sedda & Benacquista (2019) Arca Sedda, M., & Benacquista, M. 2019, MNRAS, 482, 2991
  * Arca-Sedda & Capuzzo-Dolcetta (2017) Arca-Sedda, M., & Capuzzo-Dolcetta, R. 2017, ArXiv e-prints
  * Askar et al. (2017) Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36
  * Banerjee (2017) Banerjee, S. 2017, MNRAS, 467, 524
  * Banerjee (2018a) —. 2018a, MNRAS, 473, 909
  * Banerjee (2018b) —. 2018b, MNRAS, 481, 5123
  * Banerjee (2020) Banerjee, S. 2020, Phys. Rev. D, 102, 103002
  * Banerjee (2021a) Banerjee, S. 2021a, arXiv e-prints, arXiv:2109.14612
  * Banerjee (2021b) —. 2021b, MNRAS, 500, 3002
  * Banerjee (2021c) —. 2021c, MNRAS, 503, 3371
  * Banerjee et al. (2010) Banerjee, S., Baumgardt, H., & Kroupa, P. 2010, MNRAS, 402, 371
  * Banerjee et al. (2020) Banerjee, S., Belczynski, K., Fryer, C. L., Berczik, P., Hurley, J. R., Spurzem, R., & Wang, L. 2020, A&A, 639, A41
  * Bavera et al. (2021) Bavera, S. S., et al. 2021, A&A, 647, A153
  * Belczynski (2020) Belczynski, K. 2020, ApJ, 905, L15
  * Belczynski et al. (2010a) Belczynski, K., Bulik, T., Fryer, C. L., Ruiter, A., Valsecchi, F., Vink, J. S., & Hurley, J. R. 2010a, ApJ, 714, 1217
  * Belczynski et al. (2010b) Belczynski, K., Dominik, M., Bulik, T., O’Shaughnessy, R., Fryer, C. L., & Holz, D. E. 2010b, ApJ, 715, L138
  * Belczynski et al. (2016a) Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016a, Nature, 534, 512
  * Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407
  * Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., Taam, R. E., Zezas, A., Bulik, T., Maccarone, T. J., & Ivanova, N. 2008, ApJS, 174, 223
  * Belczynski et al. (2007) Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504
  * Belczynski et al. (2012) Belczynski, K., Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012, ApJ, 757, 91
  * Belczynski et al. (2016b) Belczynski, K., et al. 2016b, A&A, 594, A97
  * Belczynski et al. (2020) —. 2020, A&A, 636, A104
  * Belczynski et al. (2022) —. 2022, ApJ, 925, 69
  * Björklund et al. (2022) Björklund, R., Sundqvist, J. O., Singh, S. M., Puls, J., & Najarro, F. 2022, arXiv e-prints, arXiv:2203.08218
  * Broekgaarden et al. (2021) Broekgaarden, F. S., et al. 2021, arXiv e-prints, arXiv:2112.05763
  * Chatterjee et al. (2017) Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2017, ApJ, 834, 68
  * Chattopadhyay et al. (2022) Chattopadhyay, D., Hurley, J., Stevenson, S., & Raidani, A. 2022, arXiv e-prints, arXiv:2202.08924
  * Chen et al. (2021) Chen, H.-Y., Holz, D. E., Miller, J., Evans, M., Vitale, S., & Creighton, J. 2021, Classical and Quantum Gravity, 38, 055010
  * Chruslinska & Nelemans (2019) Chruslinska, M., & Nelemans, G. 2019, MNRAS, 1974
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833
  * Costa et al. (2021) Costa, G., Bressan, A., Mapelli, M., Marigo, P., Iorio, G., & Spera, M. 2021, MNRAS, 501, 4514
  * de Mink & Belczynski (2015) de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58
  * de Mink & Mandel (2016) de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545
  * Di Carlo et al. (2020) Di Carlo, U. N., et al. 2020, MNRAS, 498, 495
  * Dominik et al. (2012a) Dominik, M., Belczynski, K., Fryer, C., Holz, D., Berti, B., Bulik, T., Mandel, I., & O’Shaughnessy, R. 2012a, ApJ, 759, 52
  * Dominik et al. (2012b) Dominik, M., Belczynski, K., Fryer, C., Holz, D. E., Berti, E., Bulik, T., Mand el, I., & O’Shaughnessy, R. 2012b, ApJ, 759, 52
  * Edelman et al. (2022) Edelman, B., Doctor, Z., Godfrey, J., & Farr, B. 2022, ApJ, 924, 101
  * El-Badry et al. (2019) El-Badry, K., Quataert, E., Weisz, D. R., Choksi, N., & Boylan-Kolchin, M. 2019, Monthly Notices of the Royal Astronomical Society, 482, 4528
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36
  * Farrell et al. (2021) Farrell, E., Groh, J. H., Hirschi, R., Murphy, L., Kaiser, E., Ekström, S., Georgy, C., & Meynet, G. 2021, MNRAS, 502, L40
  * Fields et al. (2018) Fields, C. E., Timmes, F. X., Farmer, R., Petermann, I., Wolf, W. M., & Couch, S. M. 2018, ApJS, 234, 19
  * Fishbach et al. (2021) Fishbach, M., et al. 2021, Astrophys. J., 912, 98
  * Ford & McKernan (2021) Ford, K. E. S., & McKernan, B. 2021, arXiv e-prints, arXiv:2109.03212
  * Fragione & Banerjee (2020) Fragione, G., & Banerjee, S. 2020, ApJ, 901, L16
  * Fragione et al. (2021) Fragione, G., Kocsis, B., Rasio, F. A., & Silk, J. 2021, arXiv e-prints, arXiv:2107.04639
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., Dominik, M., Kalogera, V., & Holz, D. E. 2012, ApJ, 749, 91
  * Fuller & Ma (2019) Fuller, J., & Ma, L. 2019, ApJ, 881, L1
  * Gaburov et al. (2008) Gaburov, E., Lombardi, J. C., & Portegies Zwart, S. 2008, MNRAS, 383, L5
  * Gallegos-Garcia et al. (2021) Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, ApJ, 922, 110
  * García et al. (2021) García, F., Simaz Bunzel, A., Chaty, S., Porter, E., & Chassande-Mottin, E. 2021, A&A, 649, A114
  * Gerosa & Kesden (2016) Gerosa, D., & Kesden, M. 2016, PRECESSION: Python toolbox for dynamics of spinning black-hole binaries, Astrophysics Source Code Library, record ascl:1611.004
  * Gessner & Janka (2018) Gessner, A., & Janka, H.-T. 2018, ApJ, 865, 61
  * Giacobbo & Mapelli (2018) Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011
  * Heggie & Hut (2003) Heggie, D., & Hut, P. 2003, Classical and Quantum Gravity, 20, 4504
  * Hénon (1971) Hénon, M. 1971, Astrophysics and Space Science, 13, 284
  * Hénon (1971) Hénon, M. H. 1971, Ap&SS, 14, 151
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
  * J. C. Lombardi et al. (2002) J. C. Lombardi, J., Warren, J. S., Rasio, F. A., Sills, A., & Warren, A. R. 2002, The Astrophysical Journal, 568, 939
  * Joshi et al. (2000) Joshi, K. J., Rasio, F. A., & Portegies Zwart, S. P. 2000, The Astrophysical Journal, 540, 969
  * King et al. (2001) King, A. R., Davies, M. B., Ward, M. J., Fabbiano, G., & Elvis, M. 2001, ApJ, 552, L109
  * Kinugawa et al. (2020) Kinugawa, T., Nakamura, T., & Nakano, H. 2020, arXiv e-prints, arXiv:2009.06922
  * Klencki et al. (2018) Klencki, J., Moe, M., Gladysz, W., Chruslinska, M., Holz, D. E., & Belczynski, K. 2018, Astron. Astrophys., 619, A77
  * Klencki et al. (2021) Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, A&A, 645, A54
  * Kremer et al. (2020) Kremer, K., et al. 2020, ApJ, 903, 45
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Leung et al. (2019) Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, ApJ, 887, 72
  * Lipunov et al. (1997) Lipunov, V. M., Postnov, K. A., & Prokhorov, M. E. 1997, Astronomy Letters, 23, 492
  * MacLeod et al. (2017) MacLeod, M., Antoni, A., Murguia-Berthier, A., Macias, P., & Ramirez-Ruiz, E. 2017, ApJ, 838, 56
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, ApJ, 840, 39
  * Mandel & Broekgaarden (2021) Mandel, I., & Broekgaarden, F. S. 2021, arXiv e-prints, arXiv:2107.14239
  * Mandel & de Mink (2016) Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634
  * Mapelli (2016) Mapelli, M. 2016, MNRAS, 459, 3432
  * Mapelli et al. (2021) Mapelli, M., et al. 2021, MNRAS, 505, 339
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50
  * Mikkola & Merritt (2008) Mikkola, S., & Merritt, D. 2008, AJ, 135, 2398
  * Mikkola & Tanikawa (1999) Mikkola, S., & Tanikawa, K. 1999, Monthly Notices of the Royal Astronomical Society, 310, 745
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15
  * Mondal et al. (2020) Mondal, S., Belczyński, K., Wiktorowicz, G., Lasota, J.-P., & King, A. R. 2020, MNRAS, 491, 2747
  * Neijssel et al. (2019a) Neijssel, C. J., et al. 2019a, MNRAS, 490, 3740
  * Neijssel et al. (2019b) —. 2019b, MNRAS, 490, 3740
  * Nitadori & Aarseth (2012) Nitadori, K., & Aarseth, S. J. 2012, Monthly Notices of the Royal Astronomical Society, 424, 545
  * O’Connor & Ott (2011) O’Connor, E., & Ott, C. D. 2011, ApJ, 730, 70
  * Olejak & Belczynski (2021) Olejak, A., & Belczynski, K. 2021, ApJ, 921, L2
  * Olejak et al. (2021) Olejak, A., Belczynski, K., & Ivanova, N. 2021, A&A, 651, A100
  * Olejak et al. (2022) Olejak, A., Fryer, C. L., Belczynski, K., & Baibhav, V. 2022, arXiv e-prints, arXiv:2204.09061
  * Pattabiraman et al. (2013) Pattabiraman, B., Umbreit, S., Liao, W. K., Choudhary, A., Kalogera, V., Memik, G., & Rasio, F. A. 2013, The Astrophysical Journal, Supplement Series, 204, 16
  * Pavlovskii et al. (2017) Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
  * Podsiadlowski et al. (2004) Podsiadlowski, P., Langer, N., Poelarends, A. J. T., Rappaport, S., Heger, A., & Pfahl, E. 2004, ApJ, 612, 1044
  * Podsiadlowski et al. (2003) Podsiadlowski, P., Rappaport, S., & Han, Z. 2003, MNRAS, 341, 385
  * Portegies Zwart & McMillan (2000) Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17
  * Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431
  * Qin et al. (2018) Qin, Y., Fragos, T., Meynet, G., Andrews, J., Sørensen, M., & Song, H. F. 2018, A&A, 616, A28
  * Rastello et al. (2019) Rastello, S., Amaro-Seoane, P., Arca-Sedda, M., Capuzzo-Dolcetta, R., Fragione, G., & Tosta e Melo, I. 2019, MNRAS, 483, 1233
  * Rastello et al. (2021) Rastello, S., Mapelli, M., Di Carlo, U. N., Iorio, G., Ballone, A., Giacobbo, N., Santoliquido, F., & Torniamenti, S. 2021, MNRAS, 507, 3612
  * Rodriguez et al. (2018a) Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., Kremer, K., Rasio, F. A., Samsing, J., Ye, C. S., & Zevin, M. 2018a, Phys. Rev. D, 98, 123005
  * Rodriguez et al. (2018b) Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., & Rasio, F. A. 2018b, Phys. Rev. Lett., 120, 151101
  * Rodriguez et al. (2016a) Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016a, Phys. Rev. D, 93, 084029
  * Rodriguez et al. (2016b) Rodriguez, C. L., Haster, C.-J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016b, ApJ, 824, L8
  * Rodriguez & Loeb (2018) Rodriguez, C. L., & Loeb, A. 2018, The Astrophysical Journal Letters, 866, L5
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., Chatterjee, S., Kremer, K., Rasio, F. A., & Ye, C. S. 2019, Phys. Rev. D, 100, 043027
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., Chatterjee, S., Kremer, K., Rasio, F. A., & Ye, C. S. 2019, Physical Review D, 100, 43027
  * Roulet et al. (2021) Roulet, J., Chia, H. S., Olsen, S., Dai, L., Venumadhav, T., Zackay, B., & Zaldarriaga, M. 2021, Phys. Rev. D, 104, 083010
  * Safarzadeh et al. (2020) Safarzadeh, M., Hamers, A. S., Loeb, A., & Berger, E. 2020, ApJ, 888, L3
  * Samsing et al. (2020) Samsing, J., et al. 2020, arXiv e-prints, arXiv:2010.09765
  * Sana & Evans (2011) Sana, H., & Evans, C. J. 2011, in IAU Symposium, Vol. 272, Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits, ed. C. Neiner, G. Wade, G. Meynet, & G. Peters, 474–485
  * Sana et al. (2012) Sana, H., et al. 2012, Science, 337, 444
  * Sana et al. (2013) —. 2013, A&A, 550, A107
  * Santoliquido et al. (2021) Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021, MNRAS, 502, 4877
  * Schneider et al. (2019) Schneider, A. S., Roberts, L. F., Ott, C. D., & O’Connor, E. 2019, Phys. Rev. C, 100, 055802
  * Sharda et al. (2021) Sharda, P., Krumholz, M. R., Wisnioski, E., Acharyya, A., Federrath, C., & Forbes, J. C. 2021, MNRAS, 504, 53
  * Spera et al. (2019) Spera, M., Mapelli, M., Giacobbo, N., Trani, A. A., Bressan, A., & Costa, G. 2019, MNRAS, 485, 889
  * Spitzer (1987) Spitzer, L. 1987, Dynamical evolution of globular clusters, Princeton University Press, Princeton, NJ, 191 p.
  * Stevenson et al. (2019) Stevenson, S., Sampson, M., Powell, J., Vigna-Gómez, A., Neijssel, C. J., Szécsi, D., & Mandel, I. 2019, arXiv e-prints
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration et al. 2021a, arXiv e-prints, arXiv:2111.03606
  * The LIGO Scientific Collaboration et al. (2021b) —. 2021b, ApJ, 913, L7
  * The LIGO Scientific Collaboration et al. (2021c) —. 2021c, arXiv e-prints, arXiv:2111.03634
  * Timmes et al. (1996) Timmes, F. X., Woosley, S. E., & Weaver, T. A. 1996, ApJ, 457, 834
  * Tiwari (2022) Tiwari, V. 2022, ApJ, 928, 155
  * Tiwari & Fairhurst (2021) Tiwari, V., & Fairhurst, S. 2021, ApJ, 913, L19
  * van den Heuvel et al. (2017) van den Heuvel, E. P. J., Portegies Zwart, S. F., & de Mink, S. E. 2017, MNRAS, 471, 4256
  * van Son et al. (2021) van Son, L. A. C., et al. 2021, arXiv e-prints, arXiv:2110.01634
  * van Son et al. (2022) —. 2022, ApJ, 931, 17
  * Vink & de Koter (2005) Vink, J. S., & de Koter, A. 2005, A&A, 442, 587
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Voss & Tauris (2003) Voss, R., & Tauris, T. M. 2003, MNRAS, 342, 1169
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
  * Zampieri & Roberts (2009) Zampieri, L., & Roberts, T. P. 2009, MNRAS, 400, 677
  * Zevin (2020) Zevin, M. 2020, Constraining the Origins of Binary Black Holes using Multiple Formation Pathways
  * Zevin et al. (2020) Zevin, M., Spera, M., Berry, C. P. L., & Kalogera, V. 2020, ApJ, 899, L1
  * Zevin et al. (2021) Zevin, M., et al. 2021, The Astrophysical Journal, 910, 152
  * Zhang (2013) Zhang, Q. S. 2013, ApJS, 205, 18

##  5\. APPENDIX A

All four binary models show a general increase in the average intrinsic total
BH-BH merger mass with redshift (see Fig. 1). However, there is an exception
from this trend for the CE BH-BH formation channel at low redshifts
($z\lesssim 2-3$). Here, we provide an explanation for this behavior.

For this channel, BH-BH merger formation depends sensitively on radial
expansion of BH progenitor stars and thus on metallicity (Belczynski et al.,
2010b). For our CE model, we only allow for CEs to develop and successfully
form close binaries for late donor evolutionary stages (e.g., during core He-
burning). Significant radial expansion during core He-burning is found at low
metallicity, as high metallicity H-rich stellar envelopes are usually severely
depleted by stellar winds for massive stars. As metallicity of stellar
populations increases with decreasing redshift, the CE channel BH-BH merger
formation efficiency significantly drops as seen in Figure 3. For low
redshifts ($z\lesssim 2-3$) the average metallicity of stellar populations
becomes significant ($Z>0.5{\rm~{}Z}_{\odot}$, see Fig.7 of Belczynski et al.
(2020)) and the formation efficiency of BH-BH mergers is very small. The low
mass BH-BH mergers that form from high metallicity stars do not contribute
significantly to overall population of BH-BH mergers at these redshifts. The
majority of mergers at these redshifts originate from low-metallicity stars
that still form in the tail of metallicity distribution even to $z=0$. Massive
mergers that form at low metallicity at high redshifts do not contribute
significantly at low redshifts as merger rates fall off steeply with delay
time (Dominik et al., 2012a). This sharp drop in BH-BH merger formation
efficiency with metallicity (Fig. 3) is responsible for the flat behavior of
the average BH-BH merger mass at low redshifts. This does not affect RLOF
channels as they depend on expansion of stars at much earlier evolutionary
stages (e.g., during the Hertzsprung gap) which are not affected as much by
stellar winds and metallicity.

Figure 3.—  BH-BH merger formation efficiency per unit mass as a function of
metallicity for the CE and three RLOF isolated binary evolution channels.
Figure 4.—  Top: average birth redshift of binary progenitors of BH-BH mergers
that merge at a given redshift for four isolated binary formation channels.
Note that black and cyan lines are right under blue line. Bottom: average
metallicity of binaries that produce BH-BH mergers at a given redshift.

##  6\. APPENDIX B

Here, we explain the difference in average total mass of BH-BH mergers between
the CE formation channel and RLOF channels. It is clearly visible that RLOF
channels produce much more massive BH-BH mergers at all redshifts (see Fig.
1).

In Figures 5 and 6 we present formation scenarios for typical BH-BH mergers at
redshift $z=6$ in the CE and stable RLOF (standard IMF, delayed SN engine)
isolated binary evolution channels, respectively. The two chosen examples of
BH-BH systems have their total mass and metallicity close to the average mass
$\overline{M}_{\rm tot,int}=29.3,\ 47.7{\rm~{}M}_{\odot}$ (CE, RLOF) and the
average metallicity $\overline{Z}=0.0013,\ 0.0012$ (CE, RLOF) at $z=6$ for the
two channels (see Fig. 4). For the CE channel we show example with $M_{\rm
tot,int}=31{\rm~{}M}_{\odot}$, and for the standard RLOF channel with $M_{\rm
tot,int}=47{\rm~{}M}_{\odot}$. Note the very clear total mass difference of
BH-BH mergers from the two formation scenarios. Below we explain the origin of
this difference.

The two examples, despite having rather different primary stars ($\sim
38{\rm~{}M}_{\odot}$ and $\sim 65{\rm~{}M}_{\odot}$ at ZAMS: $t=0$ Myr),
produce the similar mass BHs ($\sim 14{\rm~{}M}_{\odot}$) out of the primary
stars. This comes from the fact that both primaries are stripped of their
H-rich envelopes (first RLOF phase in each case) at different evolutionary
stages and their exposed helium cores, Wolf-Rayet (WR) stars, have almost the
same mass. In the CE scenario, the first RLOF begins when the primary is on
Hertzsprung gap and its helium core is developed after the entirety of main-
sequence H-burning. In the RLOF scenario, the first RLOF begins when the
primary is still on the main sequence, and due to mass loss the star is de-
rejuvenetad, producing a much smaller core than would be expected for a single
star with the same initial mass. Since both first-formed BHs have the same
mass, it means that mass difference between the two BH-BH mergers is directly
connected to the secondary star, its mass, and its evolution.

In the CE channel, the secondary mass is $34{\rm~{}M}_{\odot}$ at $t=0$ Myr.
It increases to $45{\rm~{}M}_{\odot}$ after first stable RLOF event. Then the
secondary loses mass in stellar winds to initiate the second RLOF at a mass of
$27{\rm~{}M}_{\odot}$ when it is core He-burning. This second RLOF is a CE
event that leaves the secondary as a naked helium core (WR star) with a mass
of $17{\rm~{}M}_{\odot}$, which quickly collapses to a BH with similar mass.
Note the significant orbital separation decay during the CE phase, from
$a=482{\rm~{}R}_{\odot}$ to $a=9.5{\rm~{}R}_{\odot}$. If we made the binary at
this stage instead proceed through stable RLOF, then there would be not enough
orbital contraction to produce a merging BH-BH system. Simply, the secondary
star would not have enough mass in its envelope to allow for enough associated
angular momentum loss from the binary to provide significant orbital decay.
The orbital decay would be rather small: from $a=482{\rm~{}R}_{\odot}$ to
$a=440{\rm~{}R}_{\odot}$, producing a BH-BH system with an inspiral time much
larger than the Hubble time.

In the RLOF channel, the secondary starts with a mass of
$64{\rm~{}M}_{\odot}$. It increases to $87{\rm~{}M}_{\odot}$ after the first
stable RLOF event. Then the secondary, with almost no mass loss:
$85{\rm~{}M}_{\odot}$ (mostly main-sequence evolution with a relatively small
stellar radius), initiates the second RLOF event right after the main sequence
when it is on the Hertzsprung gap. The second RLOF is stable thermal-timescale
mass transfer. The RLOF strips the secondary of its envelope, leaving a WR
star with a mass of $38{\rm~{}M}_{\odot}$. After losing about
$5{\rm~{}M}_{\odot}$ to winds, the WR star collapses to a massive BH of
$33{\rm~{}M}_{\odot}$. This second RLOF event reduces orbital separation from
$a=295{\rm~{}R}_{\odot}$ to $a=6.2{\rm~{}R}_{\odot}$, producing a close BH-BH
system that can merge in the Hubble time. The high mass of the donor star
(secondary) in this event is required to allow for significant decay of the
orbit. During thermal-timescale mass transfer, almost the entirety of the
donor’s envelope is lost from the binary (highly over-Eddington mass transfer
rate), taking away binary angular momentum and reducing the size of the orbit
(as long as donor mass is larger than accretor mass). In order to form a BH-BH
binary system close enough to merge in the Hubble time through two stable RLOF
episodes, the donor star needs to be few times more massive than its BH
companion at the onset of the second RLOF event (van den Heuvel et al., 2017;
Olejak et al., 2021). Only such unequal mass binaries ($q_{\rm
RLOF}=\frac{M_{\rm don}}{{M_{\rm acc}}}\gtrsim 3$) are able to eject a large
enough amount of mass and orbital angular momentum to reduce the orbital
separation enough to produce merging BH-BH systems.

If at onset of the second RLOF event we applied CE evolution instead of stable
mass transfer, then the ejection of very massive envelope of the secondary
($47{\rm~{}M}_{\odot}$) would lead to post-CE separation of only
$a=0.7{\rm~{}R}_{\odot}$. This is smaller than the radius of a secondary core
(a WR star: $>1-2{\rm~{}R}_{\odot}$; Hurley et al. (2000)) and thus CE
evolution would lead to a merger of the BH with the massive secondary and no
formation of a BH-BH system. Note that in our standard simulations we would
not even attempt to perform CE evolution for such a case, as the donor
(secondary) is on Hertzsprung gap at the onset of RLOF and has a radiative
envelope that would not allow for the development/survival of a CE (Belczynski
et al., 2007; Klencki et al., 2021).

Of course, to alleviate this issue one may require the formation of a much
wider binary at the onset of the second RLOF event in this case. Then, not
only may there be enough orbital energy to successfully expel the CE, but also
the donor star would be core He-burning and much more extended, as well as
having a convective envelope. However, such star has already evolved through
most its life and lost a significant fraction of its H-rich envelope through
stellar winds. In such a case, even very efficient CE evolution may not be
enough to decay the orbit below the threshold of formation of merging BH-BH
systems. It is possible, but it requires fine tuning, and for a given model of
CE efficiency and adopted stellar winds only few systems form massive BH-BH
mergers in the CE channel. This was already proposed and explained, although
in somewhat different terms, by van Son et al. (2021), and our detailed
description is provided only to serve as a replication study.

Figure 5.—  The formation scenario for a typical BH-BH merger at redshift
$z\approx 6$ in the classical (CE) isolated binary evolution channel. The
shaded part of the evolution shows the most important part of evolution that
leads to orbital decay, allowing for the formation of close BH-BH binary that
can merge in the Hubble time.  Figure 6.—  The formation scenario for a
typical BH-BH merger at redshift $z\approx 6$ in the stable RLOF channel
(standard IMF, delayed SN engine). The shaded part of the evolution shows the
most important part of evolution that leads to orbital decay, allowing for the
formation of close BH-BH binary that can merge in the Hubble time.


Nlc}(X,\Delta)}$, which we think of an $\mathbb{R}$-line bundle on ${\rm
Nlc}(X,\Delta)$, is semi-ample over $Z$. Let

$(X,\Delta+A)=:(X_{1},\Delta_{1}+A_{1})\dashrightarrow\cdots\dashrightarrow(X_{i},\Delta_{i}+A_{i})\dashrightarrow\cdots$

be a sequence of steps of a $(K_{X}+\Delta+A)$-MMP over $Z$ with scaling of
$A$ such that if we put

$\lambda_{i}:=\\{\mu\in\mathbb{R}_{\geq
0}\,|\,\text{$(K_{X_{i}}+\Delta_{i}+A_{i})+\mu A_{i}$ is nef over $Z$}\\}$

for each $i\geq 1$, then ${\rm lim}_{i\to\infty}\lambda_{i}=0$. Then
$\lambda_{m}=0$ for some $m$ and $(X_{m},\Delta_{m}+A_{m})$ is a good minimal
model of $(X,\Delta+A)$ over $Z$.

###### Proof.

The proof of Theorem 5.3 works without any change. ∎

###  5.2. Proofs of corollaries

In this subsection we prove corollaries.

######  Corollary 5.5.

Let $\pi\colon X\to Z$ be a projective morphism of normal quasi-projective
varieties. Let $(X,\Delta)$ be a normal pair and let $A$ be a $\pi$-ample
$\mathbb{R}$-divisor on $X$ such that $K_{X}+\Delta+A$ is $\pi$-pseudo-
effective. Suppose that ${\rm NNef}(K_{X}+\Delta+A/Z)\cap{\rm
Nlc}(X,\Delta)=\emptyset$. Suppose in addition that $(K_{X}+\Delta+A)|_{{\rm
Nlc}(X,\Delta)}$, which we think of an $\mathbb{R}$-line bundle on ${\rm
Nlc}(X,\Delta)$, is semi-ample over $Z$. Then ${\rm
Bs}|K_{X}+\Delta+A/Z|_{\mathbb{R}}\cap{\rm Nlc}(X,\Delta)=\emptyset$.

###### Proof.

This immediately follows from Theorem 5.3. ∎

######  Corollary 5.6.

Let $\pi\colon X\to Z$ be a projective morphism of normal quasi-projective
varieties. Let $(X,\Delta)$ be a normal pair such that $K_{X}+\Delta$ is
$\pi$-pseudo-effective. Suppose that ${\rm NNef}(K_{X}+\Delta/Z)\cap{\rm
Nlc}(X,\Delta)=\emptyset$. Let $A$ be an ample $\mathbb{R}$-divisor on $X$.
Then there exists a sequence of steps of a $(K_{X}+\Delta)$-MMP over $Z$ with
scaling of $A$

|
---|---
$\textstyle{(X,\Delta)=:(X_{1},\Delta_{1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{i},\Delta_{i})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\varphi_{i}}$$\textstyle{(X_{i+1},\Delta_{i+1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{\varphi^{\prime}_{i}}$$\textstyle{\cdots,}$$\textstyle{W_{i}}$

where $(X_{i},\Delta_{i})\to W_{i}\leftarrow(X_{i+1},\Delta_{i+1})$ is a step
of the $(K_{X}+\Delta)$-MMP over $Z$, such that

  * •

the non-isomorphic locus of the MMP is disjoint from ${\rm Nlc}(X,\Delta)$,

  * •

$\rho(X_{i}/W_{i})=1$ for every $\geq 1$, and

  * •

if we put

$\lambda_{i}:={\rm inf}\\{\mu\in\mathbb{R}_{\geq
0}\,|\,\text{$K_{X_{i}}+\Delta_{i}+\mu A_{i}$ is nef over $Z$}\\}$

for each $i\geq 1$, then ${\rm lim}_{i\to\infty}\lambda_{i}=0$.

Moreover, if $X$ is $\mathbb{Q}$-factorial, then all $X_{i}$ in the MMP are
also $\mathbb{Q}$-factorial.

###### Proof.

By Corollary 3.12, there exists a sequence of steps of a $(K_{X}+\Delta)$-MMP
over $Z$ with scaling of $A$

|
---|---
$\textstyle{(X,\Delta)=:(X_{1},\Delta_{1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{i},\Delta_{i})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{i+1},\Delta_{i+1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots,}$$\textstyle{W_{i}}$

where $(X_{i},\Delta_{i})\to W_{i}\leftarrow(X_{i+1},\Delta_{i+1})$ is a step
of a $(K_{X_{i}}+\Delta_{i})$-MMP over $Z$, satisfying the following.

  * •

the non-isomorphic locus of the MMP is disjoint from ${\rm Nlc}(X,\Delta)$,

  * •

$\rho(X_{i}/W_{i})=1$ for every $\geq 1$,

  * •

if we put

$\lambda_{i}:={\rm inf}\\{\mu\in\mathbb{R}_{\geq
0}\,|\,\text{$K_{X_{i}}+\Delta_{i}+\mu A_{i}$ is nef over $Z$}\\}$

for each $i\geq 1$ and $\lambda:={\rm lim}_{i\to\infty}\lambda_{i}$, then the
MMP terminates after finitely many steps or otherwise $\lambda\neq\lambda_{i}$
for every $i\geq 1$, and

  * •

if $X$ is $\mathbb{Q}$-factorial, then all $X_{i}$ in the MMP are also
$\mathbb{Q}$-factorial.

Then the MMP is a $(K_{X}+\Delta+\lambda A)$-over $Z$ with scaling of $A$. If
$\lambda>0$, then we have ${\rm NNef}(K_{X}+\Delta+\lambda A/Z)\cap{\rm
Nlc}(X,\Delta)=\emptyset$, and $(K_{X}+\Delta+\lambda A)|_{{\rm
Nlc}(X,\Delta)}$ is ample over $Z$ because $(K_{X}+\Delta)|_{{\rm
Nlc}(X,\Delta)}$ is nef over $Z$ and we may use [KM98, Proof of Theorem 1.38]
and Nakai–Moishezon’s criterion for the ampleness. This contradicts Theorem
5.4, and therefore $\lambda=0$. This MMP satisfies the conditions of Corollary
5.6. ∎

######  Corollary 5.7 (cf. [G11]).

Let $(X,\Delta)$ be a projective normal pair such that $K_{X}+\Delta$ is
pseudo-effective and the non-nef locus of $K_{X}+\Delta$ is disjoint from
${\rm Nlc}(X,\Delta)$. Suppose that the numerical dimension
$\kappa_{\sigma}(K_{X}+\Delta)$ (see [N04, V, 2.5. Definition]) is zero. Put
$X_{1}=X$ and $\Delta_{1}=\Delta$. Then there exists a sequence of steps of a
$(K_{X}+\Delta)$-MMP

|  |
---|---|---
$\textstyle{(X_{1},\Delta_{1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{i},\Delta_{i})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{i+1},\Delta_{i+1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{m},\Delta_{m}),}$$\textstyle{W_{i}}$

where $(X_{i},\Delta_{i})\to W_{i}\leftarrow(X_{i+1},\Delta_{i+1})$ is a step
of the $(K_{X}+\Delta)$-MMP, such that

  * •

the non-isomorphic locus of the MMP is disjoint from ${\rm Nlc}(X,\Delta)$,

  * •

$\rho(X_{i}/W_{i})=1$ for every $\geq 1$, and

  * •

$K_{X_{m}}+\Delta_{m}\equiv 0$.

Moreover, if $X$ is $\mathbb{Q}$-factorial, then all $X_{i}$ in the MMP are
also $\mathbb{Q}$-factorial.

###### Proof.

By Corollary 5.6, there exists a sequence of steps of a $(K_{X}+\Delta)$-MMP
with scaling of an ample $\mathbb{R}$-divisor $A$

|  |
---|---|---
$\textstyle{(X_{1},\Delta_{1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{i},\Delta_{i})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{i+1},\Delta_{i+1})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{\cdots\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\textstyle{(X_{m},\Delta_{m}),}$$\textstyle{W_{i}}$

where $(X_{i},\Delta_{i})\to W_{i}\leftarrow(X_{i+1},\Delta_{i+1})$ is a step
of the $(K_{X}+\Delta)$-MMP over $Z$, such that

  * •

the non-isomorphic locus of the MMP is disjoint from ${\rm Nlc}(X,\Delta)$,

  * •

$\rho(X_{i}/W_{i})=1$ for every $\geq 1$, and

  * •

if we put

$\lambda_{i}:={\rm inf}\\{\mu\in\mathbb{R}_{\geq
0}\,|\,\text{$K_{X_{i}}+\Delta_{i}+\mu A_{i}$ is nef over $Z$}\\}$

for each $i\geq 1$, then ${\rm lim}_{i\to\infty}\lambda_{i}=0$.

Moreover, if $X$ is $\mathbb{Q}$-factorial, then all $X_{i}$ in the MMP are
$\mathbb{Q}$-factorial. Prime divisors contracted by the MMP are components of
the negative part of the Nakayama–Zariski decomposition of $K_{X}+\Delta$.
Hence, there exists $m$ such that $K_{X_{m}}+\Delta_{m}$ is the limit of
movable divisors. By the same argument as in [G11, Proof of Theorem 5.1], we
see that $K_{X_{m}}+\Delta_{m}\equiv 0$. The MMP satisfies the properties of
Corollary 5.7. ∎

######  Corollary 5.8 (cf. [BBP13], [TX23]).

Let $(X,\Delta)$ be a projective normal pair such that $K_{X}+\Delta$ is
pseudo-effective and the non-nef locus of $K_{X}+\Delta$ is disjoint from
${\rm Nlc}(X,\Delta)$. Then

${\rm NNef}(K_{X}+\Delta)=\boldsymbol{\rm B}_{-}(K_{X}+\Delta),$

where $\boldsymbol{\rm B}_{-}(K_{X}+\Delta)$ is the diminished base locus
(also called the restricted base locus) of $K_{X}+\Delta$ as in [BBP13] and
[TX23]. Furthermore, every irreducible component of the non-nef locus of
$K_{X}+\Delta$ is uniruled.

###### Proof.

The argument in [BBP13] or [TX23] works with no changes because we may use
Corollary 5.6 and [F21, Theorem 1.12]. ∎

######  Corollary 5.9.

Let $\pi\colon X\to Z$ be a projective morphism of normal quasi-projective
varieties. Let $(X,\Delta)$ be a normal pair such that $\Delta$ is a
$\mathbb{Q}$-divisor on $X$. Let $A$ be a $\pi$-ample $\mathbb{Q}$-divisor on
$X$ such that $K_{X}+\Delta+A$ is $\pi$-pseudo-effective. Suppose that ${\rm
NNef}(K_{X}+\Delta+A/Z)\cap{\rm Nlc}(X,\Delta)=\emptyset$ and that
$(K_{X}+\Delta+A)|_{{\rm Nlc}(X,\Delta)}$, which we think of an
$\mathbb{Q}$-line bundle on ${\rm Nlc}(X,\Delta)$, is semi-ample over $Z$.
Then the sheaf of graded $\pi_{*}\mathcal{O}_{X}$-algebra

$\underset{m\in\mathbb{Z}_{\geq 0}}{\bigoplus}\pi_{*}\mathcal{O}_{X}(\lfloor
m(K_{X}+\Delta+A)\rfloor)$

is finitely generated.

###### Proof.

This immediately follows from Theorem 5.3. ∎

##  6\. Examples

In this section we collect some examples.

The first example shows that a step of an MMP for normal pairs does not always
exist and the abundance conjecture for normal pairs does not hold in general.

######  Example 6.1.

Let $E$ be an elliptic curve with a very ample divisor $H$ and define

$X:=\mathbb{P}_{E}(\mathcal{O}_{E}\oplus\mathcal{O}_{E}(-H))\overset{f}{\longrightarrow}E.$

Let $S$ be the unique section of $\mathcal{O}_{X}(1)$. Then
$K_{X}+S+(S+f^{*}H)\sim 0$ and $S+f^{*}H$ is base point free. The contraction

$\pi\colon X\to Z$

induced by $S+f^{*}H$ is birational and we have ${\rm Ex}(\pi)=S$. Thus $-S$
is $\pi$-ample. Let $D$ be a non-torsion Cartier divisor on $X$ of degree
zero. Since $D+S+f^{*}H$ is nef and big, we can write
$D+S+f^{*}H\sim_{\mathbb{Q}}A+B$ such that $A$ is an effective ample
$\mathbb{Q}$-divisor and $B$ is an effective $\mathbb{Q}$-divisor on $X$. Now
we have

$K_{X}+S+B+A\sim_{\mathbb{Q}}f^{*}D$

and $D$ is not semi-ample.

  * •

With notation as above, we consider $\pi\colon X\to Z$ and the pair
$(X,2S+B+A)$. Since $-S$ is $\pi$-ample, $\pi$ is a $(K_{X}+2S+B+A)$-negative
extremal contraction. However, the divisor
$\pi_{*}(K_{X}+2S+B+A)\sim_{\mathbb{Q}}\pi_{*}f^{*}D$ is not
$\mathbb{R}$-Cartier since $rH$ is not $\mathbb{R}$-linearly equivalent to $D$
for any $r\in\mathbb{R}$ (cf. [F17, Proposition 7.2.8]). This example implies
that we cannot always construct a step of an MMP even in the case of surfaces.

  * •

With notation as above, consider the pair $(X,S+B+A)$. By construction,
$K_{X}+S+B+A$ is nef but not semi-ample. This pair shows that the abundance
conjecture does not always hold for normal pairs with polarizations.

The next example shows that the non-nef locus of the log canonical divisor of
a pair is not necessarily Zariski closed. In particular, minimal models of
normal pairs does not always exist.

######  Example 6.2.

For a blow up $V\to\mathbb{P}^{3}$ at nine very general points and

$X:=\mathbb{P}_{V}(\mathcal{O}_{V}\oplus\mathcal{O}_{V}(1))\overset{f}{\longrightarrow}V,$

where $\mathcal{O}_{V}(1)$ is very ample line bundle on $V$, Lesieutre [L14,
Theorem 1.1] constructed a big $\mathbb{R}$-divisor $D$ on $X$ such that the
diminished base locus

$\boldsymbol{\rm B}_{-}(D):=\bigcup_{\begin{subarray}{c}\text{$A$: ample}\\\
\text{$D+A$: $\mathbb{Q}$-Cartier}\end{subarray}}{\rm Bs}|D+A|_{\mathbb{Q}}$

is a countable union of curves. By construction, there is an effective
$\mathbb{R}$-divisor $\Delta_{V}$ on $V$ such that
$K_{V}+\Delta_{V}\sim_{\mathbb{R}}0$, and therefore
$K_{X}+\Delta\sim_{\mathbb{R}}0$ for some effective $\mathbb{R}$-divisor
$\Delta$ on $X$. We can write $D\sim_{\mathbb{R}}A+B$ such that $A$ is an
effective ample $\mathbb{R}$-divisor and $B$ is an effective
$\mathbb{R}$-divisor on $X$. Then

$K_{X}+\Delta+B+A\sim_{\mathbb{R}}D.$

Moreover, by [N04, V, 1.3 Theorem], it follows that $\boldsymbol{\rm
B}_{-}(D)$ coincides with the non-nef locus of $D$. Therefore, the pair
$(X,\Delta+B+A)$ satisfies the condition that the non-nef locus of
$K_{X}+\Delta+B+A$ is not Zariski closed. If $(X,\Delta+B+A)$ has a minimal
model $(X^{\prime},\Gamma^{\prime})$, taking a common resolution $g\colon Y\to
X$ and $g^{\prime}\colon Y\to X^{\prime}$ of $X\dashrightarrow X^{\prime}$,
then we can write

$g^{*}(K_{X}+\Delta+B+A)=g^{\prime*}(K_{X^{\prime}}+\Gamma^{\prime})+E$

for some effective $g^{\prime}$-exceptional $\mathbb{R}$-divisor $E$ on $Y$.
Since $K_{X^{\prime}}+\Gamma^{\prime}$ is nef, we have

${\rm NNef}(g^{*}(K_{X}+\Delta+B+A))={\rm
NNef}(g^{\prime*}(K_{X^{\prime}}+\Gamma^{\prime})+E)={\rm Supp}\,E.$

Then ${\rm NNef}(K_{X}+\Delta+B+A)=g({\rm
NNef}(g^{*}(K_{X}+\Delta+B+A)))=g({\rm Supp}\,E)$ is Zariski closed, a
contradiction. Therefore, $(X,\Delta+B+A)$ does not have a minimal model.

## References

  * [A04] F. Ambro, Shokurov’s boundary property, J. Differential Geom. 67 (2004), no. 2, 229–255.
  * [AK17] F. Ambro, J. Kollár, Minimal models of semi-log-canonical pairs. In Moduli of K-stable varieties, Springer INdAM Ser., 31, Springer, Cham (2019), 1–13.
  * [B12] C. Birkar, Existence of log canonical flips and a special LMMP, Publ. Math. Inst. Hautes Études Sci. 115 (2012), no. 1, 325–368.
  * [BCHM10] C. Birkar, P. Cascini, C. D. Hacon, J. McKernan, Existence of minimal models for varieties of log general type, J. Amer. Math. Soc. 23 (2010), no. 2, 405–468.
  * [BZ16] C. Birkar, D. Q. Zhang, Effectivity of Iitaka fibrations and pluricanonical systems of polarized pairs, Publ. Math. Inst. Hautes Études Sci. 123 (2016), no. 1, 283–331.
  * [BBP13] S. Boucksom, A. Broustet, G. Pacienza, Uniruledness of stable base loci of adjoint linear systems via Mori theory, Math. Z. 275 (2013), no. 1–2, 499–507.
  * [F07a] O. Fujino, What is log terminal? In Flips for $3$-folds and $4$-folds, Oxford University Press (2007).
  * [F07b] O. Fujino, Special termination and reduction to pl flips. In Flips for $3$-folds and $4$-folds, Oxford University Press (2007).
  * [F11] O. Fujino, Fundamental theorems for the log minimal model program, Publ. Res. Inst. Math. Sci. 47 (2011), no. 3, 727–789.
  * [F14] O. Fujino, Fundamental theorems for semi log canonical pairs, Algebraic Geom. 1 (2014), no. 2, 194–228.
  * [F17] O. Fujino, Foundations of the minimal model program, MSJ Mem. 35, Mathematical Society of Japan, Tokyo, 2017.
  * [F21] O. Fujino, Cone theorem and Mori hyperbolicity. preprint (2021), arXiv:2102.11986v3.
  * [FH23] O. Fujino, K. Hashizume, Adjunction and inversion of adjunction, Nagoya Math. J. 249 (2023), 119–147.
  * [FST11] O. Fujino, K. Schwede, S. Takagi, Supplements to non-lc ideal sheaves. RIMS Kokyuroku Bessatsu, B24, Res. Inst. Math. Sci. (RIMS), Kyoto, 2011, 1–46.
  * [G11] Y. Gongyo, On the minimal model theory for dlt pairs of numerical log kodaira dimension zero, Math. Res. Lett., 18 (2011) , no. 5, 991–1000.
  * [HL21] C. D. Hacon, J. Liu, Existence of flips for generalized lc pairs, to appear in Camb. J. Math.
  * [HX13] C. D. Hacon, C. Xu, Existence of log canonical closures, Invent. Math. 192 (2013), no. 1, 161–195.
  * [H19] K. Hashizume, Remarks on special kinds of the relative log minimal model program, Manuscripta Math. 160 (2019), no. 3, 285–314.
  * [HH20] K. Hashizume, Z. Hu, On minimal model theory for log abundant lc pairs, J. Reine Angew. Math., 767 (2020), 109–159.
  * [K13] J. Kollár, Singularities of the Minimal Model Program, Cambridge Tracts in Mathematics 200. Cambridge University Press, Cambridge, 2013.
  * [KM98] J. Kollár, S. Mori, Birational geometry of algebraic varieties. With the collaboration of C. H. Clemens and A. Corti. Translated from the 1998 Japanese original. Cambridge Tracts in Mathematics 134. Cambridge University Press, Cambridge, 1998.
  * [LP20] V. Lazić, T. Peternell, On generalized abundance, I, Publ. Res. Inst. Math. Sci. 56 (2020), no. 2, 353–389.
  * [LT22] V. Lazić, N. Tsakanikas, Special MMP for log canonical generalised pairs, with an appendix joint with Xiaowei Jiang, Selecta Math. (N.S.) (2022), Article no. 89.
  * [L14] J. Lesieutre, The diminished base locus is not always closed, Compositio Math. 150 (2014), 1729–1741.
  * [N04] N. Nakayama, Zariski-decomposition and abundance, MSJ Mem. 14, Mathematical Society of Japan, Tokyo, 2004.
  * [TX23] N. Tsakanikas, Z. Xie, Comparison and uniruledness of asymptotic base loci, preprint (2023), arXiv: 230901031v1.