Issue 
A&A
Volume 636, April 2020



Article Number  A104  
Number of page(s)  40  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201936528  
Published online  24 April 2020 
Evolutionary roads leading to low effective spins, high black hole masses, and O1/O2 rates for LIGO/Virgo binary black holes^{⋆}
^{1}
Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00716 Warsaw, Poland
email: chrisbelczynski@gmail.com
^{2}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{3}
Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA
^{4}
Astronomical Observatory, Warsaw University, Al. Ujazdowskie 4, 00478 Warsaw, Poland
^{5}
Center for Theoretical Physics, Polish Academy of Sciences, Al. Lotnikow 32/46, 02668 Warsaw, Poland
^{6}
Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
^{7}
Department of Astronomy, University of Geneva, Chemin des Maillettes 51, 1290 Versoix, Switzerland
^{8}
CCS2, MSD409, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
^{9}
University of Chicago, Chicago, IL 60637, USA
^{10}
Center for Computational Relativity and Gravitation, Rochester Institute of Technology, Rochester, NY 14623, USA
^{11}
Department of Physics, Syracuse University, Syracuse, NY 13224, USA
^{12}
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo, Kashiwa, Chiba 2778583, Japan
^{13}
TAPIR, Walter Burke Institute for Theoretical Physics, Mailcode 35017, Caltech, Pasadena, CA 91125, USA
^{14}
Department of Astronomy & Astrophysics, University of California, 1156 High Street, Santa Cruz, CA 95064, USA
^{15}
LennardJones Laboratories, Keele University, Keele ST5 5BG, UK
^{16}
School of Physics & Astronomy & Institute for Gravitational Wave Astronomy, University of Birmingham, Birmingham, UK
^{17}
Lund Observatory, Department of Astronomy, and Theoretical Physics, Lund University, Box 43, 221 00 Lund, Sweden
^{18}
Physics Department, Kenyon College, 201 North College RD, Gambier, OH 43022, USA
^{19}
Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA
^{20}
National Astronomical Observatories & University of the Chinese Academy of Sciences, Beijing, PR China
^{21}
Department of Astronomy and Joint SpaceScience Institute University of Maryland, College Park, MD 207422421, USA
^{22}
Department of Physics, University of Oregon, Eugene, OR 97403, USA
^{23}
Institut d’Astrophysique de Paris, CNRS et Sorbonne Université, UMR 7095, 98bis Bd Arago, 75014 Paris, France
Received:
20
August
2019
Accepted:
5
March
2020
All ten LIGO/Virgo binary black hole (BHBH) coalescences reported following the O1/O2 runs have nearzero effective spins. There are only three potential explanations for this. If the BH spin magnitudes are large, then: (i) either both BH spin vectors must be nearly in the orbital plane or (ii) the spin angular momenta of the BHs must be oppositely directed and similar in magnitude. Then there is also the possibility that (iii) the BH spin magnitudes are small. We consider the third hypothesis within the framework of the classical isolated binary evolution scenario of the BHBH merger formation. We test three models of angular momentum transport in massive stars: a mildly efficient transport by meridional currents (as employed in the Geneva code), an efficient transport by the TaylerSpruit magnetic dynamo (as implemented in the MESA code), and a veryefficient transport (as proposed by Fuller et al.) to calculate natal BH spins. We allow for binary evolution to increase the BH spins through accretion and account for the potential spinup of stars through tidal interactions. Additionally, we update the calculations of the stellarorigin BH masses, including revisions to the history of star formation and to the chemical evolution across cosmic time. We find that we can simultaneously match the observed BHBH merger rate density and BH masses and BHBH effective spins. Models with efficient angular momentum transport are favored. The updated stellarmass weighted gasphase metallicity evolution now used in our models appears to be key for obtaining an improved reproduction of the LIGO/Virgo merger rate estimate. Mass losses during the pairinstability pulsation supernova phase are likely to be overestimated if the merger GW170729 hosts a BH more massive than 50 M_{⊙}. We also estimate rates of black holeneutron star (BHNS) mergers from recent LIGO/Virgo observations. If, in fact. angular momentum transport in massive stars is efficient, then any (electromagnetic or gravitational wave) observation of a rapidly spinning BH would indicate either a very effective tidal spin up of the progenitor star (homogeneous evolution, highmass Xray binary formation through case A mass transfer, or a spin up of a WolfRayet star in a close binary by a close companion), significant mass accretion by the hole, or a BH formation through the merger of two or more BHs (in a dense stellar cluster).
Key words: stars: massive / black hole physics / gravitational waves
Our updated models of BHBH, BHNS and NSNS mergers are now publicly available at www.syntheticuniverse.org under the tab “Download/2020: Double Compact Objects/Belczynski et al. 2020”
© ESO 2020
1. Introduction
The LIGO/Virgo Collaboration has reported the detection of ten binary black hole (BHBH) mergers during the O1/O2 observations: GW150914, GW151012, GW151226, GW170104, GW170608, GW170729, GW170809, GW170814, GW170818, and GW170823 (Abbott et al. 2019a,b). We list the basic properties of these mergers in Table 1.
LIGO/Virgo O1/O2 BHBH mergers.
For our analysis, we only used the LIGO/Virgo–reported events; however, three additional BHBH mergers – GW170121, GW170304, and GW170727 – have been detected with a high level of confidence (P_{astro} > 0.98^{1}) by Venumadhav et al. (2019). They have low effective spins that are consistent with the triggers reported by the LIGO/Virgo collaboration. The three other events reported by Venumadhav et al. (2019) are not likely to be of astrophysical origin. In particular, GW170403, with a highly negative effective spin, has been reported as the least secure event. Two more events were reported with high effective spins: one with (Venumadhav et al. 2019) and the second with (Zackay et al. 2019); however, these triggers have only a 56% and 71% chance of being of astrophysical origin, respectively.
The majority of the reported mergers contain “heavy” BHs with component masses > 20 M_{⊙} and they are consistent with being formed by isolated binary evolution of stars with metallicities ≲10% Z_{⊙} and initial masses in the range of 40−100 M_{⊙}, while lower mass BHs may have formed from lower mass stars or at higher metallicity (Belczynski et al. 2010a, 2016a). A typical evolution involves the interaction of stars in a binary through mass transfer and a common envelope phase (Tutukov & Yungelson 1993; Lipunov et al. 1997; Belczynski et al. 2002; Voss & Tauris 2003; Dominik et al. 2012; Belczynski et al. 2016b; Eldridge & Stanway 2016; Woosley 2016; Stevenson et al. 2017; Kruckow et al. 2018; Hainich et al. 2018; Marchant et al. 2019; Spera et al. 2019) and its outcome was predicted to produce the first LIGO/Virgo sources (Belczynski et al. 2010b).
The dynamical formation scenario of BHBH mergers is an evolutionary channel alternative to isolated binary evolution that could operate in globular clusters, nuclear clusters or open, young clusters (Miller & Hamilton 2002a,b; Portegies Zwart et al. 2004; Gültekin et al. 2004, 2006; O’Leary et al. 2007; Sadowski et al. 2008; Downing et al. 2010; Benacquista & Downing 2013; Bae et al. 2014; Chatterjee et al. 2017; Mapelli 2016; Hurley et al. 2016; Rodriguez et al. 2016a; VanLandingham et al. 2016; Askar et al. 2017; ArcaSedda & CapuzzoDolcetta 2019; Samsing 2018; Morawski et al. 2018; Banerjee 2018; Di Carlo et al. 2019; Zevin et al. 2019; Rodriguez et al. 2018; Perna et al. 2019). The dynamical channel can also explain the range of BH masses observed by LIGO/Virgo.
However, these two basic scenarios, dynamical and isolated, predict different spinorbit misalignment distributions, with dynamical formation generating nearly isotropic distributions (Mandel & O’Shaughnessy 2010; Rodriguez et al. 2016b), while binary evolution favors distributions that show mostly small to moderate misalignments, with only a small fraction of mergers reaching high misalignments (Wysocki et al. 2018; Gerosa et al. 2018). A hybrid BHBH merger formation channel that involves isolated triple (dynamically interacting) star systems (Antonini et al. 2017) favors BH spin vectors that are located in the BHBH orbital plane, with ∼90° misalignment with respect to the orbital angular momentum if the tertiary dominates the angular momentum of the system (Antonini et al. 2018).
The spin orientations of the BHs can therefore provide important information about their formation (Farr et al. 2017, 2018; Vitale et al. 2017a). However, the effect of spin is subdominant in gravitational waveforms, so spins are more difficult to measure than masses. The waveform is most sensitive to the binary’s effective spin,
with a_{spin1, 2} being the BH spin magnitudes and Θ_{1, 2} being the angles between the BH spins and the orbital angular momentum. The dimensionless BH spin magnitude is defined as:
where c is the speed of light, G is the gravitational constant, and M_{BH} and J_{BH} are, respectively, the mass and angular momentum of the BH.
All ten of the BHBH binaries observed by LIGO/Virgo are consistent with effective spins near zero. Indeed, within the 90% confidence levels reported by LIGO/Virgo all mergers are consistent with −0.1 ≲ χ_{eff} ≲ 0.1 (see Table 1). If the BH spin magnitudes are large, then either (i) both black hole spin vectors are close to the orbital plane, so that cosΘ_{1, 2} ≈ 0, or (ii) the black hole angular momenta are nearly equal in magnitude but close to oppositely directed. Alternatively, small values of the effective spin parameter may be obtained if (iii) BH spin magnitudes are small.
In this study, we investigate the third possibility, namely, the approach that allows the spin magnitudes of BHs in BHBH mergers detected by LIGO/Virgo to be small. Our study is limited to the classical isolated binary evolution BHBH formation scenario. We test several models of natal BH spins to predict the effective spin parameters of BHBH mergers in the local Universe and to compare them with LIGO/Virgo observations. Compared to Belczynski et al. (2016a,b,c), we incorporated updated mass loss by pairinstability pulsation supernovae, revised our model of accretion onto BHs in close binaries, allowed for effective tidal spinup of WolfRayet stars, and adopted a new model for the evolution of metallicity and the star formation rate across cosmic time.
Here we briefly summarize the most important ingredients and results of our study. Based on single stellar evolution calculations, we introduce three BH natal spin models (see Sect. 2.1) and we argue that these models can be reasonably used in binary evolution as implemented in our calculations (see Sect. 4.2). We show that initial star rotation does not play a significant role, while angular momentum transport plays crucial role in setting the BH natal spin (see Sect. 4.1). In binary evolution, accretion and mass transfer (see Sect. 2.4) do not play significant role, while tidal interactions (see Sect. 2.5) may increase BH spins for about 20−30% of the BHBH mergers (see Sect. 3.4). Merger rates for BHBH, neutron starneutron star (NSNS), and black holeneutron star (BHNS) binaries are sensitive to assumptions regarding the common envelope and natal kicks, but they also strongly depend on the assumed cosmic chemical evolution model, while the change of merger rates with redshift depends mostly on the cosmic starformation history (see Sect. 3.2).
2. Model
2.1. Natal black hole spin
Various rotating star models differ in terms of the physics of their rotation and, in particular, the efficiency of angular momentum transport. Here we test three different models of angular momentum transport for massive stars. We adopted moderate angular momentum transport, along with Geneva stellar models (Eggenberger et al. 2008; Ekström et al. 2012) that are based on Zahn (1992) theory, in which angular momentum is mainly transported by meridional currents: the shellular model. Efficient angular momentum transport was adopted with standard MESA stellar models (calculated in this study), which use an efficient transport mechanism mediated by the action of the socalled TaylerSpruit magnetic dynamo in the radiative zone (Spruit 1999, 2002). Superefficient angular momentum transport is based on the TaylerSpruit dynamo, which is subsequently modified to include stronger magnetic field generation that allows for more efficient angular momentum transport (Fuller et al. 2019; Ma & Fuller 2019; Fuller & Ma 2019).
To test the validity of each of these models through gravitational wave astrophysics, we compared their predictions with the LIGO/Virgo effective–spin estimates. For the Geneva (moderate) and MESA (efficient) angular momentum transport, we use models based on a wide range of metallicity: Z = 0.014; 0.006; 0.002; 0.0004 (see Appendices A.1 and A.2). We assume that stars on the zero age main sequence (ZAMS) have an equatorial velocity equal to 40% of the critical velocity, defined as the velocity at which the centrifugal acceleration completely balances gravity^{2}. To test the importance of the initial rotation rates, we ran several models based on different assumptions about this initial condition; for more, see Sect. 4.1 and Figs. 23 and A.1. We find that the final core rotation rate is almost independent of the initial rotation rate but depends strongly on the angular momentum transport process included in the simulation (i.e., nonmagnetic or magnetic). For the Fuller (superefficient) model of angular momentum transport, we assume that the natal BH spin depends neither on initial stellar rotation nor on metallicity, since in each case almost all angular momentum is removed from the stellar core (Fuller et al. 2019; Ma & Fuller 2019; Fuller & Ma 2019).
The stellar models provide the angular momentum in each zone that corresponds to a spherical shell in the star. We assume that angular momentum is conserved in the collapse phase and we calculate the angular momentum of the compact remnant by summing the angular momentum of the zones with enclosed mass lower than the compact remnant mass.
From the amount of angular momentum contained in the collapsing core, we calculate the dimensionless BH spin magnitude a_{spin}. However, we limit the spin magnitude to 0.9: a_{spin} = min(a_{spin}, 0.9) to account for any potential processes that might remove some angular momentum during BH formation. We arbitrarily picked 0.9 as a maximum for the spin rate. The difficulty in predicting the fastest spin values is that they require highangular momentum, diskforming infalls. There are several spurious mechanisms which can remove angular momentum from these accretion disks, but for now we do not have means to perform predictive calculations of this effect. Disk winds are one example, yet the actual amount of material and angular momentum that is removed by such winds can vary significantly (Vlahakis & Königl 2001; Surman et al. 2006; Janiuk 2017). Therefore, we have used a_{spin} = 0.9 to denote any value between 0.9 and 1.0, depending on these accretion loss mechanisms.
In Figs. 1 and 2, we show the BH spins as a function of the progenitor’s CO core mass for the Geneva and MESA models. We approximate the natal BH spin by simple fits. These rough fits are meant to reproduce the general trends in the data and are used in our population synthesis calculations. We list the actual data in Tables A.1 and A.2, so other fits can be attempted if desired.
Fig. 1.
Magnitude of natal BH spin as a function of the CO core mass of the collapsing star for the Geneva stellar models with 40% critical initial velocity and mild angular momentum transport by meridional currents. Lines mark our adopted model for natal BH spins and its dependence on metallicity (see Eq. (3)). The star’s CO core mass may be used as a proxy for the BH mass (see Appendix A.3). 
Fig. 2.
Magnitude of natal BH spin as a function of the CO core mass of the collapsing star for the MESA stellar models with 40% critical initial velocity and the TaylerSpruit magnetic dynamo (efficient) angular momentum transport. Lines mark our adopted model for natal BH spins and its dependence on metallicity (see Eq. (4)). 
For the Geneva models the natal BH spin may be approximated by:
with a = −0.088 for all models; b = 2.258, m_{1} = 16.0 M_{⊙}, m_{2} = 24.2 M_{⊙}, a_{low} = 0.13 for Z = 0.014; b = 3.578, m_{1} = 31.0 M_{⊙}, m_{2} = 37.8 M_{⊙}, a_{low} = 0.25 for Z = 0.006; b = 2.434, m_{1} = 18.0 M_{⊙}, m_{2} = 27.7 M_{⊙}, a_{low} = 0.0 for Z = 0.002; and b = 3.666, m_{1} = 32.0 M_{⊙}, m_{2} = 38.8 M_{⊙}, a_{low} = 0.25 for Z = 0.0004. We note that progenitor stars with CO cores less massive than ∼20 M_{⊙} (low to intermediatemass BHs) tend to produce highspin BHs (a_{spin} ∼ 0.8−0.9), but higher mass stars (massive BHs) tend to produce lowspin BHs (a_{spin} ∼ 0−0.3). This general trend is easily understood. Stellar winds during the evolution of a massive star can carry away a considerable amount of angular momentum (Meynet et al. 2015). For the most massive stars, this mass loss is extensive, effectively removing angular momentum and producing lowspin BHs. The data points also show a nonmonotonic dependence on metallicity, which is the result of a complex and metallicitydependent interplay between the strength of stellar winds, the extent of the Hburning shell, and the model for the efficiency of element diffusion within the meridional current (see Appendix A.1).
For MESA models the natal BH spin may be approximated by:
with a_{1} = −0.0016, b_{1} = 0.115, m_{1} = ∞ for Z = 0.014; a_{1} = −0.0006, b_{1} = 0.105, m_{1} = ∞ for Z = 0.006; a_{1} = 0.0076, b_{1} = 0.050, a_{2} = −0.0019, b_{2} = 0.165, m_{1} = 12.09 M_{⊙} for Z = 0.002; and a_{1} = −0.0010, b_{1} = 0.125, m_{1} = ∞ for Z = 0.0004. The MESA models include the TaylerSpruit magnetic dynamo and thus models of all masses and at all metallicities end up with low BH spin magnitudes in the range of 0.05 ≲ a ≲ 0.15. There is a mild tendency for lower metallicity and lower initial mass models to end up with slightly higher BH spin magnitudes but the dependence is much weaker than in the Geneva models, which do not include magnetic field related transport of angular momentum.
Finally, for the Fuller model that employs superefficient angular momentum transport, we adopt a single value of BH natal spin for stars of all masses and all metallicities:
We note that single stellar models (of different mass and metallicity) presented by Fuller & Ma (2019) produce very low spins for BHs in the range of a = 0.003−0.035, with a typical value of a ∼ 0.01.
In our binary population synthesis calculations, we employ the abovepresented natal BH spin estimates obtained from single stellar evolution (see Sect. 2.7). The discussion of binary interactions that can affect these natal BH spins (that are and are not taken into account) is given in Sect. 4.2.
2.2. Black hole spin misalignment
To calculate the effective spin of BHBH mergers we need to know the misalignment angles of the two BH spin vectors with respect to the orbital angular momentum vector: Θ_{1} and Θ_{2} (see Eq. (1)).
In our estimation of these angles, we ignore the potential effects of tides and mass transfer and accretion; see Gerosa et al. (2013) and Wysocki et al. (2018) for alternatives and further discussion. However, in our calculations, we take into account the binary components spin precession. We assume that the two stars are born with spins that are fully aligned with the ZAMS binary orbital angular momentum (L_{0}). At each BH formation, the natal kick may change the orbit and its orientation in space. After the first BH formation, the new orbital angular momentum is L_{1}, whereas following the second BH formation it is L_{2}. After the first BH formation, if there is a natal kick, the binary component spins S_{1} and S_{2} are not aligned with the binary angular momentum vector. We assume that after the first BH formation both the BH spin and the companion star spin are small compared to the binary orbital angular momentum. With this assumption, the total angular momentum of the binary is equal to the orbital angular momentum. Therefore, both the binary component spins precess around L_{1}. The precession periods for the BH and the companion star spin are different. In this approximation the spin of the BH (or the star) can be described as a sum of two components , where is the spin component parallel to the orbital angular momentum and is the spin component perpendicular to the orbital angular momentum. During precession the former is constant and the latter has a constant value but rotates in the plane perpendicular to the direction of the angular orbital momentum. The precession period of each spin is different (Hamilton & Sarazin 1982). Thus, at the time of the second BH formation we chose random positions of each component vector on its precessing trajectory and they become and . The natal kick (if any) changes again the binary orbit orientation in space, which then becomes L_{2}. We followed all these changes to calculate the effective spin:
We note that although (after the second BH formation) misaligned BH spins are subject to precession, the value of χ_{eff} is expected to remain constant during the long inspiral towards the LIGO/Virgo band (Gerosa et al. 2015).
2.3. Black hole masses
In our calculation of BH masses, we use formulas assuming both the rapid and delayed supernova engines (Fryer et al. 2012). The rapid development of the engine naturally creates a mass gap between NSs and BHs (dearth of compact objects in the range of 2−5 M_{⊙}), while the delayed engine does produce compact objects in this range (Belczynski et al. 2012a). For singlestar evolution, the delayed model minimum BH mass is 2.5 M_{⊙}, while it is 5 M_{⊙} for the rapid model. We note that binary evolution may create light BHs (∼2.5−5 M_{⊙}) by accretion induced collapse of a NS to a BH independent of supernova model (Belczynski & Taam 2004).
The maximum mass of a BH is set by the maximum mass of a star and wind mass loss rates and is sensitive to the potential explosive massloss during the final stages of the star’s life. We adopted a rather conservative maximum initial mass for the stars: M_{ZAMS} < 150 M_{⊙}, although there appears to be evidence that stars with masses of ∼200−300 M_{⊙} may exist even in the local Universe (Crowther et al. 2010). For stellar winds we adopt formulas based on theoretical predictions of radiation driven mass loss (Vink et al. 2001) with inclusion of Luminous Blue Variable mass loss (Belczynski et al. 2010a). These wind mass loss rates may be overestimated by as much as one order of magnitude (Oskinova et al. 2011; Ramachandran et al. 2019). Therefore, we allow (conservatively) for the reduction of stellar wind mass loss rates to f_{wind} = 0.3 of their currently adopted values.
We also allow for pairinstability supernovas (PSNs) to entirely disrupt massive stars (stars with He core mass in the range of M_{He} = 65−135 M_{⊙}) and leaving no NS/BH remnant. For somewhat lowermass stars (stars with He core mass in the range of M_{He} ≈ 40−65 M_{⊙}), we allow for pairinstability pulsation supernovas (PPSNs). This process may remove the outer layers of a star, but does not lead to its disruption and allows for the BH formation in core collapse. We adopt three models for PPSN mass loss. In the first model, we adopt strong PPSNs that are assumed to always remove the entire star mass above the inner 45.0 M_{⊙} (Belczynski et al. 2016c) for stars with M_{He} = 45−65 M_{⊙}. Therefore, the post PPSN star mass is:
In the second model, we adopt recent PPSN calculations that allow for as much as 51.2 M_{⊙} of the star to remain bound after a PPSN (Leung et al. 2019). In this moderate PPSN model, the explosive mass loss depends on the He core mass. The post PPSN star mass is:
Finally, we adopt a weak PPSN model that allows only for 50% of the mass loss calculated by Leung et al. (2019). In this scheme, the post PPSN star mass may reach 55.6 M_{⊙}, and we approximate the post PPSN star mass with:
All three models are shown in Fig. 3. For the two models based on calculations from Leung et al. (2019), we note a steep decrease in post PPSN star mass for He core masses above ∼60 M_{⊙}. The masses of these He stars are very close to the boundary mass between the PPSN and the PSN. As the mass of the He star increases, the central temperature at the bounce increases and oxygen burning becomes more explosive and causes a larger amount of mass ejection. Eventually, explosive oxygen burning in the ∼65 M_{⊙} He star produces large enough nuclear energy to disrupt the star completely with no BH remnant, in other words, it induces PSN. The amount of mass ejection increases steeply as the He star mass increases from ∼60 M_{⊙} to ∼65 M_{⊙} because the oxygen burning rate is very sensitive to the temperature. Thus the remnant BH mass decreases steeply as the He star mass approaches 65 M_{⊙}.
Fig. 3.
Adopted models for pairinstability pulsation supernova mass loss. For a given He core mass we show the mass of a star after PPSN mass loss. Moderate PPSN mass loss is adopted directly from Leung et al. (2019), while its modified (50% reduced mass loss) version is presented as weak PPSN model. Strong PPSN are adopted from Belczynski et al. (2016c). 
In the PPSN mass regime, massive BHs are formed and, according to our scheme (Fryer et al. 2012), these BHs form through collapse of the entire star into a BH (direct BH formation). However, we allow for some mass loss in neutrinos. In the original formulas we have adopted a high neutrino fractional mass loss (f_{neu} = 0.1; see e.g., Belczynski et al. 2016a), whereas, here, we also allow also for much smaller mass loss in some models (f_{neu} = 0.01). The final BH mass, in case of the direct BH formation, is calculated from:
where M_{star} denotes the star’s mass just prior the core collapse.
2.4. Black hole accretion model
Mass accretion may increase the BH mass and spin after its formation. Here we test two models of accretion from a stable Roche lobe overflow (RLOF; not common envelope [CE]) masstransfer or from stellar winds. The first model (efficient BH accretion) is based on the results of global, axisymmetric simulations of accretion disks with αP viscosity, disk winds and photon trapping performed by Ohsuga (2007). Belczynski et al. (2008a) obtained a fit to the results of Ohsuga (2007):
where Ṁ_{BH} is the mass accumulation rate onto the BH, Ṁ_{crit 1} = 2.6 × 10^{−8}(M_{BH}/10M_{⊙})M_{⊙} yr^{−1} is the critical mass transfer rate (obtained from numerical simulations) above which the accretion onto the BH is not fully efficient, and Ṁ_{don} is the mass transfer rate from the donor to the accretion disk around the BH. We note that above the critical mass transfer rate some mass transferred from the donor is ejected from the system.
The second model (inefficient BH accretion) uses the analytical prescription of Shakura & Sunyaev (1973). In this model, disk winds play a significant role effectively limiting the mass accumulation for superEddington mass transfer rates by limiting the local accretion rate to its Eddington value. This model is supported by the numerical results of strong outflow from the supercritical accretion disk (Abolmasov et al. 2009; Sądowski et al. 2016), as well as by observations of many ultraluminous Xray sources in centers of unusually large and bright emission nebulae which are powered by the accumulated kinetic energy of the outflow (Pakull & Mirioni 2003; King 2004), and also supported by the fact that the neutron star in Cygnus X2 while being subject to high mass transfer rate (∼10^{−5}M_{⊙} yr^{−1}) has ejected most of the mass transferred from the donor star (King & Ritter 1999).
In this model we assume that the mass accumulation rate onto the BH is given by:
where (1 − f_{1}) is the fraction of the mass transferred to the disk from the donor which is lost in wind from the outer part of the accretion disk. It is difficult to determine the wind mass loss rate from the outer part of the disk (see Sądowski et al. 2016), and we treat f_{1} as a parameter in our model. For the current calculations we adopt f_{1} = 1, indicating that there is no wind mass loss in the outer part of the disk, which maximizes accumulation of mass onto the BH. Similarly, (1 − f_{2}) represents the fraction of mass which is lost from the inner part of the accretion disk. The value of f_{2} depends on the mass transfer rate as:
where the critical mass transfer rate is M_{⊙} yr^{−1}, where the hydrogen mass fraction X is 0.7 for Hrich donor stars and 0.0 for Hdeficient donor stars, R_{ISCO} is the innermost stable circular orbit radius, and R_{sph} is the spherisation radius, where the disk’s height becomes comparable to the radius. The R_{ISCO} depends on the spin value of the accreting BH and varies between 0.5R_{S} (for a maximally prograde spinning BH) to 4.5R_{S} (for a maximally retrograde spinning BH), where is the Schwarzchild radius of a BH. The R_{sph} is given by (Shakura & Sunyaev 1973) as:
The two equations above make up a simplified description of the change of the accretion mode onto compact objects. In particular, they lead to a jump in mass accretion rate onto a compact object when mass transfer from a donor star is equal to the critical mass transfer rate (M_{don} = Ṁ_{crit2}). This simplification does not influence the results of our evolutionary calculations.
In our calculations, we always assume a prograde BH spin, and the initial BH spin magnitude is adopted from a given stellar model (see Sect. 2.1) and then increased by accretion as detailed in Belczynski et al. (2008a). At subEddington accretion rates (Ṁ_{don} ≤ Ṁ_{Edd}), we assume that there is no mass loss from the inner part of the disk (f_{2} = 1).
The mass accumulation onto a BH in the inefficient BH accretion model is always Eddingtonlimited (Ṁ_{BH} < Ṁ_{Edd}), whereas in the efficient BH accretion model it increases monotonically with the mass transfer rate and may significantly exceed Ṁ_{Edd}. These two models produce the same accumulation on a BH for mass transfer rates below Ṁ_{crit1} = 0.06 × Ṁ_{Edd} where both models give exactly the same prescription (Ṁ_{BH} = Ṁ_{don}, noting that Ṁ_{crit1} = 0.06 Ṁ_{crit2}).
As a result, the evolution of BH binaries which go through a phase of stable mass transfer with high (supercritical) mass transfer rate may be considerably different under different assumptions about the BH accumulation efficiency. However, the dominant formation channel for BHBH mergers (also Figs. 8 and 9 in Belczynski et al. 2016a) contains no such phase. Along the way to formation for most BHBH mergers in our models, the accretion of matter by a BH takes place only during a shortlived CE event or as result of capturing a fraction of a stellar wind from the companion. Accretion during CE can be best related to the BondiHoyle accretion and has been found to be rather insignificant (less than ∼1 M_{⊙} for a typical 30 M_{⊙} BH, see Appendix A.7 for a detailed discussion of recent calculations of accretion during CE). The windfed accretion during the subsequent phase of a compact BH – WolfRayet (BHWR) binary evolution is even less significant, partially due to small windmass loss rates from low metallicity systems (which are the progenitors of most BHBH mergers).
That said, the BH accretion models presented in this section can play an important role in certain subdominant channels for the BHBH merger formation. In Fig. 4 we showcase the time evolution of the mass transfer and accretion rates as well as the BH spin and mass in a system in which the BH accretion is particularly significant. The system began its evolution as two ZAMS stars with masses of 84.6 M_{⊙} and 48 M_{⊙} formed at a very low metallicity Z = 0.0002 in a binary with separation of about 1900 R_{⊙}. At the age of 3.8 Myr the primary, now a 66 M_{⊙} heliumcore burning star with a radius of 800 R_{⊙}, goes into RLOF and initiates a CE phase. As a result, the primary becomes a 36 M_{⊙} WR star and the binary separation decreases down to 30 R_{⊙}. With the orbit already being quite compact, the companion MS star (∼47 M_{⊙}) initiates another RLOF and starts stable mass transfer back onto the WR primary. The WR star grows to about 39 M_{⊙}, before collapsing directly into a 35.54 M_{⊙} BH at the age of 4 Myr (with no natal kick, 10% of mass being lost in neutrinos, and the natal spin of a_{init} = 0.832 adopted from Geneva BH natal spin model). Soon thereafter, the companion MS initiates a RLOF again and the first phase of stable, supercritical mass transfer onto the newlyformed BH begins (Ṁ_{don} ∼ a few × 10^{−5}M_{⊙} yr^{−1}). At that moment the system may be potentially observable as an ultraluminous Xray source (see e.g., Wiktorowicz et al. 2019, for the recent analysis of these objects). The mass transfer continues up until about t = 4.8 Myr (Fig. 4), at which point the donor is at the very end of its MS evolution and contracts a bit, causing a temporary detachment. The mass transfer starts again when the companion begins to expand on its HG at the age of about t = 5.9 Myr having mass of 28 M_{⊙}. The rate is high (∼10^{−3}M_{⊙} yr^{−1}), but the phase is shortlived (7.7 Myr). Finally, the last phase of mass transfer occurs when the companion is a slowly expanding core helium burning star and terminates at about 5.25 Myr.
Fig. 4.
Comparison of the two black hole accretion models employed in the StarTrack code. We present the evolution of the mass transfer rate from a donor star (Ṁ_{don}; top panel, solid lines), BH mass accumulation rate (Ṁ_{BH}; top panel, dashed lines), BH mass (M_{BH}; middle panel), and BH spin (a_{spin}; bottom panel) during the RLOF stable mass transfer phases. Critical mass transfer rates (above which mass ejection from a system is expected) are provided for reference. During the first part of the RLOF (t ≈ 4−4.8 Myr) the donor is a MS star, during the shortduration peak it is a HG star (t ≈ 4.9 Myr), and for the remaining RLOF it is a heliumcore burning star. See Sect. 2.4 for a description of the full evolutionary path and both accretion models. 
The net result of the subsequent stable mass transfer phases shown in Fig. 4 is an increase of the binary separation which, at the point of the final detachment, is about 160 R_{⊙} (similar in both accretion models). At that point the donor star has been stripped almost down to its helium core and has a remaining mass of only about 10 M_{⊙}.
Since accumulation of mass onto the BH (dashed lines in the top panel of Fig. 4) is different in the two models for supercritical mass transfer rates, the final BH mass is also different: 42.08 M_{⊙} for the efficient BH accretion model and 35.66 M_{⊙} for the inefficient BH accretion model (Fig. 4, middle panel). We also note that the BH spin, was increased to 0.948 for the efficient model and to 0.835 for the inefficient model.
This particular binary ends its evolution at t ≈ 5.6 Myr forming a BHBH with masses of 35.7 M_{⊙} and 6.2 M_{⊙} for the inefficient model, or 42.1 M_{⊙} and 6.0 M_{⊙} for the efficient model. The delay time to merger is about 3 Myr. The second BH is formed through a supernova explosion with a natal kick. As the separation before the supernova explosion was relatively large in both accretion models (∼160 R_{⊙}), the formation of a binary that could merge within the Hubble time was only possible thanks to the preferentially oriented natal kick, which decreased the separation down to ∼90 R_{⊙} and, more importantly, induced a high eccentricity of e > 0.988 (from a e = 0.0 presupernova value). Such an influence of a natal kick on the binary orbit is very rare (Andrews & Zezas 2019), so the above example is an extreme example of a BHBH merger formation. In the vast majority of systems, the impact of the accretion model on the BHBH formation is much smaller.
2.5. Tidal interactions
In this section, we discuss the efficiency of tidal spinup in close binaries. Our standard implementation of tides follows from Hurley et al. (2002) and is based on the standard Zahn (1977, 1989) theory as updated by more recent calibrations (Claret 2007). These prescriptions, which are laid out in Belczynski et al. (2008b), result in rather weak tidal interactions between stars in binary systems. This is particularly true for very massive stars (e.g., progenitors of BHs) that evolve so fast that the tides do not affect significantly their rotation.
In our classical BHBH formation scenario (Belczynski et al. 2016a), the only evolutionary phase at which tides could spinup either of the BH progenitors happens after the CE phase when the initially wide binary (a ≳ 1000 R_{⊙}) is reduced to a rather tight orbit (a ≲ 100 R_{⊙}). Such a configuration will consist of at least one stripped stellar core (WR star): BHWR, WRBH or WRWR binaries.
Recently Kushnir et al. (2016, 2017), Zaldarriaga et al. (2018), Hotokezaka & Piran (2017), Qin et al. (2018) investigated the strength and treatment of tides and they argued that WR stars may potentially be significantly spun up if placed in very close binaries with massive companions (e.g., immediate progenitors of BHBH mergers).
To determine the upper limit to the spins produced by tidal locking, we ignore the orbital evolution and assume that the entire star is tidally locked to the orbit (see the discussion below), meaning that it is rigidly rotating. The orbital period then provides the spin period of the WR star and implementing these spin periods into our stellar models, we can calculate the angular momentum in the star and use this angular momentum to obtain the spin period of the collapsed star, assumed to be the BH spin. Figure 5 shows the resultant blackhole spin magnitudes for our Z = 0.014 and Z = 0.0004 MESA models. For most of the models with orbital periods in the range of P_{orb} = 0.1−1.3 d, the resultant BH spin magnitude can be fit by the exponential:
Fig. 5.
Black hole spin magnitudes as a function of the orbital period for our Z = 0.014 (circles) and Z = 0.0004 (triangles) MESA models. The color coding corresponds to the black hole remnants with masses: M_{BH} < 15 M_{⊙} (black/blue), 15 M_{⊙} < M_{BH} < 30 M_{⊙} (cyan/red), M_{BH} > 30 M_{⊙} (green/magenta). Binary systems with short orbital periods: 0.1−1.3 d and WR stars will produce BHs with broad range of spin magnitudes: 0.15−1. Systems with orbital periods below 0.1 d form black holes from tidally locked WR stars, and the BH spins are maximal. Binaries with an orbital period above 1.3 d produce BHs with spins below 0.2 and typically are not tidally locked at all. A set of the most massive, lowestmetallicity stars have such dense cores that tidal spinup does not dramatically increase their angular momentum and these stars require orbital periods of less than 0.1 d to have spin values above 0.1 (lower set of green triangles). However, we ignore this fact, and let these cores to be spunup to estimate the maximal effect of tidal interactions in our models (higher set of green triangles). The black curve shows our fit to the BH spin magnitude as a function of orbital period. 
where P_{orb} is the orbital period in seconds and P_{0} = 4000 s. For systems with orbital periods below 0.1 d, the resultant BH spin is maximal. For systems with orbital periods longer than ∼1.3 d, tidal locking takes longer than the duration of the relevant evolutionary phase and there is no significant spinup (in which case the tidal spinup is ignored). The fit with Eq. (15) is shown as the black curve in Fig. 5.
We note that the proposal by Kushnir et al. (2016, 2017), Zaldarriaga et al. (2018), Hotokezaka & Piran (2017), according to which every BH formed from a WR star subject to tidal interactions is spunup to maximum (a_{spin} = 1), is subject to three caveats. First, for systems that are undergoing tidal synchronization it is assumed that the BH is formed with a maximal spin without close scrutiny of the WR star structure and its evolution under strong tidal interactions. Second, the scheme ignores the fact that WR wind massloss (depending on the metallicity of the WR star) may widen the system pushing it into a regime where synchronization is not maintained. Third, the tidal locking of the WR star pumps orbital energy into the WR star spin and thus causes the orbit to decay. This may lead to a merger of a WR star with its companion barring the formation of a BHBH binary. A detailed analysis of this complex interplay of stellar structure, stellar winds, and tidal interactions in the context of long Gamma Ray Bursts and BHBH formation was presented by Detmers et al. (2008) and by Qin et al. (2018). In particular, Detmers et al. (2008) found that the majority of close BHWR systems that are subject to strong tidal interactions either evolve to long periods (for high metallicity) or undergo a component merger before the BHBH formation (for low metallicity). On the other hand, Qin et al. (2018) found that most close BHWR systems will increase their periods resulting in a wide range of secondary BH spins (a_{spin} = 0−1). In the light of these detailed calculations, our adopted simple model (see Eq. (15)) resembles the Qin et al. (2018) scheme and allows for a broad range of BH spin magnitudes if the WR star that is forming a BH is subject to strong tidal interactions. We should stress that we use this model only as an alternative to our standard assumption that the effect of tidal spinup in close binaries on the BH natal spin magnitude is ignored (as we use single stellar models to estimate natal BH spins; see Sect. 2.1).
Recently, Bavera et al. (2020) updated some of the Qin et al. (2018) calculations. We should note that these studies follow the angular momentum transport in stars only during the WR stage that will form the secondborn BH in BHBH merger. Additionally, it is assumed that the first BH is born with zero spin (a = 0) and that the WR star that will form the second BH is also born spinless. Such simplifying assumptions allow to match LIGO/Virgo BHBH effective spins which are mostly consistent with zero. However, it was only tested for initially nonspinning WR stars, and therefore did not take into account the fact that for inefficient angular momentum transport and high initial stellar rotation, WR stars can be born with high spins. Then for highspinning WR stars and for inefficient angular momentum transport a significant fraction of the secondborn BHs would have high spins independent of tidal interactions. This is demonstrated by our models that employ inefficient (Geneva) angular momentum transport, in which BHs (both firstborn and secondborn) may have high spins (see Figs. 1 and 18).
2.6. Cosmic star formation history and metallicity evolution
Since 2016 (e.g., Belczynski et al. 2016a), we have been using the cosmic star formation density (SFRD) determinations from Madau & Dickinson (2014), which are based on a number of deep UV and infrared galaxy surveys. We refer to this SFRD as the “old SFRD” formula. Here, we adopt two bestfitting comoving SFRDs from Madau & Fragos (2017), which update the previous formula by better reproducing a number of recent 4 < z < 10 results:
The correction factor 𝒦_{IMF} adjusts the SFRD for the assumed IMF to the Salpeter IMF (i.e., 𝒦_{IMF; Salpeter} = 1.0). In this work we adopt a three component broken powerlaw Kroupa IMF with α_{3} = −2.35 (see Appendix A.7 for details), for which the correction factor is 𝒦_{IMF; Kroupa} ≈ 0.66 (Madau & Fragos 2017).
We refer to the SFRD given by Eq. (16) as to the “low SFRD”, since it is based on a conversion from luminosity density to SFRD in which all published galaxy luminosity functions (LFs) have been conservatively integrated down to the same limiting luminosity of 0.03 L_{*}, where L_{*} is the characteristic luminosity of a Schechter function. For completeness, in the following we shall also provide results for a “high” model in which the comoving SFRD is increased at high redshifts to account for a steep faint end of the LF:
Hereafter, we refer to this formula as the “high SFRD”. These expressions bracket uncertainties in the contribution to the early SFRD by galaxies fainter than −16 mag, but agree at z < 2. Our adopted SFRDs are shown in Fig. 6.
Fig. 6.
Cosmic star formation histories adopted in our modeling (see Sect. 2.6 for details). We note how the different SFRDs agree at low redshifts (z ≲ 2), while they disagree by as much as a factor of three at earlier epochs. 
Our modeling of the metallicity of the starforming gas has also been updated as follows. In previous calculations, we used the evolving mean cosmic metallicity Z(z) from Madau & Dickinson (2014), which is the ratio between the total mass density of heavy elements produced over cosmic time and the cosmic baryon density. This quantity had to be increased by 0.5 dex to better account for various supernovae and GRB observations (Belczynski et al. 2016a). At each redshift, z = 0−15, we further assumed that the distribution of log(Z/Z_{⊙}) was a Gaussian centered at the average metallicity and with a dispersion of σ = 0.5 dex. Stellar populations (Population I/II stars) with various metallicities were subsequently evolved, to check whether they produce NSNS, BHNS or BHBH mergers detectable by LIGO/Virgo for a given instrument’s sensitivity (see Appendix A.8).
This approach does not properly describe the metallicity evolution of the starforming gas within galaxies, as only a small fraction of the baryons in the Universe are polluted with heavy elements and take part in the baryon cycle of galaxy evolution (e.g., Chruslinska et al. 2019). Our improved calculations follow Madau & Fragos (2017). The gasphase oxygen abundance is known to correlate strongly with the total stellar mass of starforming galaxies: this “massmetallicity relation” (MZR) has been shown to extend down to lowluminosity galaxies with stellar masses ∼10^{6} M_{⊙} (Berg et al. 2012) and out to a redshift of 3.5 (Maiolino et al. 2008). Numerical modeling by Zahid et al. (2014) suggests that the MZR originates from a more fundamental, universal relationship between the metallicity and the stellartogas mass ratio that is followed by all galaxies as they evolve. We have assumed that the Zahid et al. (2014) MZR holds at all redshifts, and integrated this relation over the evolving galaxy stellar mass function at 0 < z < 7 to compute a mean stellar massweighted gasphase metallicity (see Madau & Fragos 2017 for details and references) as:
As before, we assume the same distribution of log(Z/Z_{⊙}), but centered around this new average. We note that this is a simplification and this distribution may not be symmetric (Chruslinska et al. 2019). Recently, Chruslinska & Nelemans (2019) contributed an observationbased distribution of the star formation rate density over metallicities and redshifts and discussed the uncertainty of this distribution. We note that the metallicity distribution used in our paper peaks at similar metallicities as in the highmetallicity extreme reported by those authors at z ≲ 3 and at noticeably higher metallicities at z > 3. However, at such high redshifts, the distribution is poorly constrained by current observations.
There is no consensus on the value of solar metallicity, and no value within the range of Z_{⊙} = 0.012−0.02 can be rejected at the moment (Vagnozzi 2019). In our models we adopt either Z_{⊙} = 0.014 or Z_{⊙} = 0.02. A comparison between the adopted old and new mean gasphase metallicities versus redshift is shown in Fig. 7.
Fig. 7.
Stellar massweighted gasphase metallicity versus redshift, Z(z). At every epoch, we adopted a Gaussian distribution of log(Z/Z_{⊙}), centered at the mean metallicity and with dispersion σ = 0.5 dex. We note that the new metallicity for the star forming gas is noticeably higher than the mean cosmic metallicity adopted in the past (see Sect. 2.6). 
2.7. Calculations
Our binary evolution calculations were performed with the upgraded population synthesis code StarTrack (Belczynski et al. 2002, 2008b). Improvements to the code include updates to the treatment of the common envelope (CE) evolution, compact object mass calculations including the effect of pairinstability pulsation supernovae and pairinstability supernovae, and new BH natal spin prescriptions (Sect. 2.1), among other upgrades (see Appendix A.7).
We considered fourteen different realizations (models) of our classical isolated binary evolution to test whether it is possible to form LIGO/Virgo BHBH mergers with the observed rates, masses, and effective spins. The first two models correspond to our previous calculations with fallbackdecreased, BH massdependent (M10), and high, massindependent (M13) BH natal kicks with input physics listed in Table 2 and detailed in (Belczynski et al. 2016a,c).
Binary evolution models.
The next four models include different input physics on mass transfer, BH accretion in the CE phase, and the approximation of effects of stellar rotation on the BH mass (see Appendix A.7). These models include natal kicks which are fallbackdecreased, BH massdependent (M20); small, BH massindependent (M26; σ = 70 km s^{−1}); intermediate, BH massindependent (M25; σ = 130 km s^{−1}); and high, BH massindependent (M23; σ = 265 km s^{−1}).
The above six models employ BH natal spins adopted from the Geneva model with mild angular momentum transport in massive stars (Eq. (3)), SFRD and average metallicity evolution with redshift from Madau & Dickinson (2014), and high value of solar metallicity Z = 0.02. The next eight models represent our current update (2019) and tests of input physics.
M30 employs the rapid supernova engine model for NS/BH mass from Fryer et al. (2012) supplemented with weak PPSN (Eq. (9)), 1% neutrino mass loss at the BH formation and 10% neutrino mass loss at the NS formation (Eq. (10)), fallbackdecreased, BH/NS massdependent, natal kicks (1D σ = 265 km s^{−1}), 50% nonconservative RLOF, 5% BondiHoyle accretion during CE, inefficient accretion onto BH/NS during stable mass transfer and capture from stellar winds (Sect. 2.4), BH natal spins from MESA (Eq. (4)), SFRD and average metallicity evolution from Madau & Fragos (2017) (Eqs. (16) and (18)), and we adopt a low value for the solar metallicity Z = 0.014. In this model, we do not take into account additional effects of rotation as is done in model M20, we use the standard stellar winds from Vink et al. (2001) with addition of LBV winds from Belczynski et al. (2010a), and we employ the initial binary parameters from Sana et al. (2012) as discussed in de Mink & Belczynski (2015).
Models M33 and M35 are the same as model M30, but with mass independent NS/BH natal kicks (not affected by fallback): high natal kicks with 1D σ = 265 km s^{−1} (M33) and intermediate natal kicks with 1D σ = 130 km s^{−1} (M35).
Model M40 is the same as model M30, but with a different model for BH spins: the Fuller model (Eq. (5)). Model M43 is the same as model M40, but with massindependent NS/BH natal kicks with 1D σ = 265 km s^{−1}.
Model M50 is the same as model M30, but with mass loss rates reduced to 30% for all the stars (see Appendix A.7 for a justification).
Model M60 is the same as model M30, but with strong PPSN (Eq. (7)) with 10% neutrino mass loss for both BH and NS.
Model M70 is the same as model M30, but with moderate PPSN (Eq. (8)).
Models that allow CE with HG donors are marked as submodels MXX.A, while models that do not allow CE with HG donors are marked as as MXX.B (Belczynski et al. 2007).
In all models, we assume BH natal kicks to be randomly oriented, thus generating BH spin misalignments with respect to the orbital angular momentum. We note that our models are by no means exhaustive in terms of probing the evolutionary uncertainties. However, they allow us to test the key parameters that set BH spin magnitudes and misalignments, and thus the effective spin: angular momentum transport in massive stars, accretion onto the BH and its progenitor, tidal interactions and BH natal kicks. Table 2 gives an overview of the models, and in Appendix A.7 we describe the details of our binary population synthesis calculations.
The StarTrack population synthesis code is used to generate populations of BHBH/BHNS/NSNS systems. The star’s initial properties (mass and metallicity) are used to calculate stellar evolution. Binary interactions (mass gain and loss in RLOF and CE events) are taken into account through rejuvenation and derejuvenation in estimating the final stellar properties (mass and CO core mass) at the time of corecollapse, which are then used to obtain the NS/BH natal mass (see Appendix A.3). In binary calculations, we use nonrotating stellar models (Hurley et al. 2000). We record whether any of the binary components accreted significant amounts of mass (≳10% of its own mass). If accretion occurred during main sequence and if the accreting star was of low metallicity (Z < 0.002), then we assume in models M20, M23, M25, and M26 that the star will produce a more massive CO core (greater by 20%) and, thus, a more massive NS/BH to mimic the effects of increased mixing due to rapid rotation induced by accretion (see Appendix A.7).
At BH formation we use single stellar models to estimate the BH natal spin magnitude through the COcore massspin relations proposed in Sect. 2.1. In these estimates, we assume that a star with a given CO core mass (as estimated from a binary evolution) forms a BH with spin given by single stellar models. This scheme ignores the effects of mass accretion on the stellar spin of the BH progenitor that may increase the BH’s natal spin. Obviously, this is far from perfect, and stellar rotation models will need to be fully integrated into binary population synthesis in the future. However, this is impossible at the moment due to the lack of stellar models with consistent input physics that would appropriately sample mass, metallicity and rotation for massive stars and naked stellar cores (as stellar cores are much more often exposed in binary evolution than is single star evolution).
We note also that we use moderately high initial stellar spins that assume 40% of breakup velocity in single stellar models. Depending on a mass and stellar structure of a model, this gives a range of 250−450 km s^{−1} initial rotation speeds at the equator (see Tables A.1 and A.2). For comparison, the observed spins of massive stars show a bimodal distribution (RamírezAgudelo et al. 2013), with one large peak at ∼100 km s^{−1} and another small peak at ∼400 km s^{−1}.
A firstborn BH can accrete mass from its unevolved companion either during the RLOF/CE stages or from the companion’s winds. These effects are included in our calculations, and the BH spin is increased accordingly. The secondborn BH preserves its natal spin as it does not accrete mass. However, we note that we allow for the possibility of the tidal spin up of BH progenitors and allow for significant increase of the BH spin as compared to BH spins that result from single stellar models (Sect. 2.5).
In all cases, we assume that the stellar spins are initially aligned (Θ_{1} = Θ_{2} = 0) with the orbital angular momentum of the main–sequence star binary. At BH formation we estimate the spin vector misalignment due to natal kicks. We allow for BH spin realignment neither during the mass transfer phases, nor by tidal interaction between the stars in a binary. We note that with this approach we may be overestimating the BH spin misalignment, but only for models with high to moderate BH natal kicks (see Sect. 2.2 for details).
3. Results
We estimated the double compact object merger rate densities, merger detection rates, merger masses, and BHBH effective spins using the methods presented in Belczynski et al. (2016b,a) along with updates and revisions presented in the present study (see Sect. 2). In the following sections, we discuss some particular properties of our models. In our study, we do not exhaust the information that can be extracted from our models. Focusing on BHBH mergers, we instead show some particular examples of what can be obtained with population synthesis modeling in context of the LIGO/Virgo sources. For now, we compare our models to the LIGO/Virgo observations showing that some models fit the data better than other, but only in terms of the observed allowed ranges of rates, masses, and effective spins. We note that we do not yet attempt to match particular distributions’ shapes (for BH masses and effective spins) or the rate of increase of merger rates with redshift. However, anyone interested in such comparisons can easily perform them on our models as we make them publicly accessible through our website^{3}^{,} ^{4}.
3.1. Binary evolution of BHBH mergers
In this section, we present several examples of binary evolution leading to the formation of BHBH mergers. In the framework of the Geneva model of angular momentum transport, it is challenging to explain mergers with very low effective spin parameters as lots of BHs are formed with high or moderate spin magnitudes (see Fig. 1). Yet, it is not impossible. Therefore, we show examples of evolution that can lead to the formation of merger resembling GW170104, which has one of the lowest measured effective spin parameters: −0.24 < χ_{eff} < 0.13 (90% confidence limits). In the framework of MESA angular momentum transport, it is challenging to explain mergers with moderate and high effective spin parameters as lots of BHs are formed with low spin magnitudes (see Fig. 2). Yet, we show that we can form mergers that are consistent even with GW170729, which is the merger that has the highest effective spin yet measured: 0.11 < χ_{eff} < 0.57 (90% confidence limits).
3.1.1. Case of GW170104
Here we present a proofofprinciple scenario demonstrating that isolated binary evolution with the Geneva model of angular momentum transport can form a BHBH merger with BH masses and effective spin compatible with LIGO’s observation of GW170104; in particular, given its very low effective spin.
Within our models with the Geneva angular momentum transport (M10, M13, M20, M23, M25, M26), we search for systems with BH masses and effective spins within LIGO’s 90% credible limits: 25.4 < M_{BH1} < 38.2 M_{⊙}, 15.6 < M_{BH2} < 25.0 M_{⊙}, −0.24 < χ_{eff} < 0.13. The upper bound on χ_{eff} may actually be as high as χ_{eff} ≈ 0.2 (Appendix A.5). For example, within model M20, it is indeed possible to form a BHBH merger resembling GW170104: M_{BH1} = 33.3 M_{⊙}, M_{BH2} = 24.7 M_{⊙}, χ_{eff} = 0.09. The evolutionary history of such a merger is presented in Fig. 8. We note that model M20 is rather conservative regarding assumptions on natal kicks, which are strongly suppressed by the fallback material. Massive BH spins are mostly aligned with the binary angular momentum (cosΘ_{1} = cosΘ_{2} = 1), thus maximizing the value of the effective spin. For all the other models (with the exception of M10) the BH spins will tend to be misaligned, decreasing the value of the effective spin and making it even easier to produce systems with low effective spin values (as observed in GW170104).
Fig. 8.
Example of a possible route leading to the formation of a BHBH merger similar to GW170104. This example follows the classical isolated binary evolution channel. In this model (M20.B) we employ the Geneva BH natal spins, assume that massive BHs do not receive natal kicks and that their spins are aligned with the binary angular momentum (Θ_{1} = Θ_{2} = 0°), producing an upper limit on the effective spin parameter (χ_{eff}). Yet, this system has χ_{eff} = 0.09, within LIGO’s 90% credible limits for GW170104 [ − 0.24:0.13]. Both BH masses are also within the limits: M_{BH1} = 31.0 M_{⊙} [25.4:38.2] and M_{BH2} = 20.1 M_{⊙} [15.6:25.0]. 
Although the BH spin can be modified by accretion from a binary companion, the amount of matter accreted in our calculations is very modest, and the accretioninduced spinup of BHs is not significant. In the example shown in Fig. 8, the firstborn massive BH forms with no spin (a_{spin1} = 0), then it accretes very little mass from its MS companion wind increasing its spin only to a_{spin1} = 0.02. Most of the accretion occurs during the CE phase (0.4 M_{⊙}) and the BH increases its spin to a_{spin1} = 0.05. Finally, the BH accretes very little mass from its WolfRayet star companion wind, increasing its spin to its final value of a_{spin1} = 0.053. The secondborn BH forms with a natal spin of a_{spin1} = 0.14 and it is not spun up, as it does not accrete any material. In other words, we predict that LIGO/Virgo observations of BHBH mergers will probe the natal BH spin distribution, up to evolutionary effects of order 0.05 in dimensionless spin.
We also note that this particular system is not subject to a potential WRstar tidal spinup in a BHWR binary (the last evolutionary stage before the BHBH formation). The system separation at this stage (a ≈ 35 R_{⊙}) is too large for the tides to effectively spinup WR star (P_{orb} > 1.3 d, see Sect. 2.5). The full details of this evolutionary example are given in Appendix A.4.
3.1.2. Case of GW170729
In a subset of our binary evolution models M30.A/B, M33.A/B, M35.A/B, M50.A/B, M60.A/B, and M70.A/B, the natal BH spins are obtained from stellar models calculated with the MESA code under the assumption of efficient angular momentum transport in the stellar interiors (see Sect. 2.1 and Fig. 2). In this framework, the initial BH spins are always small (a_{spin} ≲ 0.15), mostly independent of the progenitor mass and metallicity. Small natal BH spin values can, in principle, be increased during further evolution as result of mass accretion. However, in the isolated binary evolution channel for BHBH mergers, the first formed BH can only accrete mass during the CE inspiral and through accretion of the wind from its stellar companion. We find that in our simulations neither of those two processes leads to a significant increase of the BH spin which is primarily the consequence of: (i) inefficient (5−10%) BondiHoyle accretion rate onto BH in the shortlived CE phase (e.g., MacLeod et al. 2017) and (ii) small wind mass loss rates from low metallicity stars (Vink et al. 2001). This means that the small natal BH spins in the framework of efficient angular momentum transport in stellar interiors do result in small effective spin values of BHBH mergers (typically χ_{eff} ≲ 0.25).
Here we present a proofofprinciple scenario demonstrating that even the BHBH merger event with the largest effective spin value reported to date (GW170729, ) can be reconstructed in the isolated–binary evolution channel with small natal BH spin values given by Eq. (4). Within our models with efficient (MESA) angular momentum transport we search for systems with BH masses and effective spin within LIGO’s 90% credible limits for GW170729: 40.4 < M_{BH1} < 67.2 M_{⊙}, 24.2 < M_{BH2} < 43.4 M_{⊙}, 0.11 < χ_{eff} < 0.57. For example, within the model M30.B it is indeed possible to form a BHBH merger resembling GW170729: M_{BH1} = 55.0 M_{⊙}, M_{BH2} = 32.8 M_{⊙}, χ_{eff} = 0.137. The evolutionary history of such a merger is presented in Fig. 9. Both the BHs where formed in direct collapse with small but nonzero natal spins (a_{spin1} = 0.093 and a_{spin2} = 0.073). For reasons discussed above, the BH formed first accreted only a very modest amount of mass during binary evolution (< 1.0 M_{⊙}), which led to an increase of its spin to a_{spin1} = 0.176. The final effective spin of the merging BHBH system is χ_{eff} = 0.137, and thus it is in (marginal) agreement with the 90% credible limits for GW170729.
Fig. 9.
Example of a possible route leading to the formation of a BHBH merger similar to GW170729. This example follows the classical isolated binary evolution channel. In this model (M30.B) we employ the MESA BH natal spins, assume that massive BHs do not receive natal kicks and that their spins are aligned with the binary angular momentum (Θ_{1} = Θ_{2} = 0°), producing an upper limit on the effective spin parameter (χ_{eff}). We note that given the small binary separation at the BHWR stage (a ≈ 7 R_{⊙}), the WR star is most likely going to become tidally synchronized (see Sect. 2.5). If we assume that the rapidly rotating WR star collapses into rapidly spinning BH (a_{spin2} = 1) then the effective spin of the presented system would increase from χ_{eff} = 0.137 (no tidal spinup; presented in the figure) to χ_{eff} = 0.484 (full tidal spinup). Both values are within the LIGO/Virgo 90% credible limits for GW170729 [0.11:0.57]. Both BH masses are also within the limits: M_{BH1} = 50.6 M_{⊙} [40.4:67.2] and M_{BH2} = 34.3 M_{⊙} [24.2:43.4]. 
It should be noted that the spin value for the secondary BH presented in the above example (Fig. 9; a_{spin2} = 0.073) is most likely underestimated because this example neglects the potential tidal spinup of the secondary WR star (direct progenitor of the second BH) during the BHWR stage (see Sect. 2.5). Since the orbital separation during the BHWR phase is very small (a ≈ 7 R_{⊙}), tidal interactions are expected to be efficient in spinningup the WR star. In fact, in the case of this system the timescale of WR tidal synchronization as approximated by Zaldarriaga et al. (2018, see our Eq. (15)) is only t_{sync} ≈ 11 yr, which is a few orders of magnitude shorter than the duration of BHWR stage (∼2 × 10^{5} yr). Thus, the secondary WR star has most likely become fully synchronized with the orbit by the time it collapsed to form the second BH. If we assume that such a rapidly rotating WR star collapses into an also rapidly spinning BH (a_{spin2} = 1.0) then the effective spin of the BHBH merger in Fig. 9 would become χ_{eff} = 0.484, still well within the 90% credible limits for GW170729. The large range of χ_{eff} values that spans from 0.137 to 0.484 encloses all the possible results of the uncertain process of tidal spinup in BHWR binaries and its effect on the secondary BH spin. As a word of caution, we note that the simple fact that a particular BHBH system can be reproduced in a given binary evolution channel does not in itself guarantee consistency between population synthesis results and observations. Such consistency can only be tested by comparing the entire populations. We discuss the distributions of BHBH mergers parameters obtained in our simulation in the following sections.
3.2. Merger rate density and detection rate
Table 3 summarize our predictions for the local (redshift z ∼ 0) merger rate density, the detection rate for LIGO/Virgo’s midhigh sensitivity curve (that may be taken as an approximation of O3 observing run) along with the maximum horizon redshift for the best located and oriented source in a given merger type category: NSNS, BHNS, BHBH.
Local merger rate densities and LIGO/Virgo detection rates.
The predicted BHBH merger rate densities vary between 1.24 and 1368 Gpc^{−3} yr^{−1} across our models. A number of models (M13.A, M23.A/B, M25.B, M26.B, M30.B, M33.A, M40.B, M43.A, M50.B, M60.B, M70.B) produce rates within the allowable range determined by the first 10 LIGO/Virgo detections (9.7−101 Gpc^{−3} yr^{−1}: Abbott et al. 2019b). The predicted NSNS merger rate densities vary between 49.3 and 524 Gpc^{−3} yr^{−1} for the tested models. A number of models (M10.A,M13.A, M20.A, M23.A, M25.A, M26.A, M30.A/B, M33.A, M35.A/B, M40.A/B, M43.A, M50.A, M60.A/B, M70.A/B) produce rates within allowable range determined by the first LIGO/Virgo detection (110−3840 Gpc^{−3} yr^{−1}: Abbott et al. 2019b). The predicted BHNS merger rate densities vary between 0.48 and 297 Gpc^{−3} yr^{−1} for the tested models. So far, all the models produce rates within an upper limit determined by the nondetection in the O1/O2 LIGO/Virgo runs (< 610 Gpc^{−3} yr^{−1}: Abbott et al. 2019b).
Figure 10 shows the intrinsic BHBH merger rate density (not weighted by LIGO/Virgo detection probability) evolution with redshift for models M10.B and M30.B that employ different star formation rate and cosmic evolution of metallicity. We note that model M10.B which employs Z(z) measured from stars generates higher (by factor of ∼5 at low redshifts z < 2) rates than the model M30.B that employs Z(z) as measured by the metallicity of star forming gas at any given redshift. The choice of Z(z) is one of the most important factors affecting local (low z) merger rate densities for BHBH mergers. This is because with decreasing metallicity we expect (i) an increase of the BH mass (Belczynski et al. 2010a) and (ii) an increase of the BHBH merger formation efficiency per unit mass (Belczynski et al. 2010b). Therefore, models with lower Z(z) result in higher BHBH merger rates.
Fig. 10.
Merger rate density of BHBH mergers from Population I/II stars. We employ models M10.B and M30.B to illustrate the effects of our assumptions on the star formation rate and cosmic metallicity evolution on the BHBH mergers. LIGO/Virgo O1/O2 constraint on the BHBH merger rate in local Universe is also shown. We note that the decrease of the BHBH merger rate density at lowredshift (from old to new models) is due to the average metallicity of stars at any given redshift in old models, as it is much lower than the average metallicity of stars in new models (see Sect. 3.2 for details). 
We note that the difference in SFRD(z) between the two models does not affect significantly the rates as the differences in star formation are relatively small (see Fig. 6). We also note that for the two updated (new) models of the SFRD(z), there is virtually no difference in the BHBH merger rates at low redshifts (z < 2), while differences (∼ factor of a few) begin to appear only at higher redshifts. For our updated models, we use a SFRD(z) with low star formation at high redshifts. This particular choice does not affect any of our conclusions for LIGO/Virgo as these instruments are not expected to probe sources at high redshifts (z > 2).
There are different PPSN models employed in models M10.B (strong PPSN; low maximum BH mass) and M30.B (weak PPSN; high maximum BH mass). However, this does not affect significantly the intrinsic merger rate density in a given volume that is presented in Fig. 10. This is because the number of BHBH binaries is about the same in each of the PPSN models. However, the differences in BH masses affect notably the detection rates. For example, compare the detection rate of BHBH mergers in models M30.A/B and M60.A/B that differ only by PPSN input physics (see Table 3).
Figure 11 demonstrates the effects of the CE phase and natal kicks on the intrinsic BHBH merger rate density. We use a sequence of four models with exactly the same input physics that differ only with regard to their various assumptions about the CE phase and natal kicks. With the optimistic approach to the CE and natal kicks (that allows for survival of the CE phase with HG donors and very low or zero BH natal kicks; model M30.A) the local (z = 0) merger rate density is relatively high: R_{BHBH} = 641 Gpc^{−3} yr^{−1}. When we apply a more restrictive approach to the CE phase, but keep the same natal kicks as above (CE survival not allowed for HG donors; model M30.B), then we note a significant decrease in the rate: R_{BHBH} = 43.7 Gpc^{−3} yr^{−1}. Addition of the moderate natal kicks while keeping the restrictive CEphase approach (model M35.B) decreases the rate further: R_{BHBH} = 4.69 Gpc^{−3} yr^{−1}. An additional increase of natal kicks (model M33.B) with the same restrictive CEphase approach, brings the rate to a very low value: R_{BHBH} = 1.37 Gpc^{−3} yr^{−1}. A comparison with the LIGO/Virgo rate estimate shows that models M30.A, M35.B, M33.B are excluded. It indicates that within our limited sample of evolutionary models, some combinations of CEphase approach and natal kicks can be excluded as not likely. For example, moderate to high natal kicks and a restrictive CEphase treatment are an unlikely combination. Note, however, that we cannot draw conclusions about the rates based solely on the natal kicks or solely on the CEphase approach because the results are degenerate with respect to these two major factors. For example, if we apply the optimistic CEphase approach to models with moderate to high natal kicks then both these models fit within (or very close to) the LIGO/Virgo rate estimate. In particular, we find R_{BHBH} = 109 Gpc^{−3} yr^{−1} for model M35.A (moderate natal kicks and optimistic CE=phase) and R_{BHBH} = 19.5 Gpc^{−3} yr^{−1} for model M33.A (high natal kicks and optimistic CEphase).
Figure 12 shows the intrinsic merger rate density for all types of mergers: NSNS, BHNS and BHBH. As an example, we use model M30.A/B. As discussed above, the local BHBH intrinsic merger rate density for model M30.B (R_{BHBH} = 43.7 Gpc^{−3} yr^{−1}) is within the LIGO/Virgo estimate. The predicted local BHNS intrinsic merger rate density for model M30.B (R_{BHNS} = 11.1 Gpc^{−3} yr^{−1}) is within the LIGO/Virgo upper limit. The local intrinsic merger rate density for NSNS systems is only just above LIGO/Virgo 90% level lower limit for model M30.B (R_{NSNS} = 122 Gpc^{−3} yr^{−1}), while it is well within the LIGO/Virgo estimate for model M30.A (R_{NSNS} = 426 Gpc^{−3} yr^{−1}).
Fig. 11.
Merger rate density of BHBH mergers for various assumptions about the common envelope and BH natal kicks. The possibility of development and survival of the CE phase with a Hertzsprung gap donor, increases the lowredshift BHBH merger rate density by ∼1 order of magnitude: compare models M30.A and M30.B (a similar effect is found for other models). Natal kicks tend to decrease the BHBH merger rate density by ∼1.5 order of magnitude: from fallback attenuated natal kicks (almost no BH natal kicks: model M30.B) to full scale BH natal kicks (as high as observed for Galactic single pulsars: model M33.B). We note that the degeneracy between the tested input physics; by applying high natal kicks we can bring model M30.A down to agree with the LIGO/Virgo estimate, or by allowing for HG CE we can bring models M35.B and M33.B up to also match the LIGO/Virgo constraint. 
Fig. 12.
Merger rate density of NSNS, BHNS and BHBH mergers for model M30. For comparison, we show LIGO/Virgo O1/O2 constraints on merger rate densities. While actual estimates are available for the BHBH and NSNS mergers, only an upper limit is available for the BHNS mergers from O1/O2 data (but see Appendix A.9). We note that the presented merger rate densities are consistent with the LIGO/Virgo estimates. 
Detection rates (Table 3 is based on each source merger redshift, its mass (we use mass dependent waveforms for each merger) and take into account the LIGO detector antenna (peanutshaped) pattern. For this particular estimate we employ LIGO sensitivity labeled “midhigh” that approximately corresponds to the current (O3) LIGO sensitivity (for details, see Appendix A.8). The detection rate is then a convolution of the merger rate density (and its change with redshift within the LIGO horizon for a given source type) and the mass distribution of given source type (NSNS, BHNS, BHBH). For clarity, we list the horizon redshift for the best located/oriented (directly overhead source with an orbital plane perpendicular to the line of sight) source within each merger type. This naturally corresponds to the furthest detectable redshift (distance) of the most massive source within each merger type found in our simulations.
3.3. BH masses
In this section we discuss the masses of BHs in BHBH mergers, focusing on either the individual component masses, M_{BH1} (more massive component) and M_{BH2} (less massive component), or the total merger mass, M_{tot} = M_{BH1} + M_{BH2}.
Figure 13 shows the intrinsic (source frame) distributions of the individual BH masses (M_{BH1} and M_{BH2} shown together in one distribution) and the total BHBH merger masses for our three models of PPSN: M30.B (weak PPSN), M70.B (moderate PPSN), and M60.B (strong PPSN). These distributions are weighted by the intrinsic merger rate density for all BHBH mergers (independent of mass) within a given redshift (z < 2). As clearly seen, the different PPSN treatments in the models shown here impact the maximum BH mass generated. In our calculations, the individual BH masses in merging binaries extend to M_{BH, max} ∼ 40 M_{⊙} in M60.B, M_{BH, max} ∼ 50 M_{⊙} in M70.B, and M_{BH, max} ∼ 55 M_{⊙} in M30.B. For comparison, among all the candidates reported by LIGO/Virgo (see Table 1), the largest mean mass for any individual BH is 50.6 M_{⊙} (for GW170729); even most optimistically (within the 90% credible limits), the largest BH reported in a binary is only 67.2 M_{⊙} (again for GW170729). The total BHBH merger mass reaches M_{BH, max} ∼ 80 M_{⊙} in M60.B, M_{BH, max} ∼ 100 M_{⊙} in M70.B, and M_{BH, max} ∼ 110 M_{⊙} in M30.B. For comparison, LIGO/Virgo total mean BHBH mass estimates reach 86.3 M_{⊙}, while the 90% confidence level on the most massive BHBH merger is as high as 100 M_{⊙}. The distributions approximately resemble steep powerlaws: ∼M^{−3.6} for individual BH mass, and ∼M^{−4.0} for total BHBH mass. [Unless otherwise noted, all exponents refer to a power law fit to all merging binaries, covering the full mass range]. The power law index depends weakly on the PPSN model.
Figure 14 illustrates the effect of different Z(z) assumptions on source frame BH mass distributions. We choose models M10.B and M60.B that utilize the same strong PPSN model, but employ slow Z(z) evolution based on an older estimate and faster Z(z) evolution based on a more recent estimate, respectively. For both models the range of BH masses (whether individual BH masses or total BHBH merger masses) is very similar. The same PPSN model leads to a similar cutoff BH mass at the high mass end. However, we note the change in steepness of the powerlaws which approximate these distributions. For individual BH masses the slope changes from ∼M^{−2.9} for model M10.B to ∼M^{−3.3} for model M60.B. For total BHBH masses the slope changes from ∼M^{−1.7} for model M10.B to ∼M^{−2.9} for model M60.B. This comes from the fact that in models with high numbers of low metallicity stars (e.g., M10.B) BH masses are (on average) higher and therefore the distributions are flatter than in models with small number of low metallicity stars (e.g., M60.B).
Figure 15 shows the intrinsic (source frame) distributions of the more massive BH mass (M_{BH1}) weighted by the intrinsic merger rate density for all BHBH mergers (independent of mass) within a given redshift (z < 2). Here we present two groups of models. In one, we alter the amount of PPSN massloss: M60.B (strong massloss), M70.B (moderate massloss) and M30.B (weak massloss). All these models tend to have similar and rather steep powerlaw distribution of larger BH mass (∝M^{−3.6}). In the second group we show models that tend to produce shallower powerlaw distributions (∝M^{−2.4}): M10.B (old Z(z) relation with standard stellar winds and restrictive CE treatment), M50.A (new Z(z) relation with 30% reduced stellar winds and optimistic CE), and M50.B (new Z(z) relation with 30% reduced stellar winds and restrictive CE). Clearly models that tend to produce more low metallicity stars (e.g., M10) and models with reduced stellar winds (e.g., M50) generate more massive BHs and produce shallower distributions. However, our mass distributions are not perfect power laws, particularly at low mass. For example, a power law approximation to our simulation data restricted to the mass range of reported LIGO/Virgo observations will recover a steeper power law exponent. For comparison, LIGO/Virgo infers a phenomenological pure power law distribution of ∝M^{+0.1} − M^{−3.1} (Abbott et al. 2019a based on observations spanning M_{BH1}; see their Fig. 1). When comparing our calculations to LIGO/Virgo observations in this way, the estimated power law exponent depends sensitively on the choice of the lowmass cutoff M_{BH1}, where our models predict a substantial structure in the mass distribution.
Fig. 13.
BH masses in BHBH mergers within the design advanced–LIGO horizon (z < 2) for models M60.B (strong pairinstability pulsation supernovae), M70.B (moderate PPSN), and M30.B (weak PPSN). Top panel: distribution of individual BH masses, bottom panel: total BHBH system mass. All distributions are intrinsic (neither redshifted nor weighted by detection probability). For comparison, we also indicate the range of LIGO/Virgo O1/O2 mean mass estimates and their most narrow and the widest allowed range within 90% confidence limits. Individual BH masses may reach ∼40 M_{⊙} (M60.B) and ∼55 M_{⊙} (M30.B), while the total BHBH mass may reach ∼80 M_{⊙} (M60.B) and ∼110 M_{⊙} (M30.B). Note that these distributions only very approximately resemble powerlaws: ∼M^{−3.6} (for individual BH masses) and ∼M^{−4.0} (for total BHBH mass). Powerlaw fits (dashed black lines) were performed for model M30.B in the loglog space. 
Fig. 14.
BH masses in BHBH mergers within design advanced–LIGO sensitivity (z < 2) for models M10.B (that employs old average metallicity cosmic evolution) and M60.B (new average metallicity evolution). BH masses are calculated with the same formulas in both models. Top panel: distribution of individual BH masses, bottom panel: total BHBH system mass. All distributions are intrinsic (not redshifted nor weighted by detection probability). We note that application of new metallicity evolution with redshift (less lowZ stars: less high mass BHs) results in somewhat steeper BH mass distributions as contrasted with application of old metallicity evolution (more lowZ stars: more high mass BHs). Powerlaw fits (dashed black lines) were performed for both models in the loglog space. 
Fig. 15.
Larger (more massive) BH mass distributions in BHBH mergers within design advancedLIGO sensitivity (z < 2) for models M60.B, M70.B, M30.B, M10.B, M50.A, M50.B. Powerlaw fits (dashed black lines) were performed for model M30.B in the top panel (∝M^{−3.6}), and for models M50.A (∝M^{−2.3}) and M50.B (∝M^{−2.4}) in the bottom panel in the loglog space. These distributions can be compared with other theoretical predictions (e.g., Fig. 1 of Mapelli et al. 2019, or Fig. 4 of Stevenson et al. 2019 that seem to be close to ∝M^{−2} − M^{−3}; we note that powerlaws shown in the latter work are not correct) or with the LIGO/Virgo observational estimate of ∝M^{−1.6} with large uncertainty on the powerlaw index: ∝M^{+0.1} − M^{−3.1} (Fig. 1 in Abbott et al. 2019a). Some of our models (e.g., top panel) show somewhat steeper relations, while other models (bottom panel) are in agreement with LIGO/Virgo estimate. 
Figure 16 shows the total BHBH merger mass distributions weighted by the LIGO detection probability (dependent on mass). We show models with different assumptions on PPSN (M30.B, M60.B, M70.B) and contrast them with LIGO/Virgo observations. These distributions are flatter than the intrinsic mass distributions as more massive mergers are detected in larger volumes (and thus in larger numbers) than lower mass mergers. Naturally, the maximum BHBH mass cutoffs are the same as for the intrinsic distributions. W note that within the observational uncertainties, all three of these models can reproduce the range of the detected BHBH total masses. However, a closer inspection of the distribution of observational points shows an overabundance of points near M_{tot} ∼ 50−70 M_{⊙} with respect to the models. At this point, we do not attempt to match this apparent peak in the observed distribution (this can wait for another published set of detections from LIGO/Virgo) but we have only tested a single potential change to the input physics (see below) that may affect the shape of the total mass distribution for massive BHBH mergers.
In Fig. 17, we replot the model M30.B (with weak PPSN and restrictive CEphase treatment) employing the standard theoretical estimates of wind massloss rates (rather high wind massloss rates; Vink et al. 2001). Additionally, we show the distributions for models M50.A and M50.B for which we reduce all stellar wind mass loss to 30% of the values used in all other models, but the rest of the input physics is exactly the same as in model M30.A/B. We note the appearance of a very prominent peak in the total BHBH mass distribution at M_{tot} ∼ 80−100 M_{⊙} in model M50.B (weak PPSN and restrictive CEphase treatment) as compared to the model M30.B distribution. That peak disappears in model M50.A (weak PPSN and optimistic CEphase treatment), but the distribution becomes much flatter than for model M30.B. These models seem to better reproduce the distribution of observational points. In Sect. 4.4 we present a brief discussion comparing synthetic StarTrack data with LIGO/Virgo data on BH masses.
Fig. 16.
Detectionweighted BHBH merger total intrinsic (not redshifted) mass for models with different assumptions on the pairinstability pulsation supernovae: M70.B (strong PPSN), M60.B (moderate PPSN), M30.B (weak PPSN). For comparison, we also show LIGO/Virgo O1/O2 mean total mass estimates and their 90% confidence limits. 
Fig. 17.
Detectionweighted BHBH merger total intrinsic (not redshifted) mass for models with different assumptions on wind mass loss and common envelope development: M50.B (weak stellar winds and optimistic CE), M50.B (weak stellar winds and pessimistic CE), M30.B (strong stellar winds and pessimistic CE). For comparison, we also show LIGO/Virgo O1/O2 mean total mass estimates and their 90% confidence limits. We note that although no model seems to reproduce the observed LIGO/Virgo BHBH total mass distribution, we can expect that some interplay of several parameters (winds, CE, PPSN, Z(z)) may, in the future, possibly reproduce the shape of the observed distribution. For example, stronger PPSN and higher Z(z) (not shown here) will tend to shift the highest BH masses to lower values. This calls for further parameter study calculations. 
3.4. BHBH effective spins
In this section, we present our predictions for the effective spins of BHBH mergers. Having shown that we can reproduce the effective spins (whether low or high) with some of our adopted BH spin models (Geneva or MESA), we wish to see whether we can also reproduce the effective spin distribution (or, rather, the range of observed values) of all of the reported LIGO/Virgo BHBH mergers.
In Fig. 18, we show the measurements of the effective spins in the LIGO/Virgo BHBH merger observations, superimposed with the χ_{eff} predictions calculated using the BH natal spins from the Geneva evolutionary calculations. All of the LIGO/Virgo observations cluster around zero (χ_{eff} ∼ 0). For the comparison, we use two models: model M20.B, with fallbackdecreased BH natal kicks (effectively no kicks for massive BHs, and small kicks for low mass BHs); and the model M23.B, with high BH kicks that are independent of the BH mass (with 1dimensional σ = 265 km s^{−1}). Our model distributions peak at high values (χ_{eff} ∼ 0.9) and then quickly fall off toward small values (see the logarithmic scale on Fig. 18).
Mergers with high values of the effective spin host BHs with high spins. In the Geneva model, these BHs can form with either low or high masses, depending on the chemical composition of the progenitor stars (Fig. 1). Mergers with low effective spins originate predominantly from systems with highmass BHs. Strong natal kicks (as in model M23.B) produce higher misalignment of the BH spins with respect to the binary’s angular momentum vector during the BH formation, and this decreases the value of the effective spin. We note that although both distributions are similar and peak at high values, natal kicks decrease the effective spin parameter (average χ_{eff} = 0.3; M23.B) as compared to model with almost no BH kicks (χ_{eff} = 0.7; M20.B). The fraction of BHBH mergers with negative χ_{eff} is sizeable in model M23.B (27.4%), while it is negligible in model M20.B (∼0.3%).
We can recover all the reported values within both models (more easily with model M23.B than with M20.B), but our predicted χ_{eff} distributions using Geneva spins as inputs are not consistent with the current LIGO/Virgo data. For example, for each of our models and for N observations, we can ask whether we would expect at least one measurement with a maximum value of χ_{eff} above what has been observed thus far; that is, 1 − P(< χ_{eff})^{N}, where P(< χ_{eff}) is the cumulative distribution implied by the detectionweighted χ_{eff} distribution reported in Fig. 18. Even considering only subsets of events, such as the highmass or lowmass events, we find a small (1.1 × 10^{−4} for M23.B) to infinitesimal (6.7 × 10^{−10} for M20.B) probability that the predicted spin distribution is compatible with observations reported to date.
Fig. 18.
Detectionweighted distribution of effective spin parameter of BHBH mergers for model M20.B (fallback decreased BH kicks; no BH kicks for massive BHs) and M23.B (high BH natal kicks with 1D σ = 265 km s^{−1} for all BHs) with Geneva mildly efficient angular momentum transport. We note that both distributions peak at high values (χ_{eff} ∼ 0.9). Natal kicks decrease the effective spin parameter (average χ_{eff} = 0.3; M23.B) as compared to model with almost no BH kicks (average χ_{eff} = 0.7; M20.B). For comparison, we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. Although we can recover all the reported values, the predicted peak of χ_{eff} distribution is not coincident with the current LIGO/Virgo data. It indicates that BHs have typically lower spins than resulting from the Geneva model for BH natal spins, or that the detected BHBH mergers are not formed in the classical isolated binary evolution. 
We note that to reach this conclusion, we have used two extreme natal kick models. In particular, we tested a model with high BH natal kicks that tends to maximize the BH spinorbit misalignment and therefore decrease our predicted effective spins. Although BH natal kicks as high as those adopted in model M23.B (with average 3dimensional velocity of ∼400 km s^{−1}) cannot yet be observationally excluded (Belczynski et al. 2016b), it is unlikely for BHs to receive larger natal kicks than the ones that have been adopted in this model (Mandel 2016; Repetto et al. 2017). We note that in these models, we allow neither for any processes that could realign the BH spins, nor for an effective tidal spinup of stars in the binary progenitors of BHBH mergers. These processes might increase the effective spin parameter for some BHBH progenitors, thus shifting our results for Geneva model of BH natal spins even further away from the LIGO/Virgo observations.
In Fig. 19, we present the effective spin distribution for the model M30.B that employs the MESA BH natal spins. We show two versions of this model: one in which we do not allow and one in which we do allow for efficient tidal interactions in close binaries (that host WR stars; see Sect. 2.5). The version with no efficient tides shows a narrow distribution of effective spins: −0.2 ≲ χ_{eff} ≲ 0.2 that is peaked at positive values: average (χ_{eff} = 0.15). This limited range of effective spins follows directly from the underlying BH natal spin model that allows only for narrow range of BH spin magnitudes a_{spin} ∼ 0.05−0.15 (see Fig. 2). The BH spins are, at most, moderately increased by accretion from the binary companion (see Figs. 8 and 9) and this tends to slightly extend the distribution to higher positive values. The distribution is only moderately affected by smalltono BH natal kicks (fallback decreased kicks are employed in this model) producing a small population of BHBH mergers with negative effective spins. On the other hand the variation with the efficient WR tides generates a rather broad effective spin distribution: −0.5 ≲ χ_{eff} ≲ 1.0 with a peak at: χ_{eff} ∼ 0.15 (∼80%) and a tail with χ_{eff} ≳ 0.25 (27%). This is because in this variation about 1/3 of systems are subject to significant tidal spinup of at least one binary component and they produce large BH spins. The tail shows a significant drop beyond χ_{eff} ≳ 0.6. Systems with 0.25 ≲ χ_{eff} ≲ 0.6 are these in which only one binary component was subject to tidal spinup, while systems with χ_{eff} ≳ 0.6 are those with both binary components being subject to tidal spinup in close WRWR binaries.
Fig. 19.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M30.B with MESA efficient angular momentum transport and fallback decreased BH kicks (no BH kicks for massive BHs). We either do not allow for efficient tidal spinup of WR stars that are the most common immediate progenitors of BHs in our models (natal BH spin is calculated directly from MESA stellar models) or we take it into account (natal BH spin is then calculated as described in Sect. 2.5 if the WR star progenitor was subject to an efficient tidal spinup). For the “no WR tides” approach we find a rather narrow distribution of effective spins (−0.2 ≲ χ_{eff} ≲ 0.2) that is peaked at positive values (average χ_{eff} = 0.15). For efficient “WR tides” the distribution is broad (−0.5 ≲ χ_{eff} ≲ 1.0) with a peak at χ_{eff} ∼ 0.15 (∼73%) and a tail with χ_{eff} ≳ 0.25 (27%). For comparison we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. 
In Fig. 20, we present the effective spin distribution for model M40.B that employs Fuller BH natal spins. We also show two versions of this model: one without efficient WR tides and one with the tides. The distribution is very narrow: −0.1 ≲ χ_{eff} ≲ 0.1 and is peaked at positive values: average χ_{eff} = 0.05. This comes directly from the assumption of very low natal BH spins: a_{spin} = 0.01 (see Eq. (5)). For such a low value of natal BH spins, the effective spin of BHBH mergers is χ_{eff} ∼ 0. Binary accretion onto the secondborn BH spreads effective spins to χ_{eff} ∼ 0.1, while small natal kicks applied to some (lowmass) BHs in this model create a tail extending to low negative values χ_{eff} ∼ −0.1. For efficient “WR tides” the distribution is broad: −0.5 ≲ χ_{eff} ≲ 1.0 with a peak at χ_{eff} ∼ 0.05 (∼78%) and a tail with χ_{eff} ≳ 0.25 (22%).
Fig. 20.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M40.B with Fuller superefficient angular momentum transport and fallback decreased BH kicks (no BH kicks for massive BHs). We either do not allow for the efficient tidal spinup of WR stars that are the most common immediate progenitors of BHs in our models (natal BH spin is calculated directly from MESA stellar models) or we take it into account (natal BH spin is then calculated as described in Sect. 2.5 if the WR star progenitor was subject to an efficient tidal spinup). For the “no WR tides” approach we find very narrow distribution of effective spins (−0.1 ≲ χ_{eff} ≲ 0.1) that is peaked at positive values (average χ_{eff} = 0.05). For efficient “WR tides” the distribution is broad (−0.5 ≲ χ_{eff} ≲ 1.0) with a peak at χ_{eff} ∼ 0.05 (∼78%) and a tail with χ_{eff} ≳ 0.25 (22%). For comparison, we show 90% credible limits (blue arrows) and the most likely values (blue stars) of effective spin parameter for ten LIGO/Virgo BHBH mergers. 
In Figs. 21 and 22, we present the effect of high BH natal kicks on the effective spin distributions in models that employ MESA and Fuller BH natal spins. For MESA BH natal kicks we use two models: M30.B (fallbackdecreased natal kicks: lowtono BH kicks), and M33.B (high natal kicks, independent of BH mass, drawn from a 1D Maxwellian distribution with σ = 265 km s^{−1}). We note that the average effective spin decreases from model M30.B: χ_{eff} = 0.15 to model M33.B: χ_{eff} = 0.04. This is an effect of the natal kicks that tend to misalign BH spins and lower the effective spin. For Fuller BH natal kicks, we use two models: M40.B (fallbackdecreased natal kicks) and M43.B (high natal kicks). We also note that the average effective spin decreases with increasing natal kicks: from model M40.B: χ_{eff} = 0.05 to model M43.B: χ_{eff} = 0.004.
Fig. 21.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M30.B with fallback decreased natal kicks (low natal BH kicks for lowmass BHs and no natal kicks for highmass BHs) and for model M33.B with high natal kicks (all BHs, independent of mass, are subject to natal kicks drawn from a 1D Maxwellian distribution with σ = 265 km s^{−1}). We note that the average effective spin decreases from model M30.B (χ_{eff} = 0.15) to model M33.B (χ_{eff} = 0.04) as an effect of natal kicks that tend to misalign BH spins and lower the effective spin (see Eq. (1)). For comparison, we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. 
Fig. 22.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M40.B with fallback decreased natal kicks (low natal BH kicks for lowmass BHs and no natal kicks for highmass BHs) and for model M43.B with high natal kicks (all BHs, independent of mass, are subject to natal kicks drawn from a 1D Maxwellian distribution with σ = 265 km s^{−1}). We note that the average effective spin decreases from model M40.B (χ_{eff} = 0.05) to model M43.B (χ_{eff} = 0.004) as an effect of the natal kicks that tend to misalign BH spins and lower the effective spin (see Eq. (1)). For comparison, we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. 
4. Discussion
4.1. Angular momentum transport in massive stars
In this section, we discuss the dependence of the final angular momentum of the star at core collapse on angular momentum transport prescriptions and initial (ZAMS) rotation rate.
In order to evaluate the dependence of the final angular momentum on angular momentum transport prescriptions and initial rotation rate, we ran three additional 32 M_{⊙} models at a metallicity of Z = 0.002: a slow initial rotation (V_{ini} = 100 km s^{−1}) nonmagnetic (“noTS”, where TS stands for the TaylerSpruit dynamo) model with both the Geneva and MESA code as well as a slow initial rotation (V_{ini} = 100 km s^{−1}) magnetic (“TS”) model with the MESA code. The specific angular momentum profile of these models at the end of core Heburning are shown in Fig. 23. The helium core, which will form the bulk of the black hole, extends from the center to about 12 M_{⊙} in all the models plotted. For comparison, we also show the specific angular momentum profiles for our two fast initial rotation models (V_{ini}/V_{crit} = 40%: V_{ini} ≈ 340 km s^{−1}) for 32 M_{⊙} star with Z = 0.002 calculated with MESATS and Geneva noTS assumptions. We note that these two fast models are used to calculate BH spin magnitudes employed in our population synthesis calculations (and that are referred to as MESA BH spins and Geneva BH spins).
Fig. 23.
Specific angular momentum profile at the end of core Heburning for the 32 M_{⊙} models at Z = 0.002 calculated with a range of physical ingredients and initial rotation rates: the Geneva fast rotation model (model listed in Table A.1 with V_{ini}/V_{crit} = 40%, V_{ini} = 341.9 km s^{−1}, dotted blue line), the Geneva slow rotation model (V_{ini} = 100 km s^{−1}, solid orange line). Both these Geneva models do not include the TaylerSpruit magnetic dynamo (“noTS”). We also show: the MESA slow rotation nonmagnetic (“noTS”) model (V_{ini} = 100 km s^{−1}, solid green line), the MESA slow rotation magnetic (“TS”) model (V_{ini} = 100 km s^{−1}, dashed red line), and the MESA fast rotation magnetic (“TS”) model (model listed in Table A.2, dotdashed purple line). The brown thin dotdashed line corresponds to the specific angular momentum of the MESA fast rotation model multiplied by a factor of 10 to facilitate comparisons. 
In comparing various models, we observe the following dependencies. First, comparing the magnetic slow and fast rotation MESA models (dashed red and purple lines, respectively), we see that they end with a very similar final angular momentum (a_{spin} = 0.084 and a_{spin} = 0.087 for the slow and fast rotation models, respectively^{5}). There is, thus, very little dependence on the initial velocity in magnetic models. This is due to a weaker magnetic instability (and thus less efficient angular momentum transport) in slower rotation models. We would thus expect slower rotation magnetic models to have a very similar final angular momentum content as the fast models listed in Table A.2. Second, in comparing the magnetic (MESA slow and fast rotation; dashed red and purple lines) and nonmagnetic models (MESA and GENEVA noTS slow rotation and GENEVA fast rotation noTS models; blue, green and orange lines), we see that the nonmagnetic models have final angular momenta that are more than a factor of 10 higher than the magnetic models (the angular momentum profile of the MESA fast rotation TS model multiplied by a factor of 10 is shown to facilitate the comparison; brown line). This means that for the “lower” end of the massive star range, nonmagnetic models predict BHs rotating near or above critical rotation (a_{spin} ≳ 0.9). Third, comparing the slow rotation, nonmagnetic Geneva and MESA models (solid orange and green lines, respectively), the final angular momentum is very similar. The small differences in angular momentum profiles arise from small differences of the core sizes between the two codes (a result of the different treatment of the convective boundary mixing). We thus see that the two codes give consistent results when using a similar treatment of the angular momentum transport (without magnetic fields).
The angular momentum correlations listed above apply directly to black hole spins. Thus, we expect slow and fastrotating magnetic models to end up in low BH spins (a_{spin} ∼ 0.1) and nonmagnetic slow and fastrotating models to result in BHs with large spins (a_{spin} ≳ 0.8). The change of the initial stellar rotation from fast (∼340 km s^{−1}; 40% critical) to slow (∼100 km s^{−1}) does not significantly affect our estimates of the BH natal spin. In Appendix A.2, we present more detailed information and show that even when we increase rotation to ∼680 km s^{−1} (80% critical) our conclusions remain unchanged (see also Fig. A.1). We note that the observed distribution of massive star spins is rather broad (0−600 km s^{−1}), with two distinctive peaks: one (dominant) at ∼100 km s^{−1} and one (small) at ∼400 km s^{−1} (RamírezAgudelo et al. 2013). These correlations were studied for a single mass (32 M_{⊙}: the lower end of the massive star range) at a single metallicity (Z = 0.002) and with two models of angular momentum transport (the TaylerSpruit theory for the magnetic dynamo and meridional currents which dominate angular momentum transport while TS dynamo is switched off). More indepth studies, especially for higher masses, are needed to confirm these correlations.
Our models also make predictions for the birth spin periods of neutron stars and it is worthwhile comparing these predictions with the estimates of these spin periods from observations. Stars with zeroage main sequence masses lying between ∼8−20 M_{⊙} (Fryer et al. 1999, 2012) and more massive stars whose mass loss is sufficiently extreme to dramatically alter the CO core mass collapse to form neutron stars (through normal supernovae). If the CO cores are rotating and the magnetic fields are sufficiently strong, these neutron stars will emit as pulsars. By studying the distribution of spins in these pulsars, we can place constraints on the rotation period of the stars prior to collapse. This spin distribution also could provide some indication of the role angular momentum can play in normal “corecollapse” supernovae. However, in Appendix A.6 we study this problem and arrive to the conclusion that NS spins cannot be directly used to test angular momentum transport in massive stars as various (highly uncertain) mechanisms may spin down or spin up NS after its formation.
4.2. Using single stellar models in binary simulations
Although we use single stellar models to estimate the BH natal spins (see Sect. 2.1), we use binary evolution models to estimate the BHBH merger effective spins (see Sect. 2.6). This is not fully consistent as our approach does not take into account an important effect in binary evolution that may influence the natal BH spin, that is, the spinup of stars (progenitors of BHs) during the RLOF phases. This process typically affects the secondary star of the BHBH progenitors (Figs. 8 and 9; see the second evolutionary stage in both cases). Depending on the donor mass and our assumption on how much mass is accreted during the RLOF (in most models, we allow for 50% of mass lost by the donor to be accreted by its companion; but see models M20.A/B–M26.A/B) the main sequence secondary accretes anywhere from several to several tens of solar masses. This is most likely enough to significantly spinup even a massive star (Packet 1981). However, there is a good chance that at least some of this extra angular momentum is later removed from the secondary star.
When the accretion and spinup end, the secondary star is still on the main sequence. It will be therefore subject to angular momentum transport and losses (through winds) during a significant fraction of its subsequent lifetime. Obviously, there is less time in our binary evolution sequences for secondaries to lose extra angular momentum as compared with single stellar models initiated at ZAMS. On the other hand, there are binary interactions that help the secondary stars to get rid of extra angular momentum after the RLOF/accretion phase.
This happens when the secondary star expands significantly after the RLOF and leaving the main sequence and is subject to tidal spindown as this occurs on a very wide orbit (Figs. 8 and 9: see the third and fourth evolutionary stage in both cases). Then the expansion leads to a CE phase, and the entire secondary’s Hrich envelope (with all its angular momentum content) is removed from the binary (Figs. 8 and 9: see the fifth evolutionary stage in both cases). Any extraangular momentum from accretion during the main sequence that is left in the Hecore of the secondary can contribute to an extra spin of the BH with respect to single stellar models but this would not change our conclusions significantly. If we increased the spin magnitudes of the secondary BHs, then the distributions of BHBH effective spin parameters would become somewhat wider and the effective spins would mostly shift to larger positive values. This would be qualitatively similar to the effect of efficient WR tides on secondary stars (see Figs. 19 and 20) for models with the effective angular momentum transport (MESA and Fuller models). For the model with inefficient angular momentum transport (Geneva model), the BH spins (on average) are already high so increasing stellar spin may only shift this model (in terms of effective spin parameter) even further away from LIGO/Virgo observations.
Additionally, scenario, in which the star is stripped of its envelope, does not generally affect the trend in rotation rates discussed in Sect. 4.1. Indeed, the effective angular momentum transport through a magnetic dynamo only operates when there is a shear allowing the development of a magnetohydrodynamic instability. Hence, different rotation rates between the core and the envelope are crucial since they create the shear. In models, including magnetic fields, the transport of most of the angular momentum occurs during the short phase between core hydrogen depletion and core helium ignition (during this short phase, the core contracts and spins up whereas the envelope expands and slows down, thus creating a strong shear). Therefore, any stripping through binary interaction during core helium burning or afterwards does not affect significantly the slowdown of the cores. Stripping during the MS, however, is more complicated. There are two opposite effects. The first is that mass loss removes also angular momentum. The second is that if the entire hydrogenrich layer is removed, then there is no envelope left to expand at the end of the MS and the spindown of the core by internal angular momentum transport (coreenvelope coupling) is reduced. In nonmagnetic models, spindown of the core is weak anyway and models still end up with relatively large spin values, a_{spin} : 0.25 to maximum spin (see Table A.1). In magnetic models, removal of the hydrogenrich outer layer significantly reduces the spindown of the core in the contraction phase after the MS and the core may keep its fast rotation rate (see also Fuller & Ma 2019).
The majority of binaries that produce NSNS/BHNS/BHBH mergers in our models have large initial orbital separations (≳1000 R_{⊙}; de Mink & Belczynski 2015). Therefore, binary interactions (CE or stable RLOF; see Figs. 8 and 9) that remove Hrich stellar envelopes happen either at the end of HG or during CHeB so that the values of the natal compact body spin are not significantly affected.
Therefore, one may safely conclude that neither accretion spinup of stars nor increasing the initial stellar rotation rate (see Sect. 4.1) can significantly modify our conclusions about the natal spins of merging binary compact objects resulting from isolated evolution. However,we do not include quasihomogeneous evolution scenarios discussed by Chrimes et al. (2020).
We note that other processes that can affect BH spin in binary evolution are taken into account in our calculations. Besides tidal interactions that we probe with recent models presented in literature (see Sect. 2.5), we also account for accretion onto BH that may lead to BH spinup. The spinup is calculated using the formalism presented by Belczynski et al. (2008a). There are two potential phases of accretion onto the firstformed BH in our major evolutionary scenario: during the CE phase from the Hrich envelope of the secondary and after the CE phase from the Herich wind of the WR secondary (Figs. 8 and 9; see the fifth and sixth evolutionary stage in both cases). Since accretion leads to (usually very small or only modest) increase of the first BH spin magnitude, it tends to slightly broaden the effective spin parameter distribution.
4.3. BHBH/BHNS/NSNS merger rate densities
The predicted merger rate density of coalescing double compact object binaries depends on the particular formation scenario which is assumed. In the following two sections we discuss our merger rates obtained for isolated binary evolution and contrast them with BHBH merger rates obtained from dynamical evolution in globular clusters. In summary, the typical globular cluster rates (5 Gpc^{−3} yr^{−1}) can explain about 10% of the LIGO/Virgo BHBH mergers (the empirical rate peaks at ∼50 Gpc^{−3} yr^{−1}) while the remaining 90% of the observed rate can be easily explained by isolated binary BHBH formation (see Table 3).
4.3.1. Isolated binary evolution
In this paper, we present a range of models and their rate densities. The rate density in each case is the result of several assumptions: starting from the cosmic star formation, metallicity evolution and continuing through the initial binary parameters and the binary evolution model up to the implied delay time (between the birth of a binary and the final merger of two compact objects) distribution. For each model, we calculated the local merger rate density (z = 0) as well as the predicted LIGO/Virgo detection rate in O3; see Table 3. The local rate densities of BHBH and NSNS mergers that we obtained are directly comparable with the rate limits inferred by LIGO/Virgo at the conclusion of the O1/O2 observing runs; see Fig. 24. Additionally, based on the recent LIGO/Virgo candidate (S190814bv) for a very likely (∼99% in the reported massbased source classification) BHNS merger from August 14 (LVC 2019a,b), we estimate the BHNS merger rate density of ℛ = 1.6−60 Gpc^{−3} yr^{−1} (see Appendix A.9 for more), which can also be confronted with our models.
Fig. 24.
Comparison of the local merger rate densities of BHBH and NSNS mergers from all our models (see Table 3) with the current limits inferred from the O1/O2 LIGO/Virgo observational runs: 9.7−101 Gpc^{−3} yr^{−1} for the BHBH mergers and 110−3840 Gpc^{−3} yr^{−1} for the NSNS events (Abbott et al. 2019b). Models that are consistent with the observational limits: M13.A, M23.A, M33.A, M43.A, M30.B, M40.B, M60.B, M70.B. We note that while some of the models are consistent with the observational constraints and others are not, it is not straightforward to draw conclusions about the physical ingredients of the models at this stage due to degeneracies in the impact of various assumptions on the theoretical merger rates (see Sect. 4.3.1). 
In addition to the overall rate of events, the relative rates of the three generic categories of events (BHBH, BHNS, NSNS) provide particularly valuable constraints. With this in mind, in Figs. 24 and 25, we show the rates of all three source categories. The horizontal and vertical bands represent the allowed rates consistent with current LIGO/Virgo data. Superimposed on the LIGO/Virgo bands are the rate predictions from a range of our models, with the red diamonds from submodel A (the optimistic approach to common envelope, in which CE events with Hertzsprung gap donors are allowed to survive), and the blue diamonds from submodel B (the pessimistic approach to common envelope, in which CE events with Hertzsprung gap donors are not allowed to survive). We note that in both figures, there are many models that fall in the central sweet spot of the figure.
Fig. 25.
Comparison of the local merger rate densities of BHBH and BHNS mergers from all our models (see Table 3) with the current limits: 9.7−101 Gpc^{−3} yr^{−1} for the BHBH mergers (LIGO/Virgo O1/O2) and 1.6−60 Gpc^{−3} yr^{−1} for the BHNS events (based on the LIGO/Virgo O3 candidate of the first BHNS system, S190814bv, LVC 2019a,b; see Appendix A.9 for details). Models that are consistent with the observational limits: M13.A, M23.A, M33.A, M43.A, M23.B, M25.B, M30.B, M40.B, M50.B, M60.B, M70.B. We also note that the BHNS merger rates from our models span the range between 0.48 and 297 Gpc^{−3} yr^{−1}, which is consistent with the upper limit determined by their nondetection in O1/O2 LIGO/Virgo runs (< 610 Gpc^{−3} yr^{−1}, Abbott et al. 2019b). 
Looking at Fig. 24, we see that the existing O1/O2 run bounds on the NSNS and BHBH merger rate densities already allow excluding some of the models. In the standard model group with the 2016 assumptions the only model consistent with the data is M13.A, which requires high BH/NS natal kicks. If we allow some modified physics, we can also add model M23.A. From the group of models with updated physics of 2019, we can select two: M30.B ad M33.A, that are consistent with observations. Allowing for models with some modifications we can also include the models M40.B and M43.A, as well as M60.B and M70.B. Thus, we can see that the 2016 standard model is consistent with the merger rate density O1/O2 runs data provided that the natal kicks are large. This may be in contradiction to some observations of BH in binaries that indicate small BH natal kicks, however, such a conclusion requires more detailed studies. The updated model of 2019 is consistent with the data. It also allows for certain modifications such as probing the origin of BH spins or the inclusion of various models of PPSNs.
The inclusion of the BHNS merger rate estimate (based on the O1/O2 runs and ongoing O3 run^{6}) (see Appendix A.9) in Fig. 25 offers additional insights. Of particular interest are models which fall in the central overlap region of both Figs. 24 and 25. For example, M13.A, M23.A, M33.A, M43.A, M30.B, M40.B, M60.B, and M70.B fall within the observed range for all three rates. A number of trends are evident from looking at both figures together. For example, we note that the B class of submodels, without the inclusion of HG donors in the CE evolution, brings many of the models into the center of the BHBH/BHNS plot. However, this same change brings some of the models into tension or only into marginal agreement with the data on the BHBH/NSNS plot. In other words, submodel A is almost entirely excluded in Fig. 25, while submodel B is almost entirely excluded in Fig. 24. This might be a preliminary indication that the NSNS systems do allow for CE with HG donors, while BHNS systems do not. However, there are other factors that are relevant for the relative merger rates of different source types. For example, improving the metallicity evolution (which shifts the net star formation across the cosmic history towards higher metallicities) appears to shift down the BHBH merger rates while leaving the NSNS merger rates relatively unchanged (see also Chruslinska et al. 2019).
In all cases that we consider, the merger rate density increases with increasing redshift, peaks at about z ≈ 2 and then decreases. A rough approximation of the scaling of the merger rate density with redshift in the range of z ≲ 2 can be approximated as ∝(1 + z)^{1.5}; however, the exponent is not wellconstrained. It is quite interesting to note that the assumed shape of the star formation rate history within the allowed bounds does not strongly influence the shape of the dependence of the merger rate density on redshift within z ≲ 2; see Fig. 10.
We expect that most of the remaining ambiguities should be resolved with the release of analysis of the O3 data. Assuming that these observations will yield about a hundred detections, we expect the bounds on the rates to narrow by a factor of three. This should allow us to strongly constrain the models presented in this paper based on the event rates alone. It should be stressed, however, that even then it will not be straightforward to constrain the underlying physics. Variations in different ingredients of a model can have a degenerate effect on the resulting merger rates. For instance, the inclusion of HG donors in the CE evolution (going from submodels B to A) leads to an increase in the merger rates, but so does the lowering of BH/NS natal kicks (see Sect. 3.2 and Fig. 11 for details). The CE treatment and the natal kicks are only two of the many ingredients of our models, each of those ingredients being to some extent uncertain. This introduces a degeneracy that cannot easily be resolved at the moment. In the future, we can hope that the higher number of detected mergers allow us to combine all the inferred observables (e.g., rates, masses, spins, redshifts, host galaxies) in order to break the degeneracies and put better constraints on the model’s input physics.
That said, submodels A: M13.A, M33.A, M43.A appear to be consistent with all rate estimates (although on the high end of the BHNS limits). This shows the effect of the high natal kicks that reduce high BHBH merger rates typically found in the submodels A. At the same time, these high natal kicks do no affect the NSNS merger rates significantly, as in our standard approach (fallback decreased kicks) NS natal kicks are already high. Submodels B: M30.B, M40.B, M60.B and M70.B, that appear consistent with all rate estimates, indicate that high average cosmic metallicity in unison with exclusion of CE events with HG donors may reduce high BHBH merger rates in comparison with older models and with submodels A. We note that these models include our current standard approach to input physics, with the two efficient angular momentum transport mechanisms (MESA and Fuller models) and the three prescriptions for PPSN (weak, moderate, and strong) mass loss. We note that some models are very close to being consistent with all rate constraints. Model M50.B is one such example. In this model, we reduce stellar winds to 30% of their standard values. Typically, BHs form with higher mass in this model. Additionally, with such low wind massloss rates for massive stars it is possible to form ∼60 M_{⊙} BHs in a solar metallicity environment (Z = 0.02) while still managing to avoid PPSN mass loss (see Belczynski et al. 2020 for discussion). At such massloss rates single stars (or stars in wide noninteracting binaries) may form ∼40 M_{⊙} He cores with ∼20−30 M_{⊙} Hrich envelopes leading to massive BH formation even at high metallicity. Such massive BHs would be found either as single BHs through microlensing surveys (e.g., Wyrzykowski & Mandel 2020) or in wide binaries through radial velocity surveys (motion of the companion star). We note that such massive BHs would not necessarily add to either LIGO/Virgo or Xray binary populations that form (at least in classical isolated binary evolution) from interacting binaries in which stars lose their Hrich envelopes (e.g., this study, Wiktorowicz et al. 2014). There is, of course, a complex interdependence between the various aspects of these models and this subset of models does not fully cover the parameter space. Nonetheless, these general trends may provide important clues about the underlying stellar and binary physics.
4.3.2. Dynamical evolution in globular clusters
In dense stellar environments, even binaries that are initially too wide to inspiral via gravitational radiation in a Hubble time can be hardened and induced to merge by binarysingle and binarybinary interactions. Moreover, single BHs can segregate to the cluster center due to dynamical friction and form close binary systems that can merge due to gravitational radiation. Thus dense stellar systems such as globular clusters can be highly efficient per stellar mass at producing BHBH mergers. However, only a small fraction of stars are in globular clusters or similarly dense systems (∼10^{−4} − 10^{−3}, depending on the type of galaxy). This turns out to mean that the plausible rate density of mergers in dense stellar systems is low. A strong upper bound can be obtained (see Mandel & Farmer 2018) by noting that there is roughly one globular cluster per Mpc^{3} in the local Universe. If each globular has ∼10^{6} stars (the actual average is a few times lower than this), ∼10^{−3} of stars become black holes, and all black holes pair up and merge in a Hubble time of ∼10^{10} years, then at most the rate will be ∼10^{6} × 10^{−3}/(10^{10} yr × 1 Mpc^{3}), or ∼100 Gpc^{−3} yr^{−1}. In reality, the expected rate from clusters is ten to 100 times lower, for a discussion see Mandel & Farmer (2018).
It does seem likely that globular cluster masses at formation were a few times larger than they are today. For example, it has been estimated that ∼60%−70% of cluster mass is lost in a Hubble time due to dynamical evolution and interactions with the host galaxy (Chatterjee et al. 2010; Giersz et al. 2013; Webb & Leigh 2015). However, recent realistic estimates of BHBH merger rate densities from globular clusters are typically few Gpc^{−3} yr^{−1} (Rodriguez et al. 2016c; Askar et al. 2017; Park et al. 2017; Hong et al. 2018; Choksi et al. 2019). Recoil during binarysingle and binarybinary interactions can be sufficient to eject a binary from its host globular, with the result that the majority of such mergers are expected to occur outside the globular (Portegies Zwart & McMillan 2000 and many subsequent papers).
4.4. BH masses
An important test for population synthesis models is the direct comparison of their predictions of the mass distribution of binary black holes with the observed distribution from gravitational wave detections. The masses of black holes are among the best measured quantities of the LIGO/Virgo sources. In particular, the chirp mass is often measured to high accuracy (O(5%) or better) for lower mass mergers with a significant number of inspiral cycles detected in the LIGO/Virgo sensitivity band. For the higher mass systems, the total mass can be measured with reasonable accuracy. The individual component masses are often more poorly measured, and equivalently, the mass ratio is often relatively poorly constrained. For this discussion, we focus on the shape of the mass distribution.
As LIGO/Virgo continues to add to the sample of binary black holes, the inferred underlying mass distributions become increasingly wellconstrained. Some preliminary limits have already been presented in Abbott et al. (2019a). In particular, that paper employs a number of different fits to the mass distribution of the more massive components of the detected binaries. For the purposes of comparison, we use the Model B from this paper, which consists of a powerlaw in M_{BH1} with upper and lower mass cutoffs and a powerlaw in the binary mass ratio. The mass cutoff accounts for the dearth of highmass black holes (Fishbach & Holz 2017), as would be expected from PSN and PPSN (Belczynski et al. 2016c). LIGO/Virgo finds a power law index in the range of α = −3.1 to α = 0.1 with peak probability at α = −1.6^{7}, and with an upper mass cutoff between 36.3 M_{⊙} and 57 M_{⊙}.
In Fig. 26, we show the intrinsic rate density of mergers averaged out to z = 0.5 versus the mass of the more massive BHs in BHBH mergers (M_{BH1}; upper figure) and the mass ratio q = M_{BH2}/M_{BH1} (lower figure) for several of our models and compare them with the fits to Model B from Abbott et al. (2019a). Although Model B does not have the freedom to capture some of the features of the population synthesis distributions, the overall rate densities of the mass distribution found in the LIGO/Virgo collaboration analysis (shown in gray) are quite similar to those obtained, for example, in models M30.B, M40.B, and M50.B. We emphasize that the shapes of the primary mass distributions we obtain in these StarTrack models are better fit to exponentials than to powerlaws and that there are some interesting features in the mass distributions at low BH mass. Nevertheless, compared to our theoretical models, the LIGO/Virgo Model B fits to the 10 O1/O2 BHBHs exhibit a similar falloff of the primary mass distribution towards higher masses, lower and upper mass cutoffs, and a falloff of the mass ratio distribution at more unequal masses.
Fig. 26.
Intrinsic merger rate distribution as a function of primary black hole mass M_{BH1} (top) and binary mass ratio M_{BH2}/M_{BH1} of BHBH mergers (bottom) for different StarTrack model choices. The gray bands show the 90% confidence limits on the populations inferred in Abbott et al. (2019a) using the phenomenological Model B. 
At this point, we do not attempt any more elaborate fits to the O1/O2 data than the comparison to the LIGO/Virgo phenomenological models, as these can wait till more LIGO/Virgo data points become available. For anyone interested in such fits, all our data is available online. With O(100) BHBH mergers expected for the entirety of O3, we will be able to test the existence of some of the finer features seen in our mass distributions as well as the goodnessoffit of exponentials versus powerlaws to the primary masses.
4.5. BHBH effective spins
The rate of BHBH mergers with different effective spins provides several reliable features that we can use to corroborate and constrain our models with present and future GW observations.
First and foremost, as has been discussed in our and other previous work, our binary evolution models only rarely, if ever, produce binaries with χ_{eff} significantly below zero (see Figs. 18–21). This void provides an opportunity to identify the unique contribution from alternative formation channels as our models cannot produce a preponderance of events with negative effective spins.
Second, our predicted effective spin distributions have sharp cutofflike features in each of the one (or sometimes two) subpopulations that dominate the detection rate versus effective spin. For example, the black hole spin distribution drops sharply for χ_{eff} ≳ 0.9 in the inefficient Geneva angular momentum transport mechanism (Fig. 18), above 0.2 for the binary black holes not spun up by tides in the MESA models; above 0.6 for the subpopulation of BHBH with a tidallyspunup WR progenitor in the MESA models (Fig. 19); and near 0.1 for the Fuller models (Fig. 20). Sharp cutoffs like these are rapidly identified empirically, allowing future gravitationalwave observations to constrain the BH natal spin distribution, and the physics of tidal spinup.
Third, for models with WR tides (e.g., M30.B and M40.B shown in Figs. 19 and 20), our effective spin distribution has two wellseparated subpopulations (a peak and a tail), associated with the principal channel and the WRspinup channel. In other words, the gravitationalwave population may allow us to reliably associate specific binary physics to specific binary black hole mergers, and thus to measure the relative proportions (and properties) of BHBH systems forming through each channel.
While our binary evolution models have many parameters, because of the limited role of accretion on the BH spin, very few parameters have more impact on the χ_{eff} distribution than the physics we have described above. For example, we show with Fig. 22 that BH kicks have a relatively modest effect on the shape of the χ_{eff} distribution. Strong BH kicks can for example make the χ_{eff} distribution more symmetric and isotropic, effectively by disrupting binaries which would otherwise dominate the sharp cutoffs at the largest values of χ_{eff}. In other words, these kicks make the χ_{eff} merger rate smaller and with a less prominent (and therefore more difficult to identify) cutoff. Although this poorly constrained physics of natal kicks can reduce the rate at which χ_{eff} measurements can inform about the binary evolution and diminish the contrast of the features we have described above, they cannot erase them, particularly with regard to features such as the second population of WR tidally spunup stars or the low probability of mergers with large negative χ_{eff}.
In Fig. 27, we show the cumulative density functions (CDFs) of the effective spin χ_{eff} for a select set of our models and compare them to the χ_{eff} measured for the 10 O1/O2 BHBH mergers. For each model we show the χ_{eff} CDF of detectable events as a color solid line. In the corresponding colored shaded regions, we show the 90% range of 5000 mock observed CDFs under each model. To generate a mock observed CDF, we draw 10 χ_{eff} values from the detected population under a given model and add random Gaussian noise with σ = 0.05, which is approximately the uncertainty in χ_{eff} likelihoods from the events of GWTC1 (we resample any samples with χ_{eff} > 1 after adding noise). The black line shows the χ_{eff} CDF of median likelihoods of the O1/O2 LVC events, and the gray region bounds CDFs of the 5th and 95th percentiles of the likelihoods from LVC parameter estimation. Notably, these effective spins are clustered around zero effective spin, indicating that model M10.B’s predicted spins are disfavored by the data. Models M30.B (no WR tides), M30.B (WR tides), and M40.B are much more consistent with these observations, and with the LIGO/Virgo collaboration analysis in Abbott et al. (2019a), which finds a preference for most effective spins to be very close to zero (see Fig. 12 in Abbott et al. 2019a). That analysis finds that only a few tens of percent of the observed BHBH mergers have χ_{eff} > 0.05, and of those with χ_{eff} > 0.05, almost all are likely to have positive χ_{eff}. For the moment the observed χ_{eff} distribution does not allow to test the importance of WR star tides during the evolution towards a BH even if it seems to slightly favor the notide path. Overall, the limited BHBH sample from O1/O2 does not strongly constrain the BH spins other than to the conclusion that the effective spins tend to be near zero. However, it suggests that O3 should provide data that will allow selecting evolutionary path that lead to BHBH mergers.
Fig. 27.
Cumulative density functions (CDFs) of four StarTrack models versus effective spins χ_{eff}. Solid color lines correspond to χ_{eff} CDF of detectable events for each model. The black line shows the χ_{eff} CDF of median likelihoods of the O1/O2 LVC events, and the gray region bounds CDFs of the 5th and 95th percentiles of the likelihoods from LVC parameter estimation. In the corresponding colored shaded regions, we show the 90% range of 5000 mock observed CDFs under each model. To generate a mock observed CDF, we draw 10 χ_{eff} values from the detected population under a given model and add random Gaussian noise with σ = 0.05, which is approximately the uncertainty in χ_{eff} likelihoods from the events of GWTC1 (we resample any samples with χ_{eff} > 1 after adding noise). 
4.6. Dichotomy of LIGO/Virgo BHs and HMXB BHs
LIGO/Virgo BHs are mostly massive (see Table 1) and it seems that the spins of these BHs are small. Unfortunately, the spin magnitudes of massive BHs with M_{BH} > 15 M_{⊙} are not constrained by other observational means. For example, the three most massive BHs in windfed high mass Xray binaries (HMXB) for which we have spin estimates are at roughly 1/3 or 1/2 the mass of the largermass LIGO/Virgo BHs. The estimated masses and spins of these BHs can be found online at^{8}: LMC X1, M_{BH} = 10.9 ± 1.6 M_{⊙} (a_{spin} = 0.92); Cyg X1, M_{BH} = 14.8 ± 0.1 M_{⊙} (a_{spin} > 0.983); M33 X7, M_{BH} = 15.7 ± 1.5 M_{⊙} (a_{spin} = 0.84). We note that BH spins can be measured by two methods: disk reflection and disk continuum. For LMC1 and Cyg X1 BH spins from both methods are consistent, while for M33 X7 BH spin was estimated only through disk continuum (Miller & Miller 2015). There is an apparent tension between the spin estimates of LIGO/Virgo BHs and BHs in these HMXBs. However, this tension can be possibly avoided if HMXBs and LIGO/Virgo systems form through different evolutionary scenarios.
First, none of these three HMXBs is expected to produce a BHBH merger. Future evolution of these systems was studied and it was shown that none of these systems will form a BHBH merger. They will either end up as single objects (in which BH merges with its massive companion star), or at best they may form BHNS systems (Belczynski et al. 2012b). We note that this conclusion was reached in the framework of the classical binary evolution that does not take into account homogeneous evolution of rapidly rotating stars that allows for the formation of BHBH mergers in alternative ways (Marchant et al. 2016; Mandel & de Mink 2016; de Mink & Mandel 2016).
Second, it was argued that HMXBs and LIGO/Virgo systems may form through different evolutionary scenarios. As explained in the past, the classical isolated binary evolution channel forms BHBH mergers from (initially) very wide binaries (a ≳ 1000 R_{⊙}; see Figs. 8 or 9 or de Mink & Belczynski 2015; Belczynski et al. 2016a). Therefore, if the tidal spinup operates, it only acts at the very end of the evolution during the BHWR stage (potentially) spinning up some fraction of the secondborn BHs in BHBH mergers (see Sects. 2.5 and 3.4). This makes the majority of BHs in BHBH mergers to have low spins if efficient angular momentum transport is adopted for massive stars (as in our MESA or Fuller models). On the other hand, HMXBs were proposed to form from initially very close binaries (a ∼ 100 R_{⊙}, e.g., Valsecchi et al. 2010; Qin et al. 2019). In such scenario a system with a massive main sequence primary is tidally locked so, through its evolution, its spin is kept high at the expense of the orbital angular momentum. At such small orbital separations the binary undergoes case A mass transfer or homogeneous evolution that keeps reducing Hrich envelope of the donor primary and keeps it from expanding. After main sequence evolution, the primary is a compact WR star that quickly collapses to a BH. The BH is found in the near proximity of its companion star (if natal kick is small), naturally producing a (windfed) HMXB. Possibly, tidal interactions can keep the primary spin high until the end of its nuclear evolution which would lead to the formation of a rapidly spinning BH. However, the detailed evolutionary calculations have shown that if efficient (TaylerSpruit dynamo) angular momentum transport is assumed for the primary star, then properties of the three HMXBs considered here cannot be reproduced (Qin et al. 2019). Therefore, the tension still exists, albeit it seems like HMXBs and LIGO/Virgo BHBH mergers originate from different initial populations of binaries (close versus wide) but specific details of evolution leading to the formation of HMXBs need to be worked out.
We note that formation scenarios that require rapidly spinning stars (whether they lead to formation of BHBH mergers or BH HMXBs) are not available in our study. Rapid rotation induces efficient mixing and reduces radial expansion of stars in homogeneous or semihomogeneous evolution. Our calculations can only be applied to slow or moderatelyrotating stars in classical binary evolution. Both binary channels (homogeneous and classical) do not exclude one another, but rather they both operate in Universe. Finally, we should stress that estimates of BH spins through properties of their accretion disks are not direct spin measurements but they are modeldependent, so there is still a possibility that the tension is only apparent.
4.7. Origin of LIGO/Virgo BHs
At the moment (only O1/O2: data available), it seems that the classical isolated binary evolution formation channel is capable of explaining all the basic properties of the LIGO/Virgo BHBH mergers. This does not mean than other channels do not contribute to LIGO/Virgo detections. For example, if we compare the average rate estimates of the isolated binary channels (see Table 3) and the dynamical formation channels in globular clusters (see Sect. 4.3.2), it appears that 1 in 10 BHBH mergers might come from a globular cluster. There may be yet another dynamical contribution from young open clusters to BHBH formations (Ziosi et al. 2014; Di Carlo et al. 2019). Even the homogeneous isolated binary channel is not excluded by the observed LIGO/Virgo low BHBH effective spins, as apparently BHs with a broad range of spins (a_{spin} ∼ 0.2−1) can be produced also in this scenario (Ph. Podsiadlowski, priv. comm.). In the published models with rather sparse metallicity sampling and conservative evolutionary assumptions (e.g., no pairinstability pulsation supernova mass loss^{9}), BH spins are found in somewhat narrower range: a_{spin} ∼ 0.4−1 for M_{BH} ≲ 60 M_{⊙} (e.g., see Fig. 9 of Marchant et al. 2016).
If a BHBH merger is discovered with either of the binary components in the mass range 70/80 ≲ M_{BH} ≲ 135 M_{⊙} with a high spin of a_{spin} ≈ 0.7, this would strongly suggest a dynamical formation scenario. The lowermass limit corresponds to the pairinstability pulsation supernova effects on the presupernova star and on the remnant BH mass (Woosley 2017; Limongi & Chieffi 2018; Belczynski et al. 2020), while the upper mass limit corresponds to the end of the pairinstability supernova process, that is believed to disrupt the entire star without BH formation (Fryer et al. 2001; Heger & Woosley 2002). The lack of BHs in this mass range is referred to as the “second mass gap” (Belczynski et al. 2014; Spera et al. 2015; Marchant et al. 2016; Fishbach & Holz 2017). It seems unlikely that isolated binaries can fill this gap, but repeated BHBH mergers in dense environments could produce such heavy BHs, and it is expected that they would have moderately high spins a_{spin} ≈ 0.7 (Gerosa & Berti 2017; Fishbach et al. 2017). However, it cannot be excluded that a BH with high mass and low spin is formed by a merger of two BHs (Belczynski & Banerjee 2020). Alternatively, the detection of a lowspin BH with mass within the second mass gap may point to either (i) inconsistencies in pairinstability supernova theory or (ii) a primordial BH origin (Green 2017). There is also an issue of the second mass gap width. The lower bound may be as high as 70−80 M_{⊙} depending on details of input physics in stellar models (Limongi & Chieffi 2018; Belczynski et al. 2020). The upper bound may change if fusion reaction rates of heavy elements that are involved in pairinstability supernovae change. We note that the reaction rates are uncertain (Fields et al. 2018). Possibly, the second mass gap is narrower than currently believed.
In the classical isolated binary evolution channel, we show that if LIGO/Virgo observations of BHBH mergers remain consistent with small effective spins χ_{eff} ≃ 0, this would indicate that low natal spins are common in BHs. In this case, our work shows that stellar models with mild rotational coupling between the stellar interior (core) and the outer zones (envelope) would be disfavored. Our main finding is that angular momentum transport in massive stars (so far unconstrained by the electromagnetic observations) is more efficient than predicted by the Geneva shellular model. This demonstrates how LIGO/Virgo observations can be used to make clear astrophysical inferences and guide stellar evolution astrophysics. Our conclusion is not subject to the main known population synthesis uncertainties, since we have explored a broad range of key parameters: BH natal kicks, initial star rotation, tides, spinup by accretion onto BHs.
We note that adjusting the angular momentum transport to produce low spinning BHs in BHBH mergers, to fit LIGO/Virgo observations leads to several astrophysical inferences and constraints on massive binary evolution. (i) RLOF between two massive stars (see second evolutionary stage in BHBH formation; Fig. 8 or 9) cannot effectively spin up the core of the accreting MS star (since this would produce a BH with large spin). This can be avoided if the RLOF is highly nonconservative and not much mass is accreted onto MS star (this is not excluded by any EM observations). We note that we have assumed that 50−80% of the mass is lost in RLOF, but in our models this fraction can be easily increased. Alternatively, if accretion is significant and if the entire star is spunup; then we obtain an additional constraint on the angular momentum transport. It needs to be effective enough to remove most of the angular momentum from the highly spinning core to the envelope in about 1 Myr (time between RLOF and CE which expels the envelope from the binary). (ii) Spinup of the firstformed BH by accretion in the CE phase must be negligible as predicted by recent studies (MacLeod et al. 2017; MurguiaBerthier et al. 2017; Holgado et al. 2018, see the fifth evolutionary stage in BHBH formation; Figs. 8 and 9). (iii) Tidal torques are not effective in BHWR binaries that form BHBH mergers as the WR star spinup would lead to the formation of a highly spinning BH (see the sixth evolutionary stage in BHBH formation; Figs. 8 and 9). In our evolution, a small, but significant fraction (∼20−30%; see Figs. 19 and 20) of the BHWR/WRBH/WRWR binaries are found on orbits smaller than ∼10−20 R_{⊙} (orbital periods P_{orb} < 1.3 d) for which the tides are expected to be effective for WR stars (Kushnir et al. 2016; Zaldarriaga et al. 2018; Qin et al. 2018). Therefore if there are no detections of highly spinning BHs in LIGO/Virgo observations it will imply that tides are not as effective as argued by recent work or that BHBH mergers are not produced by the isolated classical binary evolution (Hotokezaka & Piran 2017). We note, however, that already one event so far (GW170729) may have an effective spin as high as χ_{eff} = 0.57 (90% credible limits; see Table 1). On the other hand, among ten events, it can be expected that one outlier will appear within the 90% credible limits. Therefore, at the moment the comparison of models with observations remains inconclusive with respect to tides. Yet, the comparison of models predicting tidally spunup stars, producing highly spinning BHs, with the number of LIGO/Virgo high effective spin BHBH mergers may help to constrain the strength of tides in near future.
Finally, we note that our main conclusion, which states that angular momentum transport needs to be more effective than predicted by shellular model for massive stars, lends support to a similar conclusion that was reached for lowmass stars based on Kepler asteroseismology data (Cantiello et al. 2014).
5. Conclusions
In this work, we updated our method of population synthesis calculations with revised and extended input physics that is important for the formation of double compact object mergers: BHBH, BHNS, and NSNS. We introduced new models of pairinstability pulsation supernovae that are crucial in mass estimates of heavy black holes. We employed the recent estimates of starformation rate density and cosmic metallicity evolution in our calculations. These two factors play an important role in setting the double compact object merger rates (Belczynski et al. 2010b; Dominik et al. 2013; Chruslinska et al. 2019; Neijssel et al. 2019). We introduced new models for BH natal spin that allow for a direct comparison with LIGO/Virgo estimates of the effective spin parameter of BHBH mergers. These models connect gravitationalwave observations with detailed stellar evolution calculations of angular momentum transport in massive stars that is otherwise hidden from electromagnetic observations. We also allowed for significantly decreased stellar winds with respect to the standard wind prescriptions (Vink et al. 2001) used in modeling. Winds are an important factor that sets the shape of black hole mass distribution. Finally, updated prescriptions of accretion onto compact objects in stable Rochelobe overflow, common envelope, and from stellar winds were used in our calculations.
We summarize our main findings in the following points.

1.
Our study is the first to employ detailed single stellar evolutionary calculations in binary population synthesis to compare LIGO/Virgo BHBH effective spins with those resulting from several angular momentum transport mechanisms. If LIGO/Virgo BHBH mergers originate from the classical isolated binary evolution channel, then their low effective spins inform about the effective angular momentum transport in massive stars and most likely disfavor effective tidal interactions in close binaries with WR stars. According to our evolutionary framework, the TaylerSpruit magnetic dynamo (as implemented in the MESA or in the Fuller model) reproduces the effective spin measurements very well, while meridional currents (as implemented in the Geneva model) are inconsistent with the current LIGO/Virgo data.

2.
For some of our models, the predicted merger rate densities of BHBH, BHNS, and NSNS systems are within the LIGO/Virgo empirical estimates. The rates strongly depend on the cosmic metallicity evolution, the choice of NS/BH natal kicks and the treatment of the common envelope phase. Due to the similar effects of these factors on the rates, it still is not possible to derive definitive conclusions about any of these pieces of input physics individually. However, some combinations of parameters may be already excluded. Interparameter degeneracies make up an important factor that cannot be overlooked in deriving astrophysical conclusions from gravitational wave observations.

3.
The range of the observed BH masses appears to be in agreement with all three of our adopted PPSN models: from strong PPSN with maximum BH mass of M_{BH, max} ∼ 40 M_{⊙}, to moderate PPSN with M_{BH, max} ∼ 50 M_{⊙}, and to weak PPSN with M_{BH, max} ∼ 55 M_{⊙}. If heavier BHs are observed it will indicate that, either the pairinstability pulsation supernovae and the pairinstability supernovae do not work as predicted, or that BHs with very heavy mass originate from other formation channels. The values of BH spins may distinguish these two possibilities. If a massive BH (∼100 M_{⊙}) is observed with low spin it will indicate a classical isolated binary evolution formation channel, but such an observation will cast shadow on PPSN/PSN predictions. If, however, such a BH is found to have a large spin as expected from consecutive BH mergers, then its presence will point to formation in a dense stellar environment (e.g., globular cluster).

4.
Based on our limited set of models and the initial ten O1/O2 detections, if we were to make a set of statements on the physics of massive binaries, we could venture to say that with our revised cosmic metallicity evolution, it seems that to simultaneously reproduce the observed BHBH, BHNS, and NSNS merger rates one should use submodels B (no CE allowed with HG donors) if BH kicks are low, but if these kicks are high, submodels A (CE allowed with HG donors) should be preferred. If individual BH masses are in fact as high as LIGO/Virgo reports for the most likely values from O1/O2 run (e.g., 50.6 M_{⊙} for GW170729) then we can already exclude strong mass ejection during the PPSN.
We note that similar conclusions have been reached by Bavera et al. (2020) and Neijssel et al. (2019), although these authors used somewhat different approaches to their detailed stellar evolution models and a different population synthesis code was used.
GENEC and MESA use two slightly different definitions of the critical velocity. Maeder & Meynet (2000) discusses the two definitions. MESA models use their Eq. (3.12), whereas GENEC models use the minimum between the values obtained from their Eqs. (3.14) and (3.19), respectively.
We note that the LIGO/Virgo power law index is defined with a minus sign in front (see Eq. (2) of Abbott et al. 2019a).
We note that pairinstability pulsation supernova mass loss, if taken into account, can reduce spin by about 30% (Marchant et al. 2019).
Acknowledgments
We have benefited from comments from Selma de Mink, Enrico RamirezRuiz, Ilya Mandel, Thomas Janka, Serena Repetto, Simon Stevenson, Sambaran Banerjee, Tom Maccarone, Craig Heinke, Phil Charles, Doron Kushnir, Tsvi Piran, Miguel Holgado, Antonio Claret, Tassos Fragos, Philipp Podsiadlowski, and Matt Benacquista. We would like to thank thousands of Universe@home users that have provided their personal computers and phones for our simulations, and in particular to Krzysztof Piszczek (program IT managers). KB acknowledges support from the Polish National Science Center (NCN) grant Maestro (2018/30/A/ST9/00050) and KB, SM and JPL from the NCN grant OPUS (2015/19/B/ST9/01099). JK and MC acknowledge support from the Netherlands Organisation for Scientific Research (NWO). C.E.F. acknowledges support from a Predoctoral Fellowship administered by the National Academies of Sciences, Engineering, and Medicine on behalf of the Ford Foundation, an Edward J Petry Graduate Fellowship from Michigan State University, and the National Science Foundation Graduate Research Fellowship Program under grant number DGE1424871. MG ans AA were partially supported by the National Science Centre, Poland, through the grant UMO2016/23/B/ST9/02732. This work has been supported by the EU COST Action CA16117 (ChETEC). RH, KN, KB acknowledge support from the World Premier International Research Centre Initiative (WPI Initiative), MEXT, Japan. D.G. is supported by Leverhulme Trust Grant No. RPG2019350 and by NASA through Einstein Postdoctoral Fellowship Grant No. PF6170152 awarded by the Chandra Xray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS803060. TB acknowledges support from the Polish National Science Center (NCN) Grant No. UMO2014/14/M/ST9/00707. EB is supported by NSF Grants No. PHY1607130, AST1716715, and by FCT contract IF/00797/2014/CP1214/CT0012 under the IF2014 Programme. C.G., G.M., and S.E. acknowledge support from the Swiss National Science Foundation (project number 200020160119). ROS is supported by NSF AST1412449, PHY1505629, and PHY1607520. DAB acknowledges support from National Science Foundation Grant No. PHY1404395. G.W. is partly supported by the President’s International Fellowship Initiative (PIFI) of the Chinese Academy of Sciences under grant no.2018PM0017 and by the Strategic Priority Research Program of the Chinese Academy of Science Multiwaveband Gravitational Wave Universe (Grant No. XDB23040000). AA is supported by the Carl Tryggers Foundation for Scientific Research through the grant CTS 17:113. MCM acknowledges support by NASA grant 80NSSC18K0527. This research was supported in part by the National Science Foundation under Grant No. NSF PHY1748958 to KITP UC Santa Barbara. Several authors (KB, EB, CLF, DAB, DEH, RO, KN, MCM, GM, JPL) were supported by KITP UC Santa Barbara to work on this project. JPL was supported by a grant from the French Space Agency CNES. S.C.L acknowledges support by the funding HSTAR15021.001A. DEH was partially supported by NSF grant PHY1708081, and also by the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation. He also gratefully acknowledges support from a Marion and Stuart Rice Award. We thank the Niels Bohr Institute for its hospitality while part of this work was completed, and acknowledge the Kavli Foundation and the DNRF for supporting the 2017 Kavli Summer Program.
References
 Abbott, B. P., et al. (LIGO Scientific and Virgo Collaboration) 2017, Phys. Rev. Lett., 118, 221101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, ApJ, 882, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, Phys. Rev. X, 9, 031040 [Google Scholar]
 Abolmasov, P., Karpov, S., & Kotani, T. 2009, PASJ, 61, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Abt, H. A., Gomez, A. E., & Levy, S. G. 1990, ApJS, 74, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, J. J., & Zezas, A. 2019, MNRAS, 486, 3213 [NASA ADS] [CrossRef] [Google Scholar]
 Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Antonini, F., Rodriguez, C. L., Petrovich, C., & Fischer, C. L. 2018, MNRAS, 480, L58 [Google Scholar]
 ArcaSedda, M., & CapuzzoDolcetta, R. 2019, MNRAS, 483, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Askar, A., Szkudlarek, M., GondekRosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36 [NASA ADS] [CrossRef] [Google Scholar]
 Bae, Y.B., Kim, C., & Lee, H. M. 2014, MNRAS, 440, 2714 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, S. 2018, MNRAS, 473, 909 [NASA ADS] [CrossRef] [Google Scholar]
 Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Bavera, S. S., Fragos, T., Qin, Y., et al. 2020, A&A, 635, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belczynski, K., & Banerjee, S. 2020, ArXiv eprints [arXiv:2002.08050] [Google Scholar]
 Belczynski, K., & Taam, R. E. 2004, ApJ, 603, 690 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Taam, R. E., Rantsiou, E., & van der Sluys, M. 2008a, ApJ, 682, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008b, ApJS, 174, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010a, ApJ, 714, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Dominik, M., Bulik, T., et al. 2010b, ApJ, 715, L138 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012a, ApJ, 757, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Bulik, T., & Fryer, C. L. 2012b, ArXiv eprints [arXiv:1208.2422] [Google Scholar]
 Belczynski, K., Buonanno, A., Cantiello, M., et al. 2014, ApJ, 789, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016a, Nature, 534, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Repetto, S., Holz, D. E., et al. 2016b, ApJ, 819, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Heger, A., Gladysz, W., et al. 2016c, A&A, 594, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020, ApJ, 890, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Benacquista, M. J., & Downing, J. M. B. 2013, Liv. Rev. Relativ., 16, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Berg, D. A., Skillman, E. D., Marble, A. R., et al. 2012, ApJ, 754, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., & Mezzacappa, A. 2007, Nature, 445, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Bouret, J. C., Lanz, T., Hillier, D. J., et al. 2015, MNRAS, 449, 1545 [NASA ADS] [CrossRef] [Google Scholar]
 Cantiello, M., Mankovich, C., Bildsten, L., ChristensenDalsgaard, J., & Paxton, B. 2014, ApJ, 788, 93 [Google Scholar]
 Chatterjee, S., Fregeau, J. M., Umbreit, S., & Rasio, F. A. 2010, ApJ, 719, 915 [NASA ADS] [CrossRef] [Google Scholar]
 Chatterjee, S., Rodriguez, C. L., Kalogera, V., & Rasio, F. A. 2017, ApJ, 836, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Chini, R., Hoffmeister, V. H., Nasseri, A., Stahl, O., & Zinnecker, H. 2012, MNRAS, 424, 1925 [NASA ADS] [CrossRef] [Google Scholar]
 Choksi, N., Volonteri, M., Colpi, M., Gnedin, O. Y., & Li, H. 2019, ApJ, 873, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Chrimes, A. A., Stanway, E. R., & Eldridge, J. J. 2020, MNRAS, 491, 3479 [NASA ADS] [CrossRef] [Google Scholar]
 Chruslinska, M., & Nelemans, G. 2019, MNRAS, 488, 5300 [NASA ADS] [Google Scholar]
 Chruslinska, M., Nelemans, G., & Belczynski, K. 2019, MNRAS, 482, 5012 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A. 2007, A&A, 467, 1389 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545 [Google Scholar]
 Detmers, R. G., Langer, N., Podsiadlowski, P., & Izzard, R. G. 2008, A&A, 484, 831 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, M., Belczynski, K., Fryer, C., et al. 2013, ApJ, 779, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Downing, J. M. B., Benacquista, M. J., Giersz, M., & Spurzem, R. 2010, MNRAS, 407, 1946 [NASA ADS] [CrossRef] [Google Scholar]
 Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269 [Google Scholar]
 Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eldridge, J. J., & Stanway, E. R. 2016, MNRAS, 462, 3302 [NASA ADS] [CrossRef] [Google Scholar]
 Farmer, R., Fields, C. E., Petermann, I., et al. 2016, ApJS, 227, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Farr, W. M., Stevenson, S., Miller, M. C., et al. 2017, Nature, 548, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Farr, B., Holz, D. E., & Farr, W. M. 2018, ApJ, 854, L9 [NASA ADS] [CrossRef] [Google Scholar]
 FaucherGiguère, C.A., & Kaspi, V. M. 2006, ApJ, 643, 332 [NASA ADS] [CrossRef] [Google Scholar]
 Fields, C. E., Timmes, F. X., Farmer, R., et al. 2018, ApJS, 234, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Fishbach, M., & Holz, D. E. 2017, ApJ, 851, L25 [Google Scholar]
 Fishbach, M., Holz, D. E., & Farr, B. 2017, ApJ, 840, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Foglizzo, T., Guilet, J., & Sato, J. 2009, in SF2A2009: Proceedings of the Annual Meeting of the French Society of Astronomy and Astrophysics, eds. M. HeydariMalayeri, C. Reyl’E, & R. Samadi, 143 [Google Scholar]
 Freytag, B., Ludwig, H. G., & Steffen, M. 1996, A&A, 313, 497 [NASA ADS] [Google Scholar]
 Fryer, C. L., & Heger, A. 2000, ApJ, 541, 1033 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., & Kusenko, A. 2006, ApJS, 163, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., & Young, P. A. 2007, ApJ, 659, 1438 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Woosley, S. E., & Hartmann, D. H. 1999, ApJ, 526, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Woosley, S. E., & Heger, A. 2001, ApJ, 550, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Fryer, C. L., LloydRonning, N., Wollaeger, R., et al. 2019, Eur. Phys. J. A, 55, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuller, J., & Ma, L. 2019, ApJ, 881, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661 [NASA ADS] [Google Scholar]
 Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gerosa, D., & Berti, E. 2017, Phys. Rev. D, 95, 124046 [NASA ADS] [CrossRef] [Google Scholar]
 Gerosa, D., Kesden, M., Berti, E., O’Shaughnessy, R., & Sperhake, U. 2013, Phys. Rev. D, 87, 104028 [NASA ADS] [CrossRef] [Google Scholar]
 Gerosa, D., Kesden, M., Sperhake, U., Berti, E., & O’Shaughnessy, R. 2015, Phys. Rev. D, 92, 064016 [NASA ADS] [CrossRef] [Google Scholar]
 Gerosa, D., Berti, E., O’Shaughnessy, R., et al. 2018, Phys. Rev. D, 98, 084036 [NASA ADS] [CrossRef] [Google Scholar]
 Giersz, M., Heggie, D. C., Hurley, J. R., & Hypki, A. 2013, MNRAS, 431, 2184 [NASA ADS] [CrossRef] [Google Scholar]
 Glebbeek, E., Gaburov, E., de Mink, S. E., Pols, O. R., & Portegies Zwart, S. F. 2009, A&A, 497, 255 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Green, A. M. 2017, Phys. Rev. D, 96, 043020 [NASA ADS] [CrossRef] [Google Scholar]
 Groh, J. H., Ekström, S., Georgy, C., et al. 2019, A&A, 627, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gültekin, K., Miller, M. C., & Hamilton, D. P. 2004, ApJ, 616, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Gültekin, K., Miller, M. C., & Hamilton, D. P. 2006, ApJ, 640, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Hainich, R., Oskinova, L. M., Shenar, T., et al. 2018, A&A, 609, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamilton, A. J. S., & Sarazin, G. L. 1982, MNRAS, 198, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532 [NASA ADS] [CrossRef] [Google Scholar]
 Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Herwig, F. 2000, A&A, 360, 952 [NASA ADS] [Google Scholar]
 Holgado, A. M., Ricker, P. M., & Huerta, E. A. 2018, ApJ, 857, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Hong, J., Vesperini, E., Askar, A., et al. 2018, MNRAS, 480, 5645 [NASA ADS] [CrossRef] [Google Scholar]
 Hotokezaka, K., & Piran, T. 2017, ApJ, 842, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Sippel, A. C., Tout, C. A., & Aarseth, S. J. 2016, PASA, 33, e036 [NASA ADS] [CrossRef] [Google Scholar]
 Igoshev, A. P., & Popov, S. B. 2013, MNRAS, 432, 967 [NASA ADS] [CrossRef] [Google Scholar]
 Janiuk, A. 2017, ApJ, 837, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, S., Hirschi, R., Pignatari, M., et al. 2015, MNRAS, 447, 3115 [NASA ADS] [CrossRef] [Google Scholar]
 Kazeroni, R., Guilet, J., & Foglizzo, T. 2016, MNRAS, 456, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, S., Husa, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044007 [NASA ADS] [CrossRef] [Google Scholar]
 King, A. R. 2004, Nucl. Phys. B Proc. Suppl., 132, 376 [NASA ADS] [CrossRef] [Google Scholar]
 King, A. R., & Ritter, H. 1999, MNRAS, 309, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Klencki, J., Moe, M., Gladysz, W., et al. 2018, A&A, 619, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., Weidner, C., PflammAltenburg, J., et al. 2013, The Stellar and SubStellar Initial Mass Function of Simple and Composite Populations, 115 [Google Scholar]
 Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908 [NASA ADS] [CrossRef] [Google Scholar]
 Kusenko, A., & Segrè, G. 1996, Phys. Rev. Lett., 77, 4872 [NASA ADS] [CrossRef] [Google Scholar]
 Kushnir, D., Zaldarriaga, M., Kollmeier, J. A., & Waldman, R. 2016, MNRAS, 462, 844 [NASA ADS] [CrossRef] [Google Scholar]
 Kushnir, D., Zaldarriaga, M., Kollmeier, J. A., & Waldman, R. 2017, MNRAS, 467, 2146 [NASA ADS] [CrossRef] [Google Scholar]
 Lai, D., & Qian, Y.Z. 1998, ApJ, 495, L103 [NASA ADS] [CrossRef] [Google Scholar]
 Leung, S.C., Nomoto, K., & Blinnikov, S. 2019, ApJ, 887, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Lipunov, V. M., Postnov, K. A., & Prokhorov, M. E. 1997, Astron. Lett., 23, 492 [NASA ADS] [Google Scholar]
 LVC 2019a, GCN Circ., 25324 [Google Scholar]
 LVC 2019b, GCN Circ., 25333 [Google Scholar]
 Ma, L., & Fuller, J. 2019, MNRAS, 488, 4338 [NASA ADS] [CrossRef] [Google Scholar]
 MacLeod, M., & RamirezRuiz, E. 2015, ApJ, 803, 41 [NASA ADS] [CrossRef] [Google Scholar]
 MacLeod, M., Antoni, A., MurguiaBerthier, A., Macias, P., & RamirezRuiz, E. 2017, ApJ, 838, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Madau, P., & Fragos, T. 2017, ApJ, 840, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A., & Meynet, G. 2000, A&A, 361, 159 [NASA ADS] [Google Scholar]
 Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manchester, R. N., & Taylor, J. H. 1977, Pulsars (San Francisco: W. H. Freeman) [Google Scholar]
 Mandel, I. 2016, MNRAS, 456, 578 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [Google Scholar]
 Mandel, I., & Farmer, A. 2018, ArXiv eprints [arXiv:1806.05820] [Google Scholar]
 Mandel, I., & O’Shaughnessy, R. 2010, Classical Quantum Gravity, 27, 114007 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M. 2016, MNRAS, 459, 3432 [NASA ADS] [CrossRef] [Google Scholar]
 Mapelli, M., Giacobbo, N., Santoliquido, F., & Artale, M. C. 2019, MNRAS, 487, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36 [Google Scholar]
 Marks, M., & Kroupa, P. 2012, A&A, 543, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marks, M., Kroupa, P., Dabringhausen, J., & Pawlowski, M. S. 2012, MNRAS, 422, 2246 [NASA ADS] [CrossRef] [Google Scholar]
 Meurs, E. J. A., & van den Heuvel, E. P. J. 1989, A&A, 226, 88 [NASA ADS] [Google Scholar]
 Meynet, G., & Maeder, A. 1997, A&A, 321, 465 [NASA ADS] [Google Scholar]
 Meynet, G., Chomienne, V., Ekström, S., et al. 2015, A&A, 575, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miller, M. C., & Hamilton, D. P. 2002a, MNRAS, 330, 232 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Miller, M. C., & Hamilton, D. P. 2002b, ApJ, 576, 894 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Miller, J. M. 2015, Phys. Rep., 548, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Moe, M., & Di Stefano, R. 2013, ApJ, 778, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15 [Google Scholar]
 Morawski, J., Giersz, M., Askar, A., & Belczynski, K. 2018, MNRAS, 481, 2168 [Google Scholar]
 MurguiaBerthier, A., MacLeod, M., RamirezRuiz, E., Antoni, A., & Macias, P. 2017, ApJ, 845, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Neijssel, C. J., VignaGómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740 [NASA ADS] [CrossRef] [Google Scholar]
 Nieuwenhuijzen, H., & de Jager, C. 1990, A&A, 231, 134 [NASA ADS] [Google Scholar]
 Noutsos, A., Schnitzeler, D. H. F. M., Keane, E. F., Kramer, M., & Johnston, S. 2013, MNRAS, 430, 2281 [NASA ADS] [CrossRef] [Google Scholar]
 Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227 [NASA ADS] [Google Scholar]
 Ohsuga, K. 2007, ApJ, 659, 205 [NASA ADS] [CrossRef] [Google Scholar]
 O’Leary, R. M., O’Shaughnessy, R., & Rasio, F. A. 2007, Phys. Rev. D, 76, 061504 [NASA ADS] [CrossRef] [Google Scholar]
 Oskinova, L. M., Todt, H., Ignace, R., et al. 2011, MNRAS, 416, 1456 [NASA ADS] [CrossRef] [Google Scholar]
 Oskinova, L. M., Sun, W., Evans, C. J., et al. 2013, ApJ, 765, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Packet, W. 1981, A&A, 102, 17 [NASA ADS] [Google Scholar]
 Pakull, M. W., & Mirioni, L. 2003, in Revista Mexicana de Astronomia y Astrofisica Conference Series, eds. J. Arthur, & W. J. Henney, 15, 197 [NASA ADS] [Google Scholar]
 Park, D., Kim, C., Lee, H. M., Bae, Y.B., & Belczynski, K. 2017, MNRAS, 469, 4665 [NASA ADS] [CrossRef] [Google Scholar]
 Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paxton, B., Smolec, R., Gautschy, A., et al. 2019, ApJS, 243, 10 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Perna, R., Wang, Y.H., Farr, W. M., Leigh, N., & Cantiello, M. 2019, ApJ, 878, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246 [NASA ADS] [CrossRef] [Google Scholar]
 Popov, S. B., & Turolla, R. 2012, Ap&SS, 341, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJ, 528, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Qin, Y., Marchant, P., Fragos, T., Meynet, G., & Kalogera, V. 2019, ApJ, 870, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramachandran, V., Hamann, W. R., Oskinova, L. M., et al. 2019, A&A, 625, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RamírezAgudelo, O. H., SimónDíaz, S., Sana, H., et al. 2013, A&A, 560, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rantsiou, E., Burrows, A., Nordhaus, J., & Almgren, A. 2011, ApJ, 732, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Repetto, S., & Nelemans, G. 2015, MNRAS, 453, 3341 [NASA ADS] [CrossRef] [Google Scholar]
 Repetto, S., Igoshev, A. P., & Nelemans, G. 2017, MNRAS, 467, 298 [NASA ADS] [Google Scholar]
 Ricker, P. M., & Taam, R. E. 2008, ApJ, 672, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Haster, C.J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016a, ApJ, 824, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., Zevin, M., Pankow, C., Kalogera, V., & Rasio, F. A. 2016b, ApJ, 832, L2 [Google Scholar]
 Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016c, Phys. Rev. D, 93, 084029 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, C. L., AmaroSeoane, P., Chatterjee, S., et al. 2018, Phys. Rev. D, 98, 123005 [NASA ADS] [CrossRef] [Google Scholar]
 Sadowski, A., Belczynski, K., Bulik, T., et al. 2008, ApJ, 676, 1162 [NASA ADS] [CrossRef] [Google Scholar]
 Samsing, J. 2018, Phys. Rev. D, 97, 103014 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Sądowski, A., Lasota, J.P., Abramowicz, M. A., & Narayan, R. 2016, MNRAS, 456, 3915 [NASA ADS] [CrossRef] [Google Scholar]
 Socrates, A., Blaes, O., Hungerford, A., & Fryer, C. L. 2005, ApJ, 632, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086 [NASA ADS] [CrossRef] [Google Scholar]
 Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C. 1999, A&A, 349, 189 [NASA ADS] [Google Scholar]
 Spruit, H. C. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stevenson, S., VignaGómez, A., Mandel, I., et al. 2017, Nat. Commun., 8, 14906 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, S., Sampson, M., Powell, J., et al. 2019, ApJ, 882, 121 [Google Scholar]
 Surman, R., McLaughlin, G. C., & Hix, W. R. 2006, ApJ, 643, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Tamborra, I., Hanke, F., Janka, H.T., et al. 2014, ApJ, 792, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Tutukov, A. V., & Yungelson, L. R. 1993, MNRAS, 260, 675 [NASA ADS] [CrossRef] [Google Scholar]
 Vagnozzi, S. 2019, Atoms, 7, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Valsecchi, F., Glebbeek, E., Farr, W. M., et al. 2010, Nature, 468, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Vanbeveren, D., Mennekens, N., Shara, M. M., & Moffat, A. F. J. 2018, A&A, 615, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 VanLandingham, J. H., Miller, M. C., Hamilton, D. P., & Richardson, D. C. 2016, ApJ, 828, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Venumadhav, T., Zackay, B., Roulet, J., Dai, L., & Zaldarriaga, M. 2019, ArXiv eprints [arXiv:1904.07214] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vitale, S., Lynch, R., Sturani, R., & Graff, P. 2017a, Classical Quantum Gravity, 34, 03LT01 [CrossRef] [Google Scholar]
 Vitale, S., Gerosa, D., Haster, C.J., Chatziioannou, K., & Zimmerman, A. 2017b, Phys. Rev. Lett., 119, 251103 [NASA ADS] [CrossRef] [Google Scholar]
 Vlahakis, N., & Königl, A. 2001, ApJ, 563, L129 [NASA ADS] [CrossRef] [Google Scholar]
 Voss, R., & Tauris, T. M. 2003, MNRAS, 342, 1169 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, J. J., & Leigh, N. W. C. 2015, MNRAS, 453, 3278 [NASA ADS] [CrossRef] [Google Scholar]
 Wiktorowicz, G., Belczynski, K., & Maccarone, T. 2014, Binary Systems, their Evolution and Environments, 37 [Google Scholar]
 Wiktorowicz, G., Lasota, J.P., Middleton, M., & Belczynski, K. 2019, ApJ, 875, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Williamson, A. R., Lange, J., O’Shaughnessy, R., et al. 2017, Phys. Rev. D, 96, 124041 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E. 2016, ApJ, 824, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E. 2017, ApJ, 836, 244 [NASA ADS] [CrossRef] [Google Scholar]
 Worley, A., Krastev, P. G., & Li, B.A. 2008, ApJ, 685, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Wyrzykowski, Ł., & Mandel, I. 2020, A&A, 363, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wysocki, D., Gerosa, D., O’Shaughnessy, R., et al. 2018, Phys. Rev. D, 97, 043014 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, X.J., & Li, X.D. 2010, ApJ, 722, 1985 [NASA ADS] [CrossRef] [Google Scholar]
 Zackay, B., Venumadhav, T., Dai, L., Roulet, J., & Zaldarriaga, M. 2019, Phys. Rev. D, 100, 023007 [Google Scholar]
 Zahid, H. J., Dima, G. I., Kudritzki, R.P., et al. 2014, ApJ, 791, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1977, A&A, 57, 383 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1989, A&A, 220, 112 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
 Zaldarriaga, M., Kushnir, D., & Kollmeier, J. A. 2018, MNRAS, 473, 4174 [NASA ADS] [CrossRef] [Google Scholar]
 Zevin, M., Samsing, J., Rodriguez, C., Haster, C.J., & RamirezRuiz, E. 2019, ApJ, 871, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Additional information
Data availability: data generated during this project with population synthesis code StarTrack is available online^{10} under the tab “Download/2020: Double Compact Objects/Belczynski et al. 2020”. Code availability: the StarTrack code is not an open source code and is not publicly available. New models will be calculated upon requests to the first author: chrisbelczynski@gmail.com.
A.1. Geneva stellar models
The physics included in the employed Geneva rotating models are described in detail in Eggenberger et al. (2008), Ekström et al. (2012). The models have been computed with a modest coreovershooting during the core H and Heburning phase (the core has been extended beyond the Schwarzschild limit by a length given by 10% of the local pressure scale height). Mass loss rates by stellar winds are accounted for (see the above references for the details). The effects of rotation are included according to the theory by Zahn (1992). In these models, the transport of the angular momentum and the mixing of the chemical species in the radiative zones are due to shear instabilities and meridional currents. In contrast with many stellar models in which the effects of meridional currents are accounted for by a diffusion equation, in the present models the angular momentum transport is accounted for by solving an advective equation. In essence, transport by meridional currents is an advective process. Diffusion treatments can lead not only to wrong estimates of the amplitude, but also to wrong signs (a diffusive process always tends to flatten any gradient, while an advective process can increase gradients in some circumstances). We also note that shear instabilities are the main drivers for the transport of the chemical species, while angular momentum is mainly transported by meridional currents. The angular momentum transport by meridional currents mildly couples the core rotation with that of the envelope. In the convective zones, the transport of angular momentum and chemical species is assumed to be extremely efficient: convective zones are assumed to rotate as a solid body and the chemical composition is homogenized.
In Table A.1, we present a suite of Geneva evolutionary models (Ekström et al. 2012; Georgy et al. 2013; Groh et al. 2019) and Eggenberger et al. (in prep.) with M_{zams} = 9−120 M_{⊙} for a wide range of chemical compositions: Z = 0.014−0.0004. Specifically, we are listing CO core mass for each evolutionary model along with our estimate of the associated BH spin magnitude (see Fig. 1 for a plot of the relation between the two).
The spin magnitude is very sensitive to the BH mass (; see Eq. (2)). The CO core mass may be taken as a proxy for BH mass (cf. Appendix A.3 and Fig. A.2). Below we present the behavior of the CO core mass within the Geneva models. This discussion explains the nonmonotonic dependence of BH spin on metallicity.
As naturally expected, the CO core mass increases with the initial stellar mass within models of the same metallicity. This intuitive trend is reversed only at the highest metallicity considered here (Z = 0.014), at which stellar winds are most efficient in mass removal. In particular, this is true for very massive stars (M_{zams} ≳ 100 M_{⊙}).
Another intuitive expectation is not supported by these stellar models. For a given initial stellar mass, one may naively expect that the CO core mass would increase with decreasing metallicity, as winds are getting weaker and star remains more massive. However, for example for M_{zams} = 85 M_{⊙} we find M_{CO, 2} = 26.4, 35.8, 27.4, 44.2 M_{⊙} for Z = 0.014, 0.006, 0.002, 0.0004, respectively. In the following, we explain this nonmonotonic behavior in terms of our adopted stellar evolution model.
Stellar winds are weaker at low metallicity thus this leads to increase of the CO core mass with decreasing metallicity. This general trend is mitigated by two other physical processes. Decreasing metallicity leads to the formation of extended convective Hburning shells that tend to slow down the growth of the He core (and subsequently the CO core) mass. At a high Z, massive stars lose most of their envelopes, so there is no vertical structure available for extended convective Hburning shell. The Hburning shell is compact and moves outwards (once H is totally depleted) through the envelope, adding mass to the underlying He core. At a low Z, massive stars not only retain their envelopes longer, but they are also more compact, and increased density helps to form extended convective zones. Within the extended convective shell reaching far above Hburning (which occurs only at the shell base), intensive mixing keeps bringing new fuel to the Hburning zone (keeping it in the same position), and keeps redistributing newly formed He across this large shell. Instead of adding He to the core, the newly formed He is redistributed through the material of the extended convective shell. Another process that leads to CO core mass decrease with metallicity is connected to the diffusion of elements due to the action of meridional currents and turbulence. With decreasing metallicity stars are more compact and the vertical scale of diffusion decreases, leading to less effective mixing of fresh fuel into burning zones. This, in turn, lowers the CO core mass.
For the particular set of assumptions used in the Geneva code, the nonmonotonic behavior of the CO core mass is the effect of the complex and metallicitydependent interplay of strength of stellar winds, the Hburning shell extent, and the model for the efficiency of element diffusion within the meridional current. This explains the nonmonotonic metallicity dependence of our BH natal spin model (see Fig. 1).
Geneva stellar models.
MESA stellar models with 40% critical initial rotation.
A.2. MESA stellar models
The MESA stellar models were evolved using MESA revision 10398 (Paxton et al. 2011, 2013, 2015, 2018, 2019). The models are evolved from the premain sequence to core Hedepletion. We use temporal and spatial parameters similar to those used in Farmer et al. (2016) and Fields et al. (2018) that provide convergence to the ≈10% level. The MESA models use the mesa49.net network that follows 49 isotopes from ^{1}H to ^{34}S. We include mass loss using the “Dutch” wind scheme with an efficiency value of 0.8 (Nieuwenhuijzen & de Jager 1990; Nugis & Lamers 2000; Vink et al. 2001; Glebbeek et al. 2009). We use the Ledoux criterion for convection with an efficiency parameter of α_{MLT} = 1.6, and the mlt++ approximation for convection (Paxton et al. 2013).
Additional mixing processes due to convective boundary mixing are included using the exponentially decaying diffusion coefficient framework of Herwig (2000) based on hydrodynamic simulations of Freytag et al. (1996). The following values of f were used: f = 0.014 above H and Heburning regions, f = 0.001 below H and Heburning regions, f = 0 elsewhere. We note that the additional parameter, f_{0} = 0.001 was used, where f_{0}H_{P} is the distance from the boundary inside the convective zone where the exponential decay starts. Jones et al. (2015) compare GENEC models to MESA models and find that using f = 0.022 on top of convective core H and He burning in MESA matches GENEC models with a penetrative overshoot, α_{ov} = 0.2H_{P}. Given that the GENEC models used in this paper use α_{ov} = 0.1H_{P}, using f = 0.014 in MESA yields similar core masses and lifetimes.
In MESA, rotation is implemented using the shellular approximation (Zahn 1992; Meynet & Maeder 1997). In our models, we initialize solid body, uniform rotation at the zero age mainsequence as a fraction of the Keplerian critical rotation rate, in most of our models this value is 40%. We use the suggested values from (Heger et al. 2000) for the diffusion coefficients for the transport of angular momentum and material due to various instabilities. Among the most influential of these mechanisms included in the MESA models, there is the TaylerSpruit dynamo. The modeling of this dynamo in the stellar models increases the efficiency in angular momentum transport and can lead to a significantly different angular momentum profiles and core rotation rates (Heger et al. 2005, see also Sect. 4.1).
We find that the TaylorSpruit dynamo essentially dominates the transport of angular momentum in our MESA models and is very efficient at transporting the angular momentum from the core to the upper layers, from which it is lost in winds. As a result, all our models end with a very similar (small) amount of angular momentum and thus small spin parameter a < 0.15; see Fig. 2 and also Table A.2. Given that stronger stellar winds lead to a more efficient loss of angular momentum, one would expect a dependence of the final spin on mass and metallicity. We find a hint of such a dependence in our models: the more massive as well as the higher metallicity models tend to end with systematically lower final spins (Fig. 2). However, this is not very clear in our models and the impact of mass and metallicity is subdominant to the fact that the efficient transport of angular momentum through the TaylorSpruit dynamo leads to small spins.
A star might spin up during its mainsequence evolution due to mass accretion from a companion or be born as a very fast rotator. We mimic this scenario with a single, very fast rotating MESA model (80% of critical rotation) and compare it to the slow and fastrotating 32 M_{⊙} MESA models plotted in Fig. 23 for both magnetic and nonmagnetic angular momentum transport. The specific angular momentum profile at core helium depletion of the very fast rotating models is presented in Fig. A.1. It clearly shows, as pointed out earlier in this paper, that the rotation rate at core helium depletion depends mainly on the angular momentum transport mechanism (the difference between the models using the TaylorSpruit dynamo will be further reduced during the advanced phases of stellar evolution) and not the initial rotation rate. This is a result of the evolution after the mainsequence where the TaylorSpruit dynamo extracts most of the angular momentum from the core. The core contracts due to the missing energy generation and the envelope expands as a consequence of the virial theorem and energy conservation. The contracting core spins up faster in the models with faster initial rotation, creating a stronger shear between core and envelope. The stronger shear leads to a stronger angular momentum transport of the magnetic dynamo () in a region where stratification is dominated by the chemical composition gradient (Spruit 2002), which is the case for the region above the helium core where the shear develops after the MS evolution in a massive star, resulting in similar specific angular momentum profiles. Therefore, a star that is either spun up by a binary companion or is born fast will still end its life with a slowspinning core when an efficient angular momentum transport mechanism is active.
In the nonmagnetic models presented in Fig. A.1, most of the angular momentum is transported in the short phase between core hydrogen and helium burning (as in the magnetic models) but the transport is much less efficient. Therefore, the core of models excluding magnetic fields will not be able to slow down and will end up having a high final rotation rate.
Fig. A.1.
Specific angular momentum profile at the end of core helium burning for the 32 M_{⊙} MESA models at Z = 0.002. The models were calculated with different initial rotation rates; slow (100 km s^{−1}, red), fast (40% of critical rotation, blue) and very fast (80% of critical rotation, purple). Magnetic (dashed line) and nonmagnetic models (solid line) are considered. Large dips in the curve indicate the base of a convective region, where angular momentum transport is very efficient. The maximal extent of the curves to the right indicate the mass of the star. We can see from these curves that the faster rotating models lose more mass in winds due to rotationenhanced mass loss. 
There are dips in the angular momentum profiles in Fig. A.1, one of the red curves (slowest rotation) has two dips because when helium is exhausted in the core there is a convective hydrogen zone (starting at 12 M_{⊙}) and the surface convective zone (starting at 25 M_{⊙}) present. This is not the case in the other models because (i) they evolve more to the blue (to higher temperature); (ii) faster rotating models tend to smooth temperature and chemical composition, leading to a shorter time with convective hydrogen core; and (iii) some models lose nearly their entire envelope, hence, there is no surface nor hydrogen convective zone. Models with high rotation end evolution with lower total mass as visible in Fig. A.1. Rotation enhances mass loss in general in several ways: directly by reducing the effective gravity, indirectly by rotationinduced mixing leading to increased luminosities and sometimes leading to quicker evolution to the red supergiant or WR phases. For the very fastrotating models evolving quasichemically homogeneously, the dominant effects are rotationinduced mixing leading to increased luminosities which, in turn, enhances the wind mass loss and reaching the WR phase earlier in the evolution.
We note that apart from the spin values, there are also differences between the final CO core masses between GENEC and MESA models (Tables A.1 and A.2). Some level of discrepancy is not unexpected solely due to the fact that those sets of models were run with slight differences in the efficiency of chemical mixing, criteria for convection, and treatment of convective boundary regions. More importantly, even though both sets of modes were computed using similar prescriptions for mass loss, differences in evolution in the HR diagram alone can lead to significant differences in the total amount of mass lost in winds and therefore the final CO core masses. This can be most clearly seen when comparing the most massive Solar metallicity models (Z = 0.014, initial masses 60, 85, and 125 M_{⊙}), in which case the MESA models end their evolution with noticeably higher CO core masses. We caution that in the case of such massive stars, and especially at higher metallicity, the difficulties in carrying out a numerical treatment of radiationdominated, superadiabatic envelope layers can lead to significant differences in the HR diagram evolution between different codes (Paxton et al. 2013). Given that such envelopes are likely dynamically unstable and highly turbulent, any 1D models of late evolution of massive stars should be consider highly uncertain.
All MESA inlists used to produce these models are publicly available^{11}.
A.3. CO core mass versus BH mass
The relation between the BH mass at formation and the final CO core mass of its progenitor in our simulations is set by formulae based on supernova modeling (the “Rapid” engine of Fryer et al. 2012), together with the presupernova mass of the star that is set by stellar and binary evolution. In Fig. A.2, we show the M_{CO}–M_{BH} relation for BHBH mergers detectable in the O1/O2 LIGO runs within models M30, M50, M60, and M70. The presented models are different in ways that affect the BH formation masses, that is, the PPSN/PSN prescription (see Sect. 2.3) as well as the assumed fraction of mass lost in neutrinos during a BH formation (1% in M30, M50, M70, and 10% in M60). We note that even though we plot M_{BH} at the moment of BH formation, it is very similar to the final BH mass. We also note that the CO core masses presented in Fig. A.2 are the final CO core masses (just before the BH formation), which in some cases have already been reduced due to PPSN mass loss. Those are the M_{CO} masses that we use in order to assign the newly formed BHs with natal spins (see Sect. 2.1).
Fig. A.2.
Final CO core masses of BH progenitors and their corresponding BH formation masses for the components (both primary and secondary) of detectable merging BHBH systems in four of our models: M30, M50, M60, and M70 (see Table 2 for an overview of all the model ingredients). For readability, the filled areas encapsulate 96% of all the BHs, excluding few outliers and showcasing the M_{CO} to M_{BH} relation for the representative majority. The presented models are different in ways that affect the BH formation masses, that is, the PPSN/PSN prescription (see Sect. 2.3) as well as the assumed fraction of mass lost in neutrinos during a BH formation (1% in M30, M50, M70, and 10% in M60). We note that model M70 only differs from M30 and M50 for M_{CO} ≳ 32.4 M_{⊙}, whereas model M50 only differs from M30 for M_{CO} > 37.5 M_{⊙}. Those threshold core masses correspond to the threshold helium core masses M_{He} in different models PPSN/PSN prescriptions. 
In general, the M_{CO} to M_{BH} relation is quite similar between our models and, for the most part, almost linear. The initial decrease in M_{BH} (for M_{CO} < 7.5 M_{⊙}) and the subsequent increase, followed by a change in slope at around M_{CO} = 11 M_{⊙}, are the result of changes in the fraction of material f_{fb} that falls back onto the protoNS in the “Rapid” supernovae engine (see Eq. (16) in Fryer et al. 2012). The maximum fallback (f_{fb} = 1.0, a direct collapse formation of a BH) is expected for the CO core masses M_{CO} either in range between 6 and 7 M_{⊙} or for M_{CO} > 11 M_{⊙}. Partial fallback and the ejection of some of the preSN mass for progenitors with M_{CO} within 7 and 11 M_{⊙} is what is responsible for the dip in M_{BH} in Fig. A.2 at the lower M_{CO} end.
The fact that BH masses are systematically smaller in model M60 compared to other models is a direct consequence of a higher mass fraction being lost in neutrinos during a BH formation (10% in M60 compared to 1% in other models).
The influence of PPSN/PSN kicks in above M_{CO} ≈ 32.4 M_{⊙} for models M30, M50, M70 (which corresponds to M_{He} ≈ 40 M_{⊙}), and above M_{CO} ≈ 37.5 M_{⊙} in the case of M60 (M_{He} ≈ 45 M_{⊙}). Above this core mass, PPSN removes outer part of the envelope and, therefore, reduces the final preSN mass with respect to evolution without PPSN (see Fig. 3). This leads to a smaller mass of the BH formed in direct collapse of the remaining postPPSN star. As expected, Fig. A.2 reveals that the more mass is lost in PPSN (depending on the model) the smaller the final BH mass. Finally, in the case of most massive stars (M_{He} > 65 M_{⊙}), all the models assume that a PSN disrupts the entire star and that no remnant remains.
The impact of weaker stellar winds in model M50 (30% of mass loss in winds at any evolutionary stage compared to other models) on the M_{CO} to M_{BH} relation is relatively small, only differentiating models M30 and M50 in any way for M_{CO} > 40 M_{⊙}. This is because in most cases of the BHBH merger formation, both the primary and secondary components are going to lose their entire envelope prior to core collapse anyway, due to RLOF. In fact, as many as 90% of all the BH formed from M_{CO} > 40 M_{⊙} progenitors in model M50 fall into the same area in Fig. A.2 as the BHs in model M30.
A.4. Making GW170104
Figure 8 shows the formation and evolution of a binary which results in a BHBH merger with properties similar to GW170104. This particular system was evolved within our model M20 (see Sect. 2.7 for an overview of our models and Appendix A.7 for an overview of our binary evolution calculations) under the assumptions of Genevabased BH natal spins (Fig. 1) and no WR spinup during the BHWR stage (Sect. 2.5. The progenitor binary was formed in a lowmetallicity environment Z = 0.001 (5% Z_{⊙}) as a pair of MS stars with masses of 94.6 M_{⊙} and 62.5 M_{⊙} on a wide orbit, with pericenter distance of ∼2000 R_{⊙}. As the primary ends its MS evolution and rapidly expands as a HG star, it initiates a stable caseB mass transfer. In model M20, we assume that 20% of the mass is accreted by the MS companion, while the other 80% is lost from the system. During mass transfer the primary is stripped of almost all of the Hrich envelope and it becomes a naked helium WR star with mass of 41.7 M_{⊙}. Only 0.3 Myr later it finishes its nuclear evolution and forms a BH (M_{BH1} = 32.9 M_{⊙}) in a direct collapse event (no natal kick), losing 10% of its mass in neutrino emission. The mass of the CO core of the BH progenitor is M_{CO} = 29.6 M_{⊙}, which means that the first BH is formed with zero natal spin (a_{spin1} = 0.0, see Fig. 1 for Z = 0.002 used as representative for Z = 0.001).
The companion, on the other hand, has increased its mass during the mass transfer to 71.1 M_{⊙}. It has also been rejuvenated (which we model according to Eq. (58) of Belczynski et al. 2008b) and spun up. At the offset of mass transfer, it is about 80% of its way through the MS phase. During the HG stage it expands to 1330 R_{⊙}, which is, however, not enough to cause a Roche lobe overflow (RLOF). It further increases its radius to 1650 R_{⊙} as a coreHeburning (CHeB) supergiant, at which point it initiates a dynamically unstable mass transfer and the common envelope (CE) phase. As a result, the wide binary orbit (a = 3678 R_{⊙}) decays to a = 34 R_{⊙}, and the secondary is left without any envelope as a naked WR star. We assume a 5% BondiHoyle accretion rate onto a compact object during CE, which leads the first BH to only a slight increase in mass (by 0.6 M_{⊙}) and spin (a_{spin1} = 0.05). Shortly after (0.3 Myr after the CE) the secondary forms a BH with mass M_{BH2} = 24.7 M_{⊙} in a direct collapse with no natal kick and only 10% neutrino mass loss. Because the secondary accretes a significant fraction of its mass (∼20%) while still on the MS, we assume that its increased rotation leads to a ∼20% increase of the CO core mass with respect to the nonrotating stellar models (see Appendix A.7). With the increased CO core mass of M_{CO} = 26.0 M_{⊙}, the second BH is assigned an initial spin of a_{spin2} = 0.14.
The BHBH is formed after 4.9 Myr of binary evolution on a close (a = 37.4 R_{⊙}), almost circular (e = 0.05) orbit. The time to coalescence via gravitational wave emission is 6.1 Gyr. For this particular evolution/model we assume no natal kicks at the formation of the massive BHs, so the BH spin vectors are aligned with the binary angular momentum (Θ_{1} = Θ_{2} = 0°). This produces an upper limit on the effective spin parameter (see Eq. (1)). For these particular BH masses, spins and spin tilts, we obtain a rather low effective spin χ_{eff} = 0.09. The progenitor binary forms at z = 1.2, so ∼5 Gyr after the Big Bang (close to a peak in star formation: z ≈ 2 means 3.2 Gyr after the Big Bang), and the BHBH merger takes place at z = 0.2 (∼11 Gyr after the Big Bang). The gravitational waves from the BHBH merger propagate for ∼2.5 Gyr to reach the LIGO detectors at the present time. All of the system properties are within 90% of the credible limits of GW170104 (Abbott et al. 2017, see also Appendix A.5 for a revised limits on χ_{eff}).
A.5. The effective spin parameter: χ_{eff}
The Bayesian analysis of GW170104 reported in Abbott et al. (2017) adopts prior assumptions about the relative likelihood of different spin magnitudes (uniform) and directions (isotropic). These assumptions are not suitable for comparison to the binary evolution model we adopt, which requires both individual spins to be initially aligned and then only mildly (if at all) misaligned by natal kicks. Recent detailed analysis addressing the impact of prior assumptions showed that they can indeed impact the inferred parameters (Vitale et al. 2017b; Williamson et al. 2017).
We therefore reassess the reported limit, concluding that χ_{eff} < 0.2 at approximately 90% confidence in the context of our model. We can justify this reanalysis using only the reported LIGO result on χ_{eff}, restricted to χ_{eff} > 0. Approximating the LIGO distribution as nearly normal with mean μ = −0.21 and width σ_{χ} ≃ 0.155, we construct a truncated normal distribution ∝θ(χ)exp − (χ + μ)^{2}/2σ^{2}, which has a 90% upper limit at x ≃ 0.2.
We arrive at a similar result by reanalyzing the underlying LIGO data using the same model and techniques (including the prior), then restricting to configurations with positive individual spins χ_{1, z}, χ_{2, z} > 0. Our revised upper limit is consistent with the range of plausible χ_{eff}, as reported by Abbott et al. (2017) (see Fig. 5 of the supplementary material), corresponding to approximately a 98−99% confidence limit within the strong assumptions of their original analysis.
A.6. NS spins
A.6.1. Observed pulsar spin period distribution
Although astronomers have amassed a large sample of pulsar spinperiod measurements, extrapolating from observed periods to birth periods is an open area of research (for a review, see Miller & Miller 2015). Because pulsar emission spins down the neutron star with time, any estimate of the birth period requires an estimate of the pulsar age, and an understanding of the rate of spindown from pulsar emission (including an understanding of magnetic field evolution in neutron stars). The age of a pulsar (t_{pulsar}) can be estimated assuming the angular momentum is lost through electromagnetic radiation from a pulsar with dipole magnetic fields (Manchester & Taylor 1977):
where n is the pulsar braking index (n = 3 for magnetic dipole radiation), P_{0}, P are, respectively, the initial and current pulsar spin periods. If the age of the pulsar is known by other means (e.g., the supernova remnant age), measurements of the spin period along with constraints on the braking index provide an estimate of the birth spin period. Clearly the pulsar age depends on the choice of n. In addition, other sources of angular momentum loss exist: for example, nonsphericities of the neutron star can cause gravitational wave emission, lowering the spin period.
Despite these uncertainties, astronomers have estimated the neutron star birth spin period distribution. The fastest pulsars could be born spinning less than 10 ms (the Crab pulsar is believed to have been born spinning at 17 ms). Popov & Turolla (2012) found that a Gaussian with an average spin period of 100 ms (with a 1σ deviation of 100 ms) fits the observations. FaucherGiguère & Kaspi (2006) found a slightly higher average, 300 ms (1σ deviation of 150 ms). Igoshev & Popov (2013) argue that the differences between these two studies could be explained by the choice of magnetic field evolution and either distribution could be made consistent with the data. Noutsos et al. (2013) also obtained a distribution of periods peaking below 125 ms, but found an additional set of longperiod birth spins (> 0.5 s). They found that poor age estimates limit determinations of the birth pulsar spin distribution and described methods to estimate pulsar ages kinematically.
With these uncertainties in mind, we can now compare stellar models to the observed spin distribution. As with the stars forming black holes, we can estimate the birth spins of of neutron stars from their massive star progenitors by assuming that the angular momentum of the collapsing core is accreted onto the protoneutron star along with its mass. We limit the accreted angular momentum (j_{acc}) to the centrifugally supported value:
where j_{shell} is the angular momentum of the accreting shell, r_{NS} is the neutron star radius, G is the gravitational constant and M_{encl} is the mass enclosed in the shell. For the neutron star radius, we assume a 10 km. We consider only compact remnants with masses below 2.5 M_{⊙} as this is our adopted maximum NS mass. By summing up the angular momentum from the accreted material, we obtain the total angular momentum of the compact remnant. The moment of inertia of a neutron star (I_{NS}) depends upon the equation of state (Worley et al. 2008) and can vary by roughly 50% with the choice of the equation of state but it is roughly linear with neutron star mass. For our calculations, we assume:
where M_{NS} is the neutron star mass. With the total angular momentum and moment of inertia, we can determine the birth spin period of the neutron stars from our progenitors (Fig. A.3). The spins from these models assume stars born rotating at 40% breakup velocity, producing core spins near the maximum allowed by a given angular momentum transport mechanism. We have not included any angular momentum loss mechanisms, but there are several possibilities of extracting rotational energy from the protoneutron star. For example, if the rotational energy is tapped to help driving the explosion (e.g., by interaction with a disk or when forming a magnetar), the total angular momentum of the system will be reduced. Therefore the rotation rates produced in our models are only upper limits on the real spin rates.
The results presented in Fig. A.3 demonstrate that some magnetic braking is needed to reduce the angular momentum of the core, confirming decadesold arguments for magnetic breaking (Heger et al. 2000; Fryer & Heger 2000). The Geneva models with the original distribution of angular momentum through a star produce only subms pulsars. While the spins from these models without magnetic coupling have too much angular momentum, they are ideally suited to determine the amount of coupling needed to produce the correct spin periods. We assume the angular momentum is distributed with a constant angular velocity (solidbody rotation) across the mass of different cores. If we assume the angular momentum to be conserved, we can calculate the angular velocity of the core and recalculate the spin of the pulsar produced. For example, if we sum up the angular momentum of the CO core and divide it by its moment of inertia, we get an average angular velocity. If we use this constant angular velocity to redistribute the angular momentum in the CO core, we revise our estimate of the spin of the compact remnant. For this redistribution, the spin periods remain below a ms (magenta triangles). To truly slow down the birth spin periods of neutron star, we must couple the angular momentum through the helium core (meaning that one assumes a constant angular velocity from star center to the outer boundary of the helium core). Figure A.3 shows the resultant spins if we assume coupling through the helium core with a helium core definition of X_{He} > 0.4, X_{H} < 0.01 (blue squares), and X_{He} > 0.5 (red empty circles). This produces maximum birth spin periods of between ∼10−1000 ms. Further coupling, the Hrich layers, would produce spin periods (solid green circles) that are too slow to match the observed pulsar distribution. Given that these spins are maximum spin values (recall that we are using fast rotating progenitors and assuming no angular momentum loss in the explosion mechanism), our pulsar spin observations argue against angularmomentum coupling beyond the helium core.
Fig. A.3.
Neutron star spins from our progenitors as a function of neutron star mass. We assume any remnant from our suite of models with a mass below 2.5 M_{⊙} is a neutron star and consider only these remnants. With the Geneva models, we have studied the angular momentum coupling of different burning layers and their effect on the neutron star spin. The Geneva models with the original mild coupling produce spin periods shorter than 1 ms (not shown here, but very close to magenta triangles; see below). Some angular momentum loss in the supernova engine (e.g., magnetic coupling such as a magnetar engine) would be required to slow neutron stars down for such models to match the data. The spins are nearly the same if the coupling extends through the CO core (magenta triangles). Coupling through the helium layer (blue squares and empty red circles; different helium core definitions) and hydrogen layer (green solid circles) produces slower neutron stars. The heliumcoupled models match the data (but realize that we are using rapidly rotating progenitors with no angular momentum loss during NS formation). The MESA models (not shown here) with the TaylerSpruit dynamo produce ∼10 ms pulsars; generating spins that match fastestspinning pulsars with rapidly spinning progenitor stars, an indication that the MESA models produce reasonable coupling results. The Fuller models (not shown here) do not produce rotating neutron stars, requiring spinup mechanisms in the supernova engine to match the data. For comparison we also mark (with black lines) the range of two observational estimates. More details on models and observations are given in Appendix A.6.1. 
If we compare the spin periods produced using our MESA models using a prescription similar to Heger et al. (2000), we find spin periods of roughly ∼6−10 ms, matching the fastestspinning nonrecycled pulsars and the results of Heger et al. (2000), Fryer & Heger (2000). These periods are produced using rapidly spinning progenitors, so we would expect their periods to match the fastestspinning systems and this result is a relatively good match to the observed pulsar distribution. For the Fuller models, the birth spin period would be too slow to match the data, arguing for some mechanism to spinup the neutron star. Simulations of the asymmetries in the engine do show spinup in the core. The amount of spinup, and whether it is sufficient to solely explain the pulsar spin periods, remains a matter of debate (Blondin & Mezzacappa 2007; Fryer & Young 2007; Foglizzo et al. 2009; Rantsiou et al. 2011; Kazeroni et al. 2016).
A.6.2. Rotational explosion properties
It has been argued that fast spinning magnetars can drive a subset of supernova explosions. Here we determine the role spins can play in the supernova explosion itself. With strong (magnetarstrength: ∼10^{15} G) magnetic fields, a pulsar can quickly release the rotational energy of a newly formed neutron star. To calculate the energy available for such a model, we need to estimate the rotational energy of the neutron star:
where ω is the angular velocity. If the neutron star is spinning with a ms period (e.g., original Geneva models), it can produce a 10^{51} erg explosion if it can tap 10% of the rotational energy to drive a jet. With the 6−10 ms periods (e.g., MESA models with rapidly rotating stars) we find that even if all of the rotational energy is tapped to drive an explosion, the rotation is unable to produce a normalenergy supernova. If the strong magnetic fields are formed quickly, a spinpowered magnetar engine will deposit its energy in the slowlymoving ejecta, accelerating this innermost material (and adding additional heating) but not contributing significantly to the total energy budget of the explosion.
Similarly, because the angular momentum is lowest in the cores of these stars, such models are unable to form a disk around the neutron star. Hence, engines that invoke jets produced by magnetic fields generated in a disk will not work. We note, however, that a disk can form after the formation of a 3 M_{⊙} remnant (presumably a black hole) as then more angular momentum is trapped in the remnant. It is likely that any model with TaylerSpruit efficient angular momentum transport (e.g., MESA or Fuller models) invoking these engines will have to rely upon some means to spin up the star prior to collapse (see review by Fryer et al. 1999). We further note that this statement does not depend on the initial stellar spin of models as independent of adopted initial spin, both MESA and Fuller models end up with very small angular momentum in the core (see Fig. 23). Such rare spinup events would be able to explain rare outbursts such as gammaray bursts but if these spins are correct, it is unlikely that magnetars or disks play a big role in explaining supernovae (for a review, see Fryer et al. 2019).
A.7. Binary evolution calculations
We employed the StarTrack population synthesis code (Belczynski et al. 2002, 2008b). The existing improvements relevant for massive star evolution include updates to the treatment of CE evolution (Dominik et al. 2012), the compact object masses produced by core collapse/supernovae (Fryer et al. 2012; Belczynski et al. 2012a) including the effect of pairinstability pulsation supernovae and pairinstability supernovae (Belczynski et al. 2016c), stellar binary initial conditions set by observations (de Mink & Belczynski 2015), and observationally constrained star formation rate and metallicity evolution over cosmic time (Madau & Dickinson 2014; Belczynski et al. 2016a). The code adopts by default the fallbackdecreased natal kick prescription (see below). Additionally, we explore three different models for the BH natal spins proposed in the current upgrade (Sect. 2.1) as well as three different models for PPSN/PSN (Sect. 2.3).
In our population synthesis calculations we evolve stars on a finite grid of metallicity: Z = 0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.0007, 0.0008, 0.009, 0.001, 0.0015, 0.002, 0.0025, 0.003, 0.0035, 0.004, 0.0045, 0.05, 0.006, 0.0065, 0.007, 0.0075, 0.008, 0.0085, 0.009, 0.0095, 0.01, 0.015, 0.02, 0.025, 0.03. If in our population synthesis model the metallicity is Z < 0.00089 we adopt the BH spin model corresponding to Z = 0.0004; if 0.00089 ≤ Z < 0.00346 we adopt the BH spin model corresponding to Z = 0.002; if 0.00346 ≤ Z < 0.00916 we adopt the BH spin model corresponding to Z = 0.006; and if Z ≥ 0.00916 we adopt the BH spin model corresponding to Z = 0.014. The limits are half points in decimal logarithm between the four metallicities of the BH natal spin model (Z = 0.014, 0.006, 0.002, 0.0004).
For the initial orbital period distribution of massive binaries we use f_{p}(logp/day)∝(logp/day)^{−0.5} in the range of [0.15, 5.5] (Sana et al. 2012). For the initial eccentricity distribution we use f_{e}(e)∝e^{−0.42} in the range of [0.0, 0.9]. The initial mass of the primary star is taken from a broken power law IMF: ∝M^{−1.3} for 0.08 ≤ M < 0.5 M_{⊙}, ∝M^{−2.2} for 0.5 ≤ M < 1.0 M_{⊙}, and ∝M^{−2.3} for 1.0 ≤ M ≤ 150 M_{⊙} (Kroupa 2001; Bastian et al. 2010; Duchêne & Kraus 2013). The initial secondary mass is taken from a uniform mass ratio distribution f_{q}(q)∝q^{0} in the range of q ∈ [0.1, 1].
We adopt maximum binary fraction: f_{bi} = 1.0 for stars of any mass and any metallicity. This is a sound assumption for massive stars of nonnegligible metallicity (Raghavan et al. 2010; Chini et al. 2012; Moe & Di Stefano 2013, 2017). However, the binarity may be smaller for low mass stars, for example, f_{bi} = 0.5 (2/3 of stars in binaries). For each metallicity, we evolve N = 2 × 10^{6} massive binaries (primary mass > 5 M_{⊙}, secondary mass > 3 M_{⊙}) and this corresponds to the total simulation stellar mass of M_{sim} = 1.9 × 10^{8} M_{⊙} (all stars over entire IMF) for f_{bi} = 1.0. However, had we assumed f_{bi} = 0.5 for all stars, then M_{sim} = 2.8 × 10^{8} M_{⊙}. Therefore, such a change would decrease of all our rate predictions for double compact objects by ∼30% (see Eq. (7) of Belczynski et al. 2016b).
The initial distributions described above assume that all the binary parameters are independent from each other. However, various correlations between those parameters have long been suggested by the observations (e.g., Abt et al. 1990; Duchêne & Kraus 2013). It was only recently that Moe & Di Stefano (2017), using results from more than 20 surveys of massive binary stars, were able to fit analytic functions to the correlated distributions and obtain a join probability density function f(M_{1}, q, P, e)≠f(M_{1})f(q)f(P)f(e). In some pockets of the entire parameter space they found differences larger than an order of magnitude with respect to the typically used method of combining several independent distributions. Similarly, even though a conclusive evidence for significant IMF variations with environmental conditions is still lacking, there is an increasing amount of results suggesting such departures from the IMF universality (Bastian et al. 2010; Kroupa et al. 2013). Notably, both theoretical arguments (see Klencki et al. 2018, for an overview) and observations of some GCs in the Milky Way (Marks & Kroupa 2012; Marks et al. 2012) seem to point towards a topheavy IMF^{12} in lowmetallicity galaxies, which are likely an important formation site of massive BHBH mergers. Recently, Klencki et al. (2018) analyzed the significance of the intercorrelations in initial distributions quantified by Moe & Di Stefano (2017) as well as of the possible variations in the IMF slope for the massive stars on the formation rate and properties of compact binary mergers. They found that the effect of those factors is very small compared to other uncertainties (e.g., rates affected by less than a factor of 2). Their result holds even for very significant changes in the IMF slope due to the coupling of the IMF and the cosmic SFRD. This justifies the simplified assumptions of the universal IMF and noncorrelated initial binary parameter distributions used in this study.
For old models (M10 and M13), we redid the calculations with the same input physics, but with the addition of the new distribution of natal BH spins. Additionally, we have updated the calibration for all models (decreasing rates by a factor of 0.926) to account for small inconsistencies in our previous estimates (Klencki et al. 2018). Below we comment on factors introduced in our new models (M20, M23, M25, M26).
The fraction of mass retained in the binary (f_{a}) during stable RLOF is not wellestablished, and could be fully conservative (f_{a} = 1), fully nonconservative (f_{a} = 0), or anywhere in between (e.g., Meurs & van den Heuvel 1989). Donor stars are typically the more massive components, as they are the ones that evolve (and typically expand) more quickly. Consequently, more massive donors often have a much shorter thermal timescale than their companions. For that reason, in the case of a thermaltimescale mass transfer when the mass transfer rate is related to the thermaltimescale of the donor star, the less massive companions will not have enough time to thermally readjust in order to accommodate all of the transferred mass. This is likely to lead to a large fraction of the transferred mass being ejected from the system. Even in the case of a slower, nucleartimescale mass transfer, accretion is typically expected to be small because the accretor quickly becomes spun up to its critical surface rotation velocity. We note that the picture is possibly more complicated due to the uncertain efficiency of stellar winds in carrying away the angular momentum from the accretor surface (e.g., Vanbeveren et al. 2018).
Nonetheless, we decided to revise and test our assumption on the accretion efficiency. In a previous work (as well as the remaining models in this study), we adopted f_{a} = 0.5. Recent estimates of mass transferring BHBH progenitors resulted (typically) in f_{a} < 0.5 (Stevenson et al. 2017), so in models M20M26 we adopt f_{a} = 0.2. As a consequence, the secondary stars (accretors) remain less massive than in our previous models, and this generates a wider BHBH binary mass ratio distribution than reported in our earlier studies. The efficiency of accretion in the first episode of mass transfer was noted to have possible impact on the BHNS formation (Kruckow et al. 2018). We note that Stevenson et al. (2017) did not use a single value for f_{a}, but, rather, estimated the accretion efficiency from the relative thermal timescales of the donor and accretor for a given binary.
Even a small amount of accretion during RLOF (∼ few percent of the accretor’s mass) may effectively spin up accreting stars (Packet 1981). With our adopted RLOF retention fraction of f_{a} = 0.2, accretors in BHBH progenitor binaries typically gain about 10 M_{⊙}, which is enough to spin up even very massive stars. This accretion usually happens around the middle (or shortly thereafter) of the accretor’s main sequence life, and therefore it allows for effective rotational mixing and the formation of more massive He and CO cores. Geneva stellar evolution models indicate that the CO core masses in rotating stars (40% critical velocity) are 20% more massive than in nonrotating models. So far all CO core masses calculated in our binary evolution models were obtained from nonrotating models (Hurley et al. 2000). Here we increase CO core mass of accreting lowmetallicity (Z < 0.002) MS stars by 20%. For highmetallicity stars, the effects of rotation on the CO core mass are suppressed due to angular momentum loss through stellar winds (Georgy et al. 2012). This change may increase the mass of the second BH, and also lower its spin magnitude. In all models we assume that material is lost from a binary in RLOF with specific angular momentum dJ/dt = j_{loss}[J_{orb}/(M_{don} + M_{acc})](1 − f_{a})dM_{RLOF}/dt, with j_{loss} = 1.0 (Podsiadlowski et al. 1992).
Our criteria for the mass transfer stability and the occurrence of CE are described in detail in Sect. 5 of Belczynski et al. (2008b). To treat CE evolution, we assume energy balance with fully effective conversion of orbital energy into envelope ejection (α = 1.0), while the envelope binding energy for massive stars is calibrated using a parameter λ that depends on stellar radius, mass, and metallicity for all models. For massive stars we assume λ ≈ 0.1 Xu & Li (2010). Additionally, we either allow or do not allow (namely, in submodels A and B) for HG stars to initiate and survive CE evolution Belczynski et al. (2007), Pavlovskii et al. (2017). In submodel A, only stars with well developed coreenvelope boundary (beyond HG) can successfully initiate and survive CE, depending on the energy balance. In submodel B, we also allow for HG stars to initiate and survive CE. Recent calculations show that the accretion rates onto compact objects in CE inspiral can be reduced even by a factor of 10^{−2} with respect to the rates resulting from the BondiHoyle approximation when the structure of the envelope (in particular, the density gradients around the inspiraling object) are taken into account (Ricker & Taam 2008; MacLeod et al. 2017; MurguiaBerthier et al. 2017; Holgado et al. 2018). MacLeod & RamirezRuiz (2015) argue that accretion structures forming around compact objects embedded in the CE may span a large fraction of the envelope radius, and so traverse substantial density gradients. Introducing gradients in the CE structure leads to net nonzero angular momentum of the flow around an accreting object (which is not the case in the standard Hoyle formalism) and in doing, it limits accretion: steeper density gradients correspond to smaller accretion. The typical values of the density gradients found by these authors introduce a considerable perturbation to the flow. For most density gradients considered by MacLeod et al. (2017), the accretion rate is well below 10% of BondiHoyle accretion rate. Based on these findings, we adopt f_{bond} = 5% of the BondiHoyle accretion rate onto a BH in CE in our current simulations. Therefore, massive BHs (M_{BH} ∼ 30 M_{⊙}) accrete ∼0.5 M_{⊙} in a typical CE event, as opposed to ∼1.0 M_{⊙} in our earlier calculations (f_{bond} = 10%). To assess the BondiHoyle accretion rate we follow the approach presented in Belczynski et al. (2002).
The spectroscopic analysis of OB stars in the Local Group (MW, SMC, LMC, e.g., Oskinova et al. 2013; Hainich et al. 2018; Ramachandran et al. 2019) as well as other lowmetallicity dwarf galaxies (e.g., Bouret et al. 2015) has shown a systematic offset in the wind mass loss rates between theoretical predictions (Vink et al. 2001, assumed in most our models) and the empirical log Ṁ–logL relation. Namely, the theoretical models seem to overestimate the actual wind mass loss rates of hot stars by at least a factor of a few. The most likely cause for this discrepancy is clumping of the wind, which is not accounted for in the standard models of radiatively driven smooth winds. For that reason, in model M50, we reduce all the wind mass loss rates down to 30% of the usually assumed values.
Compact remnants formed in supernovae can receive proper motions via two classes of enginedriven natal kicks: asymmetric matter ejecta or asymmetric neutrino emission. For BHs, asymmetric matter ejection mechanisms only work when matter is ejected, as opposed to prompt collapse or complete recapture of all ejected material (“fallback”). In our calculations, we expect that only a small fraction of systems eject a substantial amount of matter, enabling a substantial BH natal recoil kick. In contrast, asymmetric neutrino emission mechanisms operate even without any mass ejection, and thus they can affect any model of BH formation. Although neutrino mechanisms have been invoked to explain recoil velocities of pulsars and Xray binaries (Lai & Qian 1998; Repetto & Nelemans 2015), the proposed kick models all require strong magnetic fields. Models without strong magnetic fields are unable to produce significant neutrino kicks (Tamborra et al. 2014).
The sterile neutrino oscillation model (Kusenko & Segrè 1996; Fryer & Kusenko 2006) argues that neutrinos produced in the core could oscillate to sterile neutrinos and escape the core. Large magnetic fields align the ions and electrons, forcing both the neutrino scattering and absorption cross sections to be anisotropic. To ensure asymmetric neutrino emission, these strong magnetic fields must be at the last scattering surface for the neutrinos. If the magnetic field in the core is high enough to align the ions and electrons, the neutrinos in the core will be anisotropic. If these neutrinos oscillate into sterile neutrinos, they can escape, retaining their anisotropies and generating large natal kicks.
Alternatively, the neutrino bubble instability (Socrates et al. 2005) argues that magneticacoustic instabilities develop, transporting neutrino radiation to the photosphere. These instabilities carry neutrinos, and the luminosity escaping the neutrinosphere will be enhanced at these “bubbles”. If the magneticacoustic bubbles are globally asymmetric, the neutrino emission will also be asymmetric, producing a neutrinodriven kick. Current supernova calculations have several limitations: (i) they do not model high magnetic fields; (ii) they do not sufficiently resolve the hydrodynamics; and (iii) they do not include the neutrino oscillation physics necessary to produce these kicks. So high, neutrinodriven BH natal kicks cannot be ruled out.
In models M10 and M20, we test asymmetric mass ejection kicks, as we employ fallbackdecreased natal kicks Fryer et al. (2012). To mimic asymmetric neutrino emission mechanisms, we explore an alternative phenomenological prescription for BH natal kicks in models M13, M23, M25, M26, M33, M35, and M43, where we impart kicks which are random in direction, with magnitude drawn from a Maxwellian with a given 1dimensional σ, independent of the BH mass or its progenitor history (see Table 2).
A.8. Detectabilty of mergers in gravitational waves
All the compact object mergers are redistributed according to star formation history across cosmic time (z ≈ 0−15 for Population I and II stars; e.g., Madau & Dickinson 2014), taking into account the time delay between binary formation at Zero Age Main Sequence and the merger. For each merger at redshift z, we use phenomenological inspiral–merger–ringdown waveforms (IMRPhenomD; Khan et al. 2016) to calculate the signaltonoise ratio in the O1/O2 LIGO runs. A given merger is considered detectable (depending on its random sky location and orbital orientation with respect to the detectors) if the signaltonoise ratio (S/N) in a single detector is greater than 8. Only detectable mergers are used in our comparisons with O1/O2 data (e.g., rates in Table 3 or effective spins in Sect. 3.4).
With this method, we obtain a selfconsistent redshift distribution of mergers in the local Universe, and we also account for LIGO/Virgo detectability of our synthetic mergers (in particular, more massive mergers can be detected at larger redshifts). A more detailed description of detectability criteria can be found in Belczynski et al. (2016b,a,c).
A.9. BHNS merger rate limits from the ongoing O3 run
Here, we analyze what the constraint on the rate density of BHNS mergers that the recent LIGO/Virgo candidate (LVC 2019a,b, S190814bv)^{13} entails.
We do not know the mass of the object, but we assume that the mass of the NS M_{NS} must be in the range from 1.3 to 3 M_{⊙}, while the BH mass M_{BH} is in the range from 5 to 50 M_{⊙}. We denote the chirp mass of the system as ℳ. The O3 has lasted for about 4.5 months with the 80% uptime which given the effective observation time t_{obs} of 150 days, or 0.41 years for O1/O2, and t_{obs} of 3.6 months, or 0.30 years. For the range of the search, we assume that the sensitivity to NSNS mergers was r_{NSNS} = 80 Mpc (O1/O2) and r_{NSNS} = 135 Mpc (O3). We can estimate the range of BHNS merges as:
where we assume that the fiducial chirp mass of a NSNS system is 1.2 M_{⊙}. The observed volume is then . The rate density of BHNS mergers, estimated from this one event, can be estimated as
and the range corresponds to the limits at which we have allowed NS and BH mass to vary. This estimate will go down by approximately a factor of two if no similar objects are further detected in the remainder of O3 (0.7 yr). On the other hand, the range becomes narrower once we know the mass estimates of NS and BH. We note however, that including the effect of Poisson distribution width will broaden the result again. Thus we are confident that the above estimate is close to accurate under the assumption that the detected system was, in fact, a BHNS merger and not a BHBH merger with a very light secondary BH. We confront the merger rate limit estimated above with the results from our models in Fig. 25.
All Tables
All Figures
Fig. 1.
Magnitude of natal BH spin as a function of the CO core mass of the collapsing star for the Geneva stellar models with 40% critical initial velocity and mild angular momentum transport by meridional currents. Lines mark our adopted model for natal BH spins and its dependence on metallicity (see Eq. (3)). The star’s CO core mass may be used as a proxy for the BH mass (see Appendix A.3). 

In the text 
Fig. 2.
Magnitude of natal BH spin as a function of the CO core mass of the collapsing star for the MESA stellar models with 40% critical initial velocity and the TaylerSpruit magnetic dynamo (efficient) angular momentum transport. Lines mark our adopted model for natal BH spins and its dependence on metallicity (see Eq. (4)). 

In the text 
Fig. 3.
Adopted models for pairinstability pulsation supernova mass loss. For a given He core mass we show the mass of a star after PPSN mass loss. Moderate PPSN mass loss is adopted directly from Leung et al. (2019), while its modified (50% reduced mass loss) version is presented as weak PPSN model. Strong PPSN are adopted from Belczynski et al. (2016c). 

In the text 
Fig. 4.
Comparison of the two black hole accretion models employed in the StarTrack code. We present the evolution of the mass transfer rate from a donor star (Ṁ_{don}; top panel, solid lines), BH mass accumulation rate (Ṁ_{BH}; top panel, dashed lines), BH mass (M_{BH}; middle panel), and BH spin (a_{spin}; bottom panel) during the RLOF stable mass transfer phases. Critical mass transfer rates (above which mass ejection from a system is expected) are provided for reference. During the first part of the RLOF (t ≈ 4−4.8 Myr) the donor is a MS star, during the shortduration peak it is a HG star (t ≈ 4.9 Myr), and for the remaining RLOF it is a heliumcore burning star. See Sect. 2.4 for a description of the full evolutionary path and both accretion models. 

In the text 
Fig. 5.
Black hole spin magnitudes as a function of the orbital period for our Z = 0.014 (circles) and Z = 0.0004 (triangles) MESA models. The color coding corresponds to the black hole remnants with masses: M_{BH} < 15 M_{⊙} (black/blue), 15 M_{⊙} < M_{BH} < 30 M_{⊙} (cyan/red), M_{BH} > 30 M_{⊙} (green/magenta). Binary systems with short orbital periods: 0.1−1.3 d and WR stars will produce BHs with broad range of spin magnitudes: 0.15−1. Systems with orbital periods below 0.1 d form black holes from tidally locked WR stars, and the BH spins are maximal. Binaries with an orbital period above 1.3 d produce BHs with spins below 0.2 and typically are not tidally locked at all. A set of the most massive, lowestmetallicity stars have such dense cores that tidal spinup does not dramatically increase their angular momentum and these stars require orbital periods of less than 0.1 d to have spin values above 0.1 (lower set of green triangles). However, we ignore this fact, and let these cores to be spunup to estimate the maximal effect of tidal interactions in our models (higher set of green triangles). The black curve shows our fit to the BH spin magnitude as a function of orbital period. 

In the text 
Fig. 6.
Cosmic star formation histories adopted in our modeling (see Sect. 2.6 for details). We note how the different SFRDs agree at low redshifts (z ≲ 2), while they disagree by as much as a factor of three at earlier epochs. 

In the text 
Fig. 7.
Stellar massweighted gasphase metallicity versus redshift, Z(z). At every epoch, we adopted a Gaussian distribution of log(Z/Z_{⊙}), centered at the mean metallicity and with dispersion σ = 0.5 dex. We note that the new metallicity for the star forming gas is noticeably higher than the mean cosmic metallicity adopted in the past (see Sect. 2.6). 

In the text 
Fig. 8.
Example of a possible route leading to the formation of a BHBH merger similar to GW170104. This example follows the classical isolated binary evolution channel. In this model (M20.B) we employ the Geneva BH natal spins, assume that massive BHs do not receive natal kicks and that their spins are aligned with the binary angular momentum (Θ_{1} = Θ_{2} = 0°), producing an upper limit on the effective spin parameter (χ_{eff}). Yet, this system has χ_{eff} = 0.09, within LIGO’s 90% credible limits for GW170104 [ − 0.24:0.13]. Both BH masses are also within the limits: M_{BH1} = 31.0 M_{⊙} [25.4:38.2] and M_{BH2} = 20.1 M_{⊙} [15.6:25.0]. 

In the text 
Fig. 9.
Example of a possible route leading to the formation of a BHBH merger similar to GW170729. This example follows the classical isolated binary evolution channel. In this model (M30.B) we employ the MESA BH natal spins, assume that massive BHs do not receive natal kicks and that their spins are aligned with the binary angular momentum (Θ_{1} = Θ_{2} = 0°), producing an upper limit on the effective spin parameter (χ_{eff}). We note that given the small binary separation at the BHWR stage (a ≈ 7 R_{⊙}), the WR star is most likely going to become tidally synchronized (see Sect. 2.5). If we assume that the rapidly rotating WR star collapses into rapidly spinning BH (a_{spin2} = 1) then the effective spin of the presented system would increase from χ_{eff} = 0.137 (no tidal spinup; presented in the figure) to χ_{eff} = 0.484 (full tidal spinup). Both values are within the LIGO/Virgo 90% credible limits for GW170729 [0.11:0.57]. Both BH masses are also within the limits: M_{BH1} = 50.6 M_{⊙} [40.4:67.2] and M_{BH2} = 34.3 M_{⊙} [24.2:43.4]. 

In the text 
Fig. 10.
Merger rate density of BHBH mergers from Population I/II stars. We employ models M10.B and M30.B to illustrate the effects of our assumptions on the star formation rate and cosmic metallicity evolution on the BHBH mergers. LIGO/Virgo O1/O2 constraint on the BHBH merger rate in local Universe is also shown. We note that the decrease of the BHBH merger rate density at lowredshift (from old to new models) is due to the average metallicity of stars at any given redshift in old models, as it is much lower than the average metallicity of stars in new models (see Sect. 3.2 for details). 

In the text 
Fig. 11.
Merger rate density of BHBH mergers for various assumptions about the common envelope and BH natal kicks. The possibility of development and survival of the CE phase with a Hertzsprung gap donor, increases the lowredshift BHBH merger rate density by ∼1 order of magnitude: compare models M30.A and M30.B (a similar effect is found for other models). Natal kicks tend to decrease the BHBH merger rate density by ∼1.5 order of magnitude: from fallback attenuated natal kicks (almost no BH natal kicks: model M30.B) to full scale BH natal kicks (as high as observed for Galactic single pulsars: model M33.B). We note that the degeneracy between the tested input physics; by applying high natal kicks we can bring model M30.A down to agree with the LIGO/Virgo estimate, or by allowing for HG CE we can bring models M35.B and M33.B up to also match the LIGO/Virgo constraint. 

In the text 
Fig. 12.
Merger rate density of NSNS, BHNS and BHBH mergers for model M30. For comparison, we show LIGO/Virgo O1/O2 constraints on merger rate densities. While actual estimates are available for the BHBH and NSNS mergers, only an upper limit is available for the BHNS mergers from O1/O2 data (but see Appendix A.9). We note that the presented merger rate densities are consistent with the LIGO/Virgo estimates. 

In the text 
Fig. 13.
BH masses in BHBH mergers within the design advanced–LIGO horizon (z < 2) for models M60.B (strong pairinstability pulsation supernovae), M70.B (moderate PPSN), and M30.B (weak PPSN). Top panel: distribution of individual BH masses, bottom panel: total BHBH system mass. All distributions are intrinsic (neither redshifted nor weighted by detection probability). For comparison, we also indicate the range of LIGO/Virgo O1/O2 mean mass estimates and their most narrow and the widest allowed range within 90% confidence limits. Individual BH masses may reach ∼40 M_{⊙} (M60.B) and ∼55 M_{⊙} (M30.B), while the total BHBH mass may reach ∼80 M_{⊙} (M60.B) and ∼110 M_{⊙} (M30.B). Note that these distributions only very approximately resemble powerlaws: ∼M^{−3.6} (for individual BH masses) and ∼M^{−4.0} (for total BHBH mass). Powerlaw fits (dashed black lines) were performed for model M30.B in the loglog space. 

In the text 
Fig. 14.
BH masses in BHBH mergers within design advanced–LIGO sensitivity (z < 2) for models M10.B (that employs old average metallicity cosmic evolution) and M60.B (new average metallicity evolution). BH masses are calculated with the same formulas in both models. Top panel: distribution of individual BH masses, bottom panel: total BHBH system mass. All distributions are intrinsic (not redshifted nor weighted by detection probability). We note that application of new metallicity evolution with redshift (less lowZ stars: less high mass BHs) results in somewhat steeper BH mass distributions as contrasted with application of old metallicity evolution (more lowZ stars: more high mass BHs). Powerlaw fits (dashed black lines) were performed for both models in the loglog space. 

In the text 
Fig. 15.
Larger (more massive) BH mass distributions in BHBH mergers within design advancedLIGO sensitivity (z < 2) for models M60.B, M70.B, M30.B, M10.B, M50.A, M50.B. Powerlaw fits (dashed black lines) were performed for model M30.B in the top panel (∝M^{−3.6}), and for models M50.A (∝M^{−2.3}) and M50.B (∝M^{−2.4}) in the bottom panel in the loglog space. These distributions can be compared with other theoretical predictions (e.g., Fig. 1 of Mapelli et al. 2019, or Fig. 4 of Stevenson et al. 2019 that seem to be close to ∝M^{−2} − M^{−3}; we note that powerlaws shown in the latter work are not correct) or with the LIGO/Virgo observational estimate of ∝M^{−1.6} with large uncertainty on the powerlaw index: ∝M^{+0.1} − M^{−3.1} (Fig. 1 in Abbott et al. 2019a). Some of our models (e.g., top panel) show somewhat steeper relations, while other models (bottom panel) are in agreement with LIGO/Virgo estimate. 

In the text 
Fig. 16.
Detectionweighted BHBH merger total intrinsic (not redshifted) mass for models with different assumptions on the pairinstability pulsation supernovae: M70.B (strong PPSN), M60.B (moderate PPSN), M30.B (weak PPSN). For comparison, we also show LIGO/Virgo O1/O2 mean total mass estimates and their 90% confidence limits. 

In the text 
Fig. 17.
Detectionweighted BHBH merger total intrinsic (not redshifted) mass for models with different assumptions on wind mass loss and common envelope development: M50.B (weak stellar winds and optimistic CE), M50.B (weak stellar winds and pessimistic CE), M30.B (strong stellar winds and pessimistic CE). For comparison, we also show LIGO/Virgo O1/O2 mean total mass estimates and their 90% confidence limits. We note that although no model seems to reproduce the observed LIGO/Virgo BHBH total mass distribution, we can expect that some interplay of several parameters (winds, CE, PPSN, Z(z)) may, in the future, possibly reproduce the shape of the observed distribution. For example, stronger PPSN and higher Z(z) (not shown here) will tend to shift the highest BH masses to lower values. This calls for further parameter study calculations. 

In the text 
Fig. 18.
Detectionweighted distribution of effective spin parameter of BHBH mergers for model M20.B (fallback decreased BH kicks; no BH kicks for massive BHs) and M23.B (high BH natal kicks with 1D σ = 265 km s^{−1} for all BHs) with Geneva mildly efficient angular momentum transport. We note that both distributions peak at high values (χ_{eff} ∼ 0.9). Natal kicks decrease the effective spin parameter (average χ_{eff} = 0.3; M23.B) as compared to model with almost no BH kicks (average χ_{eff} = 0.7; M20.B). For comparison, we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. Although we can recover all the reported values, the predicted peak of χ_{eff} distribution is not coincident with the current LIGO/Virgo data. It indicates that BHs have typically lower spins than resulting from the Geneva model for BH natal spins, or that the detected BHBH mergers are not formed in the classical isolated binary evolution. 

In the text 
Fig. 19.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M30.B with MESA efficient angular momentum transport and fallback decreased BH kicks (no BH kicks for massive BHs). We either do not allow for efficient tidal spinup of WR stars that are the most common immediate progenitors of BHs in our models (natal BH spin is calculated directly from MESA stellar models) or we take it into account (natal BH spin is then calculated as described in Sect. 2.5 if the WR star progenitor was subject to an efficient tidal spinup). For the “no WR tides” approach we find a rather narrow distribution of effective spins (−0.2 ≲ χ_{eff} ≲ 0.2) that is peaked at positive values (average χ_{eff} = 0.15). For efficient “WR tides” the distribution is broad (−0.5 ≲ χ_{eff} ≲ 1.0) with a peak at χ_{eff} ∼ 0.15 (∼73%) and a tail with χ_{eff} ≳ 0.25 (27%). For comparison we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. 

In the text 
Fig. 20.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M40.B with Fuller superefficient angular momentum transport and fallback decreased BH kicks (no BH kicks for massive BHs). We either do not allow for the efficient tidal spinup of WR stars that are the most common immediate progenitors of BHs in our models (natal BH spin is calculated directly from MESA stellar models) or we take it into account (natal BH spin is then calculated as described in Sect. 2.5 if the WR star progenitor was subject to an efficient tidal spinup). For the “no WR tides” approach we find very narrow distribution of effective spins (−0.1 ≲ χ_{eff} ≲ 0.1) that is peaked at positive values (average χ_{eff} = 0.05). For efficient “WR tides” the distribution is broad (−0.5 ≲ χ_{eff} ≲ 1.0) with a peak at χ_{eff} ∼ 0.05 (∼78%) and a tail with χ_{eff} ≳ 0.25 (22%). For comparison, we show 90% credible limits (blue arrows) and the most likely values (blue stars) of effective spin parameter for ten LIGO/Virgo BHBH mergers. 

In the text 
Fig. 21.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M30.B with fallback decreased natal kicks (low natal BH kicks for lowmass BHs and no natal kicks for highmass BHs) and for model M33.B with high natal kicks (all BHs, independent of mass, are subject to natal kicks drawn from a 1D Maxwellian distribution with σ = 265 km s^{−1}). We note that the average effective spin decreases from model M30.B (χ_{eff} = 0.15) to model M33.B (χ_{eff} = 0.04) as an effect of natal kicks that tend to misalign BH spins and lower the effective spin (see Eq. (1)). For comparison, we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. 

In the text 
Fig. 22.
Detectionweighted distribution of the effective spin parameter of BHBH mergers for model M40.B with fallback decreased natal kicks (low natal BH kicks for lowmass BHs and no natal kicks for highmass BHs) and for model M43.B with high natal kicks (all BHs, independent of mass, are subject to natal kicks drawn from a 1D Maxwellian distribution with σ = 265 km s^{−1}). We note that the average effective spin decreases from model M40.B (χ_{eff} = 0.05) to model M43.B (χ_{eff} = 0.004) as an effect of the natal kicks that tend to misalign BH spins and lower the effective spin (see Eq. (1)). For comparison, we show the 90% credible limits (blue arrows) and the most likely values (blue stars) of the effective spin parameter for ten LIGO/Virgo BHBH mergers. 

In the text 
Fig. 23.
Specific angular momentum profile at the end of core Heburning for the 32 M_{⊙} models at Z = 0.002 calculated with a range of physical ingredients and initial rotation rates: the Geneva fast rotation model (model listed in Table A.1 with V_{ini}/V_{crit} = 40%, V_{ini} = 341.9 km s^{−1}, dotted blue line), the Geneva slow rotation model (V_{ini} = 100 km s^{−1}, solid orange line). Both these Geneva models do not include the TaylerSpruit magnetic dynamo (“noTS”). We also show: the MESA slow rotation nonmagnetic (“noTS”) model (V_{ini} = 100 km s^{−1}, solid green line), the MESA slow rotation magnetic (“TS”) model (V_{ini} = 100 km s^{−1}, dashed red line), and the MESA fast rotation magnetic (“TS”) model (model listed in Table A.2, dotdashed purple line). The brown thin dotdashed line corresponds to the specific angular momentum of the MESA fast rotation model multiplied by a factor of 10 to facilitate comparisons. 

In the text 
Fig. 24.
Comparison of the local merger rate densities of BHBH and NSNS mergers from all our models (see Table 3) with the current limits inferred from the O1/O2 LIGO/Virgo observational runs: 9.7−101 Gpc^{−3} yr^{−1} for the BHBH mergers and 110−3840 Gpc^{−3} yr^{−1} for the NSNS events (Abbott et al. 2019b). Models that are consistent with the observational limits: M13.A, M23.A, M33.A, M43.A, M30.B, M40.B, M60.B, M70.B. We note that while some of the models are consistent with the observational constraints and others are not, it is not straightforward to draw conclusions about the physical ingredients of the models at this stage due to degeneracies in the impact of various assumptions on the theoretical merger rates (see Sect. 4.3.1). 

In the text 
Fig. 25.
Comparison of the local merger rate densities of BHBH and BHNS mergers from all our models (see Table 3) with the current limits: 9.7−101 Gpc^{−3} yr^{−1} for the BHBH mergers (LIGO/Virgo O1/O2) and 1.6−60 Gpc^{−3} yr^{−1} for the BHNS events (based on the LIGO/Virgo O3 candidate of the first BHNS system, S190814bv, LVC 2019a,b; see Appendix A.9 for details). Models that are consistent with the observational limits: M13.A, M23.A, M33.A, M43.A, M23.B, M25.B, M30.B, M40.B, M50.B, M60.B, M70.B. We also note that the BHNS merger rates from our models span the range between 0.48 and 297 Gpc^{−3} yr^{−1}, which is consistent with the upper limit determined by their nondetection in O1/O2 LIGO/Virgo runs (< 610 Gpc^{−3} yr^{−1}, Abbott et al. 2019b). 

In the text 
Fig. 26.
Intrinsic merger rate distribution as a function of primary black hole mass M_{BH1} (top) and binary mass ratio M_{BH2}/M_{BH1} of BHBH mergers (bottom) for different StarTrack model choices. The gray bands show the 90% confidence limits on the populations inferred in Abbott et al. (2019a) using the phenomenological Model B. 

In the text 
Fig. 27.
Cumulative density functions (CDFs) of four StarTrack models versus effective spins χ_{eff}. Solid color lines correspond to χ_{eff} CDF of detectable events for each model. The black line shows the χ_{eff} CDF of median likelihoods of the O1/O2 LVC events, and the gray region bounds CDFs of the 5th and 95th percentiles of the likelihoods from LVC parameter estimation. In the corresponding colored shaded regions, we show the 90% range of 5000 mock observed CDFs under each model. To generate a mock observed CDF, we draw 10 χ_{eff} values from the detected population under a given model and add random Gaussian noise with σ = 0.05, which is approximately the uncertainty in χ_{eff} likelihoods from the events of GWTC1 (we resample any samples with χ_{eff} > 1 after adding noise). 

In the text 
Fig. A.1.
Specific angular momentum profile at the end of core helium burning for the 32 M_{⊙} MESA models at Z = 0.002. The models were calculated with different initial rotation rates; slow (100 km s^{−1}, red), fast (40% of critical rotation, blue) and very fast (80% of critical rotation, purple). Magnetic (dashed line) and nonmagnetic models (solid line) are considered. Large dips in the curve indicate the base of a convective region, where angular momentum transport is very efficient. The maximal extent of the curves to the right indicate the mass of the star. We can see from these curves that the faster rotating models lose more mass in winds due to rotationenhanced mass loss. 

In the text 
Fig. A.2.
Final CO core masses of BH progenitors and their corresponding BH formation masses for the components (both primary and secondary) of detectable merging BHBH systems in four of our models: M30, M50, M60, and M70 (see Table 2 for an overview of all the model ingredients). For readability, the filled areas encapsulate 96% of all the BHs, excluding few outliers and showcasing the M_{CO} to M_{BH} relation for the representative majority. The presented models are different in ways that affect the BH formation masses, that is, the PPSN/PSN prescription (see Sect. 2.3) as well as the assumed fraction of mass lost in neutrinos during a BH formation (1% in M30, M50, M70, and 10% in M60). We note that model M70 only differs from M30 and M50 for M_{CO} ≳ 32.4 M_{⊙}, whereas model M50 only differs from M30 for M_{CO} > 37.5 M_{⊙}. Those threshold core masses correspond to the threshold helium core masses M_{He} in different models PPSN/PSN prescriptions. 

In the text 
Fig. A.3.
Neutron star spins from our progenitors as a function of neutron star mass. We assume any remnant from our suite of models with a mass below 2.5 M_{⊙} is a neutron star and consider only these remnants. With the Geneva models, we have studied the angular momentum coupling of different burning layers and their effect on the neutron star spin. The Geneva models with the original mild coupling produce spin periods shorter than 1 ms (not shown here, but very close to magenta triangles; see below). Some angular momentum loss in the supernova engine (e.g., magnetic coupling such as a magnetar engine) would be required to slow neutron stars down for such models to match the data. The spins are nearly the same if the coupling extends through the CO core (magenta triangles). Coupling through the helium layer (blue squares and empty red circles; different helium core definitions) and hydrogen layer (green solid circles) produces slower neutron stars. The heliumcoupled models match the data (but realize that we are using rapidly rotating progenitors with no angular momentum loss during NS formation). The MESA models (not shown here) with the TaylerSpruit dynamo produce ∼10 ms pulsars; generating spins that match fastestspinning pulsars with rapidly spinning progenitor stars, an indication that the MESA models produce reasonable coupling results. The Fuller models (not shown here) do not produce rotating neutron stars, requiring spinup mechanisms in the supernova engine to match the data. For comparison we also mark (with black lines) the range of two observational estimates. More details on models and observations are given in Appendix A.6.1. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.