Most massive stars are members of a binary or a higherorder stellar system 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 compactobject binaries and to decipher the physical processes that lead to their formation. Here we present POSYDON a novel publicly available binary population synthesis code that incorporates full stellar structure and binaryevolution modeling using the MESA code throughout the whole evolution of the binaries. The use of POSYDON enables the selfconsistent treatment of physical processes in stellar and binary evolution including realistic masstransfer calculations and assessment of stability internal angularmomentum transport and tides stellar core sizes masstransfer rates and orbital periods. This paper describes the detailed methodology and implementation of POSYDON including the assumed physics of stellar and binary evolution the extensive grids of detailed single and binarystar models the postprocessing classification and interpolation methods we developed for use with the grids and the treatment of evolutionary phases that are not based on precalculated grids. The first version of POSYDON targets binaries with massive primary stars potential progenitors of neutron stars or black holes at solar metallicity.
['Binary stars', 'Close binary stars', 'Compact binary stars', 'Interacting binary stars', 'X-ray binary stars', 'Compact objects', 'Stellar remnants', 'Black holes', 'Neutron stars', 'Gravitational wave sources', 'Stellar evolutionary models', 'Stellar populations', '154', '254', '283', '801', '1811', '288', '1627', '162', '1108', '677', '2046', '1622', 'Astrophysics - Solar and Stellar Astrophysics']
POSYDON A Generalpurpose Population Synthesis Code with Detailed Binaryevolution Simulations
{'POSYDON : A General-Purpose Population Synthesis Code with Detailed Binary-Evolution Simulations': "Tassos Fragos, 1 Jeff J. Andrews, 2, 3 SimoneS.Bavera, 1 ChristopherP.L.Berry, 2, 4 ScottCoughlin, 2 AaronDotter, 2 PrabinGiri, 5 VickyKalogera, 2, 6 AggelosKatsaggelos, 2, 7 KonstantinosKovlakas, 1 ShamalLalvani, 7 DevinaMisra, 1 Philipp M. Srivastava, 2, 7 YingQin, 8 KyleA.Rocha, 2, 6 Jaime Rom'an-Garza, 9 JuanGabrielSerra, 2, 7 PetterStahle, 10 MengSun, 2 XuTeng, 5 GoceTrajcevski, 5 NamHaiTran, 11 Zepei Xing, 1 EmmanouilZapartas, 1 and Michael Zevin 12, 13  \n1 D'epartement d'Astronomie, Universit'e de Gen'eve, Chemin Pegasi 51, CH-1290 Versoix, Switzerland 2 Evanston, IL 60208, USA  \nCenter for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road, 3 Department of Physics, University of Florida, 2001 Museum Rd, Gainesville, FL 32611, USA 4 Institute for Gravitational Research, University of Glasgow, Kelvin Building, University Avenue, Glasgow, G12 8QQ, Scotland 5 Department of Electrical and Computer Engineering, Iowa State University, 2520 Osborn Dr, Ames, IA 50011, USA 6 Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA 7 Electrical and Computer Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA 8 Department of Physics, Anhui Normal University, Wuhu, Anhui 241000, China 9 Universidad de Monterrey, Ave. Morones Prieto 4500 Pte., C.P. 66283, San Pedro Garza Garc'ıa, Nuevo Le'on, M'exico. 10 D'epartement d'Informatique, Universit'e de Gen'eve, Route de Drize 7, CH-1227 Carouge, Switzerland 11 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark 12 Kavli Institute for Cosmological Physics, The University of Chicago, 5640 South Ellis Avenue, Chicago, Illinois 60637, USA 13 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 POSYDON , 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 POSYDON 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 POSYDON , 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 POSYDON targets binaries with massive primary stars (potential progenitors of neutron stars or black holes) at solar metallicity.', '1. INTRODUCTION': "Throughout their lives, stars affect their surroundings via the immense energy radiated across the electromag-  \nCorresponding author: T. Fragos  \nAnastasios.Fragkos@unige.ch  \nnetic 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 ( glyph[greaterorsimilar] 8 M glyph[circledot] ) 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  \nsupernova (SN) or a γ -ray burst that release more energy in 10 s than our Sun in 10 10 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 starformation (e.g., Nomoto et al. 2013; Hopkins et al. 2011, 2012).  \nIt 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).  \nOver 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 LEGUSsurveys; 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., PanSTARRS, 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.  \nAspects of the astrophysics of all these different types of stellar binaries can be obtained from observations and modeling of present-day properties of individual, wellstudied 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 distri-  \nbutions 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.  \nOver 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.  \nBPS codes stand in stark contrast to detailed stellarstructure 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 ∼ 10100 CPU hours for the modeling of a single system (e.g., Paxton et al. 2019).  \nAs 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 precalculated 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.  \nStudies 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,a; 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., Osglyph[suppress]lowski et al. 2011; Andrews et al. 2015; Chruslinska et al. 2017; Vigna-G'omez 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).  \nHowever, 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-tocalibrate 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 masstransfer rate, especially during thermal-timescale masstransfer 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 selfconsistent modeling of the end of a mass transfer phase (e.g., accounting for a potential partial stripping of the envelope).  \nIn this work, we build upon the combined experience gained from the large body of BPS studies to date, to create POSYDON (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. 1 With POSYDON 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), POSYDON 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 POSYDON 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-  \ncalculated 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 corecollapse 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 POSYDON code. The definitions of all the symbols used throughout this paper can be found in Table 1.  \nTable 1 . List of variables used throughout the paper.  \nTable 1 continued  \nTable 1 (continued)  \nTable 1 continued  \nTable 1 (continued)  \nTable 1 (continued)", '2. OVERVIEW OF THE STRUCTURE OF POSYDON': "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 POSYDON 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 ∼ 10 5 -10 6 binaries, far too few to accurately model a Milky Way sized population with ∼ 10 11 stellar binaries. As a further complication, codes like MESA cannot run individual binaries from start to  \nfinish; key physics, including CE phases and supernovae require the code to be stopped and restarted. 2  \nPOSYDON solves the problems associated with stellarstructure codes by using extensive, pre-calculated grids of single and binary stellar-evolution models, covering the parameter space relevant for the formation of highmass 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, POSYDON 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 postprocessed, 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 nearestneighbor matching, and in that case no classification or interpolation methods are required.  \nThe second major component of POSYDON 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-thefly calculations, where we use the single-star grids coupled to binary evolution routines, i.e. orbital evolu-  \ntion due to tides, stellar winds, magnetic breaking and gravitational-wave emission (Section 8.1).  \nIn the POSYDON 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 POSYDON 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 POSYDON 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 POSYDON , 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). 3 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  \nFigure 1. The structure of POSYDON (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 POSYDON 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.  \n<!-- image -->  \ntaken 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: 1 H, 3 He, 4 He, 12 C, 14 N, 16 O, 20 Ne, 24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 56 Cr, 52 Fe, 54 Fe, 56 Fe and 56 Ni, plus protons and neutrons (for the purpose of photo-disintegration).', '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 eff ( τ ) relation, and the photosphere pressure via P = τg/κ 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  \nwind 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 POSYDON , changing the wind prescription would require the computation of new set of single- and binary-star model grids.  \nFor stars with initial masses above 8 M glyph[circledot] we use the MESA Dutch scheme, which consists of de Jager et al. (1988) for T eff < 10 , 000 K and Vink et al. (2001) for T eff > 11 , 000 K. In cases where T eff > 11 , 000 K and the surface 1 H 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 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).  \nFor stars with initial masses below 8 M glyph[circledot] , we again use the Dutch scheme for stars with T eff hotter than 12,000 K. For stars with T eff less than 8,000 K, we use the Reimers (1975) wind with scaling factor η R = 0 . 1 for stars on the first ascent of the giant branch, and the  \nBloecker (1995) wind with scaling factor η B = 0 . 2 for stars in the thermally-pulsating phase. For the case of 8,000 K < T eff < 12,000 K, we calculate the wind for both the hot and cool schemes, and linearly interpolate between the two.  \nFor 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.  \nWefurther boost stellar winds to limit a star's rotation below its critical threshold ( ω s /ω s , crit ≤ 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),  \n˙ M w ( ω ) = ˙ M w (0) ( 1 1 -ω s /ω s , crit ) ξ , (1)  \nwhere ˙ M w is the star's wind mass-loss rate, and ω s and ω s , crit are the angular velocity and critical angular velocity at the surface, respectively. The default value of the exponent ξ = 0 . 43 is taken from Langer (1998). The critical angular velocity is given the expression ω 2 s , crit = (1 -L/L Edd ) GM/R 3 , where L 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; Bohm-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  \nparameter, α MLT = 1 . 93, based on results from the MIST project (Dotter et al., in preparation).  \nRotation is implemented in MESA as described in Paxton et al. (2013, 2019). Rotational mixing and angularmomentum 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).  \nMESA 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 α th = 17 . 5 (Paxton et al. 2013, Eq. 14), also referred to as C t (Charbonnel & Zahn 2007, Eq. 4), corresponding to an aspect ratio ∼ 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).  \nOvershoot mixing is treated in the exponential decay formalism (Herwig 2000; Paxton et al. 2011). For the parameter f ov describing the extent of the overshoot mixing in this formalism, we adopt an initialmass-dependent relation. For lower mass stars (initial masses less than 4 M glyph[circledot] ) we adopt a value taken from the MIST project, f 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 glyph[circledot] ), we adopt a value of f ov = 0 . 0415 motivated by the work of Brott et al. (2011), who used the step overshoot formalism. Both of these values of f ov are measured from a distance of 0 . 008 the local pressure scale height into the convection zone from the formal convectiveradiative 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 ∼ 10 larger than f ov (their Figure 3). For stars with initial masses between 4 M glyph[circledot]  \nand 8 M glyph[circledot] we smoothly ramp between the two values of f ov . The mass range was chosen to be roughly consistent with the ranges considered in the two studies.  \nWe 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': "POSYDON is predominantly a BPS tool that simulates the evolution of an ensemble of binary systems through various stages of their life. In the POSYDON 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.  \nWe 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  \nbinary corresponding to the difference between the orbital and the spin angular velocity ∆Ω, divided by the synchronization timescale τ sync :  \nδ Ω = ∆Ω τ sync δτ, (2)  \nwhere δτ is the timestep and δ Ω 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.  \nWe separately calculate τ sync for both stars (Hut 1981; Hurley et al. 2002; Paxton et al. 2015):  \n1 τ sync = 3 ( k T ) q 2 MR 2 I ( R a ) 6 . (3)  \nHere M , R and I are the mass, radius and moment of inertia of the star for which we calculate the tides, q = M 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, τ 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.  \nFor radiative regions in a star, we calculate k/T based on the dynamical tidal timescale from Zahn (1977, see also Hut 1981), where  \n( k T ) rad = √ GMR 2 a 5 (1 + q ) 5 / 6 E 2 , (4)  \nwith E 2 as the second order tidal coefficient and G the gravitational constant. 4 For the calculation of E 2 , we  \nadopt the latest prescriptions from Qin et al. (2018), who investigated the dependence of the parameter on the convective radius R conv for various metallicities and evolutionary stages, finding  \nE 2 =    10 -0 . 42 ( R conv /R ) 7 . 5 for hydrogen-rich stars 10 -0 . 93 ( R conv /R ) 6 . 7 for stripped-helium stars. (5)  \nFor 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):  \n( k T ) conv = 2 21 f conv τ conv M conv . reg . M , (6)  \nand use the shortest timescale among them. In Eq. (6), M conv . reg . is the mass of the convective region, f 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 τ conv is the convective timescale, which we take from Eq. (31) of Hurley et al. (2002, based on Rasio et al. 1996), adapted to also accomodate convective regions that are below the surface:  \nτ conv =0 . 431 [ M conv . reg . 3 L R t , conv . reg . + R b , conv . reg . 2 × ( R t , conv . reg . -R b , conv . reg . ) ] 1 / 3 . (7)  \nIn the equation above, R t , conv . reg . , R 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.  \nIn Figure 2, we show two example evolutionary tracks in a Hertzsprung-Russell Diagram of the primary star in an interacting binary of M 1 , initial = 56 . 46 M glyph[circledot] and M 2 , initial = 28 . 23 M glyph[circledot] , 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  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 2. Evolution of the donor star in two example close binary systems of initially M 1 , initial = 56 . 46 M glyph[circledot] and M 2 , initial = 28 . 23 M glyph[circledot] , 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.  \n<!-- image -->  \nimportant 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  \n<!-- image -->  \n<!-- image -->  \nappear 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 glyph[circledot] ), 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).  \nAs 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). 5 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.  \n4.2.2. Mass Accretion onto a Non-Degenerate Companion  \nThe 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.  \nFor 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 rotationallyenhanced 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).  \nWe calculate the Eddington-limited rate using standard formulae (Frank et al. 2002). We first calculate the Eddington luminosity L Edd for an accretor with mass M acc :  \nL Edd = 4 πGM acc c κ , (8)  \nwhere κ 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 κ = 0 . 2(1 + X ) cm 2 g -1 , where X is the hydrogen abundance of the donor. By setting L Edd equal to the radiation released by accreted matter as it falls into a CO's potential well ( L acc = η ˙ Mc 2 ), we can recover the Eddington-limited accretion rate,  \n˙ M Edd = 4 πGM acc κcη . (9)  \nThe dimensionless constant η sets how efficiently the rest mass energy of the incoming matter is converted to outgoing radiation,  \nη = GM acc R acc c 2 . (10)  \nFor BHs, R acc is set by the spin-dependent innermost stable circular orbit, while for NSs, we use a constant R 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 glyph[circledot] being classified as NSs, while those with mass greater than 2 . 5 M glyph[circledot] as BHs. Within our code, the only difference between these types of accretors is the corresponding η ; otherwise accretion proceeds identically regardless of the type of CO accretor.  \nAs 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 η through the radius of innermost stable circular orbit (ISCO; Podsiadlowski et al. 2003),  \nη = 1 -√ 1 -( M acc 3 M i acc ) 2 , (11)  \nwhere M i acc is the initial mass of the accreting BH and M acc its current one. In this equation it is also implicitly assumed that the birth spin of the BH is ∼ 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 spin , can be calculated following Thorne (1974) and King & Kolb (1999),  \na spin = ( 2 3 ) 1 / 2 M i acc M acc    4 -[ 18 ( M i acc M acc ) 2 -2 ] 1 / 2    . (12)  \nThese equations both assume that the spin-up occurs due to angular momentum accretion from a disk that is truncated at the ISCO.  \nWe do not explicitly stop our simulations if a NS accretes enough mass to cross our 2.5 M glyph[circledot] threshold thereby collapsing into a BH, but we do switch the η instantaneously.  \n4.2.4. Onset of CE  \nStars 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'nski & 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.  \nFirst, we assume that a dynamically unstable RLO phase is initiated whenever mass-transfer rate exceeds 0 . 1 M glyph[circledot] 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 glyph[circledot] 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.  \nAs 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 (L 2 ). In such cases the lost matter from the L 2 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 sphericalequivalent radius corresponding to L 2 . 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,  \ncf. Eq. 2) for the L 2 radius, which considers that both stars can contribute to the overflow of the L 2 volume together.  \nFor 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 trap depends on the Eddington limit of the accretor ˙ M Edd and the mass-transfer rate from the donor ˙ M donor , we limit the latter assuming an instability condition when (Begelman 1979):  \n˙ M donor ≥ ˙ M Edd 2 R L , accr R acc . (13)  \nAs 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.  \nAs 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 fullyconvective star with no nuclear fusion taking place and  \nadopting 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.  \nThe He-rich ZAMS (ZAHeMS) models are created in three steps. First, we create a pre-MS He star with 100% 4 He 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.  \nThe two sets of ZAMS and ZAHeMS models are used as a starting point for the five grids of single- and binarystar 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:  \n- · 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 init glyph[lessorsimilar] 0 . 8 M glyph[circledot] ). 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.\n- · A star becomes a WD, a condition we quantify by checking if the central degeneracy parameter Γ c (Coulomb coupling parameter) exceeds 10 (Choi et al. 2016).\n- · 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.\n- · A binary enters a CE phase, as described in Section 8.2.\n- · 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.  \nOur simulations may occasionally end prematurely before any of the aforementioned conditions are reached. This may happen because the minimum timestep limit", 'Fragos et al.': 'Table 7 (continued)  \nTable 8 . Post-processed variables referring to the final state of each binary in our grids. Quantities referring to singlestar quantities (i.e., all except termination flags and interpolation class) are prefixed with S1 and S2 to distinguish the corresponding star in the grid (e.g., S1 surface other or S2 direct mass ). These variables, along with the last values of the single and binary history variables (cf. Table 5 and Table 6; e.g., S1 log L ), comprise the final values tables stored in the grids.Table 8 continued  \nTable 8 (continued)', '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 init = 0 . 5 M glyph[circledot] to M init = 300 M glyph[circledot] with a logarithmic spacing of ∆log 10 ( M init /M glyph[circledot] ) = 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.  \nTo test their validity, we compare POSYDON evolutionary tracks to the widely used stellar evolution tracks from the Geneva (upper panel; Ekstrom 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 glyph[circledot] and 300 M glyph[circledot] . 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  \nthrough He-burning phases differs due to the way each set of models treats mixing across convective boundaries. The clearest differences between the POSYDON 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 massloss 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).  \nFigure 4 compares the final C/O core mass between POSYDON , MIST , and SSE as implemented by COSMIC . SSE models are calculated until central C-burning, while MIST and POSYDON are calculated through central Cexhaustion 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 ov = 0 . 0415 for stars with masses above 8 M glyph[circledot] , compared with f 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 glyph[circledot] ; larger masses than this are an extrapolation.  \nAs a final comparison, Figure 5 shows a HertzsprungRussell diagram of a subsample of POSYDON single stellar model where we indicate the different evolutionary POSYDON stellar states (Figure 17) across a range of stellar masses.  \n5.4. He-rich, single-star grid  \nFigure 3. Comparing POSYDON H-rich ZAMS evolutionary tracks (blue in all panels) with non-rotating Geneva (upper; Ekstrom 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 glyph[circledot] in all cases. The masses shown are the same as the Geneva grid of models between 1 and 120 M glyph[circledot] with the addition of 175 and 300 M glyph[circledot] for POSYDON , MIST , and COSMIC/SSE . The dotted, gray lines indicate constant radius at powers of 10 in R glyph[circledot] .  \n<!-- image -->  \nOur second grid of single-star evolutionary models consists of non-rotating He-rich stars with Y init = 1 -Z init and our adopted protosolar Z init = 0 . 0142. This grid consists of 40 masses ranging from M init =  \nFigure 4. Comparison of the final C/O core mass in SSE as implemented by COSMIC (magenta), MIST (green), and POSYDON (blue). Differences between MIST and POSYDON are due to the larger core overshoot parameter adopted by POSYDON . 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 glyph[circledot] ; more massive stars are an extrapolation.  \n<!-- image -->  \n/circledot  \n.  \n0 . 5 M glyph[circledot] to M init = 80 M glyph[circledot] with a logarithmic spacing of ∆log 10 ( M init /M glyph[circledot] ) = 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 glyph[circledot] do not ignite C-burning in the core, and therefore terminate as He-core WDs.  \nAs a test of our POSYDON 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 glyph[circledot] . These models match to within glyph[similarequal] 0 . 1 dex in log lifetime across the entire range of initial masses.  \nAs 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 final of 11 M glyph[circledot] 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 singleand binary-star model grids in POSYDON , we adopt the mass-loss prescription from Nugis & Lamers (2000) for He-rich stars. The latter predicts on average stronger  \nFigure 5. Hertzsprung-Russell diagram of a subsample of POSYDON single-stellar models where the different POSYDON stellar states are indicated according to the legend.  \n<!-- image -->  \n<!-- image -->  \n/circledot  \nFigure 6. Lifetimes of the POSYDON and Woosley (2019) single He-star evolution models match to within glyph[similarequal] 0.1 dex.  \n<!-- image -->  \n/circledot  \nFigure 7. Final masses of the POSYDON and Woosley (2019) single He-star evolution models. Differences at M initial glyph[greaterorsimilar] 20 M glyph[circledot] are due to the stronger wind mass loss prescription adopted by POSYDON .Figure 8. Evolution of radius for He stars with masses between 2 and 10 M glyph[circledot] . 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.  \n<!-- image -->  \nwind 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 initial glyph[greaterorsimilar] 20 M glyph[circledot] .  \nFinally, we show the radius evolution of the He star tracks for the mass interval of 2-10 M glyph[circledot] in Figure 8. Herich 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 Ha-  \nbets (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 glyph[circledot] , 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': "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 orb . We consider 52 values of initial primary masses, ranging from M 1 = 6 . 23 M glyph[circledot] to M 1 = 120 M glyph[circledot] with a logarithmic spacing of ∆ log 10 ( M 1 /M glyph[circledot] ) = 0 . 025 dex, and 20 values of initial binary mass ratios, ranging from q = 0 . 05 to q = 1 with a spacing of ∆ q = 0 . 05. Finally, we cover 56 values of initial orbital period, ranging from P orb = 0 . 7 days to P orb = 6105days with a logarithmic spacing of ∆ log 10 ( P orb / 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.  \nWe 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.  \nIn 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 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 sys-  \ntems 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).  \nSeparately, the color of each marker indicates that particular binary's mass transfer history. Systems with sufficiently close initial P 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).  \nComparison 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 orb glyph[lessorsimilar] 5 days tend to result in dynamically unstable mass transfer. (iii) There is a large region of binaries with initial primary mass glyph[similarequal] 40-50 M glyph[circledot] 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.  \nWe 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 massratio slices as in Figure 9, the final rotational rate of the", 'q = 0.7': "<!-- image -->  \nFigure 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.  \n<!-- image -->  \nsecondary (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 angularmomentum accretion during their evolution. Many of them have reached critical rotation, ( ω s /ω 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.  \nWe find a small subset of initially very close systems in the bottom right corner (log 10 ( M 1 /M glyph[circledot] ) > 1 . 75 and log 10 ( P orb / days) glyph[similarequal] 0 . 5) that retain a significant rotational rate even though they avoid mass transfer. In binaries with such tight orbits, tidal forces between the  \nstars 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': "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 CO , and the orbital period P orb . We consider 40 values of initial primary masses, ranging from M 1 = 0 . 5 M glyph[circledot] to M 1 = 120 M glyph[circledot] with a logarithmic spacing of ∆ log 10 ( M 1 /M glyph[circledot] ) = 0 . 06 dex, and 21 values of initial CO masses, ranging from M CO = 1 M glyph[circledot] to M CO = 35 . 88 M glyph[circledot] with a logarithmic spacing of ∆log 10 ( M CO /M glyph[circledot] ) = 0 . 074 dex. Finally, we cover 30 values of initial orbital period, ranging from P orb = 1 . 26 days to P orb = 3162days with a logarithmic spacing of ∆log 10 ( P orb / days) = 0 . 13 dex. Our choice of CO mass range covers massive WD, NS, and BH accretors.  \nFigure 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, ( ω s /ω 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.  \n<!-- image -->  \nOur 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 nondegenerate 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 selfconsistently, following the prescriptions described in Section 4.1. Finally, binaries that never reach the onset of  \nRLO 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 COat the onset of RLO, within a larger infrastructure to completely evolve binaries from ZAMS to double COs.  \nFigure 11 shows two slices of the grid with different CO masses, M CO = 1 . 43 M glyph[circledot] to represent a NS accretor and M CO = 14 . 66 M glyph[circledot] 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 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  \n- Initial RLOF\n- Not converged\n- Unstable RLOF' during MS\n- Stable RLOF during MS  \nFigure 12 shows the relative changes in the accretor masses in the same two slices in M 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 glyph[similarequal] 1 day, and pre-RLO mass ratios in the range q ∼ 1-2 (defining q = M CO , i /M 1 , i ). Despite  \n<!-- image -->", '14.66 Mco Mo': "Figure 11. View of two slices, for two different values of initial CO masses ( M CO = 1 . 43 M glyph[circledot] on the left, M CO = 14 . 66 M glyph[circledot] 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.  \n<!-- image -->  \nstars 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 glyph[lessorsimilar] 1 M glyph[circledot] ) 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.  \nExamining the stability of the mass-transfer phase, Figure 11 shows that nearly every donor star accreting onto a 14 . 66 M glyph[circledot] BH does so stably, whereas only the lower mass accretors ( M glyph[lessorsimilar] 4.5 M glyph[circledot] ) 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 L 2 overflow leading to dynamical instability as done in our work.  \nour 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).  \nThe 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 masstransfer 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.", 'Termination conditions': '- Stable RLOF during postMS\n- Unstable RLOF during stripped He star  \nFigure 12. Relative increase in the mass of the CO ( M CO , f -M CO , i ) /M CO , i due to accretion for systems where the nondegenerate star reached the end of its life. The grid slices are the same as shown in Figure 11. Although accretion is Eddingtonlimited, COs in binaries with pre-RLO mass ratios in the range q ∼ 1-2 (defining q = M CO , i /M 1 , i ) and short initial periods, which will experience long-duration mass-transfer phases, manage to accrete a significant amount of mass.  \n<!-- image -->', '5.7. Binaries consisting of a compact object and a He-rich star': "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 CO , and the orbital period P orb . We consider 40 values of initial primary masses, ranging from M 1 = 0 . 5 M glyph[circledot] to M 1 = 80 M glyph[circledot] with a logarithmic spacing of ∆ log 10 ( M 1 /M glyph[circledot] ) = 0 . 055 dex, and 21 values of initial CO masses, ranging from M CO = 1 M glyph[circledot] to M CO = 35 . 88 M glyph[circledot] with a logarithmic spacing of ∆log 10 ( M CO /M glyph[circledot] ) = 0 . 074 dex. Finally, we cover 47 values of initial orbital period, ranging from P orb = 0 . 02 days to P orb = 1117 . 2 days with a logarithmic spacing of ∆ log 10 ( P orb / days) = 0 . 09 dex. Our procedure for generating these binaries closely follows the process described in Section 5.5 for the grid of binarystar 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.  \nFigure 14 shows an example of two slices of this grid, one corresponding to a NS companion (with M CO glyph[similarequal] 1 . 43 M glyph[circledot] ) and one corresponding to a BH (with a M CO glyph[similarequal] 14 . 66 M glyph[circledot] ). 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.  \nWhen 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 1 ≤ 3 M glyph[circledot] ) mass transfer up to wide orbital periods ( P orb < 10 3 days). This occurs because low-mass Hestars expand their He-rich envelope much farther during their later He-shell and C-burning phases (Figure 8).  \nBoth slices of the grid present two islands of failed simulations, one with M 1 glyph[similarequal] 1 . 8 M glyph[circledot] and P orb of the order of days and another island with M 1 glyph[lessorsimilar] 1 M glyph[circledot] and  \nFigure 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.  \n<!-- image -->  \nP 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, shortP 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 glyph[similarequal] 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.  \nIn 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 selfconsistently 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  \nthe He-stars with the highest specific angular momenta are those with either short P orb or stable mass transfer.  \nThe 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 POSYDON . 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  \nFigure 16 shows a typical example of a binary-star model, initially composed of two H-rich ZAMS stars with masses 10.50 M glyph[circledot] and 5.25 M glyph[circledot] 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.  \n<!-- image -->  \nFigure 14. View of two grid slices for two different values of initial CO masses ( M CO = 1 . 43 M glyph[circledot] on the left, M CO = 14 . 66 M glyph[circledot] 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.  \n<!-- image -->  \nthe 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': 'After having computed our grids of single- and binarystar 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 shorttimescale evolution driven by the κ -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 ( κ max ) set to 0 . 5 cm 2 g -1 . This approximation reduces the failure rate of each binary grid from glyph[similarequal] 10 . 9%, glyph[similarequal] 8 . 0%, and glyph[similarequal] 11 . 8%, for the binary-star grids composed of two Hrich 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 glyph[similarequal] 0.9%, glyph[similarequal] 1.5%, and glyph[similarequal] 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  \ninitial parameter space and compared to our interpolation accuracy (Section 7.5).  \nFor further comparison we show an adjacent binarystar 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 glyph[circledot] , 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.  \nFigure 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.  \n<!-- image -->  \nTable 3. A list of the different stellar types we adopt in POSYDON .', '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 POSYDON , 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.  \nAs 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  \ndescribed 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).  \n<!-- image -->  \n/circledot  \nFigure 16. Typical example of the evolution of a binarystar 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 glyph[circledot] and a 5.25 M glyph[circledot] 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 cm 2 / 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 glyph[circledot] ), which shows that any inaccuracy induced by our opacity approximation is a similar magnitude to the differences between neighboring simulations.  \n<!-- image -->  \nAs 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  \ncodes rely upon the k -type stellar classification introduced in Hurley et al. (2000, based on Tout et al. 1997) In POSYDON , 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.  \nAs 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 POSYDON 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.  \nAs 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).  \nAs a final step, we calculate a number of postprocessed 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 out-  \nFigure 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 glyph[circledot] , 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 nuc ).  \n<!-- image -->  \nput 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.  \nThe 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  \nentire binary-star evolutionary tracks will be addressed in future releases of POSYDON .', '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).  \nFor 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 h . The evolution of the binary is a multivariate time-series given by h i = h ( t i ), with i = 1 , · · · , N , where N is the number of steps, each corresponding to an age t i . Before the downsampling, the independent variable (age), and nonphysical parameters (e.g., model number in MESA ) are excluded from h , while all other parameters are rescaled linearly from 0 to 1.  \nFigure 18. The evolution of a 37 . 6 M glyph[circledot] star with a 22 . 6 M glyph[circledot] companion with an initial P orb of 6 . 82 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 ∼ 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).  \n<!-- image -->  \nThe 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 glyph[epsilon1] :  \ne i = ‖ h i -ˆ h i ‖ < glyph[epsilon1], (14)  \nFigure 19. The accuracy ( P 90 ( 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 ( × ; see excluded parameters annotated with an asterisk in Table 5 and Table 6). As the downsampling threshold glyph[epsilon1] 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 glyph[epsilon1] = 10 -1 , which gives us a compression ratio of glyph[similarequal] 26, while limiting any errors that could be accrued from this process to within a few per cent.  \n<!-- image -->  \nwhere ˆ h i is the interpolated point using the age as the independent variable. We use linear interpolation,  \nˆ h i = h j + t i -t j t j +1 -t j ( h j +1 -h j ) , (15)  \nwhere j and j +1 are adjacent steps of the downsampled time-series so that t j < t i < t j +1 .  \nThe search for the steps in the downsampled data is performed as follows. Initially, we include only the first ( h 1 ) and last ( h N ) points. Then we search for the intermediate point ( h j where 1 < j < N ) with maximum interpolation error. If this error is below the threshold ( e j < glyph[epsilon1] ), then the algorithm has finished; otherwise it  \nincludes this point and continues the search in the two parts of the time-series before ( h 1 to h j ) and after ( h j to 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 < glyph[epsilon1], ∀ i .  \nAdditionally, 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.  \nTo demonstrate the validity of our method, Figure 18 shows an example of the downsampling of a track using the same interpolation error threshold glyph[epsilon1] = 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 twodimensional 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.  \nThe choice of glyph[epsilon1] 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 POSYDON , 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 ∼ 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 ∼ 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 POSYDON 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  \neach 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, { x n } N n =1 , along with its corresponding scalar { y n } N n =1 and class { z n } N n =1 targets. These sets can be columnwise stacked into the matrices X ∈ R 3 × N , Y ∈ R M × N (where M is the number of output quantities), and Z of dimension 4 × N (because we classify each run into one of four broadly defined categories; see Section 7.2). More specifically, each x n ∈ R 3 contains the initial masses and orbital period of the n -th binary in the grid, whereas y n and 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.  \nA 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.  \nFirst, we consider a non-linear transformation of the data using the logarithm: inputs can be transformed as log x i and targets as log y i or log( -y i ) if 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.  \nWe 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.  \nAs a second step, we apply a min-max scaling to the inputs so that the transformed features x t i are confined to the range [ -1 , 1], and we standardize the outputs such that they have zero-mean and unit variance,  \ny t i = y i -y i σ y i . (16)  \nFigure 20. The optimal k in our k NN classification scheme for each of our three grids. The bACC is calculated using 10fold 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.  \n<!-- image -->  \nThe 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 POSYDON approach to evolving binary systems. Therefore, we separate our binaries into four categories based on their mass-transfer histories. The categories are: sta-  \nble 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.  \nFor 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 POSYDON 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.  \nWe 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.  \nFigure 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.  \nFigure 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.  \nTo 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 ∼ 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  \nFigure 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 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.  \n<!-- image -->  \nFigure 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.  \n<!-- image -->  \nis evident that incorrectly classified validation binaries are very rare.  \nTo 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 sys-  \ntems. 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  \ngrid, 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.  \nWe 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 y i , log 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.  \nThe 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 activelearning 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  \nFigure 23. The convective envelope radius R conv . reg . and width D 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.  \n<!-- image -->  \nstar. 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.  \nTo 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 eff of the interpolated star is derived. In the case of inequalities between quantities (Type II), the quantity that must be less than another is limited by the value of the latter. Finally, there are cases where all  \nquantities ought to add to a certain value (Type III). For instance, the fractional chemical abundances of a star's core must, by definition, sum to unity. We ensure these constraints are satisfied by normalizing our interpolated outputs. In one case described in the Table 4 footnotes, a parameter is subject to two separate constraints, in which case we are careful to apply them in the correct order.  \nIn the case of the constraint involving the interpolated quantities R conv . reg . and D conv . reg . , the middle point ( ≡ ( R t , conv . reg + R b , conv . reg ) / 2) and the thickness ( ≡ R t , conv . reg -R b , conv . reg ) of the convective region for the computation of tides, respectively, a special treatment is required. Both quantities must be positive and less than the star's radius. However, constraining them independently as in other Type II constraints does not work as the inner and outer radius of the convective region must both be inside the star: 0 < R conv . reg . -D conv . reg . / 2 < R conv . reg . + D conv . reg . / 2 < R . We decompose this relationship into three inequalities: D conv . reg . / 2 ≥ 0, R conv . reg . ≥ D conv . reg . / 2, and R conv . reg . + D conv . reg . / 2 ≤ R . In the D conv . reg . -R conv . reg . plane, the constraints form a feasibility region in the shape of a triangle with the vertices (0 , 0) , (0 , R ) , and ( R,R/ 2), where R is a fixed values. If the constraints are violated, then the interpolated values lie outside of the triangle. The triangle's centroid ( R/ 3 , R/ 2) and the point corresponding to the interpolated values define a line l . The intersection between l and the border of the  \ntriangle satisfies the constraint inequality, and is used to assign new values to the parameters. Figure 23 provides a pictorial representation of the algorithm.  \nTo assess how often our constraints defined in Table 4 are violated in practice without imposing constraints, we interpolated 3000 binaries for each one of the three grids using the random initial conditions of their validation sets (Section 7.5). For each binary, we checked all the constraints (two checks per binary system, and 23 checks per non-degenerate companion star) and counted the violations. In the case of Type I constraints (equations), we consider violation a relative error of more than 0 . 001 in the inferred quantity. Overall, we found 57,035 violations in the 274,560 checks ( ∼ 20 . 8%) we performed. After applying the algorithm defined here, all violations were corrected.", '7.5. How Accurate are our Interpolation Methods?': "To assess the performance of the interpolation scheme, we use the same validation data sets that we used to evaluate our classification accuracy, described in Section 7.2. Our trained interpolators are applied to the same initial binary parameters as those of the three sets of 3000 binaries comprising our validation sets, one for each of our three binary grids. Since these binaries are not used in the training phase, the difference between this set and our interpolated predictions for them provide an ideal comparison from which we can determine the accuracy of our methods.  \nFigure 24. Interpolation scheme accuracy for ten selected parameters when applied to our grid of two H-rich stars, as calculated using our set of validation binaries. We separate our sample by their different mass-transfer histories to independently evaluate their individual accuracy. Median relative errors ( e r ) indicated by the horizontal lines in each distribution are typically 1% or lower for the stable mass-transfer and unstable mass-transfer cases (top panel) and the no mass-transfer case (green; bottom panel). Improving this accuracy will be a focus of future work.  \n<!-- image -->  \nFigure 25. The interpolation accuracy for the same 10 parameters as in Figure 24 for our grid of H-rich stars with a CO companion. Since we never use the models from this grid that do not undergo RLO, we do not evaluate the no mass-transfer binaries. In most cases, typical errors are 1% or better, but several of the distributions have tails extending towards larger e r .  \n<!-- image -->  \nFigure 24 provides the accuracies for eleven selected binary and stellar parameters for our grid of two Hrich stars evolved from ZAMS. We have split these samples by their mass-transfer histories so we can separately identify our algorithm's accuracy for the stable mass transfer (red) and unstable mass transfer (blue) cases in the top panel and no mass transfer case (green) in the bottom panel. For nearly all parameters and all classes, our median errors are below 1%. Some parameters such as age and J 1 are significantly more accurately interpolated, while others such as M C / O -core , 1 may be some-  \nwhat less accurate. The distributions are quite broad, suggesting that inaccuracies may exist when parameters show sharp variations as a function of input binary parameters. This is particularly apparent for our unstable mass transfer channel, likely a result of the relatively smaller number of simulations that enter unstable mass transfer and the varying evolutionary stages of the donor stars at the onset of the dynamical instability. For instance, the relatively large error distribution for M C/O -core , 1 for our unstable mass-transfer class is likely due to the rapid core growth during the giant  \nphase when donor stars typically enter dynamical instability. Furthermore, despite their large error distributions, some parameters, such as R 1 , have little impact on the evolution of a binary. Whether the primary star's next evolutionary phase is a core collapse (in the case of the stable mass-transfer scenario) or a CE (in the case of unstable mass transfer) the mass at the outermost part of the star has little impact on the binary's outcome. Nevertheless, we plan to improve these accuracies with future enhancements to our interpolation schemes.  \nIn Figure 25 we provide analogous results for our stable mass-transfer and unstable mass-transfer binaries for our grid of H-rich stars with a CO companion. Our models tend to show larger variations in accuracy compared with our grid of two H-rich stars. Median errors tend to range from 1% to 10% with certain parameters such as age and M C / O -core , 1 performing noticeably worse. At the same time certain parameters like M 2 are very accurately determined, as these parameters vary during the evolution of the binaries in this grid (CO companions in this grid typically accrete little mass). One ought to consider the importance of each parameter when evaluating the accuracy of our models. For instance for stars going through unstable MT M He , 1 is a much more important parameter than M C / O , 1 and Y surf has no impact on a binary's future evolution.  \nFinally, Figure 26 shows the accuracies for our grid of He-rich stars with CO companions. Any edges in the distributions are genuine representations of the underlying data. Our trained interpolators provide the most accurate predictions of our three grids; median accuracies are typically between 0.1% and 1% for the stable mass transfer and unstable mass transfer classes, and somewhat better for the no mass transfer class.  \nThe accuracies provided in Figures 24, 25, and 26 all refer to the data sets and associated interpolation objects provided in v1.0 of POSYDON . One could use the POSYDON infrastructure to evolve larger numbers of binaries than we have provided along with v1.0, which would improve our interpolation accuracy. A focus on regions where our interpolation methods are least accurate would provide the largest benefit. Using a combination of active-learning techniques, more complex machine-learning algorithms, and much more computation time, we expect that future versions of POSYDON will only exhibit substantially improved classification and interpolation accuracies (Section 11).", '7.6. Limitations of our Approach': "There are a few limitations of our approach. First, our approach first classifies the binary's type and subsequently performs interpolation. The effect of such a  \ntechnique is that by performing two optimization problems, the second of which relies on the first, it is possible to propagate error throughout the pipeline. Treating the entire problem as one optimization problem has the potential to reduce error.  \nAdditionally, we transform the grid space by logarithmic transformations before performing linear interpolation, which results in a non-linear interpolation model. Such an approach is similar to using a kernel, where a space is transformed through a kernel function to a space in which a linear model allows for accurate modeling of the behavior of the space (Theodoridis & Koutroumbas 2009). A more systematic approach is to consider a kernelized interpolation approach, by applying kernel selection techniques. In the case of a Gaussian process, for example, we may consider a whole family of functions which are specified by a kernel function, to allow for more flexibility on the prior belief of the space (MacKay 2003).  \nFinally, upon finding an interpolated value, we physically enforce the constraints, as detailed in Section 7.4. However, in principle such a technique, which does not consider the constraints in the optimization objective itself does not guarantee an optimal solution subject to the constraints. One way to incorporate the constraints in our model is to add a regularizer term in our loss function to enforce the constraints (Ivezi'c et al. 2014), i.e., the loss function balloons when constraints are violated.", '8. EVOLUTIONARY PROCESSES SEPARATE FROM SINGLE- AND BINARY-STAR MODEL GRIDS': 'Besides computing, processing, classifying, and interpolating the five separate grids of single- and binary-star models, additional steps are required to follow the complete evolution of a stellar binary from ZAMS to double CO formation (and potentially its merger). These are defined by three separate processes: orbital evolution in eccentric, detached binaries, CE evolution, and stellar core-collapse. While the latter two are standard elements of BPS codes, the need for the former requires some explanation. Binaries are intrinsically eccentric after a SN occurs, yet our pre-calculated grids of binarystar models are initiated with circular orbits. Including eccentricity as an input to our MESA models would add an additional dimension to our simulation grids, challenging our computational capabilities. Furthermore, selfconsistently modeling binary mass transfer along with stellar evolution and tides in eccentric orbits is an active area of study (Sepinsky et al. 2007, 2009, 2010; Dosopoulou & Kalogera 2016a,b; Hamers & Dosopoulou  \nFigure 26. The interpolation accuracy for the same 10 parameters as in Figure 24 and Figure 25, but for our grid of Hestars with CO companions. As with our other two grids, this grid has typical median errors below 1%; however, tails of the distribution extend towards larger e r , especially for the stable mass-transfer and unstable mass-transfer cases (top panel). Our no mass-transfer case (bottom panel) is much more accurate than the corresponding binaries in our grid of two H-rich stars in Figure 24. The apparent truncations in the distributions (e.g., M C / O -core , 1 and log 10 ( J 1 )) are genuine representations of the data.  \n<!-- image -->  \n2019), and to date no detailed binary evolution grids have included initially non-circular binaries.  \nNevertheless, in a detached binary, tidal forces cause an eccentric binary to both circularize and synchronize, an effect that must be taken into account, along with other orbital angular-momentum loss processes (e.g., wind mass loss, gravitational radiation and magnetic breaking). To specifically address this, we evolve binaries after the a SN event using a separate process described in Section 8.1. We only switch back to using the pre-calculated grid of binary-star models once RLO occurs.  \nIn the current version of POSYDON , binaries that successfully exit from a CE phase initiated by two nondegenerate stars, are also modeled following the process described in Section 8.1 (see also Figure 1). These are binaries consisted of a H-rich and a He-rich, or two He-rich stars in a close circular orbit. We follow their evolution as detached binaries until one of the two stars reaches core-collapse. In a small fraction of post-CE binaries ( ∼ 0 . 3% of the total population) typically consisted of a low-mass ( glyph[lessorsimilar] 4 M glyph[circledot] ) He-rich and a H-rich MS star, the He-rich star overflows its Roche-lobe as it expands to become a giant He star, and initiates mass transfer onto the H-rich MS star. Since we do not have a grid  \nof detailed binary-star models covering this part of the parameter space, we cannot follow further the evolution of these binaries. We plan to address this in future versions of POSYDON . Finally, the evolution of all binaries that are in orbits wider than what is covered by our grids of detailed binary-star models, and thus will never initiate RLO, is also followed as described in Section 8.1.  \nIn the subsequent sections, we provide details about how we evolve binaries through an eccentric, detached phase, a CE evolution phase, and core collapse.', '8.1. Evolution of eccentric, detached binaries 8.1.1. Matching with a single-star track': "Even though a non-degenerate star in a detached binary is influenced by its companion (for instance, due to tides), we are making the assumption that so long as RLO is not occurring, our single-star, non-rotating models provide reasonable approximations for the evolution of non-degenerate stars in these detached binaries. We first match the non-degenerate star in the detached binary system with the closest model (searching across different masses at all ages) from our single-star (both H-rich and He-rich) evolutionary tracks. The matching is achieved by minimizing the sum of the squares of key parameters describing the structure of a star. These pa-  \nFigure 27. Distributions of the difference between the matching point of a single star model in the beginning of the detached step and the values from the previous step, for various physical quantities of a non-degenerate star. We show the relative difference in mass ∆ M , He-core mass ∆ M He -core and logarithm of the radius ∆ log 10 ( R/R glyph[circledot] ), as well as the difference in the He central and surface abundances Y center and Y surf , respectively. Vertical, dashed lines from left to right delineate the 5-th and 95-th percentiles of the distributions. The distributions all suggest that the non-degenerate stars in detached binaries can be accurately matched to a single star model.  \n<!-- image -->  \n/circledot  \n-  \n/circledot  \n/circledot  \nrameters differ depending on the evolutionary phase of the star.  \nFor the matching process, we distinguish among: (i) MS stars that still have H in their core (with central H mass abundance X center > 0 . 01), (ii) stars that evolved off the MS and retain even a thin H-rich envelope (postMS, with X center < 0 . 01 and surface H mass abundance X surf > 0 . 01), and (iii) evolved stripped stars that are effectively a H-deficient core ( X surf < 0 . 01). In case (i) the parameters whose differences are minimized are the total mass of the star, its central H abundance X center and the radius of the MS star. After the first core collapse of the system occurs, most binary companions are in this state, being the initially less massive secondaries that evolved slower than the primary that has formed a CO. For case (ii), we replace the central H abundance with the He one ( Y center ) and the radius with the mass of the fully developed He core. In case (iii), for stripped He-rich stars, we use as minimized parameters the Hecore mass of the star (which is equal to its total mass), the radius, and its center He mass abundance.  \nWe normalize the chosen quantities, such that they have similar weighting in the minimization process. The normalization factors are chosen from typical ranges of each parameter for the stars that we focus on: 20 M glyph[circledot] for the total mass of MS or post-MS stars; 10 M glyph[circledot] for He-core masses of stripped He-rich stars, 2.0 for log 10 ( R/R glyph[circledot] ); and 1 . 0 for chemical mass-fraction abundances.  \nWe quantify the quality of the matching by calculating the difference of various quantities from the previous step. In Figure 27, we see that the difference between the new interpolated total mass in the be-  \nginning of the detached step and its previous value, ∆ M = M match -M prev , step , is typically better than 0 . 2 M glyph[circledot] . He-core masses are matched even more precisely, to within 0.05 M glyph[circledot] . Other parameters (we show R , Y center , Y surf in Figure 27) are also closely matched, justifying our assumption that the non-degenerate star in a detached binary can be accurately represented by a single star model, at least so long as it remains detached.  \nA CO component in the binary system is treated as a point mass and does not need a matching process. We also keep the orbital separation and eccentricity constant in this transition, and thus any small difference in the matched star's mass from the previous step results in a small relative change in the orbital period. In addition, by conserving the spin angular momentum on non-degenerate stars during the matching step, we can determine their initial angular frequency Ω at the beginning of the detached step.", '8.1.2. Further evolution of an eccentric detached binary system': "Once matched with single star models, we evolve the stars in detached binaries as essentially single stars, accounting for their effects on the binary's orbit. For a non-degenerate star, its parameters (e.g., mass, radius, moment of inertia) are evolved according to its interpolated stellar track. At the same time, its spin Ω as well as the system's a and e are evolved solving a set of coupled ordinary differential equations that describe their rate of change due to wind mass loss, tides, magnetic braking, and gravitational radiation. That we assume the back reaction of each of these effects does not significantly impact the internal structure of either star so that our single-star models are sufficiently accurate. For  \ninstance, although we follow each star's spin using the moment of inertia of the single-star models, we cannot account during this phase for the star's internal differential rotation and effects such as rotational mixing.  \nThis approach can only handle scenarios where no RLO mass transfer takes place between the two stars; as soon as a H-rich star enters RLO, we stop the binary's evolution and transition to our grid of MESA mass transfer simulations described in Section 5.6. Binaries are assumed to circularize instantaneously upon RLO with an orbital separation equal to the binary's separation at periastron. Alternatively, we also allow for a user to choose to circularize the orbit assuming angular momentum is conserved. 6 Likewise, this step of evolution also ends if a non-degenerate star reaches the end of its life, in which case the binary is sent to a step that handles core collapse (Section 8.3).  \nAthird stopping condition exists for binaries consisted of two COs, which merge due to gravitational wave radiation. Note that when modeling two COs, only effects due to gravitational wave radiation contribute to orbital evolution. We calculate the orbital decay until the merger or the maximum simulation time is reached.  \nOrbital evolution during the detached step is due to a combination of the relevant pieces of physics, which we assume have additive effect:  \n˙ a = ˙ a wind + ˙ a tides , 1 + ˙ a tides , 2 + ˙ a GR , (17)  \n˙ e = ˙ e tides , 1 + ˙ e tides , 2 + ˙ e GR , (18)  \n˙ Ω 1 = ˙ Ω wind , 1 + ˙ Ω intertia , 1 + ˙ Ω tides , 1 + ˙ Ω mb , 1 , (19)  \n˙ Ω 2 = ˙ Ω wind , 2 + ˙ Ω intertia , 2 + ˙ Ω tides , 2 + ˙ Ω mb , 2 . (20)  \nThe orbital separation, eccentricity and stellar spins are evolved using a set of self-consistent, coupled equations. We describe each of the terms in Eq. (17)-(20) below.  \nMass Loss: We ignore mass accretion onto a star (either non-degenerate or CO) from a companion star's wind. So non-degenerate stars will only lose mass due to their own stellar winds, with the mass lost carrying away the specific orbital angular momentum of the massloosing star (Jeans-mode mass loss; for a review, see Tauris & van den Heuvel 2006):  \n˙ a wind = -a ˙ M w , 1 + ˙ M w , 2 M 1 + M 2 . (21)  \nFor CO binary components, ˙ M w = 0, while in general for spherical, isotropic fast winds, the orbit-averaged ˙ e  \ndue to winds is zero. We discuss the effect of mass loss on stellar spin later in this section.  \nTides: Changes in the orbit's period and eccentricity due to tidal forces are described by a set of ordinary differential equations, according to Hut (1981). In order to be able to compute tidal spin-orbit coupling, we treat the donor star of mass M , radius R and of moment of inertia I , as a solid body rotating with angular velocity Ω. The initial angular-momentum budget of the nondegenerate star is assumed to be the same as from the end of the previous step.  \nThe change of the orbital separation due to tidal forces on the first star (subscript 1) is given by:  \n˙ a tides , 1 = -6 ( k T ) 1 M 1 ( M 1 + M 2 ) M 2 1 ( R 1 a ) 8 a (1 -e 2 ) 15 / 2 × [ f 1 ( e 2 ) -( 1 -e 2 ) 3 / 2 f 2 ( e 2 ) Ω 1 Ω orb ] , (22)  \nwhere Ω orb = 2 π/P orb is the mean orbital angular velocity. When both stars are non-degenerate, they each have their own contribution to the orbit's evolution. Therefore an analogous equation exists providing ˙ a tides , 2 , where R 1 is replaced with R 2 , the k/T term is calculated for the secondary star, and M 1 and M 2 are switched.  \nThe k/T term in Eq. (22) depends on a star's structure and the associated physical process of tidal dissipation. We calculated them separately for dynamical and equilibrium tides, in the same way as in our detailed, binarystar model grids (Section 4.1) described in Eq. (4) and Eq. (6), respectively. We apply the maximum of these two at each timestep.  \nThe change of the orbital eccentricity and the stellar spin from tidal forces is also calculated as  \n˙ e tides , 1 = -27 ( k T ) 1 M 2 ( M 1 + M 2 ) M 2 1 ( R 1 a ) 8 e (1 -e 2 ) 13 / 2 × [ f 3 ( e 2 ) -11 18 ( 1 -e 2 ) 3 / 2 f 4 ( e 2 ) Ω 1 Ω orb ] , (23)  \n˙ Ω tides , 1 =3 ( k T ) 1 ( M 2 M 1 ) 2 ( M 1 R 2 1 I 1 )( R 1 a ) 6 Ω orb (1 -e 2 ) 6 × [ f 2 ( e 2 ) -( 1 -e 2 ) 3 / 2 f 5 ( e 2 ) Ω 1 Ω orb ] . (24)  \nAs in the Eq. (22) for ˙ a tides , when the companion star is non-degenerate, ˙ e tides , 2 and ˙ Ω tides , 2 terms exist, which can be calculated by switching subscripts 1 and 2 in Eq. (23) and (24). The f i ( e 2 ), i = 1 -5 terms in Eq. (22), (23), and (24) can all be found in Hut (1981).  \nStellar Evolution: During the detached orbital evolution, we also take into account the change of stellar spin due to the evolution of the stars themselves. This includes spin down because of wind mass loss that carries away the specific angular momentum of the star's surface, as well as changes in its spin due to the evolution of its moment of inertia due to changes in its internal structure,  \n˙ Ω wind+inertia , 1 = 2 3 R 2 1 Ω 1 I 1 ˙ M 1 -Ω 1 I 1 ˙ I 1 , (25)  \nwhere ˙ I the rate of change of its moment of inertia. For binaries in which both stars are non-degenerate, an equation equivalent to Eq. (25) exists for ˙ Ω wind+inertia , 2 .  \nNon-degenerate stars tend to spin down, due to their expansion and their wind mass loss. However, they may also be spun up in phases where they contract. Therefore, for numerical-stability reasons we artificially limit the second term in Eq. (25) to +100 rad yr -2 (usually reached during a sudden contraction to form a WD). Although we take into account the effect of stellar spin on the orbit via tidal spin-orbit coupling, we do not include effects of spin on the stellar structure, such as stellar deformation, rotational mixing or rotationallyenhanced winds.  \nMagnetic Braking: In case the binary contains lowmass non-degenerate stars, spin-down due to magnetic braking can become important i.e., the loss of spin angular momentum due to ionized material ejected from the star that is trapped in its own radial magnetic field. The spin-down rate is given by  \n˙ Ω mb , 1 = τ mb , 1 I 1 , (26)  \nwhere τ mb , 1 is the torque calculated as in Eq. (36) of Rappaport et al. (1983),  \nτ mb , 1 = -6 . 82 × 10 34 dyncm ( M 1 M glyph[circledot] )( R 1 R glyph[circledot] ) γ mb ( Ω 1 / 2 π 1 / day ) 3 (27)  \nwith γ mb = 4 (Verbunt & Zwaan 1981). We apply the full torque to all non-degenerate stars below 1 . 3 M glyph[circledot] and assume no magnetic braking for stars above 1 . 5 M glyph[circledot] , with linear interpolation in-between. Again, for binaries with two low-mass, non-degenerate stars, an equation equivalent to Eq. (26) exists for ˙ Ω mb , 2 .  \nGravitational Wave Radiation: Finally, we take into account orbital changes due to gravitational radiation (Peters 1964; Junker & Schaefer 1992),  \n˙ α GR = -2 15 νc (1 -e 2 ) 9 / 2 [ G ( M 1 + M 2 ) ac 2 ] 3 [ g 1 ( e 2 ) -1 28 G ( M 1 + M 2 ) ac 2 g 2 ( e 2 ) ] , (28)  \n˙ e GR = -1 15 ν c 3 G ( M 1 + M 2 ) [ G ( M 1 + M 2 ) ac 2 ] 4 e (1 -e 2 ) 7 / 2 [ g 3 ( e 2 ) -1 56 G ( M 1 + M 2 ) ac 2 g 4 ( e 2 ) ] ,  \n(29)  \nwhere ν = M 1 M 2 / ( M 1 + M 2 ) 2 , and the functions g i ( e 2 ), i = 1 -4, are defined as:  \ng 1 ( e 2 ) = ( 96 + 292 e 2 +37 e 4 ) ( 1 -e 2 ) , (30) g 2 ( e 2 ) =(14008 + 4704 ν ) + (80124 + 21560 ν ) e 2 + (17325 + 10458 ν ) e 4 -1 (5501 1036 ν ) e 6 , (31)  \n2 -  \ng 3 ( e 2 ) = ( 304 + 121 e 2 ) ( 1 -e 2 ) , (32) g 4 ( e 2 ) =8(16705 + 4676 ν ) + 12 (9082 + 2807 ν ) e 2 -(25211 + 3388 ν ) e 4 . (33)  \nThese general-relativistic terms of orbital evolution are usually negligible apart from cases of close binaries. For binaries consisted of two COs, only Eq. (28)-(33) govern the evolution of the binary's orbit.", '8.2. Common-envelope evolution': "Binary interactions can lead to a dynamically unstable mass-transfer phase (Ivanova et al. 2013, 2020). We have described in Section 4.2.4 all the conditions that are assumed to trigger an unstable mass-transfer episode: a maximum mass-transfer rate of 0 . 1 M glyph[circledot] yr -1 , L 2 overflow, a contact phase with a post-MS star, or exceeding the threshold of the trapping radius during accretion onto a CO.  \nIf the donor star that triggered the unstable phase is in its MS phase or is a stripped He star during its He core-burning phase, we assume that the two stars promptly merge, as no distinct core has formed yet in its interior. In v1.0 of POSYDON we do not follow the further evolution of stellar merger products.  \nFor all the other donor stellar states, a trigger of unstable mass transfer is assumed to lead to a CE phase. If the donor has a H-rich envelope at the beginning of the phase, this envelope is considered to form the CE, inside of which the donor's He-rich core and its binary companion will spiral-in. For stripped donors, the He-rich envelope engulfs the companion which spirals in around the donor's C/O core. In case the companion star also has a giant-like structure with a distinct core-envelope separation, (i.e., anything but a MS star, a He star in its He-MS, or a CO), then its envelope also may contribute in a (double) CE.  \n,  \nFigure 28. Evolution of λ CE parameter of POSYDON singlestar models of different initial masses. For these calculations, the assumed core-envelope boundary is located at the point where the H mass fraction drops below 10%. The triangle, diamond, and star markers represent the start of shell H burning, the start of core He burning, and the end of core He burning, respectively.  \n<!-- image -->  \n/circledot  \nThe outcome of the CE phase is calculated using the α CE -λ CE prescription (Webbink 1984; Livio & Soker 1988), which equates a fraction α CE of the orbital energy released during the spiral-in with the binding energy of the CE. The α CE parameter is set equal to 1 in the example population runs shown in Section 10, following previous population synthesis works (e.g., Hurley et al. 2002), but is, in general, a free parameter in POSYDON .  \nThe parameter λ CE has been introduced to parametrize the binding energy of the envelope using the total stellar radius and mass (de Kool 1990). In POSYDON , λ CE values are calculated from the detailed stellar profile of the donor star at the beginning of CE (or of both stars, in case of a double CE). This is an important quantitative improvement of POSYDON , compared to pBPS codes. The latter need to adopt λ CE -value fits from the literature, based on single-star models with often inconsistent stellar-physics, and apply them to post-interacting stars. In our common CE energy calculation, we integrate both the gravitational and the internal energy of the envelope from the detailed stellar structure profiles of our binary models, subtracting the recombination energy from the internal energy. In Figure 28 we show the λ CE for a few example single-star POSYDON models. The parameter λ CE tends to decrease as the star evolves and expands as a giant. However, for initially very massive ( glyph[greaterorsimilar] 20 M glyph[circledot] ) stars that strip their H-rich layers due to their own wind mass loss, λ CE increases again. We find comparable trends and  \nvalues with other works that study the detailed stellar structures of giant stars (e.g., Kruckow et al. 2016; Klencki et al. 2021).  \nThe binding energy of the envelope (and thus the outcome of the CE phase) is sensitive to the exact assumed core-envelope boundary of the donor, as the deeper envelope layers tend to be the most tightly bound (Dewi & Tauris 2000; Ivanova et al. 2013; Fragos et al. 2019). For this reason, we allow for different core-envelope boundaries, defined for H-rich stars as the outermost layer where the H mass fraction drops below 0 . 3, 0 . 1 (default), and 0 . 01 and for stripped-stars when the sum of H and He drops below 0 . 1.  \nGiven the properties of the binary at the onset of the CE, the assumed α CE value and the estimated λ CE , one can calculate how much a binary's orbit shrinks in order for the released orbital energy to fully unbind the CE. The final post-CE orbital separation a post , CE is given by solving (Webbink 1984):  \nGM don , core M acc 2 a post , CE -GM don M acc 2 a pre , CE = GM don M don , env α CE λ CE R don , (34)  \nwhere M don , M don , core , and M don , env are the total, core and envelope masses of the donor star, R don the donor star's radius, M acc the mass of the accreting star, and a pre , CE the orbital separation at the onset of the CE. If the final estimated a post , CE is such that neither the accreting star nor the stripped core of the primary star are filling their respective Roche lobes, then the CE is consider successful and results in a detached, circular, tight binary. Alternatively, the binary is assumed to merge, and its evolution is not further followed in v1.0 of POSYDON .  \nOne complication with the flexibility we offer regarding the core-envelope boundary definition is that the post-CE stripped donor star might still contain some H in its outer layers, while in the next evolutionary steps we assume that the H envelope is fully removed. Exactly how much H remains depends on a user's choice of 0.01, 0.1, or 0.3 for a fractional H abundance when defining the core-mass boundary. We account for this inconsistency by assuming that either the remaining Hrich layers are either removed by stellar winds or these layers re-expand after the CE and are removed via stable mass transfer (e.g., Fragos et al. 2019). Both assumptions result in slight corrections to the post-CE donor masses and orbital separation. Although the former assumption is the default one in POSYDON , we find they both lead to correction at the level of only a few percent, and thus the choice between the two is in practice inconsequential.", '8.3. Core-collapse and compact-object formation': 'The end fate of stars primarily depends on their masses. The most massive stars undergo all nuclear burning phases (hydrogen, helium, carbon, neon, oxygen, silicon) up to the formation of an iron core . The iron core keeps growing by silicon shell burning to a mass of around the Chandrasekhar mass limit ∼ 1 . 44 M glyph[circledot] when electron degeneracy pressure can no longer stabilize the core and it collapses. This runaway process can lead to the explosion of a star in a SN or to a direct collapse into a BH and it is known as core-collapse SN (CCSN) (see Janka et al. 2007, for a review).  \nLower mass stars do not complete all nuclear burning phases. For stars which do not ignite oxygen but for which their He-cores masses are 1 . 4 M glyph[circledot] glyph[lessorsimilar] M He -core glyph[lessorsimilar] 2 . 5 M glyph[circledot] (Podsiadlowski et al. 2004) we assume a star collapses into a NS in an electron-capture SN (ECSN). POSYDON alternatively includes the option to determine whether a star undergoes an ECSN based on its C/O core mass: 1 . 37 M glyph[circledot] glyph[lessorsimilar] M C / O -core glyph[lessorsimilar] 1 . 43 M glyph[circledot] (Tauris et al. 2015). Stars with core masses below the lower limit for ECSN evolve into white dwarfs.  \nIn Figure 29, we show, for a slice at a fixed initial massratio q = 0 . 7 of the binary-star model grid composed of two H-rich stars, the core-collapse type as a function of initial orbital period primary star mass. The transition region between the ECSN and CCSN, occurs at ZAMS masses of glyph[similarequal] 8 M glyph[circledot] (consistent with previous studies, e.g., Nomoto 1984; Jones et al. 2014), but depends somewhat on the initial P orb .', '8.3.1. Pulsational pair-instability SN': 'During the post-carbon burning phase of massive stars (not modelled here), photons produced in the core can be energetic enough to produce electron-positron pairs, softening the equation of state and diminishing the pressure support of the core (Woosley et al. 2007, and references therein). In such stars the core rapidly contracts and the temperature increases, leading to explosive oxygen burning (e.g., Woosley & Heger 2015) that creates a series of energetic pulses which eject material from the star surface. This phenomenon of material ejection due to pulses is known as pulsational pair instability SN (PPISN) and occurs for stars with He-core masses in the range ∼ [32 , 64] M glyph[circledot] , (Yoshida et al. 2016; Woosley 2017; Marchant et al. 2019; Renzo et al. 2020). For more massive stars with He-core mass in ∼ [61 , 124] M glyph[circledot] , the first pulse is so energetic that can unbind and destroy the whole star in a so-called pair instability SN (PISN; Fowler & Hoyle 1964; Rakavy & Shaviv 1967; Barkat et al. 1967), leaving no remnant behind.  \n<!-- image -->  \nFigure 29. The core-collapse type for the q = 0 . 7 slice of our binary-star model grid composed of two H-rich stars. We distinguish between WD formation, electron caption SN (ECSN) following Podsiadlowski et al. (2005), core-collapse SN (CCSN) and pair-pulsational instability SN (PPISN) following Marchant et al. (2019). Models that did not reach the end of stellar evolution are indicated in black.  \n<!-- image -->  \nTo identify systems that will undergo PPISN and PISN, we adopt a polynomial fit, as implemented in Breivik et al. (2020, see Eq. 4 therein), to MESA singlestar simulations (at Z = 0 . 1 Z glyph[circledot] ) by Marchant et al. (2019, see Table 1 therein). This fitting formula is used to map the He core mass at carbon depletion in the range 31 . 99 M glyph[circledot] ≤ M He -core ≤ 61 . 10 M glyph[circledot] to the stellar mass collapsing to form the CO. We use PPISN models computed at 1/10-th solar metallicity, as it was shown that such a limit is independent of metallicity (Farmer et al. 2019), while highly dependent on the uncertain 12 C( α, γ ) 16 O reaction rates (Farmer et al. 2020). In our case, these reaction rates follow Cyburt et al. (2010), consistent with the rates used by Marchant et al. (2019).  \nIn Figure 29, we can identify two systems that enter the regime of PPISN as they possess a He-core with mass slightly larger than 31 . 99 M glyph[circledot] at carbon depletion. Other mass ratio slices show a few more similar systems but without becoming statistically relevant. We expect PPISNe to be more present at sub-solar metallicities, as stellar wind-mass loss at Z glyph[circledot] prevent the stars from  \nreaching carbon depletion with a He-core mass in the relevant mass range for PPISNe.', '8.3.2. Remnant baryonic mass': 'In this version of POSYDON , we calculate the mass left behind by the collapse using different models: (i) direct collapse where all the stellar mass is conserved; (ii) fits to the results of two-dimensional core-collapse models of Fryer et al. (2012); (iii) nearest neighbor interpolations of the results of the detailed one-dimensional corecollapse models of Sukhbold et al. (2016), or (iv) with the explodability criteria of Patton & Sukhbold (2020). The last is our default option.  \nFryer et al. (2012) presents two mechanisms which are known as rapid and delayed based on how quickly convective instabilities are expected to grow after core bounce. The rapid prescription produces a mass gap between BHs and NSs by assuming strong convection which allows instabilities to grow quickly after core bounce, producing a more energetic SN explosion. In contrast, the delayed mechanism produces a continuous spectrum of compact remnant masses. Both prescriptions determine the baryonic mass of the compact remnant M rembar given the pre-SN C/O core mass, M CO -core . More precisely, M C / O -core determines whether the star explodes into a SN, and what fraction, f fb , of the ejected mass falls back onto the CO. In the case that the star directly collapses to form a BH, f fb = 1. For the rapid prescription, direct collapse occurs for M C / O -core ≥ 7 . 6 M glyph[circledot] while for the delayed prescriptions direct collapse occurs for M C / O -core ≥ 11 M glyph[circledot] .  \nIn Sukhbold et al. (2016), the outcome of the collapse of their pre-SN models has been calibrated against the well-studied SN 1987A progenitor. We have implemented several of their SN engine calibrations, namely N20, S19.8, W15, W20 and W18, although for the training of the initial-final interpolation we only consider the default N20 option, which is the most optimistic option for successful explosions. In contrast to Fryer et al. (2012) results, Sukhbold et al. (2016) finds sharply varying behavior between the initial star mass and the final core properties, linked to convective carbon-burning episodes occurring in the later evolutionary phases. This results in a region of the parameter space where the outcome of the collapse, i.e., NS and BH formation, appears stochastic in its nature. To mitigate interpolation errors, we determine the remnant baryonic mass of a collapsing star using a nearest neighbor technique on the He-core mass at carbon depletion to map our stars to the Sukhbold et al. (2016) simulation results.  \nIn the Patton & Sukhbold (2020) prescription, the C/O core mass and the average carbon abundance of the  \ncore at carbon ignition are used to determine the explodability of the pre-SN core. For every single- and binarystar model in our grids, we store these two values, and by applying a k -nearest neighbour interpolation, with k = 5, we map to the explodability parameters M 4 and µ 4 from Ertl et al. (2016), as described in Patton & Sukhbold (2020). These two explodability parameters allow us to infer whether a SN is successful, and, if so, we estimate the resulting NS mass to be approximately equal to M 4 . We assume that BHs are produced only from failed explosions which result in a direct collapse. Finally, for the Patton & Sukhbold (2020) prescription, we have implemented the same SN engine options as for (Sukhbold et al. 2016) with N20 as our default option, using the updated calibration from Ertl et al. (2020).  \nIn Figure 30, we show a comparison between the final CO state for the same grid slice as Figure 29, as predicted by Fryer et al. (2012) delayed prescription compared and the Patton & Sukhbold (2020) prescription, based on the N20 engine. The differences between our choice of SN prescription are slight, but noticable when focusing on the NS/BH boundary. The Fryer et al. (2012) delayed prescription produces BHs for somewhat less massive stars, while the Patton & Sukhbold (2020) prescription shows a more variable boundary between NSs and BHs.  \nIn both Sukhbold et al. (2016) and Patton & Sukhbold (2020) prescriptions we assume fallback fractions of f fb = 1 for BHs and f fb = 0 for NSs. For Fryer et al. (2012) prescriptions the fallback fractions are computed explicitly, with the exception of NS ECSN where we assume f fb = 0.', '8.3.3. CO gravitational mass': 'To convert the remnant baryonic mass to gravitational mass, we use the prescription by Zevin et al. (2020), which is an updated version of the one by Lattimer & Yahil (1989) based on the neutrino observations of SN1987A. This new conversion caps the maximum neutrino mass loss to 0 . 5 M glyph[circledot] (C. Fryer, private communication) and removes any artificial discontinuity in the mass spectrum between NS and BH formation (in the case of direct collapse or the Fryer et al. (2012) delayed mechanism) as  \nM grav =            20 3 (√ 1 + 0 . 3 M rembar M glyph[circledot] -1 ) M glyph[circledot] , M rembar -M grav < 0 . 5 M glyph[circledot] M rembar -0 . 5 M glyph[circledot] otherwise . (35)  \nIf M grav < 2 . 5 M glyph[circledot] we classify the CO as a NS otherwise as a BH. There is a large uncertainty in the ex-  \nFigure 30. The CO state for the q = 0 . 7 grid slice of HMS-HMS MESA simulations in the initial primary mass-orbital period plane. We distinguish between the COs: white dwarf (WD), neutron star (NS) and black hole (BH) according to the legend. Models that did not reach the end of the stellar evolution are indicated in black. The two panels compare the Fryer et al. (2012) delayed core-collapse mechanism (left) with the outcome of Patton & Sukhbold (2020) N20 core-collapse engine (right).  \n<!-- image -->  \nact maximum NS mass and this range spans 2 . 0-2 . 7 M glyph[circledot] (Lattimer & Prakash 2010; Margalit & Metzger 2017; Rezzolla et al. 2018; Ai et al. 2020; Shao et al. 2020; Lim et al. 2021; Miller et al. 2021; Raaijmakers et al. 2021). In POSYDON , the maximum neutron star mass is set to M max NS = 2 . 5 M glyph[circledot] (see discussion in Abbott et al. 2020b, and references therein).  \nIn the left panel of Figure 31, we show the gravitational mass of the CO as predicted by Patton & Sukhbold (2020) N20 prescription for the same grid slice as Figure 30.  \nFinally, in the case of BH formation, POSYDON also allows us to take into account the detailed internal structure of the star at the moment of collapse. As illustrated in the next section, we can then make an estimate of the spin of the CO taking into account the angular momentum profile of collapsing star.', '8.3.4. Birth spins of COs': "We estimate the spin of the resulting BH following the collapse of the stellar profile as presented in Bavera et al. (2021, Appendix D). For convenience we summarize here the key assumption of this procedure. The final mass and spin of the BH resulting from the collapse is calculated by following the accretion history of M rembar soon after the direct collapse of the central core of the  \nstar which forms a proto-BH of mass 2 . 5 M glyph[circledot] . The mass lost in neutrinos during the formation of this proto-BH, also carry away specific angular momentum equal of that collapsing central part of the core that forms the protoBH. If M rembar < M star we assume the ejected mass takes away the outer layers of the star.  \nThe angular momentum content of the infalling material can in principle support the formation of an accretion disk. We consider a collapsing star to be a collection of shells with radius r , mass m shell and angular velocity Ω shell that falls one by one onto the central BH. A shell of mass is accreted by the BH once it reaches the BH's event horizon. The specific angular momentum of the in-falling material, j ( r, θ ) = Ω shell ( r ) r 2 sin( θ ), where θ is the polar angle, determines the properties of the accretion flow. Disk formation occurs when the specific angular momentum of the shell j ( r, θ ) exceeds the specific angular momentum of the ISCO j ISCO . This condition can be redefined as the polar angle at which disk formation occurs as  \nθ disk = arcsin ( j ISCO Ω shell ( r ) r 2 ) 1 / 2 . (36)  \nThe portion of the shell with θ < θ disk will collapse directly onto the BH on a dynamical timescale transferring j ( r, θ ) to the hole, while the portion of the shell  \nFigure 31. The CO mass M CO1 , and spin a CO1 for the q = 0 . 7 grid slice of HMS-HMS MESA simulations in the initial primary mass-orbital period plane as predicted by the Patton & Sukhbold (2020) N20 engine. We assume that both the neutrino mass loss up to 0 . 5 M glyph[circledot] and the ejected mass during core collapse carries away the corresponding angular momentum. Spinning BHs are formed in binary systems avoiding mass transfer or undergoing stable mass transfer during contact phase or Case A mass transfer, see Figure 9.  \n<!-- image -->  \nwith θ ≥ θ disk will form a disk and transfer only j ISCO to the BH. The disk will be accreted on a viscous timescale which is assumed to be much smaller than the fallback timescale of the following shell (Batta & Ramirez-Ruiz 2019).  \nTherefore each collapsing shell contributes to the angular momentum of the BH by  \nJ shell ≡ J direct + J disk = = ∫ θ disk 0 M shell Ω shell ( r ) r 2 sin 3 ( θ ) d θ + + ∫ π/ 2 θ disk M shell j ISCO sin( θ ) d θ . (37)  \nThe mass-energy accreted onto the BH from the disk is M disk = ε M shell cos( θ disk ) while the fraction η ≡ 1 -ε = 1 -[1 -2 GM BH / (3 c 2 r ISCO )] 1 / 2 is radiated away (Bardeen 1970; Thorne 1974). Here, r ISCO is the radius of the ISCO of the accreting BH. This means that the resultant BH will have mass smaller than M grav as a fraction of the disk will be radiated away. The dimen-  \nss spin parameter of the BH is updated after each shell is accreted onto the BH with the following relation  \na = cJ BH GM 2 BH , (38)  \nwhere J BH is the angular momentum of the BH and M BH its mass after accreting the directly infalling part of the shell and, if formed, the thin disk.  \nThe presented treatment is applicable only to the case of BH formation. For simplicity and the lack of firm alternatives in this version of POSYDON , we assume a zero spin for NSs.  \nIn the right panel of Figure 31, we show the CO spin for the q = 0 . 7 mass ratio slice of the binary-star grid composed of two H-rich stars, as predicted by the stellar profile collapse assuming the remnant baryonic mass is determined with Patton & Sukhbold (2020) N20 prescription.  \n8.3.5. SN kicks  \nDuring a SN, the binary system experiences abrupt mass loss, away from the center of mass, affecting its orbital parameters (Blaauw 1961; Boersma 1961). Furthermore, asymmetric ejection of matter (Janka & Mueller 1994; Burrows & Hayes 1996; Janka 2013) or asymmetric emission of neutrinos (Bisnovatyi-Kogan 1993; Socrates et al. 2005) can provide a momentum kick to the newly formed CO. Here we assume that the magnitudes of the asymmetric kicks ( v k ) are drawn from a Maxwellian distribution with dispersion σ :  \nf ( v k ) = √ 2 π v 2 k σ 3 exp ( -v 2 k 2 σ 2 ) . (39)  \nAs our fiducial assumption we take σ CCSN = 265kms -1 (Hobbs et al. 2005) and σ ECSN = 20kms -1 (Giacobbo & Mapelli 2019) for CCSN and ECSN, respectively. However, these velocities are free parameters. POSYDON supports multiple kicks rescaling options, e.g., if the prescription used to calculate the remnant baryonic mass assume mass loss, i.e., the fallback mass fraction f fb < 1, the kick is then rescaled by 1 -f fb (Fryer et al. 2012). Alternatively BH kicks are rescaled by a factor 1 . 4 M glyph[circledot] /M BH (using the gravitational mass) while NS kicks are not rescaled (our default option). Finally a user can opt to either not rescale any kicks or turn off SN kicks altogether.  \nThese kicks can tilt the orbit of the binaries, add eccentricity or disrupt it. We take into account all these orbital changes including orbital changes for eccentric binaries following the analytical calculations of Kalogera (1996); Wong et al. (2012).  \nWe assume the collapsing star to lie on the origin of the coordinate system moving in direction of positive y -axis. The companion lies on the negative x -axis and x -axis completes the right-handed coordinate system (cf. Kalogera 1996, Figure 1). The semi-major axis after the kick a f is computed given the instantaneous orbital separation r i = a i [1 -e i cos( E )] pre-SN, where E is the eccentric anomaly, as  \na f = ( 2 r i -v 2 k + v 2 r +2 v y k v r GM f tot ) -1 , (40)  \nwhere M f tot is the binary total stellar mass after the core collapse, v r is the pre-SN velocity of the collapsing star relative to the companion directed along the positive y -axis and v y k the y -axis component of the kick. The eccentricity after the kick is then  \ne f = { 1 -( v z k ) 2 +[sin( ψ )( v r + v y k ) -cos( ψ ) v x k ] 2 GM i tot a f r 2 i } 1 / 2 , (41)  \nwhere M i tot is the binary total stellar mass before the core collapse, ψ is the polar angle of the position vector of the collapsed star with respect to its pre-SN orbital velocity in the companion's reference frame and  \nsin( ψ ) = √ GM i tot (1 -e 2 i ) a i v r r i . (42)  \nIn the above equations, in the case the CO receive no natal kick v k = 0kms -1 but the star loses some mass, the orbit is still readjusted to conserve Kepler's third law.  \nWe consider a binary is disrupted if it does not satisfy the condition that demands the post-SN orbit passes through the pre-SN position (Willems et al. 2005; Flannery & van den Heuvel 1975),  \n1 -e f ≤ r i a f ≤ 1 + e f , (43)  \nor if it is outside the limits of the amount of orbital contraction or expansion that can take place for a given amount of mass loss and a given magnitude of the kick velocity (Kalogera & Lorimer 2000; Willems et al. 2005)  \n2 -M i tot M f tot ( v k v r +1 ) 2 < r i A f < 2 -M i tot M f tot ( v k v r -1 ) 2 . (44)  \nFinally, we also verify that e f does not exceed 1 or that the argument of the square root in Eq. (41) does not become negative, if this is the case the binary is considered to be disrupted.", '9. HOW POSYDON EVOLVES AN INDIVIDUAL BINARY SYSTEM': "To evolve a single binary within POSYDON , we use a hierarchy of classes. Every binary system is represented as a BinaryStar class containing two SingleStar classes, each with attributes that define their current state. To evolve the binary through each step, we have implemented a Pythonic flow, which takes the combination of a binary's state and event, and each stars' state to direct a particular binary to its next step. We have a complete flow set as default, which can self-consistently track the evolution of binaries from their ZAMS state comprised of two H-stars, through all parts of the evolutionary tree shown in Figure 1.  \nAll steps in POSYDON are Python classes that update a binary via the user-defined call method. Steps are implemented based on on-the-fly calculations for evolutionary phases such as the CE or the core collapse, or based on pre-calculated grids of detailed binary-star models, which is a novel component of POSYDON . In the latter  \ncase, to estimate the evolution of systems for which no detailed model exists for the exact initial binary properties, we use initial-final classification and interpolation algorithms, trained on the grids of detailed binary-star models (Section 7) or alternatively a nearest neighbor matching scheme. As the binary and its component stars evolve through these steps, the BinaryStar and SingleStar characteristics are appended to the objects, so every binary maintains a historical record of its evolution.  \nIf any system enters a phase that does not require further evolution, we use an end step to halt the evolution. This is used for binaries that either merge, disrupt, or reach the maximum physical time. For binaries whose evolution time ends in the middle of a step based on pre-calculated grid of models (e.g., if our star formation history randomly generates an initialization time for a binary within a few megayear of the end run time), we can no longer use our pre-trained classification and interpolation algorithms as these only apply to the end state of the binary. In these cases, we instead use the system's nearest-neighbor pre-computed track directly (interpolating full binary tracks is non-trivial and is being investigated for future versions of POSYDON ). This is the default behavior of all MESA grid steps, with various classification and interpolation algorithms ready to use.  \nFigure 32 depicts the complete evolution of one particular binary from ZAMS to the formation of a BBH system that merges within a Hubble time. Each vertical, colored band indicates an extended stage of a binary's evolution, whereas the two core-collapse events and the CE phase are essentially instantaneous processes, occurring in between the other, extended phases.  \nThe initial masses of the system are M 1 = 97 . 07 M glyph[circledot] , M 2 = 39 . 86 M glyph[circledot] in a P orb = 39 . 45 days circular orbit. The first part of the evolution of the system is based on the HMS-HMS binary grid (Section 5.5), and its subsequent evolution is followed either through the nearest neighbor interpolation (orange lines) or through initialfinal interpolation (black dots; using our classification and interpolation methods) for the same initial configuration. In both cases we adopt the same SN kick after the two core-collapse episodes in order to compare them as close as possible.  \nThe system does not experiences mass transfer before the first supernova, however the primary experience strong stellar wind mass loss during the Wolf-Rayet phase which widens the orbit . Shortly after the primary collapses into a ∼ 14 . 96 M glyph[circledot] BH with a low ( α spin ) 1 . The subsequent detached evolution (Section 8.1) of the mildly eccentric system (due to the BH natal kick), after matching the companion of the BH to a single star grid,  \nleads to a small increase of the period predominantly due to winds. Eventually the secondary star fills its Roche lobe at periastron, where we assume that the system circularizes. Mass transfer onto the BH is interpolated through the binary grid of COs with H-star companions (Section 5.6) and lasts for a few thousand years, becoming unstable and leading to a CE episode (Section 8.2). The system survives the process, forming a tight binary comprised of a BH with a stripped He star on a ∼ 0 . 2 day orbit. Tidal forces become important in this tight orbit, spinning up the He star (Section 5.7), which eventually also forms a mildly spinning ( a spin , 2 ∼ 0 . 48) BH of 7 . 90 M glyph[circledot] . The two BHs merge after 183 Myr from birth, due to gravitational wave radiation.  \nWe have specifically chosen a binary where the differences between the nearest neighbor and initial-final interpolation schemes are relatively small, so that we can accurately display how the binary evolves through each step. Differences between the two evolution options for binaries in general are significantly larger, and the initial-final interpolation scheme is our default choice.  \nThe binary shown in Figure 32 was evolved using our default configuration, although we have purpose-built POSYDON to be modular. Throughout the previous sections we have described possible changes to physical prescriptions that a user can make. However, a user can also easily supplement their own functions for specific steps, or even define the entire binary flow.  \nFinally, for debugging we keep track of any errors or warnings raised throughout a binary's evolution. This allows us to isolate the problematic step for a binary, or a stellar and binary state-event combination that our flow structure cannot handle. This error tracking can be especially useful for user-defined steps and flow structures.", '10. HOW POSYDON EVOLVES A BINARY POPULATION': 'Generating a model binary population for comparison to observations requires two separate steps: initializing individual binaries and then evolving those binaries, which we describe in Section 10.1 and Section 10.2 below. Finally, in Section 10.3 we describe a sample binary population evolved with POSYDON .', '10.1. Initialization': "The primary function of POSYDON is to produce synthetic populations of binaries which requires evolving a random distribution of binaries from ZAMS. We generate an initial population by sampling binary parameters from standard distributions.  \nComponent Masses: For the primary component mass of a binary, we implement initial mass functions (IMF) from Salpeter (1955); Kroupa et al. (1993);  \nFigure 32. Evolution over time, from ZAMS to binary BH formation, of one example binary system with initial properties M 1 = 97 . 07 M glyph[circledot] , M 2 = 39 . 86 M glyph[circledot] in a P orb = 39 . 45 days. The top two rows of panels show the evolution of the binary's P orb and e for our nearest neighbor matching scheme (orange lines) and our initial-final interpolation method (black, circle markers). In the bottom three rows, we show the primary star's (1; solid lines) and secondary star's (2; dashed-dotted lines) properties using the nearest neighbor interpolation method and compare them against the same binary using the initial-final interpolation method (circle markers for primary, star markers for secondary). The binary's evolution is followed across its different evolutionary steps (note that the time scale varies for each step), through both MESA grids (Section 5) and on-the-fly calculations (Section 8).  \n<!-- image -->  \nKroupa (2001). By default we use the Kroupa (2001) IMF in the range M 1 ∈ [7 , 150] M glyph[circledot] with α = 2 . 3. Then the secondary component mass is given by drawing the mass ratio q = M 2 /M 1 from a flat distribution [0 , 1] with M 2 ∈ [0 . 35 , 150] M glyph[circledot] .  \nOrbital Parameters: For a binary population's initial orbital period, Sana et al. (2013) describe a powerlaw in log-space, while Opik's Law describes a log-flat distribution in orbital separation rather than orbital period space (Abt 1983).  \nIn POSYDON we allow for both models to be adopted, with our default being Sana et al. (2013). The minimum orbital period is set by systems that undergo RLO at ZAMS, and a default maximum orbital period to be  \n6000 days. The orbital period distribution from Sana et al. (2013) is undefined for P orb < 1 day; therefore, we extend this distribution so that it is uniform in log 10 P orb space down to a P orb of 0.35 days. We further set the maximum P orb to 10 3 . 5 days. Although binaries may form at wider separations, we choose this upper limit to account for all interacting binary models (e.g., Figure 9).  \nSince we only model circular binaries in our detailed binary-star models with the MESA code, the binary eccentricity e at ZAMS is set to zero (see Section 4 and Section 8.1). Since we plan to generalize this assumption in future work, we also add the option to generate e from a thermal distribution (Duquennoy & Mayor 1991).  \nFigure 33. Example of the double compact object populations produced by POSYDON . Quantities shown are the final distributions, as the binary populations appear today. We separately indicate the parameters for NS-NS, BH-NS, and BH-BH systems. Results from a population synthesis using linear interpolation is indicated with solid lines, while results using the first nearest neighbor approach with dashed lines. Comparison with COSMIC (Fig. 3 from Breivik et al. 2020) shows morphological similarities, but important quantitative differences between the same binary CO populations.  \n<!-- image -->  \nStar-formation History: To account for the starformation history (SFH) of a binary population, we assign to each binary in the population the same maximum age. Then different star-formation rates (SFRs) can be modeled by modifying the distribution of birth times. In POSYDON we offer two options for the SFR: a burst of star formation or a constant SFR. For the burst model, all binaries have identical birth times some number of years prior to the end of the simulation. For a constant SFR, we randomly generate birth times from a uniform distribution within a user-defined range. The average metallicity of the Galaxy and the greater Universe evolve over time, but this is something that we cannot currently model accurately, as our grids of models presented in v1.0 of POSYDON are only calculated for stars at solar metallicity.", '10.2. Evolution': "Evolving a population of ZAMS binaries initialized using the distributions described in Section 10.1 requires implementing the procedure outlined in Section 9 for each binary. In POSYDON we have created an overarching BinaryPopulation class which is a container for a list of individual BinaryStar instances. Each BinaryStar instance is then iteratively evolved until the entire population has been processed.  \nThe BinaryPopulation class contains a number of additional capabilities to efficiently and easily evolve populations of binaries. First, evolved populations are automatically saved in an efficient hdf5 file format with two datasets: one that contains each binary as a single line providing both its initial and final states, and a second dataset that contains the entire evolutionary history of each binary. Second, populations can be evolved either serially or in parallel, so that large ( > 10 6 binaries)  \npopulations can be run quickly on a high-performance computing cluster. Third, we have implemented routines that catch and keep warnings and errors from each binary, so that the code does not crash when a single binary fails to complete its evolution due to a bug. We have found this to be useful for identifying and resolving coding bugs when implementing new physical prescriptions. Last, we have implemented various routines that allow a user to easily select only certain types of binaries (e.g., only BBHs or only double COs).  \nWe envision that a typical user will interact primarily with the BinaryPopulation class and its associated SimulationProperties class, which when combined provide the interface for customizing a particular user's binary population synthesis needs.", '10.3. Example Population': 'To demonstrate the results of a binary population synthesis run with POSYDON , we construct a basic population of 10 6 binaries generated with the default initial conditions described in Section 10.1. We choose a constant star formation history over the past 10 Gyr. The simulation on our high performance computing cluster, Trident , takes approximately 2 hours of wall time using 5 nodes, each with 20 cores, i.e., less than 1 s of CPU time per binary.  \nFrom the resulting binaries, we select those that evolve into bound NS-NS, NS-BH, or BH-BH systems. The present-day properties of these binaries are provided in Figure 33; this can be compared with the results of other studies, e.g., Figure 3 from Breivik et al. (2020) with similar initial conditions at solar metallicity. Overall we find good agreement reaffirming that the code is producing reasonable results. In the left two panels, we show distributions of the component compact object masses.  \nOur NS masses are in very close agreement with those of Breivik et al. (2020), although the BHs we produce extend to larger masses. The reason is that current stellar models with stellar structure and evolution parameters calibrated to the latest observations (e.g., overshoorting, etc) produce cores more massive compared to those in the late-nineties models used in rapid BPS codes. The third panel from the left shows the semi-latus rectum of the population. Again, the distribution results are qualitatively similar, with a preference for NS-NS systems at smaller a (1 -e 2 ), and NN-BH and BH-BH systems having larger a (1 -e 2 ). Finally, the rightmost panel of Figure 33 shows an increasing formation timescale (from ZAMS to the second SN) as we move from BH-BH to NS-BH to NS-NS systems. This is expected since BHs tend to form from more massive systems that complete their evolution more quickly. Our binary CO populations appear morphologically similar to those produced by COSMIC and described in Breivik et al. (2020). We will undertake detailed descriptions of the specific binary populations of interest to separate, science studies.', '11. SUMMARY AND FUTURE WORK': "Here we present POSYDON , a new, next-generation computational tool for general population synthesis of single and binary stars. POSYDON incorporates full stellar structure and evolution sequences for interacting binaries, using the MESA code. Compared to other existing, binary population synthesis code there are significant advances: (i) binary evolution is treated self-consistently without analytical fits of single-star evolutionary tracks and the need for simplified or artificial recipes to emulate the behavior of stars in interacting binaries; (ii) initial-final classification and interpolation methods trained on the pre-calculated grids of binary evolution models, allowing general synthetic simulations of binary populations. The code base along with the existing evolutionary-track grids are publicly available through the POSYDON collaboration's web portal (https://posydon.org) along with full documentation and tutorials for how to use the code. An advanced query system is also available for users to be able to mine the grids of single- and binarystar evolutionary tracks and download relevant data using pre-programmed and customized queries (e.g., Teng et al. 2021a,b). Finally, we provide a user-friendly webapplication that allows a user to perform small-scale simulations with POSYDON online, without the need of code installation and configuration.  \nCompared to current rapid population synthesis codes POSYDON has a smaller set of free parameters, for many of which there are already multiple options for the user to choose from. The code structure is modular and an  \nadvanced user is able to implement their own choices of evolutionary parameters; from as simple as changing the initial properties of the binaries to as complex as incorporating their own custom-made evolutionary-track grids. In this first instrument POSYDON paper we describe in detail the first version of the code, but technical and astrophysical advancements are ongoing and improved code versions will be released in the near and long-term future.  \nOur focus on the technical front is on classification and interpolation methods. Our current process of first classifying the grids and then performing interpolation can lead to errors propagating throughout the pipeline. Instead these two could be combined into a joint treatment to reduce errors, making additional use of covariances between the different grid types (Singh et al. 2016). We will also explore adopting kernelized-interpolation approaches (i.e., Wilson & Nickisch 2015; Gardner et al. 2018; Narayan et al. 2021). To illustrate the use of such an approach, it is observed in Figure 22 that the confusion matrix for our grid of He-rich stars with compact object companions has a lower accuracy for unstable mass transfer. This is perhaps due to the nonlinearly separable decision boundary, as visible in Figure 21. Adopting a kernelized-interpolation technique has the potential to increase the class-specific interpolation accuracy. In particular, kernel-based approaches such as support vector machines and Gaussian processes may prove fruitful. In principle, we may also use a neural networks. Also, improvements can be made by adopting non-Euclidean metrics when defining our distance functions in our k -nearest neighbors classifiers in Section 7.2.  \nApart from methods exploration for the existing classification and interpolation processes, we will focus on the next step of interpolations critical for astrophysical studies: interpolation between whole evolutionary tracks along time. This is a challenging problem that is of particular interest for any study that requires tracking binary properties as a function of their age (e.g., X-ray binary luminosities). In parallel, we are already working on increasing the computational efficiency of building the pre-computed grids necessary for future POSYDON versions. Specifically, for any future grid development, we will take advantage of a new active-learning method developed by our team (Rocha et al. 2022) that allows us to achieve the same classification and interpolation accuracies with a significantly smaller MESA tracks, by dynamically placing them at class boundaries in the parameter space. Such dynamic placement, informed by active learning, leads to grids with non-regular, but smart placement of evolutionary tracks. This work becomes more critical as we expand to more metallicities  \nand add eccentricity as another dimension in order to keep the POSYDON package sizes appropriate for downloads.  \nThe first version of POSYDON is fully functional as an astrophysical tool in the sense that it can be used for complete simulations of binary populations from ZAMS to either formation of a binary with two compact objects or binary distractions (stellar merger or binary disruption), but it is still limited in two ways: our current pre-computed grids are for primaries massive enough to likely form a NS or a BH and are calculated at solar metallicity only. Our next version will expand to a grid of metallicities appropriate for populations across the Universe instead of just the Milky Way. Challenges related to the convergence of single- and binary-star evolutionary models often depend on metallicity, e.g. because the different opacities will results in massive stars reaching the Eddinghton limit within their interiors at different mass ranges and evolutionary phases. Typically, however, such challenges become less severe with decreasing metallicity. Interpolation across metallicities will be a follow-up step as well, once grids at a sufficient number of metallicities have been computed. All the interpolation and classification methods we have developed scale naturally to additional dimensions in the initial conditions parameter space. Further we will expand our grids to low-mass primaries so populations with WDs can also be modeled. Expansion of the grids will inevitably increase data set sizes proportionally. The development of interpolation methods for not only the final properties of each grid but also the entire evolutionary tracks, which is one of our primary future objectives, will allow future POSYDON releases to come only with with the pre-trained interpolation objects. The latter  \nare expected to have a significantly smaller data footprint compared to the downsampled grids we currently ship with POSYDON .  \nWe will also continue to improve the physics treatment of binary evolution. Specifically we are already working on three improvements: (i) in this first version our treatment of the CE phase, once dynamical instability is recognized (taking into account the full stellar structure of the RLO star), is similar to what is done in pBPS codes, apart from the self-consistent calculation of the CE's binding energy. However, since we model binaries with MESA , we are able to treat the phase in a more physical way, either by following the CE inspiral selfconsistently using one-dimensional hydrodynamic simulations (e.g., Fragos et al. 2019) or by following the long-term response of the RLO star to losing its envelope on a very high rate and have its Roche lobe shrinking rapidly (e.g., Marchant et al. 2021; Gallegos-Garcia et al. 2021). One of the objectives of the next version of POSYDON will be to improve the physics of our CE treatment. (ii) like all binary population modeling to date, we assume that binaries circularize instantly upon RLO and for this reason we also assume that all ZAMS binaries are circular. However, the physics of secular binary evolution through mass transfer in eccentric orbits has been fully developed recently (Sepinsky et al. 2007, 2009, 2010; Dosopoulou & Kalogera 2016a,b; Hamers et al. 2021) and will be implemented in future POSYDON versions. (iii) Most recently more physical models for magnetic braking have been developed and calibrated against single-star rotational-velocity data (e.g., Van & Ivanova 2019, 2021; Gossage et al. 2021) and we will use them to update the options for magnetic braking evolution in binaries.", 'APPENDIX': 'Table 5 . The variables in history tables in our grids, taken from the MESA output. The star age is not included since we provide the age variable in the binary history (cf. Table 6) which refers to both the systems, and its components. Since these quantities are undefined for compact objects, history tables are not provided for NSs and BHs; however, the evolution of the mass is given in the binary history tables.  \nglyph[circledot]', 'Table 5 (continued)': 'Table 5 continued  \nTable 5 (continued)  \nTable 6 . The variables in binary history tables, taken from the MESA output.  \nTable 7 . Quantities of final profiles of the stars, taken from the MESA output.Table 7 continued', 'ACKNOWLEDGMENTS': "We thank Corinne Charbonnel, Alex de Koter, Ilya Mandel, Pablo Marchant, Georges Meynet, Fred Rasio, Yorick Vink, Andreas Zezas, for valuable discussion on several aspects stellar- and binary-evolution physics; Aldo Batta, Monica Gallegos-Garcia, Samuel Imperato, Chase Kimball, and Maxime Rambosson for contributing to the code base of the project, and Margaret Lazzarini and Mathieu Renzo for testing early development versions of the code and providing feedback.  \nThe POSYDON project is supported primarily by two sources: a Swiss National Science Foundation Professorship grant (PI Fragos, project number PP00P2 176868) and the Gordon and Betty Moore Foundation (PI Kalogera, grant award GBMF8477).  \nThe collaboration was also supported by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie RISE action, grant agreements No 691164 (ASTROSTAT) and No 873089 (ASTROSTAT-II). Individual team members were supported by additional sources: JJA acknowledges funding  \nglyph[circledot]  \nfrom Northwestern University through a CIERA Postdoctoral Fellowship, CPLB acknowledges support by the CIERA Board of Visitors Research Professorship, and SC through CIERA as a Computational Specialist. VK was partially supported through a CIFAR Senior Fellowship and a Guggenheim Fellowship. KK and EZ were partially supported by the Federal Commission for Scholarships for Foreign Students for the Swiss Government Excellence Scholarship (ESKAS No. 2021.0277 and ESKAS No. 2019.0091, respectively). YQ acknowledges funding from the Swiss National Science Foundation (grant P2GEP2 188242). DM and KR thank the LSSTC Data Science Fellowship Program, which is funded by LSSTCorporation, NSF Cybertraining Grant No. 1829740, the Brinson Foundation, and the Gordon and Betty Moore Foundation; their participation in the program has benefited this work. ZX was supported by the Chinese Scholarship Council (CSC). MZ was supported as an IDEAS Fellow, through the NRT IDEAS program, a research traineeship program supported by the National Science Foundation (PI Kalogera, award DGE-1450006).  \nThe computations were performed at Northwestern University on the Trident computer cluster (funded by the GBMF8477 award) and at the University of Geneva on the Baobab and Yggdrasil computer clusters. This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Of-  \nfice of the Provost, the Office for Research, and Northwestern University Information Technology.  \nSoftware: 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).", 'REFERENCES': "Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020,  \nNature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2  \nWebbink, R. F. 1984, ApJ, 277, 355, doi: 10.1086/161701  \nWiktorowicz, G., Lasota, J.-P., Middleton, M., &  \nBelczynski, K. 2019, ApJ, 875, 53,  \ndoi: 10.3847/1538-4357/ab0f27  \nWillems, B., Henninger, M., Levin, T., et al. 2005, ApJ, 625, 324, doi: 10.1086/429557  \nWilson, A. G., & Nickisch, H. 2015, in Proceedings of the 32nd International Conference on International (JMLR.org), 1775-1784.  \nConference on Machine Learning - Volume 37, ICML'15 https://arxiv.org/abs/1503.01057  \nWong, T.-W., Valsecchi, F., Fragos, T., & Kalogera, V. 2012, ApJ, 747, 111, doi: 10.1088/0004-637X/747/2/111 Woosley, S. E. 2017, ApJ, 836, 244,  \ndoi: 10.3847/1538-4357/836/2/244  \n- -. 2019, ApJ, 878, 49, doi: 10.3847/1538-4357/ab1b41  \nWoosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390, doi: 10.1038/nature06333  \n- Woosley, S. E., & Bloom, J. S. 2006, ARA&A, 44, 507,  \ndoi: 10.1146/annurev.astro.43.072103.150558 Woosley, S. E., & Heger, A. 2015, Astrophysics and Space Science Library, Vol. 412, The Deaths of Very Massive Stars, ed. J. S. Vink, 199, doi: 10.1007/978-3-319-09596-7 7 Yoon, S.-C. 2017, MNRAS, 470, 3970, doi: 10.1093/mnras/stx1496 Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS, 457, 351, doi: 10.1093/mnras/stv3002 Zahn, J. P. 1977, A&A, 500, 121 Zapartas, E., de Mink, S. E., Justham, S., et al. 2021, A&A, 645, A6, doi: 10.1051/0004-6361/202037744 Zapartas, E., de Mink, S. E., Izzard, R. G., et al. 2017, A&A, 601, A29, doi: 10.1051/0004-6361/201629685 Zapartas, E., de Mink, S. E., Justham, S., et al. 2019, A&A, 631, A5, doi: 10.1051/0004-6361/201935854 Zevin, M., Spera, M., Berry, C. P. L., & Kalogera, V. 2020, ApJL, 899, L1, doi: 10.3847/2041-8213/aba74e Zuo, Z.-Y., & Li, X.-D. 2014, MNRAS, 442, 1980, doi: 10.1093/mnras/stu993 Zuo, Z.-Y., Li, X.-D., & Gu, Q.-S. 2014, MNRAS, 437, 1187, doi: 10.1093/mnras/stt1918"}