#  POSYDON: A General-Purpose Population Synthesis Code with Detailed Binary-
Evolution Simulations

Tassos Fragos Département d’Astronomie, Université de Genève, Chemin Pegasi
51, CH-1290 Versoix, Switzerland  Jeff J. Andrews Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA), Northwestern University,
2145 Sheridan Road, Evanston, IL 60208, USA  Department of Physics, University
of Florida, 2001 Museum Rd, Gainesville, FL 32611, USA  Simone S. Bavera
Département d’Astronomie, Université de Genève, Chemin Pegasi 51, CH-1290
Versoix, Switzerland  Christopher P.L. Berry Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA), Northwestern University,
2145 Sheridan Road, Evanston, IL 60208, USA  Institute for Gravitational
Research, University of Glasgow, Kelvin Building, University Avenue, Glasgow,
G12 8QQ, Scotland  Scott Coughlin Center for Interdisciplinary Exploration and
Research in Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road,
Evanston, IL 60208, USA  Aaron Dotter Center for Interdisciplinary Exploration
and Research in Astrophysics (CIERA), Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, USA  Prabin Giri Department of Electrical and
Computer Engineering, Iowa State University, 2520 Osborn Dr, Ames, IA 50011,
USA  Vicky Kalogera Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road, Evanston,
IL 60208, USA  Department of Physics and Astronomy, Northwestern University,
2145 Sheridan Road, Evanston, IL 60208, USA  Aggelos Katsaggelos Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
Electrical and Computer Engineering, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, USA  Konstantinos Kovlakas Département d’Astronomie,
Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland  Shamal
Lalvani  Electrical and Computer Engineering, Northwestern University, 2145
Sheridan Road, Evanston, IL 60208, USA  Devina Misra Département d’Astronomie,
Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland  Philipp
M. Srivastava  Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road, Evanston,
IL 60208, USA  Electrical and Computer Engineering, Northwestern University,
2145 Sheridan Road, Evanston, IL 60208, USA  Ying Qin Department of Physics,
Anhui Normal University, Wuhu, Anhui 241000, China  Kyle A. Rocha Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
Department of Physics and Astronomy, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, USA  Jaime Román-Garza Universidad de Monterrey,
Ave. Morones Prieto 4500 Pte., C.P. 66283, San Pedro Garza García, Nuevo León,
México.  Juan Gabriel Serra  Center for Interdisciplinary Exploration and
Research in Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road,
Evanston, IL 60208, USA  Electrical and Computer Engineering, Northwestern
University, 2145 Sheridan Road, Evanston, IL 60208, USA  Petter Stahle
Département d’Informatique, Université de Genève, Route de Drize 7, CH-1227
Carouge, Switzerland  Meng Sun Center for Interdisciplinary Exploration and
Research in Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road,
Evanston, IL 60208, USA  Xu Teng Department of Electrical and Computer
Engineering, Iowa State University, 2520 Osborn Dr, Ames, IA 50011, USA  Goce
Trajcevski Department of Electrical and Computer Engineering, Iowa State
University, 2520 Osborn Dr, Ames, IA 50011, USA  Nam Hai Tran DARK, Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark
Zepei Xing Département d’Astronomie, Université de Genève, Chemin Pegasi 51,
CH-1290 Versoix, Switzerland  Emmanouil Zapartas Département d’Astronomie,
Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland  Michael
Zevin Kavli Institute for Cosmological Physics, The University of Chicago,
5640 South Ellis Avenue, Chicago, Illinois 60637, USA  Enrico Fermi Institute,
The University of Chicago, 933 East 56th Street, Chicago, Illinois 60637, USA

###### Abstract

Most massive stars are members of a binary or a higher-order stellar systems,
where the presence of a binary companion can decisively alter their evolution
via binary interactions. Interacting binaries are also important astrophysical
laboratories for the study of compact objects. Binary population synthesis
studies have been used extensively over the last two decades to interpret
observations of compact-object binaries and to decipher the physical processes
that lead to their formation. Here, we present , a novel, publicly available,
binary population synthesis code that incorporates full stellar-structure and
binary-evolution modeling, using the MESA code, throughout the whole evolution
of the binaries. The use of enables the self-consistent treatment of physical
processes in stellar and binary evolution, including: realistic mass-transfer
calculations and assessment of stability, internal angular-momentum transport
and tides, stellar core sizes, mass-transfer rates and orbital periods. This
paper describes the detailed methodology and implementation of , including the
assumed physics of stellar- and binary-evolution, the extensive grids of
detailed single- and binary-star models, the post-processing, classification
and interpolation methods we developed for use with the grids, and the
treatment of evolutionary phases that are not based on pre-calculated grids.
The first version of targets binaries with massive primary stars (potential
progenitors of neutron stars or black holes) at solar metallicity.

††journal: AAS Journals††software: This manuscript has made use of the
following Python modules: numpy (Harris et al., 2020), scipy (Virtanen et al.,
2020), pandas (McKinney et al., 2010), matplotlib (Hunter, 2007), astropy
(Robitaille et al., 2013; Price-Whelan et al., 2018), scikit-learn (Pedregosa
et al., 2011).

##  1 Introduction

Throughout their lives, stars affect their surroundings via the immense energy
radiated across the electromagnetic spectrum (e.g., Conroy, 2013; Eldridge &
Stanway, 2022) and the nuclear-processed material emitted as a stellar wind
(e.g., Kudritzki & Puls, 2000; Smith, 2014). The deaths of massive ($\gtrsim
8\,\rm M_{\odot}$) stars, even more than their lives, transform their
environments as their cores run out of nuclear fuel and collapse to form
neutron stars (NSs) and black holes (BHs). The formation of these compact
objects (COs) is often accompanied by a supernova (SN) or a $\gamma$-ray burst
that release more energy in $10\,\mathrm{s}$ than our Sun in
$10^{10}\,\mathrm{yr}$ (e.g., Janka, 2012; Burrows & Vartanyan, 2021; Woosley
& Bloom, 2006). These explosive events enrich their environments with heavier
elements while also regulating any ongoing star-formation (e.g., Nomoto et
al., 2013; Hopkins et al., 2011, 2012).

It is now established that most massive stars are members of a binary or a
higher-order stellar system (Sana et al., 2013; Moe & Di Stefano, 2017). More
often than not, the presence of a binary companion decisively alters the
evolution and final fate of both binary components via binary-interaction
processes such as tidal dissipation, mass-transfer phases, and stellar mergers
(Sana et al., 2012; De Marco & Izzard, 2017). Furthermore, interacting
binaries are arguably some of the most important astrophysical laboratories
available for the study of COs. Accretion of matter from a binary companion
gives rise to X-ray emission, bringing the system to the X-ray binary (XRB)
phase (Bhattacharya & van den Heuvel, 1991; Podsiadlowski et al., 1992), while
gravitational waves (GWs) enable us to witness the last moments of the lives
of coalescing binary COs (Abbott et al., 2016, 2021).

Over the last two decades, multi-wavelength surveys of the Milky Way and its
neighborhood, as well as numerous nearby and more distant galaxies, have
amassed large datasets of binary stellar systems. These datasets range from
targeted, high-resolution observations of galaxies in the local Universe
(e.g., the PHAT and LEGUS surveys; Dalcanton et al., 2012; Calzetti et al.,
2015), to serendipitous (e.g., the Chandra Source Catalogue 2.0; Evans et al.,
2019), all-sky (e.g., Gaia; Gaia Collaboration et al., 2018) and transient
surveys (e.g., Pan-STARRS, Zwicky Transient Facility; Kaiser et al., 2002;
Graham et al., 2019). In addition to electromagnetic surveys, there is also
the global GW observatory network of LIGO (Aasi et al., 2015), Virgo (Acernese
et al., 2015) and KAGRA (Akutsu et al., 2019), which have detected nearly a
hundred binary CO mergers (Abbott et al., 2021). Combined, these surveys are
revolutionizing our view of binary stellar systems, including CO binaries, and
their environments.

Aspects of the astrophysics of all these different types of stellar binaries
can be obtained from observations and modeling of present-day properties of
individual, well-studied systems. However, more comprehensive insight requires
understanding the statistical properties of their entire populations. For
these studies, binary population synthesis (BPS) modeling is often employed.
BPS modeling first generates initial binary populations, whose properties are
randomly sampled from probability distributions that can be observationally
constrained. Then, this initial population is evolved with a computationally
efficient simulation tool using our best understanding of the physics
dictating binary star interactions, to produce observable properties of the
target population. If the number of binaries evolved is large enough to
provide a statistically significant description of a population of interest,
then BPS can provide valuable insights about the expected rate and
distribution of the target population’s properties, the different evolutionary
pathways that lead to formation of these systems, and the effect that
different physical processes have on their evolution.

Over the last two decades, many general purpose BPS codes have been developed,
e.g., binary_c (Izzard et al., 2004, 2006, 2009), BPASS (Eldridge et al.,
2017), the Brussels code (Vanbeveren et al., 1998a, b), BSE (Hurley et al.,
2002), ComBinE (Kruckow et al., 2018), COMPAS (Stevenson et al., 2017; Riley
et al., 2022), COSMIC (Breivik et al., 2020), MOBSE (Giacobbo et al., 2018),
the Scenario Machine (Lipunov et al., 1996, 2009), SEVN (Spera et al., 2015),
SeBa (Portegies Zwart & Verbunt, 1996; Toonen et al., 2012a), StarTrack
(Belczynski et al., 2002, 2008), and TRES (Toonen et al., 2016). These have
been used in studies of a wide variety of binary populations. A hard
requirement for BPS is computational efficiency, as for most studies one would
need to model the evolution of many millions of binaries in a reasonable
computational time.

BPS codes stand in stark contrast to detailed stellar-structure and binary-
evolution codes, e.g., BEC (Heger et al., 2000; Heger & Langer, 2000), BINSTAR
(Siess et al., 2013), the Cambridge STARS code (Eggleton, 1971; Pols et al.,
1995; Eldridge & Tout, 2004; Stancliffe & Eldridge, 2009), MESA (Paxton et
al., 2015), and the TWIN code (Nelson & Eggleton, 2001; Eggleton & Kiseleva-
Eggleton, 2002), which self-consistently solve the stellar structure equations
of a binary’s component stars along with the orbital evolution. Many studies
have used detailed binary-evolution calculations (e.g., Nelson & Eggleton,
2001; Podsiadlowski et al., 2002; de Mink et al., 2007; Marchant et al., 2017;
Qin et al., 2018, 2019; Langer et al., 2020; Misra et al., 2020; Laplace et
al., 2020, 2021) to generate grids of models, varying the masses of the two
stars and the binary’s orbital period. However, in all those cases, the grids
of detailed binary tracks either cover a limited part of the initial parameter
space, or focus on a specific evolutionary phase. This limitation is
principally caused by the computational demands of detailed grids; each
simulation typically requires $\sim 10$–$100$ CPU hours for the modeling of a
single system (e.g., Paxton et al., 2019).

As a result of this computational expense, a common thread among the vast
majority of current BPS codes is that they approximate each star’s evolution,
employing either fitting formulae (e.g., SSE; Hurley et al., 2000) or look-up
tables (e.g., COMBINE; Kruckow et al., 2018) for the properties of _single
stars_ , based on grids of pre-calculated detailed, single-star models. Then,
the effects of binary interactions (e.g., Roche-lobe overflow or tides) are
modeled using approximate prescriptions and parametrizations. This modeling
approach is often called rapid or parametric BPS; throughout the remainder of
this work, we choose to use the term _parametric BPS_ (pBPS) modeling, to make
a distinction between computational efficiency and modeling accuracy. A
notable exception among BPS codes is BPASS (Eldridge et al., 2017), which uses
extensive grids of detailed binary evolution models computed with a custom
version of the Cambridge STARS binary evolution code (Stancliffe & Eldridge,
2009). In the grids of binary-star models employed in BPASS, both the primary
and the secondary stars are followed in detail, but only one at a time (for
computational-cost reasons). During the primary’s evolution, the properties of
the secondary star are approximated by formulae based on single-star models
(Hurley et al., 2000). Subsequently, once the modeling of the primary’s
evolution is completed, the secondary star’s evolution is re-computed,
accounting for mass-transfer and rejuvenation effects.

Studies using pBPS techniques have allowed us to make advances in our
understanding of the formation pathways leading to different types of binary
systems, and to interpret observations of binary populations. Such examples
are studies on white-dwarf binaries (e.g., Nelemans et al., 2001; Ruiter et
al., 2009, 2010; Toonen et al., 2012b; Breivik et al., 2018; Korol et al.,
2020), XRBs (e.g., Van Bever & Vanbeveren, 2000; Belczynski et al., 2004;
Fragos et al., 2008, 2013b, 2013a; Luo et al., 2012; Tzanavaris et al., 2013;
Tremmel et al., 2013; Zuo et al., 2014; Zuo & Li, 2014; van Haaften et al.,
2015; Artale et al., 2019; Wiktorowicz et al., 2019; Shao & Li, 2020), the
Galactic population of double NSs (DNSs; e.g., Osłowski et al., 2011; Andrews
et al., 2015; Chruslinska et al., 2017; Vigna-Gómez et al., 2018;
Chattopadhyay et al., 2020), GW sources observable by ground-based
observatories (e.g., Dominik et al., 2012, 2013, 2015; Mennekens & Vanbeveren,
2014, 2016; Belczynski et al., 2016; Klencki et al., 2018; Giacobbo et al.,
2018; Mapelli & Giacobbo, 2018; Neijssel et al., 2019; Spera et al., 2019;
Kinugawa et al., 2020; Breivik et al., 2020; Zevin et al., 2020; Broekgaarden
et al., 2021) and SNe in binary systems (e.g., De Donder & Vanbeveren, 2003;
Vanbeveren et al., 2013; Claeys et al., 2014; Zapartas et al., 2017, 2019,
2021).

However, the implicit assumption in pPBS codes that the binary components have
properties identical to single stars of the same mass in thermal equilibrium
(e.g., abundance profiles, core sizes, mass-radius relations and response to
mass-loss), as well as the lack of information about the star’s internal
structure at different critical evolutionary phases (e.g., the onset of a
dynamically unstable mass-transfer, the end of stable and unstable mass-
transfer phases or the core-collapse), may introduce systematic uncertainties
and inaccuracies. Current pBPS codes therefore rely on approximate
prescriptions for modeling binary interactions and difficult-to-calibrate
additional model parameters. These complications could be avoided by instead
employing detailed stellar-structure and binary-evolution simulations
(hereafter _detailed models_). Focusing on aspects that are relevant to the
formation of CO binaries, detailed models (i) allow for a self-consistent
estimation of the mass-transfer rate, especially during thermal-timescale
mass-transfer phases, and therefore an accurate assessment of mass-transfer
stability; (ii) allow for a more accurate description of the type and
properties of the formed CO as well as any potential associated transient
events since the internal structure of pre-core-collapse stars is known, (iii)
account for the transport of angular momentum between and within the binary
components, including its back-reaction on the structure and evolution of each
star (e.g., rotational mixing), and (iv) allow for the self-consistent
modeling of the end of a mass transfer phase (e.g., accounting for a potential
partial stripping of the envelope).

In this work, we build upon the combined experience gained from the large body
of BPS studies to date, to create (POpulation SYnthesis with Detailed binary-
evolution simulatiONs), a general-purpose code that can generate entire
populations of binaries, underpinned by detailed, self-consistent models of
stellar binaries.111 will become is publicly available at https://posydon.org
upon publication of the paper . With we aim to address many of the caveats of
pBPS codes, while at the same time maintaining much of their flexibility. In
its first release (v1.0), is limited to stars of solar metallicity, and
binaries where the primary star is massive enough to form a BH or a NS. Future
releases, which are already in development, will lift these limitations. In
Section 2, we introduce and the approach it takes to modeling binary
populations. In Sections 3 and 4 we provide the physics adopted for our
detailed models of single and binary stars, respectively. We describe the pre-
calculated grids of single and binary stellar evolution models in Section 5,
the way they are post-processed in Section 6, and our classification and
interpolation methods for their optimal use in Section 7. In Section 8 we
detail our treatment of evolutionary phases which are not based on pre-
calculated grids, such as the core-collapse and the common-envelope (CE)
phase, while in Section 9 we describe how all the aforementioned pieces come
together to model the entire evolution of a binary. In Section 10 we outline
our assumptions and methods in modeling populations of binary systems and
present some example results. We conclude in Section 11, where we present an
outlook of future development directions of the code. The definitions of all
the symbols used throughout this paper can be found in Table 1.

Table 1: List of variables used throughout the paper. Name | Description | First appears
---|---|---
$a$ |  Orbital separation |  4.1
$a_{i}$ |  Orbital separation before orbital kick |  8.3.5
$a_{f}$ |  Orbital separation after orbital kick |  8.3.5
$a_{\mathrm{pre,CE}}$ |  Orbital separation pre common envelope |  8.2
$a_{\mathrm{post,CE}}$ |  Orbital separation post common envelope |  8.2
$\dot{a}$ |  Rate of change of orbital separation |  8.1.2
$\dot{a}_{\mathrm{wind}}$ |  Rate of change of orbital separation due to wind mass loss |  8.1.2
$\dot{a}_{\mathrm{tides}}$ |  Rate of change of orbital separation due to tides |  8.1.2
$\dot{a}_{\mathrm{GR}}$ |  Rate of change of orbital separation due to gravitational-wave radiation |  8.1.2
$a_{\rm spin}$ |  Non-dimensional spin |  4.2.3
$c$ |  Speed of light |  4.2.2
$D_{\rm conv.reg.}$ |  Depth of a convective region |  7.4
$e$ |  Orbital eccentricity |  8.1.2
$e_{f}$ |  Orbital eccentricity after orbital kick |  8.3.5
$\dot{e}$ |  Rate of change of orbital eccentricity |  8.1.2
$\dot{e}_{\mathrm{tides}}$ |  Rate of change of orbital eccentricity due to tides |  8.1.2
$\dot{e}_{\mathrm{GR}}$ |  Rate of change of orbital eccentricity due to gravitational-wave radiation |  8.1.2
$E$ |  Eccentric anomaly |  8.3.5
$E_{\rm 2}$ |  Second-order tidal torque coefficient |  4.1
$f_{\rm conv}$ |  Dimensionless factor accounting for slow convective shells that cannot contribute to the tidal viscosity within an orbital timescale |  4.1
$f_{\mathrm{fb}}$ |  Fallback mass fraction |  8.3.2
$f_{\rm ov}$ |  Convective exponential overshooting parameter |  3.2.3
$g$ |  Local gravitational acceleration |  3.2.1
$G$ |  Gravitational constant |  3.2.2
$I$ |  Moment of inertia |  4.1
$\dot{I}$ |  Moment of inertia rate of change |  8.1.2
$j$ |  Specific angular momentum |  5.7
$j_{\mathrm{ISCO}}$ |  Specific angular momentum of the ISCO |  8.3.4
$J$ |  Angular momentum |  5.7
$J_{\mathrm{shell}}$ |  Stellar shell’s angular momentum |  8.3.4
$J_{\mathrm{direct}}$ |  Stellar shell’s angular momentum of directly collapsing material |  8.3.4
$J_{\mathrm{disk}}$ |  Stellar shell’s angular momentum of disk forming material |  8.3.4
$J_{\mathrm{BH}}$ |  Black hole angular momentum |  8.3.4
$k$ |  Apsidal motion constant |  4.1
$L$ |  Star’s luminosity |  3.2.2
$L_{\rm Edd}$ |  Eddington luminosity |  3.2.2
$L_{2}$ |  Second lagrange point |  8.2
$M_{1}$ |  Mass of the initially more massive star |  5.5
$M_{2}$ |  Mass of the initially less massive star |  5.5
$M_{\rm acc}$ |  Mass of the accretor |  4.2.2
$M_{\rm don}$ |  Mass of the donor |  8.2
$M_{\rm CO}$ |  Mass of the compact object |  5.6
$M_{\rm{C/O\mbox{-}core}}$ |  Mass of the C/O core |  7.4
$M_{\rm conv.reg.}$ |  Mass of a convective region |  4.1
$M_{\rm comp}$ |  Mass of the binary companion star |  4.1
$M_{\rm disk}$ |  Mass of accretion disk |  8.3.4
$M_{\rm{env}}$ |  Mass of the stellar envelope |  7.4
$M_{\rm{grav}}$ |  Remnant’s gravitational mass |  8.3.3
$M_{\rm{He\mbox{-}core}}$ |  Mass of the He core |  7.4
$M_{\rm{rembar}}$ |  Remnant’s baryonic mass |  7.4
$M^{\mathrm{max}}_{\rm{NS}}$ |  Maximum neutron star mass |  8.3.3
$\dot{M}_{\rm Edd}$ |  Mass-accretion rate corresponding to the Eddington limit |  4.2.2
$\dot{M}_{\rm w}$ |  Wind mass-loss rate |  3.2.2
$M_{\mathrm{tot}}^{\mathrm{i}}$ |  Binary stellar mass before core collapse |  8.3.5
$M_{\mathrm{tot}}^{\mathrm{f}}$ |  Binary stellar mass after core collapse |  8.3.5
$m_{\rm{shell}}$ |  Stellar shell’s mass |  8.3.4
$P$ |  Pressure |  3.2.1
$P_{\rm orb}$ |  Orbital period |  5.5
$q$ |  Binary mass ratio |  4.1
$R$ |  Stellar radius |  3.2.2
$R_{\rm acc}$ |  Radius of the accretor |  4.2.2
$R_{\rm b,conv.reg.}$ |  Radial coordinate of a convective region’s bottom boundary |  4.1
$R_{\rm{core}}$ |  Radius of the stellar core |  7.4
$R_{\rm{C/O\mbox{-}core}}$ |  Radius of the C/O core |  7.4
$R_{\rm conv}$ |  Radius of the convective core |  4.1
$R_{\rm conv.reg.}$ |  Radial coordinate of a convective region’s center |  7.4
$R_{\rm don}$ |  Radius of the donor |  8.2
$R_{\rm{He-core}}$ |  Radius of the He core |  7.4
$R_{\rm L}$ |  Roche lobe radius |  4.2.2
$R_{\rm L,acc}$ |  Roche lobe radius of the accretor |  4.2.2
$R_{\rm t,conv.reg.}$ |  Radial coordinate of a convective region’s top boundary |  4.1
$r$ |  Stellar shell’s radius |  8.3.4
$r_{i}$ |  Instantaneous orbital separation before orbital kick |  8.3.5
$T_{\rm eff}$ |  Effective temperature |  3.2.1
$T$ |  Timescale for orbital changes due to tides |  4.1
$v_{k}$ |  Magnitude of velocity kick |  8.3.5
$v_{r}$ |  Orbital velocity of the collapsing star |  8.3.5
$X$ |  Hydrogen mass function |  4.2.2
$X_{\mathrm{center}}$ |  Center hydrogen mass function |  8.1.1
$X_{\mathrm{surf}}$ |  Surface hydrogen mass function |  8.1.1
$Y$ |  Helium mass fraction |  3.1
$Y_{\mathrm{center}}$ |  Center helium mass fraction |  8.1.1
$Y_{\mathrm{surf}}$ |  Surface helium mass fraction |  7.5
$Z$ |  Metallicity (mass fraction of elements heavier than 4He) |  3.1
$\alpha_{\rm CE}$ |  Fraction of the orbital energy that contributes to the unbinding of the CE |  8.2
$\alpha_{\rm MLT}$ |  Convective mixing length parameter |  3.2.3
$\alpha_{\rm th}$ |  Thermohaline mixing parameter |  3.2.3
$\eta$ |  Dimensionless factor denoting the radiative efficiency of the accretion process |  4.2.2
$\theta$ |  Stellar profile’s polar angle |  8.3.4
$\theta_{\mathrm{disk}}$ |  Stellar profile’s polar angle of disk formation |  8.3.4
$\kappa$ |  Opacity |  3.2.1
$\lambda_{\rm CE}$ |  Parametrization of the CE’s binding energy |  8.2
$\tau$ |  Optical depth |  3.2.1
$\tau_{\rm sync}$ |  Tidal synchronization timescale |  4.1
$\tau_{\rm conv}$ |  Convective timescale |  4.1
$\tau_{\rm mb}$ |  Magnetic braking torque |  8.1.2
$\nu$ |  Reduced mass |  8.1.2
$\sigma$ |  Maxwellian distribution dispersion |  8.3.5
$\sigma_{\mathrm{CCSN}}$ |  Maxwellian distribution dispersion for CCSN kicks |  8.3.5
$\sigma_{\mathrm{ECSN}}$ |  Maxwellian distribution dispersion for ECSN kicks |  8.3.5
$\psi$ |  Binary orbital inclination with respect to before the kick |  8.3.5
$\omega_{\rm s}$ |  Surface angular velocity |  3.2.2
$\omega_{\rm s,crit}$ |  Critical surface angular velocity |
$\Omega_{\mathrm{orb}}$ |  Orbital angular velocity |  8.1.2
$\Omega_{\rm{shell}}$ |  Stellar shell’s angular velocity |  8.3.4
$\Omega$ |  Stellar angular velocity |  8.1.2
$\dot{\Omega}$ |  Stellar angular velocity rate of change |  8.1.2
$\dot{\Omega}_{\mathrm{wind}}$ |  Stellar angular velocity rate of change due to winds |  8.1.2
$\dot{\Omega}_{\mathrm{inertia}}$ |  Stellar angular velocity rate of change due to changes of star’s moment of inertia |  8.1.2
$\dot{\Omega}_{\mathrm{tides}}$ |  Stellar angular velocity rate of change due to tides |  8.1.2
$\dot{\Omega}_{\mathrm{mb}}$ |  Stellar angular velocity rate of change due to magnetic breaking |  8.1.2 3.2.2

##  2 Overview of the structure of

At its core, a BPS code requires two elements: a method to generate random
binaries at zero-age main sequence (ZAMS) and a mechanism to evolve those
binaries. The former is described in Section 10.1. Regarding the evolution of
each binary, our primary goal with is to self-consistently evolve the internal
structures of the two stars comprising a binary along with the binary’s orbit.
To achieve this goal, we have opted to employ the stellar-structure code
Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al., 2011,
2013, 2015, 2018, 2019). However, calculating the evolutionary track of a
binary using MESA can take in excess of $100$ CPU hours, a six-orders-of-
magnitude increase in computational cost compared to a typical pBPS code such
as COSMIC (Breivik et al., 2020). This means that, even with modern computing
resources, we cannot feasibly run more than $\sim 10^{5}$–$10^{6}$ binaries,
far too few to accurately model a Milky Way sized population with $\sim
10^{11}$ stellar binaries. As a further complication, codes like MESA cannot
run individual binaries from start to finish; key physics, including CE phases
and supernovae require the code to be stopped and restarted.222In principle
these phases could be run within a stellar-structure code like MESA; for
instance, more updated versions of MESA than the one we use can handle the
evolution of a binary through a CE (Marchant et al., 2021).

Figure 1: The structure of (v1.0) for modeling the evolution of a binary star.
Rectangles represent the initial and possible final outcomes of the evolution,
and black circles represent events in the evolution of a binary. Colored lines
showcase the different evolutionary steps that follows. Evolutionary steps
that are based on pre-calculated grids of detailed binary-evolution tracks are
designated with solid lines, while those that are based on computations
performed on-the-fly for each modeled binary are shown with dashed lines.

solves the problems associated with stellar-structure codes by using
extensive, pre-calculated grids of single and binary stellar-evolution models,
covering the parameter space relevant for the formation of high-mass binary
stars, with separate grids being calculated for each phase of binary
evolution. In v1.0 our grids contain a combined total of nearly 120,000
separate detailed binary simulations. To compute these grids, has an
infrastructure specifically designed for high-performance computing
environments that streamlines the process of producing large grids with
consistent physics inputs. Using this infrastructure, we have generated five
separate grids of single and binary stellar evolution models. We computed
three grids of interacting binary stars initially composed of two hydrogen
(H)-rich ZAMS stars (Section 5.5); a CO and a H-rich star at the onset of
Roche-lobe overflow (RLO; Section 5.6), and a helium (He)-rich ZAMS star with
a CO companion (Section 5.7). We further computed two grids of single H-rich
and He-rich stars (Sections 5.3 and 5.4, respectively), which we use for the
modeling of detached, non-interacting binaries. These five grids are then
post-processed, so that their data size is reduced (Section 6). We
additionally apply classification and interpolation algorithms on the outputs
of these extensive grids (Section 7), allowing us to effectively interpolate
between MESA simulations to estimate the evolution of any arbitrary binary
within some bounded region of the parameter space. As a simpler alternative,
we also provide functionality to evolve individual binaries using nearest-
neighbor matching, and in that case no classification or interpolation methods
are required.

The second major component of is the code infrastructure to follow the entire
evolution of a binary from start to end. To achieve this, we combined the
aforementioned grids (and classification and interpolation methods) with
physics dictating a binary’s evolution through key phases, including core
collapse and CO formation (Section 8.3) and CE (Section 8.2). These latter
phases are not modeled based on pre-calculated grids, but rather with on-the-
fly calculations. Similarly, the evolution of detached, potentially eccentric,
post-core-collapse binaries is also modeled with on-the-fly calculations,
where we use the single-star grids coupled to binary evolution routines, i.e.
orbital evolution due to tides, stellar winds, magnetic breaking and
gravitational-wave emission (Section 8.1).

In the approach, each separate evolutionary phase (step) has its own dedicated
function which determines the binary’s state resulting from that step, the
quantitative values characterizing that binary (e.g., masses of the two
stars), and the event describing how that state ended (e.g., onset of RLO).
Once a step is completed, the framework uses the resulting binary state and
event as well as each component stars’ states to determine an individual
binary’s next evolutionary step. The process is repeated until a binary’s
evolution is complete, resulting in a disrupted binary, a binary merger, or a
double CO. At this point, the next binary is run. The modular nature of allows
a user to also provide their own prescriptions to model each phase of
evolution, or even their own breakdown of the binary-evolution tree. As
default in , we provide a complete set of evolutionary steps which we visually
summarize in Figure 1 and present in detail in the following sections.

##  3 Adopted stellar physics

All stellar-evolution models described in this paper were computed with the
state-of-the-art, open source stellar-structure and evolution code MESA
(Paxton et al., 2011, 2013, 2015, 2018, 2019) revision 11701 together with the
20190503 version of the MESA software development kit (SDK; Townsend,
2020).333We made one minor bug fix in the MESA source code, which involves
replacing the mass of the proton with the atomic mass unit where it appears in
the code that evaluates the Potekhin & Chabrier (2010) equation of state (E.
Bauer, 2020, private communication). This change is included in later MESA
releases. MESA solves the one-dimensional stellar-structure and composition
equations. Mixing and burning processes are solved simultaneously; mixing is
treated as a diffusive process. Discussion of specific elements of stellar
physics are described in the following subsections, which are split into
microphysical and macrophysical processes. We implement all physics that are
not readily available in MESA using the functionality provided by
run_star_extras and run_binary_extras.

###  3.1 Microphysics

We adopt the Asplund et al. (2009) protosolar abundances as our initial
composition, with $Z=0.0142$ and $Y=0.2703$. The equation of state is the
standard MESA amalgamation of the SCVH (Saumon et al., 1995), OPAL (Rogers &
Nayfonov, 2002), HELM (Timmes & Swesty, 2000) and PC (Potekhin & Chabrier,
2010) equations of state (Paxton et al., 2019). Radiative opacities are taken
from Ferguson et al. (2005) and Iglesias & Rogers (1996) for the Asplund et
al. (2009) mixture, along with electron conduction opacities from Cassisi et
al. (2007). Nuclear reaction rates are drawn from the JINA Reaclib database
(Cyburt et al., 2010). All models were computed using the approx21 nuclear
reaction network that consists of 21 species: 1H, 3He, 4He, 12C, 14N, 16O,
20Ne, 24Mg, 28Si, 32S, 36Ar, 40Ca, 44Ti, 48Cr, 56Cr, 52Fe, 54Fe, 56Fe and
56Ni, plus protons and neutrons (for the purpose of photo-disintegration).

###  3.2 Macrophysics

####  3.2.1 Surface Boundary Conditions

The stellar surface boundary condition is satisfied by the simple_photosphere
option, which sets the photosphere temperature using the Eddington (1926)
$T_{\rm eff}(\tau)$ relation, and the photosphere pressure via $P=\tau
g/\kappa$ with enhancement due to radiation pressure at the photosphere (e.g.,
Paxton et al., 2011).

####  3.2.2 Stellar Winds

Stellar winds are a complex subject due to the varied physical mechanisms
(known and unknown) that drive them and their dependence on the evolutionary
state of their parent star. With this in mind, we have kept the wind
prescription as simple as possible, while still capturing the key
phenomenology of massive stellar evolution, but avoiding fine-tuning to
reproduce any single subset of observations. In , changing the wind
prescription would require the computation of new set of single- and binary-
star model grids.

For stars with initial masses above $8\,M_{\odot}$ we use the MESA Dutch
scheme, which consists of de Jager et al. (1988) for ${T_{\rm eff}}<10,000$ K
and Vink et al. (2001) for ${T_{\rm eff}}>11,000$ K. In cases where ${T_{\rm
eff}}>11,000$ K and the surface 1H mass fraction is below $0.4$, the Vink et
al. (2001) wind is replaced with the Wolf–Rayet wind of Nugis & Lamers (2000).
Between 10,000 K and 11,000 K there is a linear ramp (as a function of
${T_{\rm eff}}$) between the two wind prescriptions. We do not explicitly
include any luminous blue variable (LBV) type winds; however, our stellar
models at solar metallicity (e.g., Figure 3) do not enter the regime where
LBV-type winds are typically applied in other studies (e.g., Belczynski et
al., 2010).

For stars with initial masses below $8\,M_{\odot}$, we again use the Dutch
scheme for stars with ${T_{\rm eff}}$ hotter than 12,000 K. For stars with
${T_{\rm eff}}$ less than 8,000 K, we use the Reimers (1975) wind with scaling
factor $\eta_{R}=0.1$ for stars on the first ascent of the giant branch, and
the Bloecker (1995) wind with scaling factor $\eta_{B}=0.2$ for stars in the
thermally-pulsating phase. For the case of 8,000 K $<{T_{\rm eff}}<$ 12,000 K,
we calculate the wind for both the hot and cool schemes, and linearly
interpolate between the two.

For mass-loss rates that have an explicit dependence on metallicity $Z$, we
rescale wind mass loss based on the _initial_ metallicity, not the current
surface $Z$, as winds are driven predominantly by iron-group elements that
remain almost constant throughout stellar evolution (e.g., Vink & de Koter,
2005). The primary motivation for this approach is to avoid the dredge-up of
carbon and oxygen to the surface layers in the later phases of evolution,
which can cause surface $Z$ to approach $1$, from unduly influencing the mass-
loss rate. The only exception here is the wind prescription by Nugis & Lamers
(2000) for Wolf–Rayet stars, which is specifically calibrated to the total
surface metal content, including carbon and oxygen: in this case, we use the
current, surface $Z$ value of the He-rich star.

We further boost stellar winds to limit a star’s rotation below its critical
threshold $(\omega_{\rm s}/\omega_{\rm s,crit}\leq 1)$, so that the sum of the
centrifugal force and the photon pressure never exceeds gravity on the surface
of the star. The impact of rotation on the mass loss rate is considered as
indicated in (Heger & Langer, 1998; Langer, 1998),

$\centering\dot{M}_{\rm w}(\omega)=\dot{M}_{\rm
w}(0)\left(\frac{1}{1-\omega_{\rm s}/\omega_{\rm
s,crit}}\right)^{\xi},\@add@centering$ (1)

where $\dot{M}_{\rm w}$ is the star’s wind mass-loss rate, and $\omega_{\rm
s}$ and $\omega_{\rm s,crit}$ are the angular velocity and critical angular
velocity at the surface, respectively. The default value of the exponent
$\xi=0.43$ is taken from Langer (1998). The critical angular velocity is given
the expression $\omega_{\rm s,crit}^{2}=(1-L/L_{\rm Edd})GM/R^{3}$, where
$L_{\rm Edd}$ is the Eddington luminosity and its expression is given in Eq.
(8). This explicit boost to the wind is supplemented by an implicit numerical
scheme implemented in MESA which ensures that the rotation of a star never
exceeds its critical value.

####  3.2.3 Convection, Rotation, and Mixing Processes

Convective energy transport is modeled using mixing length theory (MLT; Böhm-
Vitense, 1958) except in superadiabatic, radiation-dominated regions where we
employ the MLT++ modifications introduced in MESA (Paxton et al., 2013) that
reduce the superadiabaticity in radiation-dominated convective regions, to
improve numerical convergence. For the condition of convective neutrality we
use the Ledoux criterion, and we use the convective premixing scheme as
described by Paxton et al. (2019). We adopt a solar-calibrated mixing length
parameter, $\alpha_{\rm MLT}=1.93$, based on results from the MIST project
(Dotter et al., in preparation).

Rotation is implemented in MESA as described in Paxton et al. (2013, 2019).
Rotational mixing and angular-momentum transport follow the MIST project (Choi
et al., 2016). It has been suggested that magnetic angular-momentum transport
processes are main candidates for efficient coupling between the stellar core
and its envelope during the post-MS (post-main sequence). Here, we adopt the
Spruit–Tayler (ST) dynamo (Spruit, 2002) that can be produced by differential
rotation in the radiative layers and amplify a seed magnetic field. Stellar
models with ST dynamo can reproduce the flat profile of the Sun (Eggenberger
et al., 2005) and observations of the final spins of both white dwarfs and
neutron stars (Heger et al., 2005; Suijs et al., 2008), but struggles to
explain the slow rotation rates of cores in red giants (Eggenberger et al.,
2012; Cantiello et al., 2014; Fuller et al., 2019).

MESA treats mixing processes in the diffusive approximation with MLT providing
the basic description. In addition to MLT convection, we consider thermohaline
mixing with the parameter $\alpha_{\rm th}=17.5$ (Paxton et al., 2013, Eq.
14), also referred to as $C_{\rm t}$ (Charbonnel & Zahn, 2007, Eq. 4),
corresponding to an aspect ratio $\sim 1$ of the instability fingers
(Kippenhahn et al., 1980). Thermohaline mixing is important during mass
accretion from an evolved primary star onto an unevolved secondary star
because the accreting material typically has a higher mean molecular weight
than the material near the surface of the secondary star (e.g., Kippenhahn et
al., 1980).

Overshoot mixing is treated in the exponential decay formalism (Herwig, 2000;
Paxton et al., 2011). For the parameter $f_{\mathrm{ov}}$ describing the
extent of the overshoot mixing in this formalism, we adopt an initial-mass-
dependent relation. For lower mass stars (initial masses less than 4
$M_{\odot}$) we adopt a value taken from the MIST project,
$f_{\mathrm{ov}}=0.016$, which is calibrated using the Sun, as well as open
clusters (Choi et al., 2016). In the high-mass regime (initial masses greater
than 8 $M_{\odot}$), we adopt a value of $f_{\mathrm{ov}}=0.0415$ motivated by
the work of Brott et al. (2011), who used the step overshoot formalism. Both
of these values of $f_{\mathrm{ov}}$ are measured from a distance of $0.008$
the local pressure scale height into the convection zone from the formal
convective-radiative boundary. This is the same approach adopted in the MIST
models (Choi et al., 2016). In order to translate between the step and
exponential-decay versions of overshoot mixing, we rely on the work of Claret
& Torres (2017) which shows that the free parameter in the step formalism is a
factor of $\sim 10$ larger than $f_{\mathrm{ov}}$ (their Figure 3). For stars
with initial masses between $4\,M_{\odot}$ and $8\,M_{\odot}$ we smoothly ramp
between the two values of $f_{\mathrm{ov}}$. The mass range was chosen to be
roughly consistent with the ranges considered in the two studies.

We include no extra mixing due to semiconvection (in the sense of Langer et
al., 1983), as this process is implicitly accounted for in the convective
premixing scheme (Paxton et al., 2019).

##  4 Adopted binary-star evolution physics

is predominantly a BPS tool that simulates the evolution of an ensemble of
binary systems through various stages of their life. In the framework, we base
the evolution of binary systems on three extended MESA binary grids, as shown
in Figure 1. One grid consists of initially detached binary systems of two
H-rich stars starting from ZAMS, where we follow the internal evolution of
both stars with detailed models (Section 5.5). A second grid consists of
H-rich stars in a semi-detached system with a CO companion (Section 5.6), and
a third grid consisting of naked helium stars in an initially detached system
with a CO companion (Section 5.7). For the non-CO components in these
binaries, we follow the same prescriptions for stellar structure and evolution
as described in Section 3. However, the internal structures of stars can be
affected by the presence of a companion, principally through tidal
interactions and mass transfer. In this section we describe how we use the
binary module within MESA to model each binary’s orbit, while self-
consistently accounting for the impact on each star’s structure.

###  4.1 Tides

Tidal forces take place in binary systems, as each star tends to be deformed
by the gravitational pull of its companion. Invoked by this gravitational
deformation, frictional forces inside a star drive a binary toward
circularization and stellar spin - orbit synchronization. In our MESA binary
grids, we assume that the initial orbit is circularized and the stellar spins
are synchronized with the orbit. This assumption should be valid especially
for close orbits of massive stars, where tides are strong (Portegies Zwart &
Verbunt, 1996; Hurley et al., 2002). As binaries evolve, their orbital periods
change as well as the individual stars’ rotation periods, potentially driving
them out of synchronization. Therefore, it is only the process of spin-orbit
coupling that is relevant for our grids of detailed binary-star models. In
Section 8.1 we discuss our treatment of eccentric, detached binaries.

We follow the linear approach to tides, which defines a timescale for
synchronization (Hut, 1981). In this approach, a torque is applied to non-
degenerate stars in a binary corresponding to the difference between the
orbital and the spin angular velocity $\Delta\Omega$, divided by the
synchronization timescale $\tau_{\rm sync}$:

$\delta\Omega=\frac{\Delta\Omega}{\tau_{\rm sync}}\delta\tau,$ (2)

where $\delta\tau$ is the timestep and $\delta\Omega$ is the change in spin
angular velocity over a particular timestep. Since every layer of a star
rotates with its own angular frequency, Eq. (2) is separately applied to every
layer. The torque applied to each layer of the star is added up and an
opposite torque is applied to the binary’s orbit to ensure angular momentum
conservation. Winds and mass transfer somewhat complicate the picture, and
Paxton et al. (2015) gives a detailed description of how these effects are
accounted for.

We separately calculate $\tau_{\rm sync}$ for both stars (Hut, 1981; Hurley et
al., 2002; Paxton et al., 2015):

$\frac{1}{\tau_{\rm
sync}}=3\left(\frac{k}{T}\right)q^{2}\frac{MR^{2}}{I}\left(\frac{R}{a}\right)^{6}.$
(3)

Here $M$, $R$ and $I$ are the mass, radius and moment of inertia of the star
for which we calculate the tides, $q=M_{\rm comp}/M$ is the binary mass ratio,
and $a$ is the orbital separation. $k$ is a dimensionless apsidal motion
constant characterizing the central condensation of the star, and $T$ is the
characteristic timescale for the orbital evolution due to tides. As Eq. (3)
shows, $\tau_{\rm sync}$ is strongly dependent on the ratio of the stellar
radius to the binary orbital separation. In practice, the quantity $k/T$ also
varies significantly, depending on whether tidal dissipation occurs
principally within convective regions (due to turbulent friction) or radiative
regions (due to dynamical tides interacting with stellar oscillations). As
stars may have both convective and radiative regions during their lives, at
every timestep taken by MESA we separately calculate the dynamical and
equilibrium tidal timescales, layer by layer, and apply the shorter of the
two.

For radiative regions in a star, we calculate $k/T$ based on the dynamical
tidal timescale from Zahn (1977, see also ), where

$\left(\frac{k}{T}\right)_{\mathrm{rad}}=\sqrt{\frac{GMR^{2}}{a^{5}}}(1+q)^{5/6}E_{2},$
(4)

with $E_{2}$ as the second order tidal coefficient and $G$ the gravitational
constant.444This equation is equivalent to Eq. (42) of Hurley et al. (2002),
apart from the typo correction of a square root, found in Sepinsky et al.
(2007). For the calculation of $E_{2}$, we adopt the latest prescriptions from
Qin et al. (2018), who investigated the dependence of the parameter on the
convective radius $R_{\rm conv}$ for various metallicities and evolutionary
stages, finding

$E_{2}=\begin{cases}10^{-0.42}\left(R_{\rm conv}/R\right)^{7.5}&\text{for
hydrogen-rich stars}\\\ 10^{-0.93}\left(R_{\rm conv}/R\right)^{6.7}&\text{for
stripped-helium stars.}\\\ \end{cases}$ (5)

For the equilibrium tidal timescale, we calculate the synchronization
timescale for each convective region in the stellar envelope using Eq. (3) and
following Hurley et al. (2002, cf. Eq. 30):

$\left(\frac{k}{T}\right)_{\mathrm{conv}}=\frac{2}{21}\frac{f_{\rm
conv}}{\tau_{\rm conv}}\frac{M_{\rm conv.reg.}}{M},$ (6)

and use the shortest timescale among them. In Eq. (6), $M_{\rm conv.reg.}$ is
the mass of the convective region, $f_{\rm conv}$ is a non-dimensional
numerical factor less than unity that takes into account slow convective
shells that cannot contribute to the tidal viscosity within an orbital
timescale, and $\tau_{\rm conv}$ is the convective timescale, which we take
from Eq. (31) of Hurley et al. (2002, based on ), adapted to also accomodate
convective regions that are below the surface:

$\displaystyle\tau_{\rm conv}$ $\displaystyle=$ $\displaystyle
0.431\left[\frac{M_{\rm conv.reg.}}{3L}\frac{R_{\rm t,conv.reg.}+R_{\rm
b,conv.reg.}}{2}\right.$ (7) $\displaystyle\qquad\times\left.(R_{\rm
t,conv.reg.}-R_{\rm b,conv.reg.})\vphantom{\frac{M_{\rm
conv.reg.}}{3L}}\right]^{1/3}.$

In the equation above, $R_{\rm t,conv.reg.}$, $R_{\rm b,conv.reg.}$, and $L$
are the radii of the top and bottom boundaries of the region, and the stellar
luminosity respectively, in solar units. We use the surface luminosity in all
calculations, as it is approximately constant throughout the envelope.
Typically, the shortest equilibrium tidal timescale corresponds to the
outermost convective region. In order to avoid tides being dominated by
potential artificial convective shells that may appear during the numerical
calculation of a star’s evolution, we only take into account regions that
consist of at least $10$ consecutive shells in our models.

Figure 2: Evolution of the donor star in two example close binary systems of
initially $M_{1,\rm{initial}}=56.46\,M_{\odot}$ and
$M_{2,\rm{initial}}=28.23\,M_{\odot}$, for two different initial orbital
periods: $3.16$ days (top) and $31.62$ days (bottom). The colors show which
tidal timescale is shortest and dominates in the tidal process: equilibrium
(blue) or dynamical (orange). Star symbols depict the main evolutionary
points.

In Figure 2, we show two example evolutionary tracks in a Hertzsprung–Russell
Diagram of the primary star in an interacting binary of
$M_{1,\rm{initial}}=56.46\,M_{\odot}$ and
$M_{2,\rm{initial}}=28.23\,M_{\odot}$, for two different initial orbital
periods ($3.16$ days, top panel; $31.62$ days, bottom panel). We show which
term of the tidal timescale dominates the tidal forces: the dynamical tidal
timescale with orange, from Eq. (3) and Eq. (4), assuming the whole star is
radiative, or the equilibrium tidal timescale with blue, from Eq. (3) and Eq.
(6), according to the most important convective region of the star. We see
that in the beginning of the evolution, the dynamical tidal timescale
dominates, as expected for the MS of massive stars that have a radiative
envelope. The $3.16$ day period system (top) initiates early RLO and does not
form convective layers massive enough for equilibrium tides to dominate, until
after the end of its MS. For the wider binary, even during the MS, the
equilibrium tidal timescale tends to become comparable to the dynamical tidal
timescale, due to convective regions that appear close to the surface of the
star. These regions include a small part of the total mass of the star (as low
as $10^{-5}M_{\odot}$), but have a significant radial thickness. The
equilibrium tidal timescale dominates in all the remaining parts of the
evolution, apart from the He core burning phase of the stripped primary in the
3 days period system.

###  4.2 Mass-transfer

Over the course of a binary’s evolution, the outermost layers of one of the
binary’s stellar components may be removed, due to the gravitational pull of
its companion. As a consequence of either the binary’s orbit decay or the
expansion of a star’s envelope, this transfer of mass is dictated by the
geometry of the Roche potential, namely the gravitational potential
constructed in the co-rotating reference frame of the binary system.

####  4.2.1 Mass Loss Rates from a Star Overflowing its Roche Lobe

To calculate mass-loss rates (due to mass transfer only) from main-sequence
(MS) stars that overfill their Roche lobes, we use the contact scheme within
MESA. This prescription is a numerical approximation; for stars overfilling
their Roche lobes at the beginning of each timestep, mass is removed such that
by the end of the timestep, the star remains confined to within its Roche
lobe. For MS stars, this approximation is consistent with more accurate
methods, as MS stars are compact, with relatively small pressure scale
heights. We choose this prescription, as it allows us to evolve binary systems
in which both stars overfill their Roche lobes simultaneously (Marchant et
al., 2016).

As stars evolve off the MS, however, they tend to expand, forming less dense
envelopes as they become giant stars. The large pressure scale heights of
giant stars cause the contact scheme to become inaccurate, and a prescription
is required that can more accurately treat these stars’ extended envelopes.
For stars with a central H abundance less than $10^{-6}$, we switch to the
Kolb scheme (Kolb & Ritter, 1990).555The current version of MESA does not
allow the reversing of the donor star in this scheme. Occasionally, an once
accreting star evolves off the MS, expands, and itself overfills its Roche
lobe. In these cases mass transfer is not calculated, and the system
predominantly leads to L2 overflow.  This prescription allows the star to
expand beyond its Roche lobe, and self-consistently calculates the rate that
mass can flow through the inner Lagrangian point based on the local fluid
conditions.

####  4.2.2 Mass Accretion onto a Non-Degenerate Companion

The evolution of a binary during a mass transfer phase depends not only on the
mass-losing star but also on the mass-gaining star. Based on the nature of the
accretor, the process of accretion is treated differently.

For binaries with a non-degenerate accretor (those in our grid of two H-rich
stars), initially all the mass lost by the donor through RLO is accepted by
the accretor. We assume that material being accreted carries the specific
angular momentum according to de Mink et al. (2013, Appendix A.3.3). This
prescription allows for the distinct treatment of accretion via direct impact
of the incoming stream on the stellar surface, or, in case of accretion onto a
more compact star, the formation of a Keplerian disk around the accretor. The
accreted angular momentum spins up the accretor, and mass accretion is
restricted when the accretor reaches critical rotation. Mass falling within a
critically rotating accretor’s gravitational potential will be ejected from
the binary with the specific angular momentum of the accretor (Paxton et al.,
2015) in the form of rotationally-enhanced stellar winds, following Eq. (1).
At the same time that accretion is spinning up the outer layers of a non-
degenerate star, internal mixing processes transport the surface angular
momentum toward deeper layers, slowing the star’s rotation rate.

####  4.2.3 Accretion onto a Degenerate Companion

Mass transfer onto a degenerate star proceeds similarly as that onto a non-
degenerate star, with a few notable exceptions. The primary exception is that
mass transfer is capped at the Eddington limited rate. For sub-Eddington mass-
transfer rates onto a CO, mass transfer is assumed to be conservative.
However, for super-Eddington rates, the excess matter is lost from the
vicinity of the accretor as an isotropic wind (i.e., with the specific angular
momentum of the accretor).

We calculate the Eddington-limited rate using standard formulae (Frank et al.,
2002). We first calculate the Eddington luminosity $L_{\rm Edd}$ for an
accretor with mass $M_{\rm acc}$:

$L_{\rm Edd}=\frac{4\pi GM_{\rm acc}c}{\kappa},$ (8)

where $\kappa$ is the opacity of the incoming material and $c$ is the speed of
light. For a fully ionized gas, Thompson scattering dominates the opacity, so
$\kappa=0.2(1+X)\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$, where $X$ is the hydrogen
abundance of the donor. By setting $L_{\rm Edd}$ equal to the radiation
released by accreted matter as it falls into a CO’s potential well ($L_{\rm
acc}=\eta\dot{M}c^{2}$), we can recover the Eddington-limited accretion rate,

$\dot{M}_{\rm Edd}=\frac{4\pi GM_{\rm acc}}{\kappa c\eta}.$ (9)

The dimensionless constant $\eta$ sets how efficiently the rest mass energy of
the incoming matter is converted to outgoing radiation,

$\eta=\frac{GM_{\rm acc}}{R_{\rm acc}c^{2}}.$ (10)

For BHs, $R_{\rm acc}$ is set by the spin-dependent innermost stable circular
orbit, while for NSs, we use a constant $R_{\rm acc}$ of $12.5$ km (Most et
al., 2018; Miller et al., 2019; Riley et al., 2019; Landry et al., 2020;
Abbott et al., 2020a; Kim et al., 2021; Biswas, 2021; Raaijmakers et al.,
2021). Our grids containing a CO components, are currently focused on NS and
BH accretors, for which we simulate a range of masses. The type of a CO is
determined solely on its mass, with COs having gravitational mass less than
$2.5\,M_{\odot}$ being classified as NSs, while those with mass greater than
$2.5\,M_{\odot}$ as BHs. Within our code, the only difference between these
types of accretors is the corresponding $\eta$; otherwise accretion proceeds
identically regardless of the type of CO accretor.

As these COs accrete material, they ought to accrete angular momentum. In our
current version of the grids, we ignore any corresponding increase in the spin
rate of NSs. For BHs, on the other hand, we self-consistently incorporate the
increase in spin frequency as well as its effect on $\eta$ through the radius
of innermost stable circular orbit (ISCO; Podsiadlowski et al., 2003),

$\eta=1-\sqrt{1-\bigg{(}\frac{M_{\rm acc}}{3M^{i}_{\rm acc}}\bigg{)}^{2}},$
(11)

where $M^{i}_{\rm acc}$ is the initial mass of the accreting BH and $M_{\rm
acc}$ its current one. In this equation it is also implicitly assumed that the
birth spin of the BH is $\sim$0, as it has been suggested by several studies
(Fragos & McClintock, 2015; Qin et al., 2018; Fuller & Ma, 2019). The
corresponding increase in the BH’s non-dimensional spin rate, $a_{\rm spin}$,
can be calculated following Thorne (1974) and King & Kolb (1999),

$a_{\rm spin}=\left(\frac{2}{3}\right)^{1/2}\frac{M^{i}_{\rm acc}}{M_{\rm
acc}}\left\\{4-\left[18\left(\frac{M^{i}_{\rm acc}}{M_{\rm
acc}}\right)^{2}-2\right]^{1/2}\right\\}.$ (12)

These equations both assume that the spin-up occurs due to angular momentum
accretion from a disk that is truncated at the ISCO.

We do not explicitly stop our simulations if a NS accretes enough mass to
cross our 2.5 $M_{\odot}$ threshold thereby collapsing into a BH, but we do
switch the $\eta$ instantaneously.

####  4.2.4 Onset of CE

Stars with radiative envelopes entering RLO respond to mass loss by shrinking
(Hjellming & Webbink, 1987); mass transfer then reaches a natural equilibrium
set by the strength of some driving force (nuclear evolution, tides, thermal
expansion or some other effect), and how quickly the orbital separation, and
thus the Roche lobe radius, change due to mass transfer through the inner
Lagrangian point. However, in certain circumstances mass transfer increases in
a runaway process, either because stars expand due to mass loss (e.g., stars
with deep convective envelopes) or because the binary’s orbit shrinks faster
than a donor star’s radius (e.g., Paczyński & Sienkiewicz, 1972). These phases
of binary evolution are notoriously difficult to model as they are
intrinsically three-dimensional processes, and they span many orders of
magnitude in spatial and temporal scales (Ivanova et al., 2013). We therefore
stop our MESA models when binaries enter dynamically unstable mass transfer;
we provide a description of how we address this phase in Section 8.2. Here we
focus on the conditions we use to identify when a binary enters dynamically
unstable mass transfer.

First, we assume that a dynamically unstable RLO phase is initiated whenever
mass-transfer rate exceeds $0.1\,M_{\odot}$ yr-1. It is expected that binaries
reaching this limit will only further increase their mass-transfer rates, as
this corresponds to a dynamical limit on the mass-loss rate for giant stars
(with dynamical timescales of years). As a check, we carried out a calibration
test where we followed the evolution of a test binary to mass transfer rates
even larger than $0.1\,M_{\odot}$ yr-1. In every test we ran, we found that
the mass-transfer rate increases to arbitrarily high rates, confirming the
validity of our limit. Assigning a limit to the mass-transfer rate also avoids
numerical issues caused by the effort of stellar models to converge with such
extreme mass loss.

As a second condition, we assume dynamically unstable RLO occurs when the
stellar radius of the expanding star extends beyond the gravitational
equipotential surface, passing through the second Lagrangian point (L2). In
such cases the lost matter from the L2 point carries substantial angular
momentum, rapidly shrinking the orbit and leading to a runaway process in
which the two stars spiral in and trigger a CE (Tylenda et al., 2011; Nandez
et al., 2014). We use the prescription from Misra et al. (2020, cf. Eq. 15–19)
to define the spherical-equivalent radius corresponding to L2. This condition
cannot occur for MS donors, since the contact scheme for RLO forces a star’s
radius to be contained within its Roche lobe. For cases where two MS stars
overfill both their Roche lobes, in an over-contact binary, we alternatively
use the prescription from Marchant et al. (2016, cf. Eq. 2) for the L2 radius,
which considers that both stars can contribute to the overflow of the L2
volume together.

Table 2:  Summary of the five detailed single- and binary-star model grids.  Initial state |  | Parameters’ range and resolution |  |
---|---|---|---|---
Star 1 | Star 2 |  | $M_{1}\,[M_{\odot}]$ | $\Delta\log_{10}M_{1}$ | $M_{2}\,[M_{\odot}]$ | $\Delta\log_{10}M_{2}$ | $q$ | $\Delta q$ | $P_{\rm orb}\,[{\rm day}]$ | $\Delta\log_{10}P_{\rm orb}$ |  $N$ aaTotal number of models in this grid. | Failures bbPercentage of models that stopped due to numerical-convergence errors before reaching one of our stopping conditions. These rates describe the finalized grids, after a series of re-runs have occurred; see Section 6.1 for details.
ZAMS | - |  | 0.5–300 | 0.014 | - | - | - | - | - | - | 200 | 1.5%
ZAHeMSccZero-age He Main Sequence stars. | - |  | 0.5–80 | 0.055 | - | - | - | - | - | - | 40 | 0%
ZAMS | ZAMS |  |  7 6.2–120 | 0.025 | - | - | 0.05–0.95 | 0.05 | 0.72–6105 | 0.07 |  56000 58240 |  0.9 1.5%
Evolved, H-rich ddAlthough this grid is initialized with H-rich stars at ZAMS, we ignore the portion of each simulated binary’s evolution prior to the onset of RLO. The initial state of Star 1 in this grid are therefore somewhat evolved.  | CO |  | 0.5–120 | 0.06 | 1–35.88 | 0.074 | - | - | 1.26-3162 | 0.13 | 25200 |  0.8 0.9%
ZAHeMS | CO |  | 0.5–80 | 0.055 | 1–35.88 | 0.074 | - | - | 0.02–1117.2 | 0.09 | 39480 |  4.7 4.8%

For CO accretors, we set a third condition for unstable mass transfer based on
the photon trapping radius (Begelman, 1979; King & Begelman, 1999). Inside
that radius photons are advected inward along with accreted matter onto the
accretor, while outside that radius, photons diffuse away. For stable mass
accretion, the photon trapping radius occurs close to the accretor; however,
as the accretion rate increases, the photon trapping radius expands. Once the
photon trapping radius reaches the Roche-lobe radius of the accretor it is
assumed to lead to a CE phase. Since the radius of the photon trapping
envelope $R_{\rm trap}$ depends on the Eddington limit of the accretor
$\dot{M}_{\rm Edd}$ and the mass-transfer rate from the donor $\dot{M}_{\rm
donor}$, we limit the latter assuming an instability condition when (Begelman,
1979):

$\dot{M}_{\rm donor}\geq{\dot{M}_{\rm Edd}}\frac{2R_{\rm L,accr}}{R_{\rm
acc}}.$ (13)

As a final condition, occasionally two stars will both overfill their Roche
lobe while one of those stars has evolved off the MS. Since the contact scheme
in MESA can only evolve binaries in which both stars are on the MS, we assume
these binaries automatically enter a CE.

As a test we compared the first three dynamical instability conditions
separately to investigate their effect and found that the limiting mass
accretion rates are all similar: as soon as a binary reaches any one of them,
the other two are close to their limits as well. Therefore, for a particular
binary in our MESA simulation, the binary is considered to enter a CE if any
one of them occurs.

##  5 Grids of detailed single- and binary-star evolution models

While in Section 3 and Section 4 we describe the physics we adopt in our
simulations, here we provide numerical details about how we produce each of
our five MESA grids. This includes our procedure for producing initial stellar
models for each of our grids (Section 5.1), our termination conditions common
across all of our grids (Section 5.2), and a description for each of our five
grids of binary simulations (Sections 5.3–5.7). We summarize the basic
properties of each grid in Table 2.

###  5.1 Zero Age Main Sequence Models

We create our own library of ZAMS models for both H-rich and He-rich stars.
For the creation of the H-rich ZAMS models we use the MESA revision 11701
template create_zams. The process begins with creating a fully-convective star
with no nuclear fusion taking place and adopting our protosolar abundances
(Section 3.1). This model is then evolved with our adopted nuclear reaction
network until the H-burning luminosity exceeds 99% of the total luminosity.

The He-rich ZAMS (ZAHeMS) models are created in three steps. First, we create
a pre-MS He star with 100% 4He in the same way that we create a H-rich pre-MS
star. In the second step we adjust the initial metallicity. In the third step
we evolve the model until the He-burning luminosity exceeds 99% of the total
luminosity.

The two sets of ZAMS and ZAHeMS models are used as a starting point for the
five grids of single- and binary-star models.

###  5.2 Termination conditions

We set conditions for the termination of our evolutionary models based on both
single-star properties and binary-star properties. If any one of these
conditions are met by an individual simulation, it is terminated at that
timestep. Our termination conditions are:

  * •

A star’s age exceeds the age of the Universe (13.8 Gyr), a condition that is
typically only met for the lowest-mass stars we simulate ($M_{\rm
init}\lesssim 0.8\,\mathrm{M}_{\odot}$). In our single-star grids, for
numerical purposes, we allow stars to evolve beyond this condition, then
truncate their evolution afterwards at the end of the MS, which may extend
beyond the age of the Universe.

  * •

A star becomes a WD, a condition we quantify by checking if the central
degeneracy parameter $\Gamma_{c}$ (Coulomb coupling parameter) exceeds 10
(Choi et al., 2016).

  * •

A star reaches the end of core C-burning, a condition triggered when the
fractional abundances of both C and He decrease below $10^{-2}$ and $10^{-6}$,
respectively, at the star’s center.

  * •

A binary enters a CE phase, as described in Section 8.2.

  * •

A star reaches the thermally-pulsating asymptotic giant branch (TP-AGB) phase
and then reaches a point of failure (numerical non-convergence) during thermal
pulsations. Because these are not uncommon and we consider the nascent WD to
be well-formed be within the AGB star, we consider this evolution to be
successful.

Our simulations may occasionally end prematurely before any of the
aforementioned conditions are reached. This may happen because the minimum
timestep limit within MESA ($10^{-6}$ s) is reached or any individual
simulation reaches our maximum runtime on our computing cluster (set to 48
hours). We provide details describing how we approach such runs in Section
6.1, but these failures are rare, occurring a few percent or less in each
grid.

All single- and binary-star models in have a final, internal structutre
profile written to correspond precisely to the last evolutionary phase
milestone. These profiles are discussed further in Section 6.2, Section 8.2,
and Section 8.3.

###  5.3 H-rich, single-star grid

Our first grid of single-star evolutionary models contains a series of non-
rotating H-stars with our adopted protosolar composition of $Y=0.2703$ and
$Z=0.0142$. The grid consists of 200 masses, ranging from
$M_{\mathrm{init}}=0.5\,M_{\odot}$ to $M_{\mathrm{init}}=300\,M_{\odot}$ with
a logarithmic spacing of $\Delta\log_{10}(M_{\mathrm{init}}/M_{\odot})=0.014$
dex. For each star, models were initialized using the procedure described in
Section 5.1, and evolved until one of the termination conditions provided in
Section 5.2 occurs.

Figure 3: Comparing H-rich ZAMS evolutionary tracks (blue in all panels) with
non-rotating Geneva (upper; Ekström et al., 2012), MIST (middle; Choi et al.,
2016), and SSE with default wind mass loss prescription from COSMIC version
3.4.0 (lower; Pols et al., 1998; Hurley et al., 2000; Breivik et al., 2020)
tracks in the Hertzsprung–Russell diagram. The mass range shown is 1–300
$M_{\odot}$ in all cases. The masses shown are the same as the Geneva grid of
models between 1 and 120 $M_{\odot}$ with the addition of 175 and 300
$M_{\odot}$ for POSYDON, MIST, and COSMIC/SSE. The dotted, gray lines indicate
constant radius at powers of 10 in $R_{\odot}$.

To test their validity, we compare evolutionary tracks to the widely used
stellar evolution tracks from the Geneva (upper panel; Ekström et al., 2012),
MIST library (center panel; Choi et al., 2016), and BSE as implemented in
COSMIC (lower panel; Pols et al., 1998; Hurley et al., 2000; Breivik et al.,
2020) groups in Figure 3. In all cases we show non-rotating models with
initial masses between 1 $M_{\odot}$ and 300 $M_{\odot}$. The pre-MS evolution
is omitted from the MIST evolutionary tracks and the TP-AGB and post-AGB
phases are omitted for clarity. All sets of tracks show a similar location for
the ZAMS; the subsequent evolution along the MS and through He-burning phases
differs due to the way each set of models treats mixing across convective
boundaries. The clearest differences between the and other models are in the
location of the hook feature near the MS turnoff for higher masses, a result
of the different adopted core overshoot treatments, and the positions of later
phases, a result of the different wind mass-loss treatments among the
different groups. The COSMIC tracks extend to larger radii and cooler
effective temperatures, which may place them in the regime of LBVs; however,
none of the other sets of evolutionary tracks enter this regime. For a more in
depth comparison see, e.g., Agrawal et al. (2020, 2022).

Figure 4:  Comparison of the final C/O core mass in SSE as implemented by
COSMIC (magenta), MIST (green), and (blue). Differences between MIST and are
due to the larger core overshoot parameter adopted by . Disagreement with the
SSE models is expected as these models are based on simulations that were only
computed for initial masses up to 50 $M_{\odot}$; more massive stars are an
extrapolation.

.

Figure 4 compares the final C/O core mass between , MIST, and SSE as
implemented by COSMIC. SSE models are calculated until central C-burning,
while MIST and are calculated through central C-exhaustion for those stars
with sufficient mass to ignite carbon or to the WD cooling sequence for lower
masses. Differences between the core masses of MIST and POSYDON are generally
due to the different overshooting parameter (we adopt $f_{\rm ov}=0.0415$ for
stars with masses above 8 $M_{\odot}$, compared with $f_{\rm ov}=0.016$
adopted by MIST). The COSMIC/SSE models exhibit different behavior at larger
masses as these prescriptions are based on stellar models that only go up to
50 $M_{\odot}$; larger masses than this are an extrapolation.

Figure 5: Hertzsprung–Russell diagram of a subsample of single-stellar models
where the different stellar states are indicated according to the legend.

As a final comparison, Figure 5 shows a Hertzsprung–Russell diagram of a
subsample of single stellar model where we indicate the different evolutionary
stellar states (Figure 17) across a range of stellar masses.

###  5.4 He-rich, single-star grid

Our second grid of single-star evolutionary models consists of non-rotating
He-rich stars with $Y_{\rm init}=1-Z_{\rm init}$ and our adopted protosolar
$Z_{\rm init}=0.0142$. This grid consists of 40 masses ranging from $M_{\rm
init}=0.5\,M_{\odot}$ to $M_{\rm init}=80\,M_{\odot}$ with a logarithmic
spacing of $\Delta\log_{10}(M_{\rm init}/M_{\odot})=0.055$ dex. For these
masses, stellar evolution models were computed starting from ZAHeMS models
(Section 5.1) and evolved until one of the termination conditions described in
Section 5.2 occurs. For all but the lowest-mass cases, the core C depletion
condition is the relevant one; models with initial masses below
$1.1\,M_{\odot}$ do not ignite C-burning in the core, and therefore terminate
as He-core WDs.

Figure 6: Lifetimes of the and Woosley (2019) single He-star evolution models
match to within $\simeq$0.1 dex.

As a test of our He-star models, we compare their lifetimes to those of the
He-star models of Woosley (2019) in Figure 6. Only the overlapping range of
initial masses is shown here; the Woosley (2019) grid includes models with
masses from 1.8–120 $M_{\odot}$. These models match to within $\simeq 0.1$ dex
in log lifetime across the entire range of initial masses.

Figure 7: Final masses of the and Woosley (2019) single He-star evolution
models. Differences at $M_{\rm initial}\gtrsim 20\,M_{\odot}$ are due to the
stronger wind mass loss prescription adopted by .

As a second test, we compare the final masses between the same two model grids
in Figure 7. Although the lifetimes are similar, the final masses show a
significant difference, particularly at higher initial He-star masses. Woosley
(2019) notes that the change in slope of the initial-final mass relation
around $M_{\rm final}$ of 11 $M_{\odot}$ is due to the mass-loss prescription
adopted for exposed CO cores; at larger initial He-star masses, the entire He-
star mass is burned to heavier elements. For all the single- and binary-star
model grids in , we adopt the mass-loss prescription from Nugis & Lamers
(2000) for He-rich stars. The latter predicts on average stronger wind mass
loss than the prescription from Yoon (2017) adopted by Woosley (2019), leading
to the substantially different final masses between the two prescriptions at
$M_{\rm initial}\gtrsim 20\,M_{\odot}$.

Figure 8: Evolution of radius for He stars with masses between 2 and 10
$M_{\odot}$. Less massive He-stars evolve slower, but expand farther when they
become giant stars. This trend implies that more-massive He-stars in binary
systems will undergo RLO for a relatively narrow range of orbital periods, a
behavior exhibited by our binary star simulations and seen in Figure 14.

Finally, we show the radius evolution of the He star tracks for the mass
interval of 2–10 $M_{\odot}$ in Figure 8. He-rich stars exhibit a peculiar
feature where less-massive stars expand farther on the giant branch than their
more massive counterparts, in agreement with results by Habets (1986). When in
a binary system, this implies that there is a relatively narrow range of
orbital periods in which massive He stars will undergo RLO. Less massive He
stars expand to hundreds of $R_{\odot}$, leading to a wide range of orbital
periods in which these stars can interact with a putative companion. This
behavior is realized in our binary star grids, and its effects are seen
explicitly in Figure 14.

###  5.5 Binaries consisting of two hydrogen-rich main-sequence stars

Figure 9:  View of two grid slices, for two different values of initial binary
mass ratio ($q=0.3$ on the left, $q=0.7$ on the right), from our grid of
binary-star models consisting of two H-rich stars, initially at ZAMS. The
different symbols summarize the evolution of each of the models. We
distinguish between models that experienced stable or no mass transfer
(squares), reaching the end of the life of one of the stars, and the ones that
stopped during mass transfer due to one of our conditions for dynamical
instability (diamonds). Different colors distinguish the evolutionary phase of
the donor star during the latest episode of mass transfer (or no RLO at all
for grey). Small black dots at low initial periods depict systems that were in
initial RLO at birth and red diamonds represent the models that stopped
prematurely for numerical reasons.

Figure 10: For the same grid slices shown in Figure 9, and only for systems
where one of the two stars reached the end of its life, we depict the final
ratio of the angular velocity of the secondary star (the initially less
massive) divided by its critical rotation rate, $(\omega_{\rm s}/\omega_{\rm
s,crit})_{2}$. In most cases where mass transfer occurred, the secondary star
accreted mass and spun up, remaining highly spinning until the end of the life
of the initially more massive star.

For modeling the evolution of two ZAMS stars in a binary system, we run a grid
of 58,240 separate binary evolution models, varying the initial mass of the
primary star $M_{1}$, the initial binary mass ratio $q=M_{2}/M_{1}$ (where
$M_{2}$ is the mass of the companion star), and orbital period $P_{\rm orb}$.
We consider $52$ values of initial primary masses, ranging from
$M_{1}=6.23\,M_{\odot}$ to $M_{1}=120\,M_{\odot}$ with a logarithmic spacing
of $\Delta\log_{10}(M_{1}/M_{\odot})=0.025$ dex, and $20$ values of initial
binary mass ratios, ranging from $q=0.05$ to $q=1$ with a spacing of $\Delta
q=0.05$. Finally, we cover $56$ values of initial orbital period, ranging from
$P_{\rm orb}=0.7$ days to $P_{\rm orb}=6105$ days with a logarithmic spacing
of $\Delta\log_{10}(P_{\rm orb}/{\rm days})=0.07$ dex, in order to explore all
binary configurations ranging from close systems in initial RLO to wide
systems that never exchange any mass.

We simulate binaries by first separately initializing two H-rich, single stars
at ZAMS following the procedure defined in Section 5.1. We then place those
stars in a binary with a second relaxation step, where we force their their
rotation periods to be synchronized with the orbital period, implicitly
assuming that the synchronization has happened during the pre-MS phase. The
latter might not be true for wide binaries, but our assumption induces
negligible rotation to the stellar components of those systems and does not
affect their further evolution. As long as both stars in the binary are under-
filling their Roche lobes after this relaxation step, we start to evolve the
binary. Evolution continues until one of the termination conditions described
in Section 5.2 occurs.

In Figure 9 we provide two two-dimensional slices of this grid, where we show
our simulation outcomes as a function of $M_{1}$ and $P_{\rm orb}$ for fixed
$q$ values. In the left panel we show one example of a mass ratio $q=0.3$, and
on the right a more-equal mass $q=0.7$ slice. Each point in the panels
represents a separate simulation from our grid. Diamond markers represent
models that terminated in a CE, while square markers represent models that
terminated when one of the stars completed its evolution (e.g., reached core C
exhaustion). These are systems that experienced either only stable mass-
transfer episodes or no mass transfer at all, so their evolution can be
continuously modeled. At the bottom of each panel we see systems that are born
filling their Roche lobes (black dots). These systems are assumed to merge,
and therefore never produce a viable binary. Finally, a small fraction of
systems never complete their evolution, producing binary stellar models that
at some point fail to converge (red diamonds).

Separately, the color of each marker indicates that particular binary’s mass
transfer history. Systems with sufficiently close initial $P_{\rm orb}$ tend
to lead to contact phases (orange) where both stars fill their Roche lobes
simultaneously. Most, but not all, of these system end up entering a CE phase.
Sufficiently widely separated (or very massive) systems never fill their Roche
lobes, and therefore never interact (gray markers). For intermediate orbital
periods, the colors differentiate the evolutionary state of the donor when the
latest mass transfer phase was initiated, ranging from MS (blue) to post-MS
(tan), to stripped He-MS (brown). Stable mass transfer causes the donor star
to be almost completely stripped of its H-rich envelope. In the latter case
(brown) the low-mass stripped donors initiate a second mass transfer phase
(Case BB mass-transfer) when they re-expand (Delgado & Thomas, 1981; Laplace
et al., 2020).

Figure 11: View of two slices, for two different values of initial CO masses
($M_{\rm CO}=1.43\,\rm M_{\odot}$ on the left, $M_{\rm CO}=14.66\,\rm
M_{\odot}$ on the right), from our grid of binary-star models consisting of
H-rich star and a CO at the onset of RLO. The different symbols summarize the
evolution of each of the models, as in Figure 9. Binaries that never initiated
mass transfer are not shown here.

Comparison between the two panels shows that the mass ratio leads to a stark
difference in the mass-transfer outcomes. Whereas nearly all systems with an
initial $q=0.7$ result in stable mass transfer, the opposite is true for our
$q=0.3$ systems. At the same time, some features between the two mass ratios
are similar: _(i)_ The boundary between interacting and non-interacting
systems seems to be insensitive to $q$ (and therefore the secondary’s mass).
At the largest orbital periods, stars do not expand far enough to overfill
their Roche lobes. At the largest masses, stars have extremely strong winds
that widen their orbits, simultaneously stripping the primary of its H-rich
envelope, and these stars never expand enough to fill their Roche lobes.
_(ii)_ Systems with initial $P_{\rm orb}\lesssim 5$ days tend to result in
dynamically unstable mass transfer. _(iii)_ There is a large region of
binaries with initial primary mass $\simeq$40–50 $M_{\odot}$ that stably
overfill their Roche lobes as post-MS stars. These stars achieve their mass
transfer stability mainly due to their strong stellar winds, which increases
the mass ratio and the orbit of the system until the moment of overflow.

We model, and keep track of, the properties of both stars in the binary system
throughout their evolution, as well as their detailed internal structure at
the end of the models. In Figure 10 we show, for the same two mass-ratio
slices as in Figure 9, the final rotational rate of the secondary (the
initially less massive) stars for systems that avoid dynamically unstable mass
transfer. Each marker’s color is set by how close each star’s rotation rate is
to its critical rate. Highly rotating secondary stars have all experienced
substantial mass and angular-momentum accretion during their evolution. Many
of them have reached critical rotation, $(\omega_{\rm s}/\omega_{\rm
s,crit})_{2}=1$, early during mass transfer, at which point further mass
accretion becomes non-conservative (c.f. Section 4.2.2). The right-hand panel
shows that the companion’s rotation rate is closely linked with $M_{1}$, as
companion stars with lower mass primary stars also have lower masses and
therefore do not lose as much angular momentum through their own stellar
winds. This behavior is independent of the assumed initial rotation of the
stars.

We find a small subset of initially very close systems in the bottom right
corner (log${}_{10}(M_{1}/M_{\odot})>1.75$ and log${}_{10}(P_{\rm orb}/{\rm
days})\simeq 0.5$) that retain a significant rotational rate even though they
avoid mass transfer. In binaries with such tight orbits, tidal forces between
the stars are sufficiently strong to keep them fast rotating, despite their
strong winds.

###  5.6 Binaries consisting of a compact object and a hydrogen-rich star, at
the onset of Roche-lobe overflow

Figure 12: Relative increase in the mass of the CO $(M_{\rm CO,f}-M_{\rm
CO,i})/M_{\rm CO,i}$ due to accretion for systems where the non-degenerate
star reached the end of its life. The grid slices are the same as shown in
Figure 11. Although accretion is Eddington-limited, COs in binaries with pre-
RLO mass ratios in the range $q\sim 1$–$2$ (defining $q=M_{\rm CO,i}/M_{\rm
1,i}$) and short initial periods, which will experience long-duration mass-
transfer phases, manage to accrete a significant amount of mass.

Figure 13: Same as Figure 12, but now the color of the symbols depict the
maximum mass-transfer rate that occurred in the evolution of each binary. A
significant part of the parameter space leads to highly super-Eddington mass-
transfer rates, albeit in short-lived phases for most cases, and thus to the
potential formation of ultra-luminous X-ray sources. This peak mass-transfer
rate refers to the rate the donor star is losing mass through RLO; accretion
onto the accretor is limited to the Eddington rate.

Our second grid of binary star simulations consists of a H-rich star in a
binary with a CO at the onset of RLO. This grid consists of 25,200 binary
evolution models, where we vary the initial mass of the primary star $M_{1}$,
the initial mass of the CO $M_{\rm CO}$, and the orbital period $P_{\rm orb}$.
We consider $40$ values of initial primary masses, ranging from
$M_{1}=0.5\,M_{\odot}$ to $M_{1}=120\,M_{\odot}$ with a logarithmic spacing of
$\Delta\log_{10}(M_{1}/M_{\odot})=0.06$ dex, and $21$ values of initial CO
masses, ranging from $M_{\rm CO}=1\,M_{\odot}$ to $M_{\rm
CO}=35.88\,M_{\odot}$ with a logarithmic spacing of $\Delta\log_{10}(M_{\rm
CO}/M_{\odot})=0.074$ dex. Finally, we cover $30$ values of initial orbital
period, ranging from $P_{\rm orb}=1.26$ days to $P_{\rm orb}=3162$ days with a
logarithmic spacing of $\Delta\log_{10}(P_{\rm orb}/{\rm days})=0.13$ dex. Our
choice of CO mass range covers massive WD, NS, and BH accretors.

Our procedure in constructing this grid is different from what was described
in Section 5.5. We start each of the simulations with binaries composed of a
ZAMS H-rich star and a CO, which in the MESA code is approximated by a point
mass. Initially, and until each of the binary models reach the onset of RLO,
we neglect orbital angular-momentum loss mechanisms, such as tides, magnetic
breaking and gravitational radiation, while we artificially enforce the
synchronization of the non-degenerate star with the orbit at all times. We do,
however, allow for wind mass-loss from the non-degenerate star, which also
results to a widening of the orbit. Once the onset of RLO is reached, we
include the effects of all orbital angular-momentum loss mechanisms and
discard the prior evolution of the system, treating the onset of RLO as the
effective starting point of our models. Furthermore, from that point onward,
we do not artificially enforce the synchronization of the non-degenerate
star’s spin rotation with the orbit, but we instead follow the tidal
synchronization process self-consistently, following the prescriptions
described in Section 4.1. Finally, binaries that never reach the onset of RLO
are not considered further; these detached binaries are modeled separately as
described in Section 8.1. There, we also provide a full explanation of how we
use this binary-star grid, composed of a H-rich star and a CO at the onset of
RLO, within a larger infrastructure to completely evolve binaries from ZAMS to
double COs.

Figure 11 shows two slices of the grid with different CO masses, $M_{\rm
CO}=1.43$ M⊙ to represent a NS accretor and $M_{\rm CO}=14.66$ M⊙ to represent
a more-massive BH accretor. The symbols depicted in Fig. 11 have the same
meaning as in Figure 9. Although our true initial binary parameters are
regularly spaced, $M_{1}$ and $P_{\rm orb}$ on the axes shown in Figure 11 are
the binary’s quantities at the onset of RLO, the effective starting point of
the models; therefore, the grid does not appear to be regularly spaced (strong
winds exhibited by massive stars tend to expand binary orbits prior to mass
transfer). We do not show those binaries that never interact (even though we
ran these simulations). As already seen in the binary-star model grid composed
of two H-rich stars (Figure 9), binaries too widely separated will never
overfill their Roche lobes, and binaries with massive H-rich stars have winds
too strong to expand into giant phases. In this grid, Figure 11 shows an
additional region of white space at low mass ($M\lesssim 1$ $M_{\odot}$) that
occurs because these stars remain on the MS for the entirety of the
simulation, never expanding to fill their Roche lobes within the age of the
Universe.

Examining the stability of the mass-transfer phase, Figure 11 shows that
nearly every donor star accreting onto a $14.66$ M⊙ BH does so stably, whereas
only the lower mass accretors ($M\lesssim$4.5 $M_{\odot}$) do so for NS
accretors. This difference is because the stability of a mass transfer in a
binary primarily depends on the mass ratio, with a higher accretor mass
allowing for higher donor masses. Our findings, at least for the case of NSs,
are consistent with recent results from Misra et al. (2020), who use the same
criteria to define the onset of L2 overflow leading to dynamical instability
as done in our work.

Figure 14:  View of two grid slices for two different values of initial CO
masses ($M_{\rm CO}=1.43\,\rm M_{\odot}$ on the left, $M_{\rm CO}=14.66\,\rm
M_{\odot}$ on the right), from our grid of binary-star models consisting of
He-rich stars and CO companions. The different symbols summarize the evolution
of each of the models, as in Figure 9.

Figure 12 shows the relative changes in the accretor masses in the same two
slices in $M_{\rm CO}$ as Figure 14. High amounts of accretion mainly depends
on two factors: a sufficiently high-mass accretion rate and a long-lasting RLO
phase. In both panels, this happens for binaries with short periods $\simeq 1$
day, and pre-RLO mass ratios in the range $q\sim 1$–$2$ (defining $q=M_{\rm
CO,i}/M_{\rm 1,i}$). Despite our assumption of Eddington-limited accretion,
for these binaries, stable accretion occurs for over a long time, and in both
cases the binaries transition to low-mass X-ray binaries. These findings are
in agreement with earlier works by Podsiadlowski et al. (2003); Fragos &
McClintock (2015); Misra et al. (2020).

The high mass-transfer rates achieved by most initial binary configurations
are explicitly shown in Figure 13, where each marker’s color corresponds to
the peak mass-transfer rate for each binary. These rates refer to the mass
being lost by the donor star due to RLO; accretion onto the accretor is still
Eddington-limited. In both panels, super-Eddington mass-transfer rates occur
in most binaries with higher peak mass-transfer rates encounters in binaries
with higher periods and larger donor star masses. However, since the larger
orbital separation of these binaries implies the donors in these systems would
be more evolved at RLO onset, compared with initially shorter-period binaries,
these mass transfer phases tend to be short-lived. Therefore, binaries with
short orbital periods (but not so short that they overfill their Roche lobes
initially) will lead to the most accretion onto a CO.

###  5.7 Binaries consisting of a compact object and a He-rich star

Figure 15: The final specific angular momentum $j_{1}=J_{1}/M_{1}$ where
$J_{1}$ is the He-star AM and $M_{1}$ its mass, at carbon depletion for our
grid of He-stars with CO companions. We only show $j_{1}$ for systems where
the non-degenerate star reached the end of its life. The grid slices are the
same as shown in Figure 14.

Our final grid of detailed binary-star simulations consists of 39,480 models
of He-rich stars with CO companions, where we vary the initial mass of the
primary star $M_{1}$, the initial mass of the CO $M_{\rm CO}$, and the orbital
period $P_{\rm orb}$. We consider $40$ values of initial primary masses,
ranging from $M_{1}=0.5\,M_{\odot}$ to $M_{1}=80\,M_{\odot}$ with a
logarithmic spacing of $\Delta\log_{10}(M_{1}/M_{\odot})=0.055$ dex, and 21
values of initial CO masses, ranging from $M_{\rm CO}=1\,M_{\odot}$ to $M_{\rm
CO}=35.88\,M_{\odot}$ with a logarithmic spacing of $\Delta\log_{10}(M_{\rm
CO}/M_{\odot})=0.074$ dex. Finally, we cover $47$ values of initial orbital
period, ranging from $P_{\rm orb}=0.02$ days to $P_{\rm orb}=1117.2$ days with
a logarithmic spacing of $\Delta\log_{10}(P_{\rm orb}/{\rm days})=0.09$ dex.
Our procedure for generating these binaries closely follows the process
described in Section 5.5 for the grid of binary-star modes composed of two
H-rich stars. Here, we replace the initial primary-star models with He-rich
stars at ZAHeMS, while the companion COs are modeled as point masses.

Figure 14 shows an example of two slices of this grid, one corresponding to a
NS companion (with $M_{\mathrm{CO}}\simeq 1.43\,M_{\odot}$) and one
corresponding to a BH (with a $M_{\mathrm{CO}}\simeq 14.66\,M_{\odot}$).
Marker shapes and color scheme follow the same convention as in Figure 11, but
since these simulations are initialized with He-stars, the symbol key is
simplified in Figure 14.

When comparing the two panels, the most apparent difference occurs at large
$M_{1}$ and short orbital period: whereas accreting NSs enter unstable mass
transfer (these systems typically end up merging in a CE, cf. Section 10), the
corresponding accreting BHs typically either overfill their Roche lobes at
ZAMS or avoid mass transfer altogether. In contrast, we find that
independently of the CO mass, systems with low He-star masses
($M_{\mathrm{1}}\leq 3\,M_{\odot}$) mass transfer up to wide orbital periods
($P_{\rm orb}<10^{3}\,\mathrm{days}$). This occurs because low-mass He-stars
expand their He-rich envelope much farther during their later He-shell and
C-burning phases (Figure 8).

Both slices of the grid present two islands of failed simulations, one with
$M_{\mathrm{1}}\simeq 1.8\,M_{\odot}$ and $P_{\rm orb}$ of the order of days
and another island with $M_{\mathrm{1}}\lesssim 1\,M_{\odot}$ and $P_{\rm
orb}$ of the order of hours. MESA has difficulty modeling the envelope’s
structure as it expands to large radii in the first island, whereas the
second, short-$P_{\rm orb}$ island is due to MESA having difficulty following
a star’s evolution into a He WD after it has been spun up due to tides and
mass-transfer. Combined, failed runs account for $\simeq 5$% of the models in
this grid. In practice we find these failed runs do not bias our population
synthesis results of merging NSs and BHs as these portions of the parameter
space predominantly lead to the formation of WDs.

In Figure 15 we show the same two grid slices, but now the marker color
corresponds to the specific angular momentum of the He-star $j_{1}$, at the
end of the simulation. MESA allows us to track this quantity, as it self-
consistently models the interplay between tides (which spin up the star),
stellar winds (which spin down the star and widen the binary), mass transfer
(which alters the orbital period), and internal angular momentum transport.
Comparing Figure 14 and Figure 15, we find that the He-stars with the highest
specific angular momenta are those with either short $P_{\rm orb}$ or stable
mass transfer.

The binary-star grid, composed of a He-rich star and a CO companion, presented
in this section closely agree with those of Qin et al. (2018) and Bavera et
al. (2020, 2021). In contrast to these previous works, the present grid
further expands the parameter space coverage to lower He-star masses and to
larger orbital periods.

##  6 Grid Post-Processing

Each single- or binary-star evolution simulation produces a series of data
files which must be parsed, analyzed, and collated before we can use them
within . Our process includes: (1) re-running any failed simulations; (2)
adding post-processed quantities to our data grids; (3) a post-processing
procedure used exclusively on our single, H-rich and He-rich star grids, which
allows for an efficient interpolation among tracks of different masses; (4)
the downsampling of our grids to reduce data size; (5) classifying each model
within our grids based on the different resulting stellar and binary types,
and (6) fitting classifiers and interpolators over the stellar and binary
parameters in each grid. We describe the first 4 steps next, while the steps
of classification and interpolation are discussed in Section 7.

###  6.1 Re-running Failed Models

Figure 16: Typical example of the evolution of a binary-star model that failed
to reach the end of the simulation due to over-resolved stellar oscillations
that eventually lead MESA to convergence problems. We show the evolution of
the primary’s radius, as a function of its mass (top panel) and its track in
the Hertzsprung–Russell diagram (bottom panel). The binary initially consists
of a 10.50 M⊙ and a 5.25 M⊙ H-rich ZAMS stars at an orbital period of 43.94
days. Comparison between our original, failed simulation (dashed, black line)
and a successful simulation where we artificially limit the radiative opacity
to $0.5\,\mathrm{cm}^{2}/\mathrm{g}$ (orange line) shows that our
approximation is typically accurate to within 0.1 dex for all stellar
parameters. We additionally show the evolution of an adjacent simulation in
our binary grid (blue line; same mass ratio and orbital period but a primary
mass of 9.91 $M_{\odot}$), which shows that any inaccuracy induced by our
opacity approximation is a similar magnitude to the differences between
neighboring simulations.

After having computed our grids of single- and binary-star models, we first
identify those runs that did not reach our desired end point (cf. Section
5.2). This can happen for a variety of reasons, many of which we have not yet
been able to eliminate. For example, one source of problematic runs appears to
deal with stellar oscillations; in certain cases, MESA tries to resolve short-
timescale evolution driven by the $\kappa$-mechanism, which dramatically
shortens the size of successive MESA steps. We address this problem by re-
running our failed binary simulations with a maximum radiative opacity
($\kappa_{\mathrm{max}}$) set to $0.5\,\mathrm{cm}^{2}\,\mathrm{g^{-1}}$. This
approximation reduces the failure rate of each binary grid from $\simeq
10.9\%$, $\simeq 8.0\%$, and $\simeq 11.8\%$, for the binary-star grids
composed of two H-rich stars, a H-rich star with a CO companion at the onset
of RLO, and a He-rich star with a CO companion, respectively, to $\simeq$0.9%,
$\simeq$1.5%, and $\simeq$4.8%. The differences in the resulting evolutionary
tracks with and without the opacity limit are generally small when compared to
differences in tracks of adjacent points in our initial parameter space and
compared to our interpolation accuracy (Section 7.5).

Figure 16 shows a typical example of a binary-star model, initially composed
of two H-rich ZAMS stars with masses 10.50 M⊙ and 5.25 M⊙ and an orbital
period of 43.94 days. This binary initially failed to reach the end of the
simulation (dashed, black line; MESA exceeded its minimum timestep limit), but
did so successfully when re-run with an upper limit to the radiative opacity
(orange line). The top panel shows that the stellar radius evolves similarly
between the two simulations as the donor star loses mass. For the radius and
effective temperature (bottom panel), the two properties most affected by an
opacity limit, differences between the two tracks are typically less than 0.1
dex.

For further comparison we show an adjacent binary-star model in the same grid
with stars of the same mass ratio and orbital period but slightly less massive
primary star of $M_{1}=9.9\,M_{\odot}$, that successfully reached the end of
the simulation without the need to limit the opacity. Although the neighboring
simulation is better able to match our failed run when the luminosity dips to
low values, comparison between all three tracks in Figure 16 suggests that any
inaccuracies accrued by our opacity limit are of a similar magnitude to any
differences between adjacent simulations in our model grids.

Table 3: A list of the different stellar types we adopt in .  Compact | Non-degenerate star states
---|---
`WD` | `H-rich_Core_H_burning` | `H-rich_Core_He_burning` | `H-rich_Shell_H_burning`
`NS` | `H-rich_Central_He_depleted` | `H-rich_Central_C_depletion` | `H-rich_non_burning`
`BH` | `stripped_He_Core_H_burning` | `stripped_He_Core_He_burning` | `stripped_He_Shell_H_burning`
| `stripped_He_Central_He_depleted` | `stripped_He_Central_C_depletion` | `stripped_He_non_burning`
Figure 17: The state of compact objects is determined during the core-collapse
step: if the mass of the pre-SN C/O core is less than $1.37\,M_{\odot}$, then
we form a WD. Otherwise a SN will occur, and the state of the CO is determined
based on its mass. For non-degenerate stars, the state is a combination of a
surface-composition state, and a nuclear burning state. The former is decided
solely by the presence or absence of hydrogen at the surface, whereas the
burning state depends on whether a species (H, He or C) has been depleted from
the core, and the burning of which species still contributes to the nuclear
luminosity ($L_{\mathrm{nuc}}$).

###  6.2 Post-Processed Quantities

Our single- and binary-star MESA simulations result in two types of files:
history files which contain the time evolution of the binary’s and its
component stars’ properties and profile files which define each non-degenerate
star’s structure. In this first version of , we save the profiles of stars
only at the end of the simulations. Combined with the MESA terminal output, we
have all the information necessary to analyze each simulation.

As a first step, we analyze the final binary properties and terminal output to
broadly determine how and why each binary simulation ended. We first identify
the small subset of binaries in which, despite the process described in
Section 6.1, the MESA simulation failed to converge; these are ignored
throughout the remainder of this work. For the successful binary-star
simulations, we have four separate conditions (e.g., Figure 15): (1) binaries
that were already in RLO when initialized at ZA(He)MS; (2) binaries in which
one star reached the end of its lifetime (i.e., one of the first three
termination conditions described in Section 5.2) and went through stable mass
transfer; (3) binaries which avoided mass transfer and the two component stars
essentially evolved in isolation; and (4) binaries that entered unstable mass
transfer (described in Section 8.2).

As a second step, we separately analyze each star and assign it a stellar type
at the end of each simulation. For COs, this is a straightforward task, as the
type of CO is dependent only on its mass. Non-degenerate stars present a more
difficult object to typecast. Most pBPS codes rely upon the $k$-type stellar
classification introduced in Hurley et al. (2000, based on ) In , we use
instead a two-term classification system, dependent upon both what part of the
star (if any) is undergoing nuclear burning and the envelope’s composition.
Table 3 provides a list of the possible stellar type combinations, while in
Figure 17 we show the algorithm for determining them.

As a third step, we assign a designation to the resulting binary configuration
in each of our simulations. These include detached for binaries in which both
stars are confined within their respective Roche lobes, RLO1 or RLO2 for
binaries in which the primary or secondary star, respectively, is overfilling
its Roche lobe; contact for binaries in which both stars are overfilling their
respective Roche lobes; not_converged for systems where the binary-star
simulations ran into numerical convergence problems, and initial_MT for
binaries initialized in RLO. These designations of stars and binary states are
used throughout and are updated by each evolutionary step. Therefore
evolutionary phases modeled with on-the-fly calculations (see Section 8) also
affect the star and binary states, which result in two additional possible
designations: merger for those binaries that have merged, and disrupted for
those binaries that have become unbound due to some process.

As a fourth step, we analyze the mass-transfer history of the binaries we
simulate, identifying the donor star’s state when RLO initiated, and whether
or not that mass transfer phase was stable or unstable (e.g.,
caseA_from_star1). Note that we use the canonical definitions for Case A, Case
B, and Case C mass transfer (for a review, see e.g., Iben, 1991). Specifically
these labels (which can be identified by the differing symbol types and colors
shown in, e.g., Figure 11) identify whether the donor star was on the MS, on
the post-MS, or a stripped He star. In cases where systems evolve through
multiple phases of mass transfer, all phases are included in the label (e.g.,
caseA/B_from_star1 if a Case A mass-transfer phase is followed by Case B one).

As a final step, we calculate a number of post-processed quantities, ranging
from parameters related to specific core-collapse mechanisms (Section 8.3) to
different CE prescriptions (Section 8.2). These are typically parameters that
require integrals over all or part of a star’s structure, which for efficiency
we pre-compute. All post-processed quantities are summarized in Table 8.

###  6.3 Resampling of single-star grids using Equivalent Evolutionary Phases

For our single-star grids, we perform an additional post-processing step,
which resamples the history output of the MESA code in a way that facilitates
the interpolation of entire evolutionary tracks. This is necessary for the
computations described in Section 8.1. Using the method from Dotter (2016), we
assign _equivalent evolutionary phases_ (EEPs) throughout the evolution of a
star. This method designates _primary_ EEPs to major structural changes to a
star (e.g., He ignition), and regularly spaced _secondary_ EEPs in between.
Primary EEPs are extracted directly from the computed stellar tracks. We then
interpolate between the time-steps to identify every quantity of a star that
we track at each secondary EEP. By applying this method to each of our single
star tracks (both H-rich and He-rich) we can more easily interpolate within
our single-star grids to find the quantities (e.g., radius, core mass)
characterizing a star of any mass, and at any point throughout its evolution.

The methodology described above is unfortunately not directly applicable to
binary-star evolutionary tracks. For this interpolation method to work, the
defined EEPs must be strictly ordered a priori. However, binary interactions,
can happen at any point during the lifetime of a binary, often more than once,
changing the order of EEPs in a non-tractable way. The interpolation of entire
binary-star evolutionary tracks will be addressed in future releases of .

###  6.4 Downsampling of binary-star grids

The evolutionary timesteps taken by MESA are typically small, producing high-
resolution binary histories and final profiles of the individual stars. To
reduce the memory footprint of the data, and decrease computation times when
modeling binary populations, we downsample the binary tracks (i.e., keep a
subset of the steps).

For each individual run in a binary-star grid, we obtain from the MESA
simulation the evolution of the binary and individual stars’ parameters, as
well as the post-processed quantities described in Section 6.2. For a total
number of parameters $m$, the state of the binary is encoded in the
$m$-dimensional vector ${\mathbf{h}}$. The evolution of the binary is a
multivariate time-series given by ${\mathbf{h}}_{i}={\mathbf{h}}(t_{i})$, with
$i=1,\cdots,N$, where $N$ is the number of steps, each corresponding to an age
$t_{i}$. Before the downsampling, the independent variable (age), and non-
physical parameters (e.g., model number in MESA) are excluded from
${\mathbf{h}}$, while all other parameters are rescaled linearly from $0$ to
$1$.

Figure 18: The evolution of a $37.6\,M_{\odot}$ star with a $22.6\,M_{\odot}$
companion with an initial $P_{\rm orb}$ of $6.82\,\rm days$. We compare the
complete track provided by MESA (orange) comprised of 3412 steps to our
downsampled track (black dots) containing on 122 steps for the binary’s
orbital period (top panel) and each stars’ radius (bottom two panels). In this
particular case, the compression ratio is ${\sim}155$, but the ratio varies
from star to star and depends on which parameters are accounted for by the
algorithm. The downsampling algorithm captures even the rapid variations seen
in the two stars’ radii between 4.8 Myr and 5.3 Myr (shown in the insets).
Figure 19: The accuracy ($P_{90}(\mathbf{r})$; top panel) and compression
ratio (bottom panel) of our downsampling algorithm as applied to the HMS–HMS
binary grid, as a function of the downsampling threshold. We compare the
performance of the algorithm when it is applied to all output columns in the
data (circle markers) or only a selected list of columns ($\times$; see
excluded parameters annotated with an asterisk in Table 5 and Table 6). As the
downsampling threshold $\epsilon$ increases from $10^{-3}$ to $10^{-1}$, the
compression ratio dramatically improves, but at the cost of the accuracy
($P_{90\%}$ represents the average of the 90-th percentile of the
interpolation relative errors across all runs and either all or selected
parameters). For our grids we use only selected columns with
$\epsilon=10^{-1}$, which gives us a compression ratio of $\simeq$26, while
limiting any errors that could be accrued from this process to within a few
per cent.

The downsampling algorithm selects a subset of the original steps, so that if
interpolated at the original timesteps $t_{i}$ the interpolation absolute
scaled error is below a chosen threshold $\epsilon$:

$e_{i}=\|{\mathbf{h}}_{i}-\hat{{\mathbf{h}}}_{i}\|<\epsilon,$ (14)

where $\hat{{\mathbf{h}}}_{i}$ is the interpolated point using the age as the
independent variable. We use linear interpolation,

$\hat{{\mathbf{h}}}_{i}={\mathbf{h}}_{j}+\frac{t_{i}-t_{j}}{t_{j+1}-t_{j}}\left({\mathbf{h}}_{j+1}-{\mathbf{h}}_{j}\right),$
(15)

where $j$ and $j{+}1$ are adjacent steps of the downsampled time-series so
that $t_{j}<t_{i}<t_{j+1}$.

The search for the steps in the downsampled data is performed as follows.
Initially, we include only the first (${\mathbf{h}}_{1}$) and last
(${\mathbf{h}}_{N}$) points. Then we search for the intermediate point
(${\mathbf{h}}_{j}$ where $1<j<N$) with maximum interpolation error. If this
error is below the threshold ($e_{j}<\epsilon$), then the algorithm has
finished; otherwise it includes this point and continues the search in the two
parts of the time-series before (${\mathbf{h}}_{1}$ to ${\mathbf{h}}_{j}$) and
after (${\mathbf{h}}_{j}$ to ${\mathbf{h}}_{N}$) the intermediate point. The
process continues until all original points are well approximated by the
interpolation of the selected subset of steps, i.e., $e_{i}<\epsilon,\forall
i$.

Additionally, we apply this downsampling method to the stellar profile data,
following the exact procedure outlined above with the mass coordinate as the
independent variable. A shell is kept not only when the interpolation error
exceeds the predefined threshold, but also when the adjacent shells that are
kept have difference in mass is larger then 0.5% of the total stellar mass.

To demonstrate the validity of our method, Figure 18 shows an example of the
downsampling of a track using the same interpolation error threshold
$\epsilon=0.1$ we use for our grids. However, here we apply it to only two
parameters (orbital period and radius of secondary star) for visualization
purposes; when the algorithm operates in a higher-dimensional space, it
retains a large fraction of the initial points to capture the overall shape,
making it hard to inspect its performance through two-dimensional plots. The
downsampled version of the track is able to follow even the rapid oscillations
occurring during the late-stage evolution of this particular binary.

The choice of $\epsilon$ is a balance between data compression and
interpolation accuracy, a trade-off we demonstrate explicitly for our model
grid composed of two H-rich stars in Figure 19. For our grids in , we set the
downsampling threshold to $0.1$ and enforce it only for a list of 22 columns
from the simulation output (see excluded parameters annotated with an asterisk
in Table 5, Table 6, and Table 7). This results to a compression factor of
${\sim}26$ with respect to original simulation data, but still sufficiently
high accuracy with respect to the original grid. The final size of the three
binary-star grids, after downsampling, is $\sim 9.3$ GB.

##  7 Our Classification and Interpolation Approach

Even after our various stages of post-processing, we cannot use the grids of
binary-star simulations within as is for modeling populations (except if we
follow a nearest-neighbor matching approach). While our binary-star
simulations have only been run for a select, finite combination of initial
masses and orbital periods, BPS requires us to have the capability to evolve a
binary anywhere within the domain of interest. We solve this problem in two
steps. First, we apply a classification method to each of our grids to
identify regions that undergo qualitatively different classes of evolution.
Then we separately apply an interpolation method to each class to calculate
stellar and binary properties. We describe the details of those methods in
Section 7.2 and Section 7.3, respectively.

###  7.1 Transformation and Rescaling of Grid Data

For classification and interpolation purposes, we can interpret each of our
binary grids as a data set that comprises $N$ input binaries,
$\\{{\mathbf{x}}_{n}\\}_{n=1}^{N}$, along with its corresponding scalar
$\\{{\mathbf{y}}_{n}\\}_{n=1}^{N}$ and class
$\\{{\mathbf{z}}_{n}\\}_{n=1}^{N}$ targets. These sets can be columnwise
stacked into the matrices ${\mathbf{X}}\in{\mathbb{R}}^{3\times N}$,
${\mathbf{Y}}\in{\mathbb{R}}^{M\times N}$ (where $M$ is the number of output
quantities), and ${\mathbf{Z}}$ of dimension ${4\times N}$ (because we
classify each run into one of four broadly defined categories; see Section
7.2). More specifically, each ${\mathbf{x}}_{n}\in{\mathbb{R}}^{3}$ contains
the initial masses and orbital period of the $n$-th binary in the grid,
whereas ${\mathbf{y}}_{n}$ and ${\mathbf{z}}_{n}$ denote the collection of
final binary and single-star quantities, and their associated classification,
respectively. The $N$ runs are distributed in a uniform three-dimensional
mesh, either on a linear or logarithmic scale following the description
provided in Section 5. This uniform grid constitutes an initial and naive way
of thoroughly covering the parameter space, which is feasible due to the low
dimensionality of the data.

Figure 20: The optimal $k$ in our $k$NN classification scheme for each of our
three grids. The bACC is calculated using 10-fold MC cross-validation, and
averages the statistical recall of each of our four classes described in
Section 7.2. The highest bACC for each grid is provided at the top of each
panel and is indicated by the vertical lines; we use these $k$s when
classifying our grids for running populations. Although we use the optimal $k$
for each grid, our results are relatively insensitive to that choice.  Figure
21: Decision boundaries of the $k$NN classifiers for a single slice in each of
the three grids as a function of the primary star’s mass on the horizontal
axis and the orbital period on the vertical axis (the choice of $q$ or $M_{\rm
CO}$ for each slice is indicated in each panel’s title). Shaded gray regions
overlaid onto class regions represent the confidence of the classifier in that
region. Points on top of the decision boundaries represent the validation
data, where the edge color of each point shows the ground truth of the given
point, and the fill color shows the classifier’s prediction. Only in rare
circumstances and only near classification boundaries does our classifier make
incorrect predictions for our validation set.

A convenient preprocessing of the data is crucial for both interpolation and
classification. We apply a series of non-linear and linear transformations to
numeric data. Choosing the optimal transformation depends on the task
(interpolation or classification), the method used for each task, and whether
we are dealing with an input or an output quantity.

First, we consider a non-linear transformation of the data using the
logarithm: inputs can be transformed as $\log{\mathbf{x}}_{i}$ and targets as
$\log{\mathbf{y}}_{i}$ or $\log(-{\mathbf{y}}_{i})$ if ${\mathbf{y}}_{i}<0$.
Classification accuracy will improve when our algorithm uses the logarithm of
the inputs for data sampled evenly in log-space. The effect on interpolation
is different: e.g., a linear interpolation in the log-space results in a non-
linear interpolation on the untransformed space. This is similar to an
approach where a non-linear space is transformed through a kernel to a space
in which a linear model allows for modeling behavior appropriately.

We automatically choose whether to apply a logarithmic scaling using a cross-
validation scheme. The optimal scaling for both inputs and outputs is chosen
using the lowest relative error out of all the feasible scalings that could be
applied to the given variable.

As a second step, we apply a min–max scaling to the inputs so that the
transformed features ${\mathbf{x}}_{i}^{\mathrm{t}}$ are confined to the range
$[-1,1]$, and we standardize the outputs such that they have zero-mean and
unit variance,

${\mathbf{y}}_{i}^{\mathrm{t}}=\frac{{\mathbf{y}}_{i}-\overline{{\mathbf{y}}}_{i}}{\sigma_{{\mathbf{y}}_{i}}}.$
(16)

The choice of scaling for the inputs derives from the uniform nature of the
input grid data. Although it is possible that we sampled our data in a non-
optimal way, in practice we find the best results occur when our data scaling
follows our grid sampling. In the case of the interpolated quantities,
standarization produces improved metrics, particularly because it is less
sensitive to outliers.

###  7.2 Classification of our Grids

Accurate classification is a critical aspect of the approach to evolving
binary systems. Therefore, we separate our binaries into four categories based
on their mass-transfer histories. The categories are: stable mass transfer,
unstable mass transfer, binaries that never interact, and those in RLO at ZAMS
(Section 6.2). In addition to using their mass-transfer history, we could
further segregate binaries into more refined classes; however, we find this to
be currently unnecessary, and we can accurately interpolate our binaries given
these four broad classes.

Figure 22: Confusion matrices for each of our three binary grids. Each value
at grid cell $c_{ij}$ represents the fraction of binaries that belong to class
$j$ (vertical axis) and were classified as class $i$ (horizontal axis). Each
row is normalized so that the sum of each row is $1$, and the color of each
cell indicates the magnitude of the value in the cell. Accuracies are all
above 90%, with the exception of unstable mass transfer for the He-rich star
with a CO companion.

For each of our three binary-star grids, we generate a classification object
that determines which of the four previously defined outcomes will be the
result of a binary with any particular combination of two masses and orbital
period. In this first version of we use a $k$-nearest neighbors ($k$NN)
classifier, a simple and robust classifier that achieves high precision in
this task. We use the Euclidean distance as distance metric for the
transformed input grid and weight each neighbor in the neighborhood
proportionally to their inverse distance.

We optimize the number of nearest neighbors we use by applying a Monte Carlo
(MC) cross-validation scheme and selecting the $k$ that produces a higher
balanced accuracy (bACC). The bACC metric averages the statistical recall for
each class (recall is the number of true positives divided by the combined
number of true positives and false negatives) to produce a metric that
accounts for any imbalances between classes.

Figure 20 shows the average cross-validation performance for our three grids
in terms of bACC as a function of the number of neighbors in the $k$NN
classifier starting from $k=1$ and highlights the location of the optimum. We
train our final classifier on our regularly spaced grids using the optimal
$k$, listed at the top of each panel in Figure 20 and indicated by a vertical
black line.

Figure 21 shows our classifier applied to one slice in each of our three
grids, with the colors indicating different regions. Overlaid grey contours,
pronounced near class boundaries, indicate classification uncertainty. We
ignore the no mass-transfer class for the grid of H-rich stars with CO
companions, as this grid only applies to interacting binaries.

To evaluate the accuracy of our classifiers, we use the validation data set
associated with each of our three binary grids. For each grid this validation
set is comprised of binaries randomly sampled with the same range and scale,
linear or logarithmic, as its training counterpart. Each of the three
validation sets contain $3000$ samples which roughly represent $\sim$5–12% of
the number of binaries in the regular grid. By applying our classifiers to the
same initial values as those of our validation binaries, we can evaluate the
accuracy of our classifiers. In Figure 21 our validation data is indicated by
points (correctly classified) and crosses (incorrectly classified). It is
evident that incorrectly classified validation binaries are very rare.

To evaluate the quantitative accuracy of our classifier, we provide a
confusion matrix for each of three grids in Figure 22. Diagonal squares
indicate the fraction of systems that were correctly classified, while off-
diagonal squares indicate the fraction of incorrectly classified systems. The
matrices are calculated such that each row sums to unity. All classes in all
grids have an accuracy in excess of 90%, often much more so, except for
unstable mass transfer for binaries with a He-rich star and a CO companion.
Examination of two slices of this grid in Figure 15 shows that the unstable
mass-transfer class comprises a relatively small portion of the overall grid,
existing at small orbital periods, small CO masses and large companion masses.
Reliable classification of small classes is difficult, but improving our
classification accuracy will be a focus of future efforts (Section 11).

###  7.3 Interpolation of our Grids

Once classified based on their mass-transfer characteristics we separately
interpolate binaries falling into each class for each of our three binary-star
simulation grids. We only interpolate quantities for three of our binary
classes, since those binaries overfilling their Roche lobe at ZA(He)MS are
dismissed.

Constraint |  Relation
---|---
Type I constraints: Equations
Kepler’s Third Law |  $a=\left[G\left(M_{1}+M_{2}\right)P_{\rm orb}^{2}/4\pi^{2}\right]^{{1}/{3}}$
Mass-transfer fraction(a) |  $x=1-\dot{M}_{\rm{sys,2}}/\dot{M}_{\rm{tr}}$
Stefan-Boltzmann Law |  $T_{\rm{eff}}=\left(L/4\pi R^{2}\sigma_{\rm{SB}}\right)^{{1}/{4}}$
Sum of nuclear luminosities |  $L_{\rm nuc}=L_{\rm H}+L_{\rm He}+L_{\rm Z}$
Type II constraints: Inequalities
Mass-loss from the system from the vicinity of a star(b) |  $\dot{M}_{\rm{sys,1}}<\dot{M}_{\rm{tr}}$ and $\dot{M}_{\rm{sys,2}}<\dot{M}_{\rm{tr}}$
Core masses and radii |  $M_{\rm{C/O-core}}<M_{\rm{He-core}}<M$ and $R_{\rm{C/O-core}}<R_{\rm{He-core}}<R$
Envelope masses and core radii(c) |  $M_{\rm{env}}{<}M$ and $R_{\rm{core}}{<}R$
Mass, thickness and middle radius of the convective region (for tides) |  $M_{\rm{conv.reg.}}<M$ and $0<R_{\rm{conv.reg.}}-(1/2)D_{\rm{conv.reg.}}<R_{\rm{conv.reg.}}+(1/2)D_{\rm{conv.reg.}}<R$
Remnant baryonic mass |  $M_{\rm{rembar}}<M$
Type III constraints: constrained sum
Central abundances |  $X_{\rm{c,H1}}+X_{\rm{c,He4}}+X_{\rm{c,C12}}+X_{\rm{c,N14}}+X_{\rm{c,O16}}+X_{\rm{c,other}}=1$
Surface abundances |  $X_{\rm{s,H1}}+X_{\rm{s,He4}}+X_{\rm{s,C12}}+X_{\rm{s,N14}}+X_{\rm{s,O16}}+X_{\rm{s,other}}=1$
Table 4: Constraints that are ensured for our interpolated quantities. Type I
relations are written such as the left-hand side indicates the quantity that
is inferred from the rest. Notes:
(a) This constraint is applied after constraint (b) since the mass-loss rate
has to be less than the mass-transfer rate. Moreover, $x$ is set to $1$ if no
mass is transferred.
(c) There are four pairs of these quantities corresponding to the radii where
the 1H fraction drops below 1%, 10% and 30%, and finally where the 4He
fraction drops below 10% for pure He stars.

We use an $N$-dimensional (where $N$ is the number of binaries) linear
interpolation: the data is divided into a set of $N$-simplices, tetrahedra in
our three-dimensional data, by means of a Delaunay triangulation (which is not
unique given the regular structure of our grids). The interpolated value for a
given point corresponds to the value at the hyperplane that passes through the
vertices of the simplex which contains the point. The choice of whether to
apply a non-linear transformation on ${\mathbf{y}}_{i}$,
$\log{\mathbf{y}}_{i}$ depends directly on that magnitude. For each output
magnitude we select the optimal scaling via MC cross-validation with $x$
iterations and $p\%$ of test data comparing the average relative error across
iterations. The final interpolator is trained using all binaries within a
particular class for each grid.

The linear interpolation method is not capable of extrapolation: the value for
any point which lies outside the convex hull defined by the constructed
Delaunay triangulation will be undetermined. Although we are not, in general,
interested in interpolating outside the training grid, there will be a small
region between the convex hull of the linear interpolation and the decision
boundary provided by the classifier where we still want to obtain system
properties. For this small sliver of parameter space, we adopt values of the
nearest point in parameter space of the same class. This is a problematic
region where the probability of belonging to the interpolable class will be
low, expressing the uncertainty we have about those binaries with the current
resolution of the grids. We are currently exploring a method of tackling this
problem by incorporating new simulations along the decision surface,
identified using an active-learning scheme (Rocha et al., 2022).

###  7.4 Ensuring Physical Congruity of Interpolated Values

The linear interpolation method described here treats each feature
independently without preserving possible physical correlations. However, the
interpolated results may produce incongruous quantities within a resultant
star. For example, the Stefan–Boltzmann law connecting the luminosity, radius,
and effective temperature of a star might not hold for an interpolated binary.
As another example, the He-core mass of a star must always be less than the
star’s total mass. We have carefully identified a number of physical
constraints within the quantities that we are interpolating that must be
satisfied by any realistic star, each of which we list in Table 4.

Figure 23:  The convective envelope radius $R_{\rm{conv.reg.}}$ and width
$D_{\rm{conv.reg.}}$ need to simultaneously agree with three separate
inequalities defined in Table 4. The feasibility region (dark grey) represents
the overlap of all three inequalities decomposed from the constraint. When our
interpolation method proposed a point outside the feasibility region, we
reassign it to a new point determined from the intersection of the border of
the feasibility region and the line drawn from the centroid of the region to
the constraint-violating, proposed point.

To address this issue, we process the interpolated quantities for a given
binary so that they respect this list of constraints, which fall into one of
three different _Types_ depending on the basis of the corrective action
required. When quantities are connected via an equation (Type I), then the
interpolated value of one parameter is ignored and inferred by solving the
equation on the interpolated values of the remaining parameters. The
Stefan–Boltzmann equation provides an example of a Type I constraint: we only
interpolate each stars’ $R$ and $L$, while $T_{\rm eff}$ of the interpolated