# The Uncertain Future of Massive Binaries Obscures the Origin of LIGO/Virgo
Sources

K. Belczynski11affiliationmark:  , A. Romagnolo11affiliationmark:  , A.
Olejak11affiliationmark:  , J. Klencki22affiliationmark:  , D.
Chattopadhyay33affiliationmark:  , S. Stevenson33affiliationmark:  , M.
Coleman Miller44affiliationmark:  , J.-P. Lasota11affiliationmark:
55affiliationmark:  , Paul A. Crowther66affiliationmark:  1 Nicolaus
Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18,
00-716 Warsaw, Poland<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
2 Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL
Nijmegen, The Netherlands<EMAIL_ADDRESS>
3 Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, John St., Hawthorn, Victoria- 3122, Australia
<EMAIL_ADDRESS><EMAIL_ADDRESS>
4 Department of Astronomy and Joint Space-Science Institute University of
Maryland, College Park, MD 20742-2421, USA<EMAIL_ADDRESS>
5 Institut d’Astrophysique de Paris, CNRS et Sorbonne Université, UMR 7095,
98bis Boulevard Arago, 75014 Paris, France<EMAIL_ADDRESS>
6 Department of Physics & Astronomy, University of Sheffield, Hounsfield Rd,
Sheffield, S3 7RH, UK<EMAIL_ADDRESS>
###### Abstract

The LIGO/Virgo gravitational–wave observatories have detected at least 50
double black hole (BH) coalescences. This sample is large enough to have
allowed several recent studies to draw conclusions about the implied branching
ratios between isolated binaries versus dense stellar clusters as the origin
of double BHs. It has also led to the exciting suggestion that the population
is highly likely to contain primordial black holes. Here we demonstrate that
such conclusions cannot yet be robust, because of the large current
uncertainties in several key aspects of binary stellar evolution. These
include the development and survival of a common envelope, the mass and
angular momentum loss during binary interactions, mixing in stellar interiors,
pair-instability mass loss and supernova outbursts. Using standard tools such
as the rapid population synthesis codes StarTrack and COMPAS and the detailed
stellar evolution code MESA, we examine as a case study the possible future
evolution of Melnick 34, the most massive known binary star system (with
initial component masses of $144{\rm~{}M}_{\odot}$ and
$131{\rm~{}M}_{\odot}$). We show that, despite its fairly well-known orbital
architecture, various assumptions regarding stellar and binary physics predict
a wide variety of outcomes: from a close BH-BH binary (which would lead to a
potentially detectable coalescence), through a wide BH-BH binary (which might
be seen in microlensing observations), or a Thorne-Żytkow object, to a
complete disruption of both objects by pair-instability supernovae. Thus since
the future of massive binaries is inherently uncertain, sound predictions
about the properties of BH-BH systems are highly challenging at this time.
Consequently, drawing conclusions about the formation channels for the
LIGO/Virgo BH-BH merger population is premature.

###### Subject headings:

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

##  1\. Introduction

The LIGO/Virgo Collaboration (LVC) has reported gravitational-wave detections
of $\sim 50$ double black hole (BH-BH) coalescences (Abbott et al., 2021). The
majority of these can be explained as originating through any of several
channels, including isolated binary evolution, dynamics in dense stellar
clusters, or primordial black holes (Mandel & Broekgaarden, 2021). Until now,
observations contain only hints of the possible origin of observed double BHs.
For example, as anticipated prior to the detections (Belczynski et al.,
2010b), many of the black holes in these binaries have masses $\sim
30{\rm~{}M}_{\odot}$ or larger, which is considerably in excess of the most
massive stellar-origin BHs known through electromagnetic observations. Another
trend is that the effective spin of the binaries (which is the mass-weighted
projection of the black hole spins onto the orbital axis) is low; this could
be an indication of random orbits from dynamical processes, or could point
toward intrinsically low spins produced by efficient angular momentum
transport in massive stars (Spruit, 2002; Farr et al., 2017; Vitale et al.,
2017; Farr et al., 2018; Fuller & Ma, 2019; Bavera et al., 2020; Belczynski et
al., 2020a).

However, there are also individual events with characteristics that may be
more challenging to explain, and which therefore hold promise for
discriminating between formation channels. One such event is GW190814 (Abbott
et al., 2020b), which is an extremely asymmetric binary consisting of a
$\approx 23{\rm~{}M}_{\odot}$ black hole and a $2.6{\rm~{}M}_{\odot}$ object
that is either the lightest black hole or the heaviest neutron star yet
detected. Another is the double black hole event GW190521 (Abbott et al.,
2020a), which has two black holes which may have masses of $\sim
85{\rm~{}M}_{\odot}$ and $\sim 65{\rm~{}M}_{\odot}$, putting them both in the
pair-instability mass gap, although it is possible that the black hole mass
ratio is farther from unity and both black holes avoid the gap (Fishbach &
Holz, 2020; Nitz & Capano, 2021).

Based on these results, several groups have recently analyzed the BH-BH
population as a whole, with special attention to outliers such as GW190814 and
GW190521, to obtain insight into the relative fraction of events from
different formation channels. For example, Zevin et al. (2021) studied a
mixture of isolated binary evolution and dynamical formation in globular
clusters and concluded that neither channel can contribute more than $70\%$ to
the LIGO/Virgo observed population of BH-BH mergers. In contrast, Ossowski
(2021) disfavored the globular cluster channel in favor of isolated binary
evolution. Franciolini et al. (2021) investigated four formation channels
(isolated binaries, globular clusters, nuclear star clusters, and primordial
BHs) and found a high likelihood that primordial BH-BH mergers are part of
LIGO/Virgo source population.

All of these studies perform proper model comparison to infer which model or
mixture of models is favored. However, model comparison requires precisely
specified models, i.e. models with precisely defined physics. Here we
emphasize that the physics uncertainties (Schootemeijer et al., 2019) in at
least one of those models, of isolated binary evolution, are sufficiently
large that (to put it in Bayesian terms) the prior dominates the conclusion.
That is, different assumptions can lead to very different outcomes, which
means that population studies are not yet at the stage allowing strong and
credible conclusions to be drawn.

We demonstrate these model weaknesses using various, different assumptions in
performing simulations whose aim is to determine the fate of the most massive
known binary system. Despite the high binary frequency of massive stars (Sana
et al., 2012), eclipsing systems with primary masses $\gg 50{\rm~{}M}_{\odot}$
are exceptionally rare. Within the Milky Way the most extreme double-lined
systems are located in young, rich star clusters: A1 within NGC 3603 (Schnurr
et al., 2008), F2 within the Arches (Lohr et al., 2018) and WR20a in
Westerlund 2 (Bonanos et al., 2004). All three are short period (days), low
eccentricity systems with main sequence Wolf-Rayet (H-rich, WNh) primaries.

The most massive double-lined eclipsing system in the Large Magellanic Cloud
(LMC) is also a short period system, Hunter 38 in the Tarantula Nebula with a
O-type primary, whose mass is $\sim 57{\rm~{}M}_{\odot}$ (Massey et al.,
2002). The most extreme LMC binaries are non eclipsing systems, also within
the Tarantula Nebula. They have minimum dynamical primary masses in excess of
$\sim 50{\rm~{}M}_{\odot}$, albeit with orbits whose periods are an order of
magnitude longer and high eccentricities. From comparison with Bonn stellar
evolutionary models at LMC metallicity (Brott et al., 2011), R139 ($P_{\rm
orb}=154$ d, $e=0.38$) has a primary O supergiant mass of $\sim
80{\rm~{}M}_{\odot}$ (Mahy et al., 2020), R144 ($P_{\rm orb}=74$ d, $e=0.51$)
has a primary WNh mass of $\sim 110{\rm~{}M}_{\odot}$ (Shenar et al., 2021)
and Melnick 34 (Mk 34, $P_{\rm orb}=155$ d, $e=0.68$) – the current record-
holder - has component WNh+WNh masses of $139^{+21}_{-18}{\rm~{}M}_{\odot}$
and $127\pm 17{\rm~{}M}_{\odot}$ (Tehrani et al., 2019) with initial masses of
$144{\rm~{}M}_{\odot}$ and $131{\rm~{}M}_{\odot}$. Pollock et al. (2018) first
established that Mk 34 is a colliding wind binary from analysis of X-ray time
series observations, while Tehrani et al. (2019) noted that the potential fate
of Mk 34 involves a double stellar mass black hole binary merger.

We predict the fate of Mk 34, using different physics assumptions, applying
two rapid population synthesis codes (StarTrack and COMPAS) and the detailed
stellar evolutionary code (MESA). We find a wide variety of possible outcomes
(see Tab. 1) in terms of the black hole masses and orbital separations, and
even in terms of whether black holes will form at all. We therefore urge
caution in drawing important and credible conclusions about the LVC BH-BH
population based on models of massive–binary evolution containing, by
necessity, uncertain physics.

Table 1Fate of Mk 34

model | $M_{\rm BH1}$ | $M_{\rm BH2}$ | $t_{\rm delay}$ | fate
---|---|---|---|---
| [${\rm~{}M}_{\odot}$] | [${\rm~{}M}_{\odot}$] | [Myr] |
StarTrack1 | 22.5 | 22.1 | 47.5 | close BH-BH?a
StarTrack2 | 35.7 | 33.3 | 10,035 | close BH-BH?a
StarTrack3 | 36.0 | 32.5 | 11,663 | close BH-BH?a
COMPAS1 | 19.8 | 20.2 | $>t_{\rm hub}$ | wide BH-BH
COMPAS2 | 31.8 | 31.7 | $>t_{\rm hub}$ | wide BH-BH
COMPAS3 | 31.8 | 31.7 | $>t_{\rm hub}$ | wide BH-BH
MESA1 | 21.9 | 51.6b | — | Thorne-Żytkow
MESA2 | 35.2 | 80.9b | — | Thorne-Żytkow
MESA3 | 35.3 | 85.4b | — | Thorne-Żytkow
Pavlovskii1 | 21.9 | 22.1 | $>t_{\rm hub}$ | wide BH-BH
Pavlovskii2 | 35.2 | 33.3 | $>t_{\rm hub}$ | wide BH-BH
Pavlovskii3 | — | — | 3.3 | stellar mergerc
QuasiSingle1 | — | — | — | PSN+PSNd
QuasiSingle2 | $\sim 60$ | $\sim 60{\rm~{}M}_{\odot}$ | $>t_{\rm hub}$ | wide BH-BH
QuasiSingle3 | $\sim 30$ | $\sim 30{\rm~{}M}_{\odot}$ | $>t_{\rm hub}$ | wide BH-BH
QuasiSingle4 | $\sim 20$ | $\sim 20{\rm~{}M}_{\odot}$ | $>t_{\rm hub}$ | wide BH-BH

a: optimistic (non-standard) StarTrack models are used to get this result
b: for MESA models we list CE donor mass in column $M_{\rm BH2}$
c: merger of post-MS star and MS star: formation of very massive single star,
fate: PSN or a single BH
d: pair-instability supernovae disrupting binary components

##  2\. Calculations

For initial properties of Mk 34 we select: $M_{\rm a}=144{\rm~{}M}_{\odot}$,
$M_{\rm b}=131{\rm~{}M}_{\odot}$, $e=0.68$, $a=760{\rm~{}R}_{\odot}$ chosen to
result in an orbital period of $P_{\rm orb}=155$ d after $0.6$ Myr (current
age, Tehrani et al., 2019) of system evolution with the StarTrack code. We
adopt the LMC metallicity of $Z=0.006$ (Rolleston et al., 2002).

###  2.1. StarTrack calculations

We use the population synthesis code StarTrack (Belczynski et al., 2020a),
which employs analytic fits to evolutionary tracks of non-rotating stellar
models (Hurley et al., 2000). We adopt standard wind losses for massive stars
from Vink et al. (2001) and LBV winds as $(dM/dt)_{\rm lbv}=f_{\rm
lbv}10^{-4}{\rm~{}M}_{\odot}{\rm~{}yr}^{-1}$ with $f_{\rm lbv}=1.5$ from
Belczynski et al. (2010a).

For stars that overfill their Roche lobes, we initiate mass transfer between
binary components and associated (if any) mass loss from binary systems. If
the binary is not circularized by tidal interactions we circularize it
($e_{\rm new}=0$) to periastron distance ($a_{\rm new}=a(1-e)$) in one
timestep and only then start Roche lobe overflow (RLOF). For nuclear-timescale
mass transfer (NTMT) and thermal-timescale mass transfer (TTMT) we use the
standard formalism, while we use a diagnostic mass ratio diagram as a
criterion for common envelope (CE) development (see Sec.5 of Belczynski et al.
(2008)). During the TTMT/NTMT the fraction of mass lost by the donor star that
is accumulated by non-degenerate companion stars is set to $f_{\rm a}=0.5$,
while the rest is lost with specific angular momentum (expressed in units of
$2\pi a^{2}/P_{\rm orb}$) of $j_{\rm loss}=1.0$ (see eq. 33 of Belczynski et
al. (2008)). The accumulation of mass on compact objects (e.g., NS/BH) is
limited by the (Eddington) critical accretion rate and mass is lost with the
specific angular momentum of the compact accretor (King et al., 2001; Mondal
et al., 2020).

We employ the delayed core-collapse supernova (SN) engine in NS/BH mass
calculation (Fryer et al., 2012) which allows for populating the lower mass
gap between NSs and BHs (Belczynski et al., 2012; Zevin et al., 2020). We
employ weak pair-instability pulsation supernova (PPSN) mass-loss and pair-
instability supernova (PSN) model that results in upper mass gap: no BHs with
mass $M_{\rm BH}\gtrsim 55{\rm~{}M}_{\odot}$ (Belczynski et al., 2020a). We
allow for the fallback decreased NS/BH natal kicks with
$\sigma=265{\rm~{}km}{\rm~{}s}^{-1}$ and no natal kicks for direct BH
formation. This is our standard input physics marked as ”StarTrack1” model in
Table 1.

The development of the CE phase is a big issue in stellar/binary astrophysics
(Ivanova et al., 2013b; Olejak et al., 2021). We are agnostic about which
systems should be sent to a CE and which should evolve through stable RLOF. In
StarTrack models we allow for the most optimistic scenario (see Sec. 3.1), and
we send nearly all systems through the CE to form (potentially) BH-BH mergers.
We do not do this on regular basis. According to our standard input physics,
donors with radiative envelopes (e.g., in the Hertzsprung gap) do not enter
the CE phase. Since it is not fully understood how exactly a CE develops, we
test various assumptions to show how these influence the future fate of binary
systems such as Mk 34. Contrasting models are being presented as well.

During CE events the entire envelope of the donor is assumed to be lost from
the binary, with the exception of compact object companions that are allowed
to accrete a small fraction of donor’s envelope at $5\%$ of the Bondi rate
(MacLeod et al., 2017a). The CE orbital decay is calculated with the standard
energy-balance formalism (Webbink, 1984) in which we adopt a $100\%$
efficiency of the orbital energy transfer ($\alpha=1.0$) into the envelope,
while the binding energy is parameterized by detailed stellar models
($\lambda$ scaling: Xu & Li (2010); Dominik et al. (2012)).

In “StarTrack2” we decrease wind mass–loss rates for LBV stars to $f_{\rm
lbv}=0.48$ and we increase the He core mass at the end of the main sequence by
a factor of $f_{\rm core}=1.5$ with respect to the original Hurley et al.
(2000) models. This model approximately reproduces the basic properties of the
$131{\rm~{}M}_{\odot}$ and $144{\rm~{}M}_{\odot}$ models at terminal-age main
sequence (TAMS) obtained in our MESA computations (see Sec. 2.3).

In “StarTrack3” we circularize massive binaries with angular momentum
conservation ($a_{\rm new}=a(1-e^{2})$), we set $j_{\rm loss}=0.1$, and
$f_{\rm a}=0.25$ while keeping the rest of the input physics as in
“StarTrack2”. This model aims to test the survival of the CE phase in an Mk
34-like future evolution (see Sec. 3). In practice, such set-up allows the
secondary star of Mk 34 to have a large radius (wide binary orbit) during the
RLOF and therefore potentially to develop a convective envelope but to survive
the CE phase.

###  2.2. COMPAS calculations

We use the population synthesis code COMPAS (Stevenson et al., 2017; Vigna-
Gómez et al., 2018; Chattopadhyay et al., 2021), which incorporates stellar
(Hurley et al., 2000) and binary evolution (Hurley et al., 2002) with updated
wind prescriptions for massive stars (Vink et al., 2001; Belczynski et al.,
2010a). The LBV wind losses are as in StarTrack (see Sec. 2.1).

In its default set-up, the mass accretion during RLOF onto degenerate stars
(i.e. black holes, neutron stars and white dwarfs) is Eddington limited in
COMPAS. For non-degenerate stars, this accretion limit is set by the ratio of
the rate at which mass is being donated to the rate of which mass can be
incorporated into the companion star. These two rates, in turn, are determined
by the thermal (Kelvin-Helmholtz) timescales of the donor and the companion
respectively (Kalogera & Webbink, 1996). The thermal timescale is an explicit
function of the mass, radius and luminosity of the star. For a star of total
mass $M$, envelope mass $M_{\rm env}$, radius $R$ and luminosity $L$, the
thermal timescale $\tau_{\rm KH}$ is given by $\tau_{\rm KH}=GMM_{\rm
env}/RL$, where $G$ is the universal gravitational constant. In binaries,
since these stellar parameters are determined by the orbital period at which
the donor overflows its Roche lobe, the thermal timescale of the donor becomes
an implicit function of the orbital properties of the binary (Schneider et
al., 2015). Thus the accretion efficiency in binaries is primarily determined
by their orbital period. The mass transfer efficiency $\beta_{\mathrm{acc}}$
(ratio of the mass gained by the donor to the mass lost from the
companion)111Equivalent of the parameter $f_{\rm a}$ in StarTrack described in
Sec. 2.1. in a non-degenerate star can be thus expressed as
$\beta_{\mathrm{acc}}=\mathrm{min}(1,10\frac{\tau_{\mathrm{KH,acc}}}{\tau_{\mathrm{KH,don}}})$,
where $\tau_{\mathrm{KH,acc}}$ and $\tau_{\mathrm{KH,don}}$ are the thermal
timescales of the accretor and donor respectively (Hurley et al., 2002;
Schneider et al., 2015; Chattopadhyay et al., 2021).

The stability of the mass loss is determined by the parameter $\zeta$
(critical mass-radius exponent for development of CE) in COMPAS using fits
from the Ge et al. (2015) simulations as described in Vigna-Gómez et al.
(2018). In nearly equal-mass, close binaries the thermal timescales of the
donor and accretor being very similar, the mass transfer is usually
conservative and remains stable. In close systems with more extreme mass
ratio, the thermal timescale of the donor being much longer than the accretor,
the mass transfer becomes non-conservative leading to a CE phase. The Ge et
al. (2015) criteria renders the mass transfer from evolved (non-main sequence,
non-degenerate) massive stars as predominantly stable (Neijssel et al., 2019),
which is very similar to the Pavlovskii et al. (2017) model (discussed in
section. 2.4). We assume an isotropic re-emission model for angular momentum
loss during non-conservative stable RLOF (Pols et al., 1998).

Unlike in StarTrack, the binary is not circularized right before RLOF (Vigna-
Gómez et al., 2018). But binaries that survive the CE events are always
circularized (Vigna-Gómez et al., 2018). All other default RLOF and CE mass
transfer specifications in COMPAS are identical to StarTrack as described in
Sec. 2.1.

Our standard model COMPAS1 utilizes the Fryer et al. (2012) “delayed’
supernovae prescription and pre-supernova core mass to post-supernova remnant
mass mapping. The (pulsational) pair-instability supernovae modeling is
implemented in COMPAS (Stevenson et al., 2019) with polynomial fitting from
the models by Marchant et al. (2018) as the default input. The natal kick
distributions (including fallback) for BHs and NSs are identical to StarTack1
model.

In the model COMPAS2 we reduce the LBV wind mass loss rate to
$f_{\mathrm{lbv}}=0.48$ (from $f_{\mathrm{lbv}}=1.5$ in COMPAS1). We also
increase the He core mass of the terminal main sequence stars by a factor of
$f_{\mathrm{core}}=1.5$ multiplied to the fitting formula from Hurley et al.
(2000) (eq. 30). COMPAS2, like StarTrack2 is adjusted to approximately
reproduce the total mass and the core mass at TAMS of $131{\rm~{}M}_{\odot}$
and $144{\rm~{}M}_{\odot}$ models computed with MESA.

The model COMPAS3 is identical to COMPAS2 but here we allow the binary to
circularize (while conserving its angular momentum) right before the onset of
RLOF as detailed for StarTrack in Sec. 2.1.

###  2.3. MESA calculations

####  2.3.1 Calibration of TAMS core masses

The fitting formula by Hurley et al. (2000) to the evolutionary tracks from
Pols et al. (1998), which are the basis for StarTrack and COMPAS population
synthesis codes, are based on stellar models computed for stars with masses up
to $40{\rm~{}M}_{\odot}$. The treatment of more massive stars in StarTrack and
COMPAS relies on extrapolation. In the mass range considered in this study
($130-145{\rm~{}M}_{\odot}$), this can lead to a significant deviation in
basic stellar properties from what detailed stellar models produce (or what is
inferred from observations). One property that is particularly inaccurate due
to the extrapolation of Hurley et al. (2000) formulae (and which is crucial
for considerations of the final fate of the Mk 34 system) is the ratio of the
helium core mass to the total star mass at TAMS ($M_{\rm core,TAMS}/M_{\rm
TAMS}$). For very massive stars this ratio is close to unity (Yusof et al.,
2013; Köhler et al., 2015). In contrast, a $144{\rm~{}M}_{\odot}$ star at
$Z=0.006$ metallicity evolved with StarTrack or COMPAS, even though
significantly stripped through MS winds ($M_{\rm TAMS}=77.4{\rm~{}M}_{\odot}$)
is far away from being a helium star with the helium core mass of only $M_{\rm
core,TAMS}=34.1{\rm~{}M}_{\odot}$.

To correct for this and calibrate the properties of the StarTrack and COMPAS
models at TAMS, we use the MESA 1D stellar-evolution code (Paxton et al.,
2011, 2013, 2015, 2018, 2019)222MESA version r15140,
http://mesa.sourceforge.net/. We compute single models of
$131{\rm~{}M}_{\odot}$ and $144{\rm~{}M}_{\odot}$ stars at $Z=0.006$
metallicity. The relative initial abundances of metals follow Grevesse et al.
(1996). We model convection by using mixing-length theory (Böhm-Vitense, 1958)
with a mixing-length parameter $\alpha=2.0$, and we adopt the Schwarzschild
criterion for convection. We used the Dutch wind setup in MESA, which combines
different prescriptions depending on the effective temperature $T_{\rm eff}$
and the fractional surface hydrogen abundance $H_{\rm sur}$. As shown in Table
2, for $T_{\rm eff}<10^{4}$ K the code uses the mass-loss rates from de Jager
et al. (1988), regardless of the hydrogen surface abundance. For $T_{\rm
eff}\geq 10^{4}$ K, MESA adopts either the Nugis & Lamers (2000) prescriptions
(if $H_{\rm sur}<0.4$), or mass-loss rates from Vink et al. (2001) (if $H_{\rm
sur}\geq 0.4$). Additionally, mass loss rates in MESA can be increased or
decreased by changing a specific scaling factor $f_{\rm wind}$. The standard
prescription of Vink et al. (2001) is known to underestimate the empirical
mass-loss rates of very massive MS stars, which increase dramatically as they
approach the Eddington limit, $\Gamma_{e}$ (Vink et al., 2011; Bestenlehner et
al., 2014; Bestenlehner, 2020). Indeed, clumping corrected mass-loss rates of
the components of Mk 34 from Tehrani et al. (2019) exceed Vink et al. (2001)
prescriptions by factors of 2–3.

Table 2Dutch Stellar Winds in MESAa

|  $T_{\rm eff}<10^{4}$ K |  $T_{\rm eff}\geq 10^{4}$ K
---|---|---
– | de Jager et al. (1988) | –
$H_{\rm sur}<0.4$ | – | Nugis & Lamers (2000)
$H_{\rm sur}\geq 0.4$ | – | Vink et al. (2001)

a: based on $T_{\rm eff}$ winds either depend or do not depend on $H_{\rm
sur}$

We account for convective overshooting above the hydrogen-burning core by
applying the step overshooting formalism, which extends the convective core by
a fraction $\delta_{ov}$ of the local pressure scale height.

We initialize our models with the initial rotation of $V_{\rm
i}=250{\rm~{}km}{\rm~{}s}^{-1}$ (guided by the analysis of Tehrani et al.,
2019). For rotational mixing, we include the effects of the Eddington-Sweet
circulation, secular shear instabilities, and the Goldreich-Schubert-Fricke
instability, with an efficiency factor $f_{\rm c}=1/30$ (Heger et al., 2000;
Brott et al., 2011).

We avoid using the MLT++ option in MESA (Paxton et al., 2013). As a result,
models that reach the red supergiant stage encounter numerical difficulties in
their superadiabatic outer envelope layers (Pavlovskii & Ivanova, 2015;
Klencki et al., 2020), which prohibits us from following their evolution to
the point of maximum radial expansion. However, for the purpose for the
current study, we are only interested in the properties of MESA models at TAMS
as well as whether or not the models expand sufficiently to lead to a RLOF in
the Mk 34 binary system. We thus stop the MESA computations when the radius of
$2000{\rm~{}R}_{\odot}$ is reached. Otherwise, we stop the simulation after
10,000 MESA steps. Such stopping conditions are sufficient for our purposes in
all the considered scenarios for the Mk 34 system.

For calibration of population synthesis models at TAMS we calculate a MESA
stellar model with initial mass $M_{\rm ZAMS}=144{\rm~{}M}_{\odot}$,
metallicity $Z=0.006$, and $V_{\rm i}=250{\rm~{}km}{\rm~{}s}^{-1}$. We apply
the standard Dutch winds ($f_{\rm wind}=1.0$, though see above) and step
overshooting of $\delta_{ov}=0.12$, to maintain consistency with the
overshooting in models by Pols et al. (1998) and the Hurley et al. (2000)
fits. At the end of the MS, this model has a mass of $M_{\rm
TAMS}=94.0{\rm~{}M}_{\odot}$ with a He core mass of $M_{\rm
core,TAMS}=66.6{\rm~{}M}_{\odot}$ (see Table 3). Post-MS expansion leads to a
maximum radius of $R_{\rm max}=1968{\rm~{}R}_{\odot}$ at the end of our
simulation, at which point the star is still expanding as a red supergiant. A
MESA stellar model with $M_{\rm ZAMS}=131{\rm~{}M}_{\odot}$, $Z=0.006$, and
$V_{\rm i}=250{\rm~{}km}{\rm~{}s}^{-1}$ results in $M_{\rm
MS}=86.6{\rm~{}M}_{\odot}$ and $M_{\rm core}=58.5{\rm~{}M}_{\odot}$ and
expands beyond $2000{\rm~{}R}_{\odot}$ (see Table 3).

The $M_{\rm ZAMS}=144{\rm~{}M}_{\odot}$ ($Z=0.006$) non-rotating StarTrack or
COMPAS model produces $M_{\rm TAMS}=77.4{\rm~{}M}_{\odot}$ and $M_{\rm
core,TAMS}=34.1{\rm~{}M}_{\odot}$. We decrease the winds during MS, keeping
the original Vink et al. (2001) prescriptions, but decreasing the LBV winds to
$f_{\rm lbv}=0.48$ to get a model with $M_{\rm TAMS}=94.2{\rm~{}M}_{\odot}$
and $M_{\rm core,TAMS}=44.5{\rm~{}M}_{\odot}$. Next, we increase the core size
by $f_{\rm core}=1.5$ to get the target values: $M_{\rm
TAMS}=94.2{\rm~{}M}_{\odot}$ and $M_{\rm core,TAMS}=66.7{\rm~{}M}_{\odot}$ in
the population synthesis codes. Note that such a massive star is already
luminous ($L\sim 3\times 10^{6}L_{\odot}$) and cold enough ($T_{\rm eff}\sim
30,000$ K) to be beyond the Humphreys-Davidson limit (Humphreys & Davidson,
1994) and subject to LBV winds on the MS. Applying the same calibration to
$M_{\rm zams}=131{\rm~{}M}_{\odot}$ ($Z=0.006$) we obtain in population
synthesis codes: $M_{\rm MS}=89.5{\rm~{}M}_{\odot}$ and $M_{\rm
core}=62.3{\rm~{}M}_{\odot}$. We apply this calibration for all metallicities.
However, one should note that some observations may be in contradiction of
metallicity-independent LBV winds (Gilkis et al., 2021).

Table 3 Properties of MESA models used to calibrate the population synthesis simulations (the two top rows) as well as a few models with increased wind mass-loss or core overshooting (see Sec. 2.3.2). We show the maximum radius reached in our simulation as well as the total mass and the He core mass at the end of the MS phase. model | $R_{\rm Max}$ | $M_{\rm TAMS}$ | $M_{\rm core,TAMS}$
---|---|---|---
| [${\rm~{}R}_{\odot}$] | [${\rm~{}M}_{\odot}$] | [${\rm~{}M}_{\odot}$]
144 ${\rm~{}M}_{\odot}$ ($\delta_{ov}=0.12$, $f_{\rm wind}=1$) |  $>$1968 | 94 | 67
131 ${\rm~{}M}_{\odot}$ ($\delta_{ov}=0.12$, $f_{\rm wind}=1$) |  $>$2000 | 87 | 58
144 ${\rm~{}M}_{\odot}$ ($\delta_{ov}=0.33$, $f_{\rm wind}=1$) | 46 | 62 | 61
131 ${\rm~{}M}_{\odot}$ ($\delta_{ov}=0.33$, $f_{\rm wind}=1$) | 41 | 59 | 58
144 ${\rm~{}M}_{\odot}$ ($\delta_{ov}=0.12$, $f_{\rm wind}=1.5$) | 40 | 57 | 55
131 ${\rm~{}M}_{\odot}$ ($\delta_{ov}=0.12$, $f_{\rm wind}=1.5$) | 40 | 63 | 53

####  2.3.2 Calibration of post-MS expansion

When considering the future fate of the Mk 34 system, a key question is
whether its very massive components will expand after the end of the MS and
initiate a mass-transfer interaction or whether they will lose their hydrogen
envelopes already during the MS and smoothly transition to become compact
helium stars, avoiding any RLOF. The two crucial aspects that affect the
degree of the post-MS expansion of very massive stars is the amount of core
overshooting and the strength of stellar winds. Here, we explore this by
computing a small grid of MESA models with different overshooting and wind
assumptions.

For the calibration of TAMS properties of $131{\rm~{}M}_{\odot}$ and
$144{\rm~{}M}_{\odot}$ models, we assumed a modest core overshooting of
$\delta_{ov}=0.12$, following the calibration to low-mass stars in open
clusters by Pols et al. (1998). More recently, Choi et al. (2016) found the
best agreement with the properties of the Sun for a MESA model with
$\delta_{ov}=0.16$. However, there is an increasing amount of evidence that
core overshooting could be significantly larger in the case of massive stars
(e.g., Brott et al., 2011; Castro et al., 2014; Claret & Torres, 2018; Scott
et al., 2021). In particular, the calibration by Brott et al. (2011) to match
the observed drop in rotational velocities of post-MS B stars (although see
Vink et al., 2010) resulted in $\delta_{ov}=0.33$, a value that has become
widely used to compute stellar models of massive stars in the recent years. On
the other hand, there is no observational calibration of core overshooting in
the case of very massive stars of masses above $100{\rm~{}M}_{\odot}$. As
such, we explore six different $\delta_{ov}$ values from a wide range between
$0.12$ and $0.5$.

The Dutch wind scheme in MESA incorporates the Vink et al. (2001) prescription
for optically thin line-driven winds of hot MS stars. However, as noted above,
very massive stars possess a sufficiently high luminosity to mass ratio on
their MS that they approach their Eddington limit, leading to high mass-loss
rates (Gräfener & Hamann, 2008; Vink et al., 2011; Bestenlehner et al., 2014;
Bestenlehner, 2020). Here, we attempt to correct for this by simply increasing
the wind scaling factor $f_{\rm wind}$ from $1.0$ to $1.5$ or to $2.0$.

We compute a grid of $131{\rm~{}M}_{\odot}$ and $144{\rm~{}M}_{\odot}$ models
with the above variations in overshooting and winds (and all the other
assumptions same as in our calibration models in Sec. 2.3.1). All the results
are shown in Appendix A and a few selected examples in Table 3. In short, we
find that any model with overshooting of $\delta_{ov}=0.33$ or higher, or a
wind multiplication factor $f_{\rm wind}\geq 1.5$ evolves to become a helium
WR star already by the end of MS, avoiding radial expansion beyond
$100{\rm~{}R}_{\odot}$ and any RLOF interaction in the Mk 34 binary. This
result is at the basis of the quasi-single evolutionary scenario for Mk 34,
see Sec. 2.5.

####  2.3.3 Calibration of envelope binding energies

Recent studies by Klencki et al. (2021) and Marchant et al. (2021) have shown
that the envelope binding energies used in StarTrack and COMPAS (i.e.
$\lambda$ scaling following Xu & Li (2010); Dominik et al. (2012)) may be
severely underestimated in the case of massive stars with outer radiative
envelopes. Note that population synthesis codes do not perform CE evolution
for Hertzsprung gap stars (radiative outer envelope) under standard
assumptions on input physics, but typically the CE is applied for core-helium
burning stars even for those with outer radiative envelopes.

To explore the effect of revised binding energies on the future fate of Mk 34,
whenever our StarTrack binary evolution calculation predicts a CE phase to
occur, we use MESA to compute a detailed stellar model of the donor star. We
then follow the method outlined in Klencki et al. (2021) to integrate through
the envelope of the MESA model and compute its binding energy. When matching
the properties of a MESA model with those from StarTrack, we ensure that the
CE donor has the same helium core mass and the same radius but allow for a
lower envelope mass in the MESA model, so that we may be under- but never
over-estimating the envelope binding energy. This allows for conservative
statements on inability of a binary to eject the donor’s envelope and the CE
survival.

###  2.4. Pavlovskii calculations

Models presented below are obtained with modified StarTrack. In particular, we
use more restrictive criteria for the CE development (Pavlovskii et al.,
2017), and allow more binaries to evolve through stable mass transfer instead
(Olejak et al., 2021). The new criteria are applied to (i) H-rich post MS
donor stars, (ii) for initial masses larger than $18{\rm~{}M}_{\odot}$, (iii)
when the mass ratio (companion to donor mass at CE onset) fulfils the
condition $q_{\rm CE}<q_{\rm crit}$ for CE to develop, where $q_{\rm
crit}=0.19-0.36$ depending on donor mass and metallicity, and (iv) when the
donor’s radius at the onset of CE fulfils specific criteria (shown at Fig. 2
and 3 of Olejak et al. (2021)) for CE to develop. These new criteria lead to
the emergence of BH-BH formation channels without CE in StarTrack simulations
(Olejak et al., 2021). This is the same channel that was proposed by van den
Heuvel et al. (2017) and that is also found in COMPAS simulations (Stevenson
et al., 2019; Neijssel et al., 2019).

Models labeled ”Pavlovskii1” and ”Pavlovskii2” correspond to models
”StarTrack1” and ”StarTrack2” but with modified CE development criteria,
respectively.

In model ”Pavlovskii3” we test different formulae for the loss of the angular
momentum during mass transfer through the L2 point given by MacLeod & Loeb
(2020): $j_{\rm loss}=j_{\rm L2}=1.2^{2}\frac{M_{\rm tot}^{2}}{M_{\rm
don}M_{\rm acc}}$ instead of our standard $j_{\rm loss}=1.0$ (Podsiadlowski et
al., 1992). We expect much higher loss of angular momentum with this
modification, which may result in another potential fate for the future
evolution of Mk 34 (see Sec. 3.4). The rest of input physics of ”Pavlovskii3”
corresponds to ”Pavlovskii2”. A change of the circularization scheme would
only increase binary orbital separation.

###  2.5. Quasi single star calculations

In this part we approximate the evolution of non-expanding components of Mk
34. This is driven by MESA models with large overshooting or increased stellar
winds that do not show significant post-MS expansion (see Tab. A). In
particular, MESA models with $\delta_{\rm ov}=0.33$ and standard winds
($f_{\rm wind}=1.0$) reach a maximum radius of $46{\rm~{}R}_{\odot}$ for
$M_{\rm ZAMS}=144{\rm~{}M}_{\odot}$ and $41{\rm~{}R}_{\odot}$ for $M_{\rm
ZAMS}=131{\rm~{}M}_{\odot}$. Additionally, important for the development of
the PPSN/PSN, stellar models that do not expand may have very different TAMS
helium core masses. For example, the $M_{\rm ZAMS}=144{\rm~{}M}_{\odot}$ model
produces $M_{\rm core,TAMS}=61{\rm~{}M}_{\odot}$ for $\delta_{\rm ov}=0.33$
and $f_{\rm wind}=1.0$ or $M_{\rm core,TAMS}=30{\rm~{}M}_{\odot}$ for
$\delta_{\rm ov}=0.4$ and $f_{\rm wind}=1.5$. The former model is possibly
subject to a PPSN/PSN while the latter is not (Woosley, 2017; Farmer et al.,
2020).

In Sec. 3.5 we explain our choice of models, showing how uncertainties can
affect the future evolution of Mk 34.

##  3\. Examples of Mk 34 Future evolution

Various predicted models of Mk 34 future evolution are illustrated in Figure
1, summarized in Table 1, and described below.

Figure 1.—  Future evolution of Melnick 34. The fate of this massive binary is
subject to a number of stellar and binary evolution uncertainties. Depending
on the adopted evolutionary model, Mk 34 may form a close or wide BH-BH
system, a Thorne-Żytkow (TZ – Thorne & Zytkow, 1977) object or end its life in
two pair-instability supernovae (PSN). RLOF: Roche lobe overflow, CE: common
envelope, TTMT: thermal-timescale mass transfer, ZAMS: Zero Age Main Sequence,
BH: black hole, GRB: gamma-ray burst, PPSN: pair-instability pulsation
supernova.
⋆: COMPAS track gives somewhat lower mass BH (formed out of an initially more
massive star) than other binary scenarios.

###  3.1. StarTrack models

A binary star resembling Mk 34 is evolved with the StarTrack2 model (Sec.
2.1). Star A (initially more massive) expands as it evolves and finally
periodically overfills its Roche lobe ($R_{\rm a}=139{\rm~{}R}_{\odot}$) at
periastron passages. At this point the orbit expanded from the initial
$a=760{\rm~{}R}_{\odot}$ to $1125{\rm~{}R}_{\odot}$ due to wind mass loss from
both binary components ($M_{\rm a}=94.2{\rm~{}M}_{\odot}$, $M_{\rm
b}=91.9{\rm~{}M}_{\odot}$) while the eccentricity remained much unchanged
($e=0.68$). The tidal circularization force is the strongest at periastron and
we assume that the orbital motion is circularized to periastron distance and
leads to normal stable RLOF on new ($a=360{\rm~{}R}_{\odot}$) circular orbit.
RLOF leads first to TTMT which subsequently transforms to a NTMT. Star A is
stripped almost entirely of its H-rich envelope ($M_{\rm
a}=67.4{\rm~{}M}_{\odot}$, $M_{\rm a,core}=67.3{\rm~{}M}_{\odot}$) while star
B accreted half of that lost envelope ($M_{\rm b}=105{\rm~{}M}_{\odot}$) while
the rest of the mass has been lost from the binary. In response the orbit
increased in size ($a=439{\rm~{}R}_{\odot}$). Star A is a Wolf-Rayet with
heavy wind mass–loss and at the end of its nuclear evolution its mass
decreases to $M_{\rm a}=35.6{\rm~{}M}_{\odot}$ ($a=1370{\rm~{}R}_{\odot}$,
$M_{\rm b}=99.1{\rm~{}M}_{\odot}$). Star A collapses directly to a BH with a
mass of $M_{\rm a}=35.2{\rm~{}M}_{\odot}$ ($1\%$ neutrino mass loss, no
baryonic mass loss, no natal kick). Then star B evolves and expands to fill
its Roche lobe in a circular orbit with $a=1515{\rm~{}R}_{\odot}$ ($M_{\rm
b}=88.5{\rm~{}M}_{\odot}$, $R_{\rm b}=700{\rm~{}R}_{\odot}$, $T_{\rm
b,eff}=8690$ K). This time in our standard approach, due to relatively high
mass ratio ($q=88.5/35.2=2.5$) RLOF is evaluated to lead to a CE phase. We
estimate the binding energy of star’s A envelope ($\lambda=0.103$) to be low
enough to be ejected at the cost of orbital energy. After envelope ejection
star B becomes a massive stripped He core ($M_{\rm b}=62.0{\rm~{}M}_{\odot}$)
and the orbit decays to $a=34.2{\rm~{}R}_{\odot}$. During CE the first-formed
BH accretes $\sim 0.5{\rm~{}M}_{\odot}$ ($M_{\rm a}=35.7{\rm~{}M}_{\odot}$).
Star B, after a Wolf-Rayet wind mass loss ($M_{\rm b}=33.7{\rm~{}M}_{\odot}$,
$a=47.8{\rm~{}R}_{\odot}$) collapses directly to a BH ($M_{\rm
b}=33.3{\rm~{}M}_{\odot}$). After $t_{\rm evol}=4.1$Myr of binary evolution a
close BH-BH binary is formed with a coalescence time of $t_{\rm
coal}=10.0$Gyr.

There are caveats in this scenario. Star B at the time of the RLOF onset has
just finished core H-fusion and is a Hertzsprung-gap star in the transition to
become a core-He burning giant. It was argued that such stars do not have a
clear core-envelope structure and that the CE phase should always lead to
merging of the donor star with its companion (Belczynski et al., 2007). This
finds some support in observations as the predicted BH-BH merger rates that
allow for such a scenario as presented above are too high to match the
empirical LIGO/Virgo estimate (see submodels A in Table 4 of Belczynski et al.
(2020a)). In addition, star B has an outer radiative envelope: at the time of
RLOF and with surface properties $R_{\rm b}=700{\rm~{}R}_{\odot}$, $T_{\rm
b,eff}=8690$K $\rightarrow L=2.5\times 10^{6}{\rm~{}L}_{\odot}$, it is well
above the effective temperature threshold below which stars at LMC metallicity
have convective envelopes ($T<3900$K; see Fig.6 of Klencki et al. (2020)).
Klencki et al. (2021) argued that massive radiative-envelope giants have
binding energies that are too high to allow for a successful CE ejection in
BH-BH merger progenitor binaries. In StarTrack models, the estimate of the
binding energy of star B with the $\lambda$ formalism is only an approximation
that is needed for the use in rapid population synthesis models. It should be
also stressed that in the standard input physics of StarTrack models we would
allow no Hertzsprung-gap star to survive a CE phase no matter what is our
estimate of the star binding energy ($\lambda$) or assumed efficiency of
orbital energy transfer to the envelope ($\alpha$). This is why we call the
StarTrack models used here as optimistic scenarios. We test this optimistic
estimate with MESA in Sec. 3.3.

In StarTrack1 model a similar scenario develops. However, since the stars and
their cores are less massive the BH masses are $\sim 10{\rm~{}M}_{\odot}$
smaller than in model StarTrack2 (see Tab. 1).

In StarTrack3 model we adjust evolutionary parameters in such a way that CE
survival is less caveated than in model StarTrack2. We alter circularization
process and we change RLOF parameters setting mass transfer/loss to obtain
wider binary than in StarTrack2 model. This allows star B to expand more
before it initiates CE phase: $R_{\rm b}=1439{\rm~{}R}_{\odot}$,
$a=3141{\rm~{}R}_{\odot}$, $M_{\rm b}=85.9{\rm~{}M}_{\odot}$, $M_{\rm
b,core}=60.1{\rm~{}M}_{\odot}$, $T_{\rm b,eff}=5986$ K, $\lambda=0.050$. This
star has almost convective envelope, but not quite so. If we perform CE energy
balance with the above parameters this system survives CE and forms close BH-
BH binary at $t_{\rm evol}=3.9$ Myr and with $t_{\rm coal}=11.7$ Gyr. Yet, the
same caveats remain as for StarTrack2 model.

###  3.2. COMPAS models

With the initial masses, orbital period and eccentricity of Mk 34 (see Sec. 2)
we evolve COMPAS1, COMPAS2 and COMPAS3 models with the individual variations
specified in Sec. 2.2. As in Sec. 3.1, we will always refer to the originally
more massive star (with $M_{\rm ZAMS}=144{\rm~{}M}_{\odot}$) as star A.

In COMPAS2 model, star A ends core-H burning with a total mass of
$94.2{\rm~{}M}_{\odot}$, while the stellar winds increase the separation to
$a=1165.3{\rm~{}R}_{\odot}$. As star A leaves the MS, its core mass is
calculated to be $M_{\rm a,core}=66.8{\rm~{}M}_{\odot}$. Star A overfills its
Roche lobe soon during post-MS evolution losing its envelope
($27.4{\rm~{}M}_{\odot}$) in a stable (fully conservative) TTMT RLOF phase.
This is the outcome of the binary being fairly wide and the donor being an
evolved (post-MS) star, as commented on in Sec. 2.2 (see also Schneider et al.
(2015)). Eccentricity remains unchanged ($e=0.68$) during this phase and the
orbit expands to $a=1374.9{\rm~{}R}_{\odot}$. Star B, which had a total mass
of $91.8{\rm~{}M}_{\odot}$ right before the mass transfer, becomes a
$119.2{\rm~{}M}_{\odot}$ MS star. Being stripped off its H-rich envelope, star
A enters the naked helium-star (WR) phase with mass $M_{\rm
a}=66.8{\rm~{}M}_{\odot}$. Evolution continues while both stars are losing
mass in winds. Star A with mass $M_{\rm a}=35.4{\rm~{}M}_{\odot}$ undergoes
direct core-collapse, forming a BH of mass $M_{\rm a}=31.8{\rm~{}M}_{\odot}$,
while star B with $M_{\rm b}=115.3{\rm~{}M}_{\odot}$ is still on the MS. The
orbital separation, right after the formation of BH A, becomes
$a=1712.1{\rm~{}R}_{\odot}$. Star B leaves the MS with mass $M_{\rm
b}=93.7{\rm~{}M}_{\odot}$ and core mass $M_{\rm b,core}=66.3{\rm~{}M}_{\odot}$
when the orbital separation is $a=2006.6{\rm~{}R}_{\odot}$. Shortly thereafter
star B fills its Roche lobe, loses most of its H-rich envelope ($\sim
27.3{\rm~{}M}_{\odot}$) and becomes a naked helium star in stable TTMT RLOF.
The post-RLOF orbital separation is decreased to $a=918.1{\rm~{}R}_{\odot}$.
Eddington limited accretion allows BH A to gain only $\sim 3.7\times
10^{-5}{\rm~{}M}_{\odot}$. The binary at the onset of this RLOF, though has
slightly smaller separation than the previous RLOF phase, is still fairly wide
and results in a stable TTMT RLOF despite rather high mass ratio
($q=93.7/31.8=2.9$). It is noted that both RLOF phases noted in COMPAS2 model
are fairly similar to Pavlovskii2 model (see Sec. 3.4). However, the binary
remains eccentric through both RLOF phases ($e=0.673$). Star B, with a mass
$M_{\rm b}=35.2{\rm~{}M}_{\odot}$, undergoes direct core-collapse and forms a
BH of $31.7{\rm~{}M}_{\odot}$. COMPAS2 model creates a BH-BH system, which at
the second BH formation has a separation $a=1372.6{\rm~{}R}_{\odot}$ and an
eccentricity of $e=0.658$. This wide BH-BH system does not merge in a Hubble
time.

The evolution within COMPAS1 model is fairly analogous to COMPAS2 model.
However, higher wind mass-loss and less internal mixing leads to formation of
a much less massive BH-BH system. BHs in this model are $\sim
20{\rm~{}M}_{\odot}$ while in COMPAS2 they are $\sim 32{\rm~{}M}_{\odot}$. The
orbital separation at BH-BH formation is $a=740.6{\rm~{}R}_{\odot}$ while the
eccentricity remains virtually unchanged ($e=0.68$) resulting in a coalescence
time longer than the Hubble time.

Results of evolution in the COMPAS3 model are also broadly similar to those of
the COMPAS2 model (see Tab. 1). The additional condition of pre-RLOF orbital
circularization, however, changes a few key points. The first RLOF (stable,
TTMT, star A to star B) decreases the orbital separation from
$11165.6{\rm~{}R}_{\odot}$ to $739.2{\rm~{}R}_{\odot}$ and the eccentricity
becomes $e=0.0$. As star B evolves, stellar winds increase this orbital
separation to $a=1097.3{\rm~{}R}_{\odot}$ right before the second RLOF
(stable, TTMT, star B to BH A). At this point, the eccentricity is $e=0.023$,
increased from the previous circularization at the formation of the first BH.
The second RLOF decreases the separation to $a=501.9{\rm~{}R}_{\odot}$, and
again the orbit is circularized. The binary orbital separation at the time of
formation of the BH-BH system (with $\sim 32{\rm~{}M}_{\odot}$ BHs) becomes
$a=779.6{\rm~{}R}_{\odot}$ while the orbital eccentricity is negligible
($e=0.054$). Though COMPAS3 model evolution decreases the binary orbital
separation substantially, we note that this change is not significant enough
to create a close double BH system that merges in a Hubble time.

We note that COMPAS3 creates the closest BH-BH binary obtained in the three
COMPAS models. We note that to create a similarly circularized double BH
system with the same masses, that merges within a Hubble time, the orbital
separation at BH-BH formation can at most be about $46.7{\rm~{}R}_{\odot}$
(Peters, 1964). A highly eccentric orbit can also decrease the merger time.
However, for a BH-BH system with same masses and orbital separation as in the
COMPAS3 model, this cut-off eccentricity should be at least $e=0.98$.

Interestingly the structure of the BH X-ray binary Cyg X-1 (though a less
massive system than Mk42) was used in StarTrack (Wiktorowicz et al., 2014) and
COMPAS (Neijssel et al., 2021) calculations to argue that the future evolution
of Cyg X-1 (Miller-Jones et al., 2021) may also lead to a wide BH-BH system
which will not merge in a Hubble time.

###  3.3. MESA models

In this section we use the MESA code to check the outcome of the CE phase
encountered in the three StarTrack models from Sec. 3.1 (see Sec. 2.3.3 and
Klencki et al., 2021, for the method).

In the model MESA2 we evolve a star with $M_{\rm zams}=133{\rm~{}M}_{\odot}$,
overshooting $\sigma_{\rm ov}=0.2$ and Dutch winds with $f_{\rm wind}=1.5$.
This model produces at some point of its post-MS evolution a star with a mass
of $M=80.9{\rm~{}M}_{\odot}$ and a He-core mass of $M_{\rm
core}=62.2{\rm~{}M}_{\odot}$, radius of $R=700{\rm~{}R}_{\odot}$ $T_{\rm
eff}\approx 9000$K (an outer radiative envelope; the envelope would not become
convective until least $T_{\rm eff}\lesssim 4500$K Klencki et al. (2020)). At
this point we calculate the envelope binding energy (obtained from integration
of the mass distribution over the entire envelope) corrected for the internal
energy of the envelope, then we subtract the BH accretion luminosity that
effectively lowers binding energy, to obtain $E_{\rm bind}=6.67\times
10^{50}{\rm~{}erg}$. This translates to $\lambda_{\rm MESA}=0.012$. This model
resembles the star B at the onset of CE in the StarTrack2 example of
evolution. The orbital energy at the onset of CE is $E_{\rm orb,i}=-0.04\times
10^{50}{\rm~{}erg}$ and the post-CE separation is $a=6.2{\rm~{}R}_{\odot}$
(corresponding to a post-CE $E_{\rm orb,f}=-6.81\times 10^{50}{\rm~{}erg}$).
This was obtained under the assumption of a $100\%$ efficiency of the orbital
energy transfer to unbind the envelope ($\alpha=1.0$). The radius of the
exposed core of star B is $2.84{\rm~{}R}_{\odot}$ ( Hurley et al. (2000)
formulae) while its new Roche lobe is only $2.64{\rm~{}R}_{\odot}$ and we
assume a CE merger in such case. In the StarTrack2 model the binding energy
($\lambda=0.103$) was underestimated as it scales with $\propto 1/\lambda$ by
factor of $\sim 9$ as compared with detailed MESA estimate ($\lambda=0.012$).

Note that since the BH accretes part of the envelope during the CE inspiral
(for the numerical treatment of this process see Appendix in Belczynski et al.
(2002)), the binary does not need to balance the entire binding energy of the
envelope ($E_{\rm bind}$) with the orbital energy ($E_{\rm orb,f}-E_{\rm
orb,i}$). However, this has no influence on our conclusion above. Accretion
onto the BH is estimated at the level of $\sim 0.5{\rm~{}M}_{\odot}$ (see Sec.
2.1), while this MESA model underestimates the mass of the stellar envelope
found in StarTrack2 simulation by $\sim 7.4{\rm~{}M}_{\odot}$ ($M_{\rm
env}=26.1{\rm~{}M}_{\odot}$ in StarTrack2 simulation, and $M_{\rm
env}=18.7{\rm~{}M}_{\odot}$ in the above MESA simulation).

The other two population synthesis models from StarTrack also produce CE
mergers (Thorne-Żytkow objects) if MESA binding energy estimate is used. In
StarTrack1 model the CE donor was estimated to have $\lambda=0.050$, while the
MESA calculation gives $\lambda=0.007$. In StarTrack3 model $\lambda=0.050$ as
contrasted with MESA estimate of $\lambda=0.008$.

###  3.4. Pavlovskii models

A binary star resembling Mk 34 is evolved with the Pavlovskii2 model (Sec.
2.4). Evolution to RLOF initiated by star B is the same as in the StarTrack2
model ($\lesssim 1\%$ differences in binary parameters are numerical).
However, here the binary undergoes stable a RLOF: TTMT instead of CE. During
mass exchange mass/loss the orbital separation changes from $a=1507\rightarrow
928{\rm~{}R}_{\odot}$ and star B is stripped from its H-rich envelope ($M_{\rm
b}=87.7\rightarrow 62.3{\rm~{}M}_{\odot}$) becoming a massive Wolf-Rayet star.
Accretion onto the BH is negligible as the mass transfer was highly supper-
Eddington ($M_{\rm a}=35.2\rightarrow 35.2{\rm~{}M}_{\odot}$) and the loss of
angular momentum (given the mass ratio at the onset of RLOF:
$q=87.7/35.2=2.5$) causes orbit to decrease in size by a factor of $1.6$; this
may be compared with the orbital decrease by a factor of $63$ during CE in the
StarTrack2 model. After Wolf-Rayet wind mass loss ($M_{\rm
b}=33.9{\rm~{}M}_{\odot}$, $a=1427{\rm~{}R}_{\odot}$) star B collapses
directly to a BH ($M_{\rm b}=33.3{\rm~{}M}_{\odot}$). After $t_{\rm
evol}=4.16$Myr of binary evolution a wide BH-BH binary is formed with a
coalescence time of $t_{\rm coal}=6.0\times 10^{6}$Gyr.

In the Pavlovskii1 model, the binary follows a similar evolution, but a lower-
mass wide BH-BH binary forms ($21.9+22.1{\rm~{}M}_{\odot}$) due to stronger
winds and lower core masses adopted in this model (see Tab. 1).

In Pavlovskii3 model, due to strongly increased loss of the orbital angular
momentum during TTMT, the system does not survive the first TTMT. It ends with
a stellar merger of the Hertzsprung gap star donor ($M_{\rm
a}=80.2{\rm~{}M}_{\odot}$) with its main-sequence star companion ($M_{\rm
b}=98.9{\rm~{}M}_{\odot}$). The actual mass and the fate of the stellar-merger
product is uncertain. Both observations and simulations of stellar mergers are
usually related to low-mass stars, which are not BH progenitors (J. C.
Lombardi et al., 2002; Tylenda & Kamiński, 2016), or are calculated for
dynamical collisions in dense stellar clusters (Glebbeek et al., 2013). It
seems that a rather low mass fraction is lost during stellar mergers (J. C.
Lombardi et al., 2002; Lombardi et al., 2006; Glebbeek et al., 2013). Assuming
that the merger product in our simulation will become a Hertzsprung-gap star
with the mass of $\sim 163{\rm~{}M}_{\odot}$ (similarly to the scheme used in
Olejak et al. (2020) with $20\%$ of the less massive star being ejected during
the merger) the single star will end its evolution either as (i) PSN leaving
no remnant if classical PSN models are used (Woosley, 2017; Leung et al.,
2019), or (ii) as a single $\sim 30-40{\rm~{}M}_{\odot}$ BH, if non-standard
PSN models are used (see Fig.1 of Belczynski, 2020, and references therein).

The amount of angular-momentum loss through the L2 Lagrangian point adopted in
the Pavlovskii3 model, $j_{\rm loss}=j_{\rm L2}=1.2^{2}\frac{M_{\rm
tot}^{2}}{M_{\rm don}M_{\rm acc}}\in[5.76,5.82]$,333This range corresponds to
changing mass of donor and accretor during RLOF. is considered to be an upper
limit whereas the standard StarTrack $j_{\rm loss}=1.0$ used in models
Pavlovskii1 and Pavlovskii2 is instead close to the lower limit as indicated
by MacLeod et al. (2018) and MacLeod & Loeb (2020). The maximal possible
$j_{\rm loss}$ which allows to avoid a stellar merger during the first TTMT
and would lead to the formation of a wide BH-BH binary from Mk 34 is $j_{\rm
loss}\in[3.74,3.92]$ ($65\%$ of $j_{\rm L2}$ of MacLeod & Loeb (2020)). This
demonstrates that even with increased angular momentum losses, it is possible
to form either a wide BH-BH binary ($t_{\rm delay}>t_{\rm hub}$) with minimal
separation of about $a=609{\rm~{}R}_{\odot}$ ($t_{\rm delay}=2.6\times 10^{5}$
Gyr $>t_{\rm hub}$) or a stellar merger but not a close BH-BH system.

###  3.5. Quasi single star evolution models

Based on MESA models (see Sec. A), current literature and simple estimates, we
follow the future evolution of the Mk 34 binary with non-expanding stars. We
put the two stars on an eccentric ($e=0.68$) and wide orbit
($a=780{\rm~{}R}_{\odot}$). These stars lose $\gtrsim 100{\rm~{}M}_{\odot}$
during their MS life in stellar winds expanding the orbital separation
($a>1000{\rm~{}R}_{\odot}$; see Sec. A). At the post-MS closest encounter of
these two stars (periastron), Roche lobe radii of both components are $R_{\rm
lobe}>100{\rm~{}R}_{\odot}$. The radii of both stars are
$R<100{\rm~{}R}_{\odot}$ for many MESA models. There is no mass exchange
between the stars.

Depending on (i) the mass and core mass of non-expanding stellar models and
(ii) the mass (in reality central temperature and density) and the range
allowed for the onset of a PPSN/PSN, we can envision several different
outcomes of Mk 34’s future evolution.

If both stars have core masses as high as $61{\rm~{}M}_{\odot}$ and
$58{\rm~{}M}_{\odot}$ at TAMS (see MESA models with $\delta_{ov}=0.33$ and
$f_{\rm wind}=1.0$ in Sec. A), these cores will reach $65{\rm~{}M}_{\odot}$ at
the time of oxygen burning, which will then become explosive leading to PSN
(Woosley, 2017). Each star gets disrupted, leaving no compact object remnant
but producing luminous PSN supernova (model: QuasiSingle1; see Tab. 1).

There is a significant caveat to the above prediction. According to recent
studies (Woosley, 2017; Limongi & Chieffi, 2018; Farmer et al., 2020; Costa et
al., 2021; Farrell et al., 2021) very low metallicity stars can produce BHs
with mass as high as $\sim 80-90{\rm~{}M}_{\odot}$. But there are also
detailed MESA stellar evolutionary models that allow for the formation of BHs
with $\sim 70{\rm~{}M}_{\odot}$ avoiding a PPSN/PSN (Belczynski et al., 2020b)
at high metallicity. If this scenario is adopted then it is expected that
$60+60{\rm~{}M}_{\odot}$ wide BH-BH binary would form ((model: QuasiSingle2).

For lower mass stars/cores at TAMS, we expect avoiding a PPSN/PSN and we
predict the formation of a wide BH-BH binary. For example, if we take MESA
models with $\delta_{ov}=0.4$ and $f_{\rm wind}=1.5$ they will produce stars
with $M=32{\rm~{}M}_{\odot}$ ($M_{\rm core,TAMS}=30{\rm~{}M}_{\odot}$) and
$M=32{\rm~{}M}_{\odot}$ ($M_{\rm core,TAMS}=30{\rm~{}M}_{\odot}$). These stars
are not subject to PPSN/PSN and depending on the post-MS stellar wind mass
loss will form $\sim 30{\rm~{}M}_{\odot}$ BHs (model: QuasiSingle3). The
formation of a wide BH-BH binary, with coalescence time exceeding the Hubble
time, is predicted.

If we push MESA models even further to higher overshooting and stronger winds
($\delta_{ov}=0.5$ and $f_{\rm wind}=2.0$) we produce stars with
$M=21{\rm~{}M}_{\odot}$ ($M_{\rm core,TAMS}=19{\rm~{}M}_{\odot}$) and
$M=21{\rm~{}M}_{\odot}$ ($M_{\rm core,TAMS}=19{\rm~{}M}_{\odot}$) at TAMS.
This will also lead to the formation of a wide BH-BH binary but with $\sim
20{\rm~{}M}_{\odot}$ BHs at most (model: QuasiSingle4).

##  4\. Discussion

We have investigated the future evolutionary tracks and fate of the most
massive known binary system Mk 34. Several interesting possibilities seem to
exist (see Tab. 1). However, it is impossible to decide with certainty (due to
various stellar and binary physics uncertainties) which predicted fate is the
correct one (if any).

If very massive stars at LMC metallicity and with moderate rotation expand
during their post-MS evolution (expected for low overshooting) then we predict
the following evolution sequence for Mk 34:

$RLOF_{\rm A}\ \rightarrow BH_{\rm A}\ \rightarrow RLOF_{\rm B}\ \rightarrow
TZ_{\rm A+B}/BH_{\rm B}$ (1)

where indices “A” and “B” mark the more- and less-massive component of Mk 34
respectively, $BH_{\rm A/B}$ denotes the BH formation from a given component,
and $TZ_{\rm A+B}$ means the formation of a Thorne-Żytkow object from both
binary components in the second RLOF. The first RLOF (initiated by star A) is
always found to be stable (TTMT/NTMT), while the second RLOF (donor: star B)
can be either stable or dynamically unstable (CE). Additionally, the binary
system may not survive the first RLOF while both stars merge forming a single
star that will be either subject to PSN (no remnant) or will form a single BH:

$RLOF_{\rm A}\ \rightarrow single\ star\ \rightarrow PSN/BH$ (2)

If such massive stars do not expand (for example, because of significant
overshooting) the future evolutionary history proceeds without any binary
interaction:

$BH_{\rm A}/PPSN_{\rm A}/PSN_{\rm A}\ \rightarrow BH_{\rm B}/PPSN_{\rm
B}/PSN_{\rm B}$ (3)

Under very optimistic conditions (development and survival of CE initiated by
a massive star with a radiative envelope and with comparable mass companion)
Mk 34 may form a heavy BH-BH merger that would be a source of gravitational-
waves. Depending on our assumptions on mass loss and mixing in stellar
interiors we find formation of a $\sim 20+20{\rm~{}M}_{\odot}$ close BH-BH
system that resembles LIGO/Virgo detection of GW190408_181802
($24.6^{+5.1}_{-3.4}+18.4^{+3.3}_{-3.6}{\rm~{}M}_{\odot}$: Abbott et al.
(2021)) or $\sim 30+30{\rm~{}M}_{\odot}$ BH-BH system that would look similar
to GW150914 ($35.6^{+4.8}_{-3.0}+30.6^{+3.0}_{-4.4}{\rm~{}M}_{\odot}$: Abbott
et al. (2016)) or to GW190828_063405
($32.1^{+5.8}_{-4.0}+26.2^{+4.6}_{-4.8}{\rm~{}M}_{\odot}$: Abbott et al.
(2021)). If this does not work, the formation of such LIGO/Virgo BH-BH mergers
can still be obtained with more realistic CE input physics in the isolated
binary evolution (Belczynski et al., 2016; Spera et al., 2019; Patton et al.,
2021).

The detailed evolutionary estimates of a very massive star envelope binding
energy do not allow for a CE survival of the Mk 34 descendant binary even
under very conservative assumptions ($100\%$ orbital energy used to eject the
envelope with the help of internal energy of the gas and accretion luminosity
from the inspiralling BH). If this is taken into account, then instead of
forming a close BH-BH binary in CE scenario, we encounter the formation of a
Thorne-Żytkow object (BH sinks into the center of the post-MS massive star).
This single object would first appear as a post-MS massive star (most likely a
classical Wolf-Rayet star), that starts to expand, cooling off and getting
redder as the envelope puffs up in response to the BH inspiralling in a H-rich
envelope. Once the BH sinks into the star’s core (the majority of star mass at
this point) the accretion of helium is extremely rapid ($\gtrsim
1{\rm~{}M}_{\odot}\ {\rm min}^{-1}$, as in a collapsar engine; Fryer &
Woosley, 1998) and the core disappears and then the rest of the star is
accreted as well and the object disappears entirely from sky. Such a transient
should be visible in optical/infrared (initial expansion of the envelope
before it collapses onto the BH), although at this moment there are no
available calculations of the light curve or spectra for such heavy mergers
($\sim 100{\rm~{}M}_{\odot}$ post-MS star and $\sim 20-30{\rm~{}M}_{\odot}$
BH). Observationally, various red novae/transients were proposed to be the
outcome of CE mergers (Tylenda et al., 2005; Ivanova et al., 2013a; Kamiński
et al., 2015; MacLeod et al., 2017b). It may be even possible that such a
merger would lead to a gamma-ray burst (GRB). If there is enough angular
momentum in the He core and BH then such a configuration may lead to formation
of jets powering a GRB (Zhang & Fryer, 2001). The angular momentum transport
in stellar interiors of massive stars is not fully constrained, although low
effective spins of LIGO/Virgo BH-BH mergers seem to indicate efficient angular
momentum transport in massive stars (Spruit, 2002; Fuller et al., 2019; Bavera
et al., 2020; Belczynski et al., 2020a). In such a case a slowly spinning He
core–BH system would have only a small chance of producing a GRB. Another
obstacle in producing a GRB in this case is $\sim 30{\rm~{}M}_{\odot}$ of
H-rich envelope for jets to punch through. Yet, there are signs that such jets
do form in massive stars in CE mergers and they try to breakout from stellar
interiors (Thöne et al., 2015).

Another possibility for the future evolution of Mk 34 is to avoid the CE phase
entirely, even if both stars in this binary expand. Since the orbital
separation is not very large for this system, any RLOF encountered in the
evolution is bound to happen when the donor star is not too large (a radiative
envelope). Additionally, since the stars in Mk 34 are of similar mass, any
RLOF is not bound to happen at extreme mass ratio. Taking this into account,
this binary may evolve through two episodes of stable RLOF (first initiated by
initially more massive star, and then by the other star). Although such RLOF
episodes may decrease the orbital separation, such orbit shrinking will be not
large enough to lead to the formation of a close BH-BH binary (i.e., with a
merger time smaller than the Hubble time) for the initial binary configuration
of Mk 34. Instead, a wide BH-BH binary forms with separation as large as
$\gtrsim 1000{\rm~{}R}_{\odot}$. Such a descendant of Mk 34 cannot be a
LIGO/Virgo source, but could possibly be detected by microlensing
observations. The magnification of the source in such a microlensing event
would last months and would include bumps typical of a binary lens and would
be potentially detectable in LMC.

Finally, it is also possible that neither of the stars in Mk 34 will
experience any significant expansion in their post-MS evolution. In such a
case, Mk 34 expands due to wind mass loss from both stars and forms a wide
system when both components end their nuclear evolution (separation of $\sim
1000{\rm~{}R}_{\odot}$). Depending on the highly uncertain mixing physics and
not fully constrained nuclear reaction rates (overshooting, rotation,
convection, carbon-fusion) that set temperature/density in stellar cores,
massive components of Mk 34 may or may not be a subject to significant mass
loss associated with PPSN during oxygen burning. We predict the formation of a
wide BH-BH system with comparable-mass BHs in the mass range $\sim
40-90{\rm~{}M}_{\odot}$. Alternatively, both stars in Mk 34 may be subject to
a full-fledged pair-instability and get disrupted in luminous PSNe’s (Higgins
et al., 2021).

To summarize, we conclude that we cannot yet predict the fate of a massive
binary such as Mk 34. The involved stellar and binary physics uncertainties
are still too overwhelming. However, our study offers several conditional
statements that shed light on the future evolution of this massive binary. If
the stars in Mk 34 expand in their post-MS evolution then they are bound to
initiate two RLOF interactions (one by each component). Although the first
interaction is always stable and does not threaten the survival of this
system, the second one is more problematic. It is not at all clear if the
second interaction will be dynamically stable. If it is not, then a CE phase
develops and the most likely fate of the system is then the merger of the two
binary components, possibly associated with a red nova or a GRB. If the second
interaction is stable then the RLOF will lead to the formation of a wide and
massive BH-BH system with a merger time much larger than the Hubble time. Such
a system is not a potential LIGO/Virgo source of high-frequency gravitational
waves, but it may produce a microlensing event. If the stars in Mk 34 do not
expand, which is also allowed by the current detailed evolutionary models,
then we predict either the formation of a wide and potentially very massive
BH-BH system, or the spectacular death of both stars in luminous pair-
instability supernovae that leave no BHs behind.

##  5\. Conclusion

For any given origin scenario (i) there is a large number of input physics
uncertainties or even unknowns (model parameters), (ii) the implementation of
the physical processes that involve these uncertainties in the numerical codes
are far from being based on first-principle physics, so even probing the full
range of a given parameter might not get the right answer, and (iii) there are
more parameters and thus uncertainties than are commonly realized (for
example, there are at least $\sim 30$ parameters in the isolated binary
evolution, even though it is usually thought that a small subset of them are
the most important in determining binary outcomes). Thus strong conclusions
are unjustified at this time.

As an example of how biases could enter the model comparisons, suppose that we
use the rapid supernova engine model of Fryer et al. (2012), which naturally
produces a $\sim 2-5{\rm~{}M}_{\odot}$ mass gap between neutron stars and
black holes. Then the discovery of the $2.6{\rm~{}M}_{\odot}$ object in
GW190814 would rule strongly against isolated binary evolution and in favor of
another channel, such as primordial black holes or multiple-generation mergers
in dense stellar systems (e.g., two neutron stars could merge to make a $\sim
2.6{\rm~{}M}_{\odot}$ black hole). But perhaps the delayed supernova engine
model of Fryer et al. (2012) is a better description; in this model the
compact objects with mass $2.6{\rm~{}M}_{\odot}$ are produced naturally and
the isolated binary formation channel is perfectly viable. Or perhaps some
other supernova model is selected by Nature, which would change the Bayes
factor between the models that are considered. Similar considerations apply to
the high-mass merger GW190521, which is consistent with binary stellar
evolution given the substantial uncertainties (Farmer et al., 2020;
Belczynski, 2020; Kinugawa et al., 2020; Vink et al., 2021; Costa et al.,
2021; Mehta et al., 2021).

We have focused on the specific system Mk 34, but our caveat extends to
analyses of the full population. For example, Olejak et al. (2021) explored
the effects on the BH-BH population that stem from different treatments of the
common envelopes. Over the range of models they studied, the population
characteristics varied drastically. For example, the BH-BH merger rate varied
from $18{\rm~{}Gpc}^{-3}{\rm~{}yr}^{-1}$ (consistent with the current LVC
estimate of $15.3-38.8{\rm~{}Gpc}^{-3}{\rm~{}yr}^{-1}$: The LIGO Scientific
Collaboration et al. (2021)) to $88{\rm~{}Gpc}^{-3}{\rm~{}yr}^{-1}$. The BH
mass distribution can be consistent with the LVC estimate ($\propto M^{-1.5}$
below $\sim 40-50{\rm~{}M}_{\odot}$ and $\propto M^{-5.3}$ at higher masses),
or very inconsistent with the estimate ($\propto M^{+2.7}$ for
$M<15{\rm~{}M}_{\odot}$ and $\propto M^{-3.3}$ for heavier black holes). The
mass ratio distribution can have one peak or two peaks. To reiterate, even
this broad range of predicted population characteristics does not include many
other possible variations of aspects of stellar and binary evolution.

What must be done to reach a stage in which we can draw firm conclusions about
the BH-BH system origins? More and better data will obviously help: for
example, if multiple events point to a compact object in the $\sim
2-5{\rm~{}M}_{\odot}$ range then this tells us that the lower mass gap is not
a major feature of the mass distribution. Rare individual events, if they have
definitively established properties, can point to particular origins. For
example, an event with many cycles that is clearly highly eccentric would
favor a dynamical origin, and a compact object with a mass
$<0.5{\rm~{}M}_{\odot}$ would signify a primordial black hole. But we
emphasize that a detailed work on the physics of BH formation in each case is
essential: statistical analyses must be grounded in both thorough and accurate
physics and astrophysics.

We thank Lukasz Wyrzykowski, Tomasz Kamiński, Chris Fryer and Daniel Holz for
useful comments on the manuscript. KB, AR, AO acknowledge support from the
Polish National Science Center grant Maestro (2018/30/A/ST9/00050). DC and SS
acknowledge the support of the Australian Research Council Centre of
Excellence for Gravitational Wave Discovery (OzGrav), through project number
CE170100004. This work made use of the OzSTAR high performance computer which
is funded by Swinburne University of Technology and the National Collaborative
Research Infrastructure Strategy (NCRIS). MCM acknowledges support from NASA
ADAP grant 80NSSC21K0649. He performed part of his work on this paper at the
Aspen Center for Physics, which is supported by National Science Foundation
grant PHY-1607611. JPL was supported in part by a grant from the French Space
Agency CNES.

## References

  * Abbott et al. (2016) Abbott, B. P., et al. 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2020a) Abbott, R., et al. 2020a, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2020b) —. 2020b, ApJ, 896, L44
  * Abbott et al. (2021) —. 2021, Physical Review X, 11, 021053
  * Bavera et al. (2020) Bavera, S. S., et al. 2020, A&A, 635, A97
  * 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. (2016) Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016, 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. (2020a) Belczynski, K., et al. 2020a, A&A, 636, A104
  * Belczynski et al. (2020b) —. 2020b, ApJ, 890, 113
  * Bestenlehner (2020) Bestenlehner, J. M. 2020, MNRAS, 493, 3938
  * Bestenlehner et al. (2014) Bestenlehner, J. M., et al. 2014, A&A, 570, A38
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Bonanos et al. (2004) Bonanos, A. Z., et al. 2004, ApJ, 611, L33
  * Brott et al. (2011) Brott, I., et al. 2011, A&A, 530, A115
  * Castro et al. (2014) Castro, N., Fossati, L., Langer, N., Simón-Díaz, S., Schneider, F. R. N., & Izzard, R. G. 2014, A&A, 570, L13
  * Chattopadhyay et al. (2021) Chattopadhyay, D., Stevenson, S., Hurley, J. R., Bailes, M., & Broekgaarden, F. 2021, MNRAS, 504, 3682
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., Cantiello, M., Paxton, B., & Johnson, B. D. 2016, ApJ, 823, 102
  * Claret & Torres (2018) Claret, A., & Torres, G. 2018, ApJ, 859, 100
  * Costa et al. (2021) Costa, G., Bressan, A., Mapelli, M., Marigo, P., Iorio, G., & Spera, M. 2021, MNRAS, 501, 4514
  * de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
  * Dominik et al. (2012) Dominik, M., Belczynski, K., Fryer, C., Holz, D., Berti, B., Bulik, T., Mandel, I., & O’Shaughnessy, R. 2012, ApJ, 759, 52
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S., Fishbach, M., & Justham, S. 2020, arXiv e-prints, arXiv:2006.06678
  * Farr et al. (2018) Farr, B., Holz, D. E., & Farr, W. M. 2018, ApJ, 854, L9
  * Farr et al. (2017) Farr, W. M., Stevenson, S., Miller, M. C., Mandel, I., Farr, B., & Vecchio, A. 2017, Nature, 548, 426
  * Farrell et al. (2021) Farrell, E., Groh, J. H., Hirschi, R., Murphy, L., Kaiser, E., Ekström, S., Georgy, C., & Meynet, G. 2021, MNRAS, 502, L40
  * Fishbach & Holz (2020) Fishbach, M., & Holz, D. E. 2020, arXiv e-prints, arXiv:2009.05472
  * Franciolini et al. (2021) Franciolini, G., et al. 2021, arXiv e-prints, arXiv:2105.03349
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., Dominik, M., Kalogera, V., & Holz, D. E. 2012, ApJ, 749, 91
  * Fryer & Woosley (1998) Fryer, C. L., & Woosley, S. E. 1998, ApJ, 502, L9
  * Fuller & Ma (2019) Fuller, J., & Ma, L. 2019, ApJ, 881, L1
  * Fuller et al. (2019) Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS
  * Ge et al. (2015) Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2015, ApJ, 812, 40
  * Gilkis et al. (2021) Gilkis, A., Shenar, T., Ramachandran, V., Jermyn, A. S., Mahy, L., Oskinova, L. M., Arcavi, I., & Sana, H. 2021, MNRAS, 503, 1884
  * Glebbeek et al. (2013) Glebbeek, E., Gaburov, E., Portegies Zwart, S., & Pols, O. R. 2013, Monthly Notices of the Royal Astronomical Society, 434, 3497–3510
  * Gräfener & Hamann (2008) Gräfener, G., & Hamann, W. R. 2008, A&A, 482, 945
  * Grevesse et al. (1996) Grevesse, N., Noels, A., & Sauval, A. J. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 99, Cosmic Abundances, ed. S. S. Holt & G. Sonneborn, 117
  * Heger et al. (2000) Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
  * Higgins et al. (2021) Higgins, E. R., Sander, A. A. C., Vink, J. S., & Hirschi, R. 2021, MNRAS, 505, 4874
  * Humphreys & Davidson (1994) Humphreys, R. M., & Davidson, K. 1994, PASP, 106, 1025
  * 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
  * Ivanova et al. (2013a) Ivanova, N., Justham, S., Avendano Nandez, J. L., & Lombardi, J. C. 2013a, Science, 339, 433
  * Ivanova et al. (2013b) Ivanova, N., et al. 2013b, A&A Rev., 21, 59
  * 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
  * Kalogera & Webbink (1996) Kalogera, V., & Webbink, R. F. 1996, Astrophys. J., 458, 301
  * Kamiński et al. (2015) Kamiński, T., Mason, E., Tylenda, R., & Schmidt, M. R. 2015, A&A, 580, A34
  * 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. (2021) Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, A&A, 645, A54
  * Klencki et al. (2020) Klencki, J., Nelemans, G., Istrate, A. G., & Pols, O. 2020, A&A, 638, A55
  * Köhler et al. (2015) Köhler, K., et al. 2015, A&A, 573, A71
  * Leung et al. (2019) Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, arXiv e-prints
  * Limongi & Chieffi (2018) Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13
  * Lohr et al. (2018) Lohr, M. E., Clark, J. S., Najarro, F., Patrick, L. R., Crowther, P. A., & Evans, C. J. 2018, A&A, 617, A66
  * Lombardi et al. (2006) Lombardi, J. C., J., Proulx, Z. F., Dooley, K. L., Theriault, E. M., Ivanova, N., & Rasio, F. A. 2006, ApJ, 640, 441
  * MacLeod et al. (2017a) MacLeod, M., Antoni, A., Murguia-Berthier, A., Macias, P., & Ramirez-Ruiz, E. 2017a, ApJ, 838, 56
  * MacLeod & Loeb (2020) MacLeod, M., & Loeb, A. 2020, ApJ, 893, 106
  * MacLeod et al. (2017b) MacLeod, M., Macias, P., Ramirez-Ruiz, E., Grindlay, J., Batta, A., & Montes, G. 2017b, ApJ, 835, 282
  * MacLeod et al. (2018) MacLeod, M., Ostriker, E. C., & Stone, J. M. 2018, The Astrophysical Journal, 863, 5
  * Mahy et al. (2020) Mahy, L., et al. 2020, A&A, 634, A118
  * Mandel & Broekgaarden (2021) Mandel, I., & Broekgaarden, F. S. 2021, arXiv e-prints, arXiv:2107.14239
  * Marchant et al. (2021) Marchant, P., Pappas, K. M. W., Gallegos-Garcia, M., Berry, C. P. L., Taam, R. E., Kalogera, V., & Podsiadlowski, P. 2021, A&A, 650, A107
  * Marchant et al. (2018) Marchant, P., Renzo, M., Farmer, R., Pappas, K. M. W., Taam, R. E., de Mink, S., & Kalogera, V. 2018, arXiv e-prints
  * Massey et al. (2002) Massey, P., Penny, L. R., & Vukovich, J. 2002, ApJ, 565, 982
  * Mehta et al. (2021) Mehta, A. K., Buonanno, A., Gair, J., Miller, M. C., Farag, E., deBoer, R. J., Wiescher, M., & Timmes, F. X. 2021, arXiv e-prints, arXiv:2105.06366
  * Miller-Jones et al. (2021) Miller-Jones, J. C. A., et al. 2021, Science, 371, 1046
  * Mondal et al. (2020) Mondal, S., Belczyński, K., Wiktorowicz, G., Lasota, J.-P., & King, A. R. 2020, MNRAS, 491, 2747
  * Neijssel et al. (2021) Neijssel, C. J., Vinciguerra, S., Vigna-Gómez, A., Hirai, R., Miller-Jones, J. C. A., Bahramian, A., Maccarone, T. J., & Mandel, I. 2021, ApJ, 908, 118
  * Neijssel et al. (2019) Neijssel, C. J., et al. 2019, MNRAS, 490, 3740
  * Nitz & Capano (2021) Nitz, A. H., & Capano, C. D. 2021, ApJ, 907, L9
  * Nugis & Lamers (2000) Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227
  * Olejak et al. (2020) Olejak, A., Belczynski, K., Bulik, T., & Sobolewska, M. 2020, A&A, 638, A94
  * Olejak et al. (2021) Olejak, A., Belczynski, K., & Ivanova, N. 2021, arXiv e-prints, arXiv:2102.05649
  * Ossowski (2021) Ossowski, M. 2021, A&A, 649, A57
  * Patton et al. (2021) Patton, R. A., Sukhbold, T., & Eldridge, J. J. 2021, arXiv e-prints, arXiv:2106.05978
  * Pavlovskii & Ivanova (2015) Pavlovskii, K., & Ivanova, N. 2015, MNRAS, 449, 4415
  * Pavlovskii et al. (2017) Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., & Timmes, F. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) —. 2015, ApJS, 220, 15
  * Paxton et al. (2018) —. 2018, ApJS, 234, 34
  * Paxton et al. (2019) —. 2019, arXiv e-prints
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
  * Podsiadlowski et al. (1992) Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246
  * Pollock et al. (2018) Pollock, A. M. T., Crowther, P. A., Tehrani, K., Broos, P. S., & Townsley, L. K. 2018, MNRAS, 474, 3228
  * Pols et al. (1998) Pols, O. R., Schröder, K.-P., Hurley, J. R., Tout, C. A., & Eggleton, P. P. 1998, MNRAS, 298, 525
  * Rolleston et al. (2002) Rolleston, W. R. J., Trundle, C., & Dufton, P. L. 2002, A&A, 396, 53
  * Sana et al. (2012) Sana, H., et al. 2012, Science, 337, 444
  * Schneider et al. (2015) Schneider, F. R. N., Izzard, R. G., Langer, N., & de Mink, S. E. 2015, ApJ, 805, 20
  * Schnurr et al. (2008) Schnurr, O., Casoli, J., Chené, A. N., Moffat, A. F. J., & St-Louis, N. 2008, MNRAS, 389, L38
  * Schootemeijer et al. (2019) Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132
  * Scott et al. (2021) Scott, L. J. A., Hirschi, R., Georgy, C., Arnett, W. D., Meakin, C., Kaiser, E. A., Ekström, S., & Yusof, N. 2021, MNRAS, 503, 4208
  * Shenar et al. (2021) Shenar, T., et al. 2021, A&A, 650, A147
  * Spera et al. (2019) Spera, M., Mapelli, M., Giacobbo, N., Trani, A. A., Bressan, A., & Costa, G. 2019, MNRAS, 485, 889
  * Spruit (2002) Spruit, H. C. 2002, A&A, 381, 923
  * 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
  * Stevenson et al. (2017) Stevenson, S., Vigna-Gómez, A., Mandel, I., Barrett, J. W., Neijssel, C. J., Perkins, D., & de Mink, S. E. 2017, Nature Communications, 8, 14906, [Nature Commun.8,14906(2017)]
  * Tehrani et al. (2019) Tehrani, K. A., Crowther, P. A., Bestenlehner, J. M., Littlefair, S. P., Pollock, A. M. T., Parker, R. J., & Schnurr, O. 2019, MNRAS, 484, 2692
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration et al. 2021, ApJ, 913, L7
  * Thöne et al. (2015) Thöne, C. C., de Ugarte Postigo, A., Fryer, C. L., & Kann, D. A. 2015, in Extragalactic Jets from Every Angle, ed. F. Massaro, C. C. Cheung, E. Lopez, & A. Siemiginowska, Vol. 313, 396–397
  * Thorne & Zytkow (1977) Thorne, K. S., & Zytkow, A. N. 1977, ApJ, 212, 832
  * Tylenda & Kamiński (2016) Tylenda, R., & Kamiński, T. 2016, A&A, 592, A134
  * Tylenda et al. (2005) Tylenda, R., Soker, N., & Szczerba, R. 2005, A&A, 441, 1099
  * van den Heuvel et al. (2017) van den Heuvel, E. P. J., Portegies Zwart, S. F., & de Mink, S. E. 2017, MNRAS, 471, 4256
  * Vigna-Gómez et al. (2018) Vigna-Gómez, A., et al. 2018, MNRAS, 481, 4009
  * Vink et al. (2010) Vink, J. S., Brott, I., Gräfener, G., Langer, N., de Koter, A., & Lennon, D. J. 2010, A&A, 512, L7
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Vink et al. (2021) Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146
  * Vink et al. (2011) Vink, J. S., Muijres, L. E., Anthonisse, B., de Koter, A., Gräfener, G., & Langer, N. 2011, A&A, 531, A132
  * Vitale et al. (2017) Vitale, S., Lynch, R., Sturani, R., & Graff, P. 2017, Classical and Quantum Gravity, 34, 03LT01
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355
  * Wiktorowicz et al. (2014) Wiktorowicz, G., Belczynski, K., & Maccarone, T. 2014, in Binary Systems, their Evolution and Environments, 37
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244
  * Xu & Li (2010) Xu, X.-J., & Li, X.-D. 2010, ApJ, 722, 1985
  * Yusof et al. (2013) Yusof, N., et al. 2013, MNRAS, 433, 1114
  * 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, ApJ, 910, 152
  * Zhang & Fryer (2001) Zhang, W., & Fryer, C. L. 2001, ApJ, 550, 357

##  Appendix A Winds-overshooting grid

Table 4Maximum radius and stellar mass, core mass and envelope mass at TAMS
from different combinations of overshoot fractions and Dutch winds scale
factors for the $144{\rm~{}M}_{\odot}$ ZAMS star (left) and the
$131{\rm~{}M}_{\odot}$ ZAMS star (right).

model | $R_{Max}$ | $M_{TAMS}$ | $M_{core,TAMS}$ | $M_{env,TAMS}$
---|---|---|---|---
| [${\rm~{}R}_{\odot}$] | [${\rm~{}M}_{\odot}$] | [${\rm~{}M}_{\odot}$] | [${\rm~{}M}_{\odot}$]
$\delta_{ov}$ 0.12 |  |  |  |
$f_{\rm wind}$ 1 |  $>$1967.640 | 94.018 | 66.578 | 27.440
$f_{\rm wind}$ 1.5 | 39.694 | 56.672 | 55.423 | 1.249
$f_{\rm wind}$ 2 | 25.870 | 30.662 | 28.014 | 2.648
$\delta_{ov}$ 0.16 |  |  |  |
$f_{\rm wind}$ 1 |  $>$2000 | 93.139 | 68.443 | 24.696
$f_{\rm wind}$ 1.5 | 34.663 | 47.202 | 45.412 | 1.790
$f_{\rm wind}$ 2 | 24.192 | 27.978 | 25.227 | 2.750
$\delta_{ov}$ 0.2 |  |  |  |
$f_{\rm wind}$ 1 |  $>$1129.641 | 91.070 | 70.021 | 21.049
$f_{\rm wind}$ 1.5 | 31.430 | 41.892 | 39.779 | 2.113
$f_{\rm wind}$ 2 | 23.346 | 26.517 | 23.794 | 2.723
$\delta_{ov}$ 0.33 |  |  |  |
$f_{\rm wind}$ 1 | 46.046 | 62.221 | 60.782 | 1.438
$f_{\rm wind}$ 1.5 | 26.474 | 34.516 | 32.304 | 2.212
$f_{\rm wind}$ 2 | 21.062 | 23.249 | 20.738 | 2.511
$\delta_{ov}$ 0.4 |  |  |  |
$f_{\rm wind}$ 1 | 36.335 | 53.797 | 52.095 | 1.702
$f_{\rm wind}$ 1.5 | 24.508 | 32.205 | 30.090 | 2.115
$f_{\rm wind}$ 2 | 20.183 | 22.135 | 19.744 | 2.391
$\delta_{ov}$ 0.5 |  |  |  |
$f_{\rm wind}$ 1 | 30.383 | 48.080 | 46.399 | 1.682
$f_{\rm wind}$ 1.5 | 22.559 | 29.984 | 28.017 | 1.968
$f_{\rm wind}$ 2 | 19.189 | 20.896 | 18.737 | 2.159

model | $R_{Max}$ | $M_{TAMS}$ | $M_{core,TAMS}$ | $M_{env,TAMS}$
---|---|---|---|---
| [${\rm~{}R}_{\odot}$] | [${\rm~{}M}_{\odot}$] | [${\rm~{}M}_{\odot}$] | [${\rm~{}M}_{\odot}$]
$\delta_{ov}$ 0.12 |  |  |  |
$f_{\rm wind}$ 1 |  $>$2000 | 86.568 | 58.490 | 28.078
$f_{\rm wind}$ 1.5 | 40.419 | 62.672 | 52.644 | 10.028
$f_{\rm wind}$ 2 | 25.070 | 32.530 | 30.080 | 2.451
$\delta_{ov}$ 0.16 |  |  |  |
$f_{\rm wind}$ 1 |  $>$2000 | 85.196 | 60.295 | 24.900
$f_{\rm wind}$ 1.5 | 33.693 | 49.404 | 47.764 | 1.640
$f_{\rm wind}$ 2 | 23.654 | 29.393 | 26.787 | 2.606
$\delta_{ov}$ 0.2 |  |  |  |
$f_{\rm wind}$ 1 |  $>$2000 | 84.293 | 61.686 | 22.607
$f_{\rm wind}$ 1.5 | 29.989 | 42.995 | 40.967 | 2.029
$f_{\rm wind}$ 2 | 22.569 | 27.366 | 24.724 | 2.642
$\delta_{ov}$ 0.33 |  |  |  |
$f_{\rm wind}$ 1 | 41.232 | 59.079 | 57.538 | 1.541
$f_{\rm wind}$ 1.5 | 24.696 | 34.222 | 31.988 | 2.234
$f_{\rm wind}$ 2 | 20.022 | 23.455 | 20.948 | 2.507
$\delta_{ov}$ 0.4 |  |  |  |
$f_{\rm wind}$ 1 | 33.053 | 51.420 | 49.583 | 1.838
$f_{\rm wind}$ 1.5 | 23.301 | 32.049 | 29.895 | 2.154
$f_{\rm wind}$ 2 | 19.042 | 22.120 | 19.747 | 2.373
$\delta_{ov}$ 0.5 |  |  |  |
$f_{\rm wind}$ 1 | 27.484 | 45.711 | 43.895 | 1.816
$f_{\rm wind}$ 1.5 | 21.311 | 29.517 | 27.527 | 1.990
$f_{\rm wind}$ 2 | 17.986 | 20.735 | 18.555 | 2.181