Importance of stable mass transfer and stellar winds for the formation of gravitational wave sources

The large number of gravitational wave (GW) detections have revealed the properties of the merging black hole binary population, but how such systems are formed is still heavily debated. Understanding the imprint of stellar physics on the observable GW population will shed light on how we can use the gravitational wave data, along with other observations, to constrain the poorly understood evolution of massive binaries. We perform a parameter study on the classical isolated binary formation channel with the population synthesis code SeBa to investigate how sensitive the properties of the coalescing binary black hole population are on the uncertainties related to first phase of mass transfer and stellar winds. We vary five assumptions: 1 and 2) the mass transfer efficiency and the angular momentum loss during the first mass transfer phase, 3) the mass transfer stability criteria for giant donors with radiative envelopes, 4) the effective temperature at which an evolved star develops a deep convective envelope, and 5) the mass loss rates of stellar winds. We find that current uncertainties related to first phase of mass transfer have a huge impact on the relative importance of different dominant channels, while the observable demographics of GW sources are not significantly affected. Our varied parameters have a complex, interrelated effect on the population properties of GW sources. Therefore, inference of massive binary physics from GW data alone remains extremely challenging, given the large uncertainties in our current models.


INTRODUCTION
Massive stars play an essential role in astrophysics.They are responsible for the chemical enrichment of the universe via stellar winds and supernovae.They are also progenitors of various interesting astrophysical phenomena, e.g.neutron stars, black holes, gamma ray bursts (e.g.Langer 2012).However, our understanding of these rare and short-lived objects is still incomplete.The population of merging compact binaries, observed via gravitational waves (GW) offer a unique but indirect way to study the evolution of these objects.Since the first detection of GWs, about a hundred merging binary black holes have been observed, which makes the inference of the population statistics of black hole-black hole binaries (BH-BH binaries) possible (Abbott et al. 2021b;Abbott et al. 2023).
The classical isolated binary channel is perhaps the most studied formation path.In these scenario, typically two main sub-channels are identified.The first one is the CEE channel, in which the key step in the formation of close BH-BH binaries is the so-called common envelope evolution (CEE, e.g.Ivanova et al. 2013).Several, earlier population synthesis studies predicted a merger rate for this channel that is broadly consistent with the currently inferred LIGO rate (see e.g.Mandel & Broekgaarden 2022), although these predictions are sensitively dependent on the highly uncertain common envelope efficiency and the binding energy of the envelope of the donor star.However, recent detailed stellar evolutionary models of Klencki et al. (2021) and Marchant et al. (2021) showed that the binding energy of evolved stars with radiative envelopes could be underestimated by prescriptions commonly used by rapid population synthesis codes.Furthermore, a deep convective envelope could potentially be developed at a significantly cooler effective temperature (Klencki et al. 2020) than previously assumed.These two developments would imply an appreciably lower predicted merger rate for this channel, possibly orders of magnitude lower than the currently inferred rate.
The second dominant channel (stable channel) involves two subsequent stable mass transfer episodes.The orbit of a binary experiencing a stable phase of mass transfer episode with black hole accretor can shrink significantly, if the mass ratio of the system is sufficiently high.This can lead to the formation of BH-BH binaries that merge due to GWs within the age of the universe (van den Heuvel et al. 2017).Earlier studies did not predict this formation path to be significant (see e.g.Dominik et al. 2012;Belczynski et al. 2016).However, the detailed simulations of Pavlovskii et al. (2017) showed that a stable mass transfer episode in a binary comprising an evolved donor star with radiative envelope and a BH accretor is more readily achieved than previously assumed.Subsequent studies, with assumptions that are in agreement with the findings of Pavlovskii et al. (2017), have shown that this formation path can be the dominant channel within the isolated binary scenario (Klencki & Nelemans 2018, Olejak et al. 2021;van Son et al. 2022a;Gallegos-Garcia et al. 2021;Bavera et al. 2021;Marchant et al. 2021;Andrews et al. 2021;Briel et al. 2023).
The properties of GW sources from the two aforementioned channels are sensitively dependent on various, highly uncertain binary evolutionary phases (e.g.mass transfer episodes).This, in principle means that observations of merging binary black holes (e.g.Abbott et al. 2020;The LIGO Scientific Collaboration et al. 2021) could be used to constrain massive binary physics.Unfortunately, there are currently too many uncertainties in massive stellar evolution to draw any meaningful conclusion (see e.g.Belczynski et al. 2022).Nevertheless, it is still essential to understand how uncertainties of binary and stellar physics affect the observable properties of the merging binary black hole population in order to correctly interpret the observed GW data in the future.
In this paper, we perform a parameter study on the classical isolated binary channel, using a rapid population synthesis code, SeBa (Portegies Zwart & Verbunt 1996, Toonen et al. 2012).In the first part of this paper, we study the uncertainties related to the first phase of mass transfer (i.e.mass transfer episodes between two hydrogen rich stars).For this, we test different assumptions regarding the angular momentum loss mode and the fraction of mass ejected during the mass transfer phase with a non-compact accretor.We also vary the mass transfer stability criteria of evolved stars with radiative envelopes and make different assumptions about the evolutionary stage at which giant stars develop convective envelopes to investigate the implications of studies such as Ge et al. (2015); Pavlovskii et al. (2017); Ge et al. (2020); Klencki et al. (2021).We make model variations using all possible combinations of parameter variations.This allows us to explore the interrelated effects of uncertainties.
In the second part of the paper, we investigate the effects of uncertainties in mass-loss rates of line-driven stellar winds.Both theoretical (Krtička & Kubát 2017;Sundqvist et al. 2019;Björklund et al. 2020;Björklund et al. 2022) and observational (e.g.Fullerton et al. 2006) studies suggest that mass loss rates of O/B stars could be overestimated by a factor of 2-3 by the prescription of Vink et al. (2001), which is commonly used in stellar evolutionary codes.As we will see, the impact of lowered mass loss rates on the demographics of GW sources sensitively depends on our assumptions regarding other, seemingly unrelated binary physics.Therefore, the second part of the paper shows an example of the importance of interrelated effects of uncertain parameters and it highlights the dangers of devising strategies to infer stellar physics directly from GW data without performing a full parameter study.
The impact of the uncertainties in binary physics on the isolated binary channel has been extensively studied with population synthesis approach in the recent years.However, most of such parameter studies typically concentrated on the episodes following the first phase of mass transfer.For example, the importance of common envelope evolution was investigated by varying parameters related to the common envelope efficiency (e.g.Vigna-Gómez et al. 2018;Bavera et al. 2021;Broekgaarden et al. 2022), the binding energy of the donor star (e.g.O'Shaughnessy et al. 2008;Vigna-Gómez et al. 2018;Dominik et al. 2012;Stevenson et al. 2015), and by making different assumptions on whether Hertzsprung gap donors can survive the common envelope phase (e.g.Dominik et al. 2012;Stevenson et al. 2015;Chruslinska et al. 2018;Vigna-Gómez et al. 2018;Broekgaarden et al. 2021Broekgaarden et al. , 2022)).The impact of core collapse was investigated by varying the magnitude of natal kicks received by the stellar remnant (e.g.O'Shaughnessy et al. 2008;Stevenson et al. 2015;Vigna-Gómez et al. 2018;Broekgaarden et al. 2021Broekgaarden et al. , 2022;;Ghodla et al. 2022;Richards et al. 2023), by applying different supernova mechanisms (e.g Dominik et al. 2012;Stevenson et al. 2015;Vigna-Gómez et al. 2018;Broekgaarden et al. 2021Broekgaarden et al. , 2022;;Román-Garza et al. 2021;van Son et al. 2022b), and by varying the maximum neutron star mass (e.g.Dominik et al. 2012;Stevenson et al. 2015;Broekgaarden et al. 2021Broekgaarden et al. , 2022)).Uncertainties regarding the second (stable) phase of mass transfer were studied by exploring the implications of super-Eddington accretion for BH accretors (e.g.Belczynski et al. 2020a;Bavera et al. 2021;Briel et al. 2023).
A few studies also considered some of those parameters, which are investigated in this paper.For example, the impact of mass transfer stability of donor stars crossing the Hertzsprung gap was studied in Vigna-Gómez et al. (2018); Olejak et al. (2021); van Son et al. (2022b).Furthermore, O' Shaughnessy et al. (2008); Dominik et al. (2012); Broekgaarden et al. (2021Broekgaarden et al. ( , 2022)); Belczynski et al. (2022); van Son et al. (2022b) investigated the impact of different accretion efficiencies for mass transfer episodes with non-compact accretors, while Chruslinska et al. (2018); Vigna-Gómez et al. (2018); Belczynski et al. (2022) tested different assumptions on the angular momentum loss mode during non-conservative mass transfer episodes.The importance of stellar winds on the properties of merging compact objects was also previously explored by O'Shaughnessy et al. (2008); Dominik et al. (2012); Stevenson et al. (2017); Renzo et al. (2017); Belczynski et al. (2020a); Belczynski et al. (2020bBelczynski et al. ( , 2022)); Broekgaarden et al. (2021).Clearly, many previous studies investigated the role of those parameters, which we also consider in this paper.However, this was done typically in a different astrophysical context (e.g. to study the formation of merging double neutron star binaries, see Vigna-Gómez et al. 2018;Chruslinska et al. 2018) or to investigate the formation of specific systems see e.g.Belczynski et al. 2020bBelczynski et al. , 2022)).More importantly, all of the above mentioned studies typically vary only one parameter with respect to their fiducial models and they never study systematically the importance of the first phase of mass transfer and how the related uncertainties can impact the population of GW sources.
The paper is organised as following.In section 1.1, we briefly review the classical isolated binary formation channel to introduce the terminology used in this paper.In section 2, we describe the code used in this study.In section 3, we show our results that highlight the importance of uncertainties in the first mass transfer phase.In section 4, we show how decreasing the line-driven winds for O/B stars and/or WR stars by a factor of 3 can change the properties of the merging binary black hole population.Finally, in section 5, we summarise our main findings.

The classical isolated binary formation channel
In this subsection, we provide a brief overview of the classical isolated formation channel for merging binary black holes.Our purpose is to introduce the terminology used in the rest of the paper.For a detailed overview on this subject, see e.g.Postnov & Yungelson (2014); Mandel & Farmer (2018); Mapelli (2020).
Black hole binaries with circularised orbits merge within the Hubble time (which we define here as 13.5 Gyrs), if their orbital separation is not larger than a few tens of solar radii (Peters 1964;Mandel & Farmer 2018).However, the radii of massive stars reaches orders of magnitude larger values than that during their evolution.Therefore, GW sources of the classical isolated binary channel must originate from interacting binaries.
We show a schematic drawing of the most common formation paths of GW sources according to our simulations in Fig. 1.At zero-age main sequence, the binaries are not interacting and their orbital separation widens due to stellar winds (stage 1).The initial primary star eventually fills its Roche-lobe and consequently the first phase of mass transfer is initiated (stage 2).In the dominant channels considered here, this mass transfer phase always occurs in a dynamically stable manner.This phase ends with the initial primary star losing its hydrogen envelope.At this stage, the donor star is a stripped helium star, and depending on the metallicity, it could launch intense stellar winds and therefore may be observed as a Wolf-Rayet star (stage 3).The mass ratio and the orbit of the binary system at this stage depend on how much matter was accreted by the secondary and how much angular momentum was lost by the binary during the first phase of mass transfer.The stripped helium star eventually forms a black hole (stage 4).It is currently uncertain whether this occurs via a supernova or a direct collapse.
After the inital primary forms a compact object, the secondary star expands as well and initiates the second mass transfer phase (stage 5).Based on the mass ratio of the system and the envelope structure of the donor, this episode can occur in a stable (stage 5a) or unstable fashion (stage 5b or 5c).In case of the latter, a common envelope phase is initiated (Ivanova et al. 2013).As the common envelope ensuing the binary exerts friction on the system, the period is expected to dramatically decrease.If orbital energy is not used efficiently to unbind the envelope, the change in the orbital separation could be of order ∼ 1000  ⊙ .This process can lead to an efficient formation of merging binary black holes (i.e.CEE channel).We distinguish two subtypes of the common envelope channel based on the evolutionary phase of the donor star during the second phase of the mass transfer.In the first type (5b), the donor star has a radiative envelope (rCEE), while in the second (5c), the donor star has a deep convective envelope (cCEE).The mass transfer stability criteria is sensitively dependent on whether the envelope of the donor is mostly radiative or convective.We also note that the binding energy of the envelope could be significantly different for these two types of evolved stars.
If the second mass transfer phase occurs in a stable manner, the orbit can shrink sufficiently, and thus lead to the formation of a GW source, if the mass ratios of the binary at the onset of the mass transfer phase are relatively high (stage 5a).For example, the orbit of a binary with a mass ratio  =  donor / accretor ∼ 3 and  accretor ≈ 30  ⊙ at the onset of the mass transfer phase, shrinks roughly by ∼ 100  ⊙ , assuming the accretion rate of the black hole is Eddington limited.This implies that binaries with initial orbital separations of a few ∼ 100  ⊙ can form GW sources efficiently via the stable channel.
In principle, GW sources could also be formed from systems, in which the first phase of mass transfer is dynamically unstable.However, we do not discuss this formation scenario in this paper.This is because the merger rates associated with this formation scenario are negligible in our models.
By the time the second phase of mass transfer occurs, the initial primary star is typically already a black hole.We note that this can occur for a very large fraction of the parameter space, especially if rejuvenation of the accretor is taken into account after the first phase of mass transfer (see e.g.Tout et al. 1997).For example, systems with an initial primary mass of  ZAMS,1 = 100  ⊙ can evolve in such a way, even, if their initial mass ratio are as close to unity as  ZAMS =  ZAMS,2 / ZAMS,1 = 0.99.We find that only those binaries form gravitational wave sources with non-negligble rates, in which the second phase of mass transfer occurs with a compact accretor.Therefore, in the rest of the paper, when we mention the second phase of mass transfer, we always refer to mass transfer episodes with BH accretors.
are especially relevant for this study, or which have been changed with respect to Toonen et al. (2012).The most relevant parameters related to stellar and binary physics used in all of our model variations are summarised in Table 1.

Treatment of binary interactions
In this section, we summarise how binary interactions are treated in SeBa.We assume tidal interactions circularise the orbit by the onset of the mass transfer (Portegies Zwart & Verbunt 1996).Change in the orbital separation and eccentricity due to gravitational waves emission are calculated according to Peters 1964.
A mass transfer episode occurs, if any of the stars in the binary fill their Roche-lobe.We calculate the Roche-lobe radius according to Eggleton (1983).The evolution of the orbital separation during a stable phase of mass transfer is determined as (e.g.Tauris & van den Heuvel 2006;van den Heuvel 1994): where   ,   are the mass of the donor and the accretor star, respectively,  is the mass accretion efficiency, ie. the amount of mass that is accreted and  is the ratio of specific angular momentum that leaves the system and the total specific angular momentum of the binary, ie./ tot = / tot , where J is the angular momentum of the binary and  tot is the total mass of the binary.Although,  and  is expected to depend on the parameters on the binary (see e.g.references in section 2.2), it is commonly assumed that these parameters are constant in rapid population synthesis codes (see e.g.Riley et al. 2022;Belczynski et al. 2008).If  is constant, Equation 1 can be integrated: where the subscript 'i' stands for initial, i.e. at onset of the mass transfer phase.
If the accretor is a black hole, we assume that the accretion is Eddington limited and that the specific angular momentum leaving the system is that of the accretor, which implies  =   /  (i.e.so-called isotropical reemission).The change in orbital separation for isotropical reemission in the limit of  → 0 (e.g.Soberman et al. 1997;van den Heuvel et al. 2017): We discuss the treatment of dynamically unstable mass transfer episodes in section 2.4.When the members of the binary lose mass via stellar winds, we assume that a fraction of it is accreted via Bondi-Hoyle accretion (Bondi & Hoyle 1944), while the rest leaves the system with a specific angular momentum of the donor, which corresponds to  =   /  .In this case, the orbit widens as:

First phase of mass transfer
If the accretor is not a remnant, we assume that  = 2.5, following Portegies Zwart & Verbunt (1996).We also test  = 1, following Podsiadlowski et al. (1992) and Belczynski et al. (2008).In this study, we test two, constant values for mass transfer efficiency when the accretor is a non-compact object;  = 0.3 and  = 0.7.When the accretor is a neutron star or black hole, we assume that the accretion is Eddington-limited and  =   /  .There are currently numerous uncertainties regarding mass transfer episodes with non-degenerate accretors (e.g.Langer 2012).The fraction of the transferred mass that is eventually ejected from the binary and the specific angular momentum that is removed from the system depends on many factors.For example on whether an accretion disk is formed during the mass transfer episode (e.g.Lubow & Shu 1975), on how efficiently the accretor star is spun up due to accretion (Packet 1981), on whether accretion is possible above the critical rotation of the accretor star (Popham & Narayan 1991) and on the response of the radius of the accretor star on thermal timescale (see e.g.Pols & Marinus 1994;Hurley et al. 2002;Riley et al. 2022).
Detailed binary evolution models of massive stars indicate a low mass transfer efficiency ( ≈ 0.1, on average) for mass transfer phases with evolved donors, if it is assumed that accretion is not possible above critical rotation (e.g.Langer et al. 2020).If the donor star is still on the main sequence, tides can be sufficiently strong to counteract the spinning up, leading to higher mass transfer efficiencies (see e.g.Sen et al. 2022).These findings are broadly consistent with a few observational studies (Petrovic et al. 2005;Shao & Li 2016 and see also de Mink et al. 2007).On the other hand, there are theoretical and observational studies, which conclude near-conservative mass transfers episodes among massive stars (e.g Schootemeĳer et al. 2018;Vinciguerra et al. 2020), therefore a consensus regarding this physical process is still missing.

Mass transfer stability criteria and treatment of mass transfer
We determine the stability of mass transfer with the use of the socalled mass-radius exponents (Soberman et al. 1997): where   is the radius of the donor star,   expresses how the Rochelobe radius reacts to mass overflow, while  ad and  th expresses how the radius of the donor reacts to mass loss during mass transfer on dynamical, and thermal timescale, respectively.Three different mass transfer modes can be distinguished: stable mass transfer on nuclear time scale (  ≤ min( ad ,  eq )), stable mass transfer on thermal timescale ( ad ≥   ≥  eq .)and unstable mass transfer (  > max( ad ,  eq )).In the first two cases, we assume that the mass transfer rate is  =   /, where  is the nuclear timescale in the first and the thermal timescale in the second case.If the mass transfer is dynamically unstable, we assume common envelope evolution (see section 2.4).
As a major simplification, we assume a constant  ad and  th for a given stellar evolutionary phase (these are summarised in Table 2).Giants with deep convective envelopes tend to have low  ad , possibly even negative.Therefore, donor stars of this type are likely to initiate unstable phases of mass transfers.At what stage the deep convective envelope develops in massive stars is still very uncertain.It is common to use effective temperature as a proxy for the evolutionary stage at which such an envelope is developed (we will note this as  eff,boundary ).We test two assumptions.First, that a deep convective envelope develops at an effective temperature logT eff = 3.73 K, following Ivanova & Taam (2004) and Belczynski et al. (2008).In the second model variation, we follow the prescription of Klencki Table 1.Summary of the most important parameters in our rapid population synthesis simulation with SeBa.The top of the table shows the standard model used in this paper, while the bottom of the table shows the model variations.We run simulations on a metallicity grid Z = 0.0001, 0.0003, 0.0005, 0.0007, 0.001, 0.003, 0.005, 0.01, 0.02 with all possible combinations of model variations unless stated otherwise in section 2. For the detailed descriptions of stellar wind models see Table 3. )  wind =  wind,WR = 1 and  LBV = 1.5 Model I 2.)  wind = 1/3,  wind,WR = 1 and  LBV = 1.5

Common envelope evolution
In this study, we model common envelope evolution by adopting the energy formalism (e.g.Webbink 1984;van den Heuvel 1976).As the common envelope engulfs and exerts friction on the binary, the orbital separation starts to shrink.It is assumed that a fraction (  ) of the energy liberated from the orbital energy is used to unbind the envelope.Then the orbital separation by the end of the CEE phase can be given as: here, the left hand side term is the binding energy of the envelope. d,core is the mass of the helium core of the donor,   is the radius of the donor,   and   are the initial and final orbital separation, respectively The term  describes the structure of the envelope (de Kool 1990).Several studies published tabulated or fitted data for  for a range of masses and evolutionary stages (e.g.Dewi & Tauris 2000;Xu & Li 2010;Loveridge et al. 2011;Claeys et al. 2014;Kruckow et al. 2016;Klencki et al. 2021).The values of these calculated  parameters can vary over orders of magnitude, depending on the radius and the mass of the donor star and on the metallicity.Klencki et al. (2021) and Kruckow et al. (2016) predict, however, that for sufficiently massive stars ( ZAMS ≳ 30  ⊙ ),  varies only by a factor of a few (∼2-5) as a function of stellar parameters and metallicity, once its radius has expanded sufficiently ( ≳ 500-1000  ⊙ ), given that the star still has mostly radiative envelopes (see e.g.Fig. 1 in Kruckow et al. 2016 andFig. C.3 in Klencki et al. 2021).
In this study, we assume a constant  = 0.05, which is in reasonable agreement with the results of Kruckow et al. (2016) for stars with  ZAMS ≳ 30  ⊙ and  ≳ 500  ⊙ .We note that we assume that Hertzsprung gap donors cannot survive CEE episodes (Dominik et al. 2012).This also implies that donor stars in GW progenitors typically have  ≳ 500  ⊙ at the onset of the CEE in our simulations.We assume   = 1.

Mass transfer episode types based on the evolutionary phase of the donor
We distinguish the following mass transfer phase types based on the evolutionary stage of the donor star.
Case A: the donor is a main sequence star (see e.g.Sen et al. 2022).If the period of the system is sufficiently short, the secondary might also fill its Roche-lobe, leading to the formation of contact systems (Pols 1994;Wellstein et al. 2001;Menon et al. 2021).The outcome of Case A mass transfer phase is expected to be very different from those which start with an evolved giant donor (i.e Case B and Case C).As opposed to giants, main sequence stars do not have fully developed helium cores.During Case A mass transfer, fusion in the developing helium core is expected to halt because of the rapidly dropping central temperatures.Consequently, the mass of the naked helium star that is left after the end of Case A mass transfer phase is lower than for binaries experiencing mass transfer episodes with evolved donor stars (e.g.see Langer et al. 2020).In populaton synthesis codes like SeBa, it is challenging to model Case A mass transfer episodes accurately.There are two major reasons for this.Firstly, the stellar tracks of Hurley et al. (2000) do not track the mass of the developing helium core on the main sequence.The core mass is only determined at the start of the Hertzsprung gap phase.Secondly, a constant mass-radius exponent is assumed for a given stellar evolutionary phase (see subsection 2.3).This means that the radius response of the donor during a Case A mass transfer phase is assumed to be the same, regardless of how much mass had already left the star (i.e.even, if in principle most of the hydrogen rich mass had already left the star).The consequence of these two points is that, the amount of mass that is transferred to the accretor can be significantly overestimated, and the mass of the black hole that the donor eventually forms can be severely underestimated by codes like SeBa.For a different approach in a binary population synthesis code, see e.g.Agrawal et al. (2023) or codes that use detailed interacting binary stellar models, such as BPASS (Eldridge et al. 2017;Stanway & Eldridge 2018), Brussels code (see e.g.Vanbeveren et al. 1998a;De Donder & Vanbeveren 2004), ComBinE (Kruckow et al. 2018) or POSYDON (Fragos et al. 2023).In such codes a more physically motivated modelling of Case A mass transfer episode is possible.To conclude, the outcome of Case A mass transfer episodes predicted by codes based on the stellar tracks of Hurley et al. (2000) should be treated with caution.Nevertheless, we still show such systems in this work for completeness.
Case B: the donor star is crossing the Hertzsprung gap.During this evolutionary phase of the donor star, a large and rapid increase in stellar radius occurs.The mass transfer phase ends with the donor star losing its hydrogen envelope, leaving a naked helium star behind (but see Laplace et al. 2020).Since the Hertzsprung gap phase lasts only for a few ∼ 10 4 years, the helium core does not have sufficient time to significantly grow and therefore the mass of the helium star (and therefore the black hole that is eventually formed) is not strongly dependent on the exact initial separation of the binary.Furthermore, these systems are not significantly affected by LBV winds (assuming steady mass loss rates).As noted in section 2.4, we follow Dominik et al. (2012) and assume that binaries with donor stars crossing the Hertzsprung gap cannot survive CEE.
Case C: the donor star is in its core helium burning phase.We distinguish two sub-categories: (i) Case Cr: core helium burning donor star with radiative envelope, (ii) Case Cc: core helium burning donor star with deep convective envelope.We assume that the mass transfer stability criteria is the same for Case B and for Case Cr mass transfer phases (see Table 2).Yet, there are important differences in the predicted outcome of these two episodes, as the core helium burning phase lasts orders of magnitude longer with a slower expansion in radius.This means that for Case Cr mass transfers, the mass of the remnant that the donor star eventually forms is sensitively dependent on the initial separation, since the mass of the helium core can grow significantly during the core-helium burning phase.Moreover, the effects of LBV winds are no longer negligible.Donor stars with deep convective envelopes have very different envelope structures and therefore different mass transfer stability criteria.As already mentioned, unstable mass transfer phases are more readily realised for these systems (see Table 2 and B1).

Supernova and natal kicks
The mass of the remnant after core collapse is computed based on the delayed supernova model from Fryer et al. (2012).This prescription determines the remnant mass as a function of CO core mass, which in SeBa is obtained from the fits of Hurley et al. (2000).The kick velocity for black holes is calculated as : Where   is the fallback,    is the canonical neutron star mass    = 1.4 ⊙ and  kick is a random velocity kick drawn from the distribution inferred by Verbunt et al. (2017) from proper motion measurements of pulsars.The distribution of Verbunt et al. (2017) is a combination of two Maxwellian functions with velocity dispersions of  = 75 km/s, and  = 315 km/s, and weights of 0.42 and 0.58, respectively.

Stellar wind prescriptions
Massive stars lose a substantial fraction of their mass via stellar winds.We can roughly group stellar wind mechanisms into three groups; line-driven winds (which also includes Wolf-Rayet winds), winds of Luminous blue variables (LBVs), and dust-driven winds.Line-driven winds can be further distinguished based on whether they are optically thin (typically stellar winds of main sequence and evolved stars with hydrogen envelopes) or they are optically thick (line-driven winds of stripped helium stars, i.e.Wolf-Rayet winds).
In the following, we briefly summarise the stellar wind prescriptions we use in this study, while Table 3 shows at what evolutionary stage these prescriptions are applied.Line driven winds of O/B stars are modelled in SeBa using the mass-loss rates from Vink et al. (2001), as long as the star is within the grid defined by Vink et al. (2001), otherwise, the empirical formula from Nieuwenhuĳzen & de Jager (1990) is used.These mass-loss rates are applied until stars reach  eff ∼ 8000 (see Table 3).
We note that Vink et al. (2001) estimates the global metallicity dependence to be  ∝  0.69 for O stars, by assuming a metallicity dependence of the final wind velocity to be  ∞ ∝  0.13 , following Leitherer et al. (1992).This is consistent with the observations as shown by Mokiem et al. (2007), although these results still depend on the findings of Leitherer et al. (1992).We note that however, Krtička (2006) and Björklund et al. (2020) find weaker metallicity dependence of escape velocity, where in the latter, the exponent can even be negative for high luminosity stars.In SeBa, we assume  ∝  0.85 (i.e.ignoring the metallicity dependence of  ∞ ), so that our choice of modelling is consistent with other population synthesis studies (e.g.Dominik et al. 2012;Mapelli 2016;Giacobbo et al. 2017;Stevenson et al. 2017, but see e.g.Eldridge & Stanway 2016 in which  ∝  0.69 is used).We assumed that solar metallicity is Z = 0.02.
If a main-sequence star is outside of the grid defined by Vink et al. (2001), we apply the empirical formula of Nieuwenhuĳzen & de Jager (1990) and assume a metallicity scaling  ∝  0.85 .We note that this is different from the what was originally suggested by Kudritzki et al. (1987), which was ∼  0.5 .We do this so that there is a consistent metallicity dependence for optically thin line-driven winds.
LBV stars, stars beyond the Humphreys-Davidson limit, experience very high mass loss rates, in the order of  LBV ∼ 10 −5 -10 −3  ⊙ yr −1 , although this is highly uncertain.Even less is known about their possible eruptions, in which huge amount of mass could be lost in a very short time (Humphreys & Davidson 1994;Vink 2012;Smith 2014).For the mass-loss rates of LBV stars, we follow the assumption of Belczynski et al. (2010), i.e. the mass loss rate is constant and has a value of  LBV = 1.5 If the star becomes a cool giant (ie  eff ≤ 8000), we calculate the mass loss rate according to Reimer's empirical formula (Reimers ).We compare it with the mass loss from Nieuwenhuĳzen & de Jager (1990) and we take the maximum value.For the mass loss rates of thermally pulsating AGB stars, we use the prescription of Vassiliadis & Wood (1993).Mass loss rates for helium stars and helium giants (Wolf-Rayet star winds) are calculated according to Sander & Vink (2020).
In section 4, we investigate the impact of uncertainties of these mass loss rates on the demographics of interacting massive binaries and GW sources.We test three stellar wind models, which are the following: 'Model I' with our standard stellar wind model (  wind =  wind,WR = 1), 'Model II' with the optically thin line driven winds scaled down by a factor 3 (  wind = 1/3,  wind,WR = 1) and 'Model III' where besides the optically thin line driven winds, Wolf-Rayetlike winds are also scaled down (  wind =  wind,WR = 1/3).For the exact description of our applied stellar winds prescriptions and scaling factors, see Table 3.With Model II, we aim to study the implications of Krtička & Kubát (2017); Sundqvist et al. (2019); Björklund et al. (2020) in a simplified way.These studies found that the prescription of Vink et al. (2001) systematically overpredicts the mass loss by a factor ∼ 2-3.With Model III.we investigate the general uncertainties in the mass loss rates of stripped helium stars (see e.g.Sander & Vink 2020).The latter becomes especially significant for interacting binaries, as envelope loss due to a mass transfer episode can significantly increase the time spent as a Wolf-Rayet star.

Initial conditions
Observations suggest that about half of the stars are in binaries (or in higher order, hierarchical systems) and this multiplicity fraction increases with increasing mass (Duchêne & Kraus 2013).In particular, Sana et al. (2012) showed that the binary fraction reaches   ≃ 0.7 for stars in the mass range  ≃ 15-60  ⊙ .The same observations showed that the orbital period distribution of these young, massive binaries favour short period systems: According to Sana et al. (2012), the mass ratio distribution of massive stars is nearly flat,   ∝  −0.1 , while the distribution of eccentricities follows   ∝  −0.42 .A similar mass ratio distribution is inferred by Kobulnicky et al. (2014), while Dunstall et al. (2015) finds   ∝  −2.9 for B stars in the Tarantula Nebula, and Sana et al. (2013) finds   ∝  −1 or O stars in the Tarantula Nebula.
The effects of uncertainties in initial conditions were investigated by de Mink & Belczynski (2015), and they found that only the changes in the initial mass function alters significantly the properties of the merging compact binary object population (but see Klencki et al. 2018).This is an important result, as there are observational and theoretical indications that the initial mass function might not be universal (see Chruślińska et al. 2020, and references therein).In our simulations, we assume these initial distributions are not correlated, although, this might not be a valid assumption (see e.g.Moe & Di Stefano 2017;Klencki et al. 2018).
In this paper, we apply the following initial conditions for our simulations: • Initial mass function: we assume a universal initial mass function of the primary star of Kroupa (2001) in a mass range of ZAMS .• Initial mass ratio distribution: we assume a uniform mass ratio distribution between 0.1 and 1, where the mass ratio is defined as  =  ZAMS,2 / ZAMS,1 , that is the ratio of the mass of the initially secondary and the mass of the initially primary star.A flat distribution is in a reasonable agreement with observations of Sana et al. (2012).
• Initial separation distribution: we assume a flat distribution in the logarithmic space of binary separation in the interval of 1 ⊙ -10 4  ⊙ ;   ≃ log(a).This is equivalent to Opik's law (Öpik 1924), ie. a uniform distribution in log(p).We note that since we sample from the distribution of semimajor axis, some of our systems have sub-day periods.However, we discard any systems that fill their Roche-lobe at zero-age main sequence and we do not take them into account for calculating event rates (see section B).

Simulation setup
In section 3, we explore the impact of uncertainties of first phase of mass transfer on the merging binary black hole population.We do this by varying four parameters in our simulations, namely: (i) , which expresses the specific angular momentum lost from the binary during a non-conservative mass transfer phase as a fraction of the total specific angular momentum of the binary.We vary  only for non-compact accretors.
(ii) , which describes the mass transfer accretion efficiency.We vary  only for non-compact accretors.
(iii)  ad,rad , which is the mass-radius exponent for giant donors with radiative envelopes.It determines the boundary between stable and unstable mass transfer episodes.
(iv)  eff,boundary , which is the effective temperature at which a deep convective envelope is expected to develop.
The details of our model variations are summarised in Table 1.Using only our standard stellar wind model, we run simulations with all possible combinations of ,   ad,rad and  eff,boundary (i.e 16 models in total, named M1..M16, see e.g.Table 4).We simulate 10 6 binaries at each value of our metallicity grid for each model variations.Our metallicity grid is defined as  = [0.02,0.01, 0.005, 0.003, 0.001, 0.0007, 0.0005, 0.0003, 0.0001].To make comparisons with other studies easier, we have converted  to the critical mass ratio for each model variation in table 𝐵1.We calculate the merger rate density in the local universe, as summarised in section B.
In section 4, in which we investigate the impact of uncertainties in mass loss rates due to stellar winds, we perform binary evolution simulations for 10 6 systems at metallicities Z = 0.02, 0.01 and 0.005.We test three stellar wind models (e.g.Model I, II, III, see section 2.7) and vary  ad,rad and , while assuming the convective envelope prescription of Ivanova & Taam (2004) and  = 0.3.
We note, however, that with the metallicity specific star formation rate model that we assume in this study (see equations B5 and B6), GW progenitors formed at  ≳ 0.005 do not contribute to the merger rate density in the local universe significantly.Consequently, in our models, different assumptions about mass loss rates of line-driven winds (which are only relevant at  ≳ 0.005) also do not affect the demographics of GW sources at  ∼ 0 significantly.However, whether merging binary black holes with masses  BH ∼ 30  ⊙ could be formed in an environment that is typical for the LMC and the SMC or even in the Milky Way (i.e.Z∼0.02-0.005)remains an important and open question (see e.g.Srinivasan et al. 2023).Furthermore, the models for metallicity-specific star formation rates, as well as the metallicity dependence of stellar winds are highly uncertain (see e.g.Chruślińska 2022).

RESULTS: THE IMPACT OF UNCERTAINTIES IN STABLE MASS TRANSFER
In Fig. 2, we show the primary mass distribution of merging binary black holes in the local universe (z = 0) for the three dominant channels: stable, rCEE, cCEE for all of our model variations.In Table 4, we show the corresponding predicted cosmic merger rate densities.
The total predicted merger rate densities ( total ) are in a broad agreement with the currently inferred rate from LIGO and Virgo observations (which is  GWTC3 = 28.3+13.9 −9.1 Gpc −3 yr −1 Abbott et al. 2021a). total of model variations with  ad,rad = 4 are within a factor of two of this inferred value, while  total of model variations with  ad,rad = 7.5 are larger than the observed rate by a factor of 3-4.In Fig. E1 and section E, we provide a comparison between the primary mass distribution of our models and inferred distribution from GWTC-3 Abbott et al. (2023).
Fig. 2 clearly demonstrates that the relative rate of each formation channel varies significantly with different model variations.Parameters  and  ad,rad have the largest impact.The stable channel dominates in 11 out of 16 model variations, and the cCEE channel dominates in the remaining 5. The rCEE channel is non-negligible only in 4 model variations.In 12 model variations, the vast majority of the most massive systems originate from the stable channel (in agreement with van Son et al. 2022a;Briel et al. 2023).
While GW sources form most efficiently via the CEE channel in models M1-M3 (see e.g. ad,rad = 4 and  = 2.5 models in Fig. E2),  total is still dominated by the stable channel in these model variations.This is due to the relatively long formation times associated with the stable channel and the monotonic increase of the cosmic star formation rate up to  ∼ 2 (see e.g.Madau & Dickinson 2014).Most of the sources of the stable channel are formed at higher redshifts than the sources of the CEE channel, and at these higher redshifts the star formation rate is also higher, which leads to an increased merger rate.(see similar conclusions in Neĳssel et al. 2019;van Son et al. 2022a;Briel et al. 2022).
Fig. 3 shows the predicted mass ratios distribution of GW sources and its dependence on primary BH mass.For relatively massive systems ( BH ≳ 20  ⊙ ), this distribution is dependent on the assumptions on the uncertain binary physics.The formation channels can yield a population of merging binary black holes with a relatively narrow range in mass ratios (0.8 ≲  ≲ 1.0) or a moderately wide mass ratio distribution (0.5 ≲  ≲ 1.0).The mass ratio distributions of the most massive GW sources of the stable channel is most sensitive to the assumed mass transfer stability parameter.In models with  ad,rad = 4, the typical mass ratios are between 0.8 ≲  ≲ 1.0, while in models with  ad,rad = 7.5, the mass ratio distribution becomes much broader, i.e. 0.5 ≲  ≲ 1.0.We notice smaller variations in the mass ratio distribution for the CEE channel with different assumptions in binary physics, though model variations with  = 1 produce somewhat broader mass ratio distributions than models with  = 2.5.The less massive systems generally have much wider mass ratio distributions in all model variations.
Overall, we do not notice significant variations in the merger rate density, BH mass range, or the shape of the mass distributions across different model variations.Our results, therefore, show that while the main GW observables do not depend sensitively on the uncertainties studied here, the relative importance of dominant channels do.In other words, the formation paths of the majority of GW progenitors can be entirely different depending on the assumptions on how the first phase of mass transfer proceeds and yet the predicted demographics of merging binary black holes are very similar (though we neglect spins in this study, see e.g.Bavera et al. 2021).This highlights why it is extremely challenging to infer physics of massive binary evolution solely from GW observations, given the huge uncertainties in the current models.
Next, we briefly summarise the most important ways how uncertainties related to the first phase of mass transfer can affect the relative importance of the dominant GW formation channels, while in sections 3.1 -3.4,we discuss these effects in detail.
(i) The assumed angular momentum loss mode has a strong impact on which formation channel dominates.Typically, the stable channel dominates in models with  = 2.5 (see e.g.first row in Fig. Table 4.A summary of our predicted rates for all of our model variation with our standard stellar wind models (i.e.Model I) 2), while the cCEE channel dominates, if  = 1, however the latter depends on the assumed  ad,rad as well (compare e.g.third and fourth row in Fig. 2).
(ii) An increased  ad,rad leads to a higher merger rate of the stable channel.However, this increase is only significant for relatively lower mass black holes ( BH,1 ≲ 20  ⊙ ).With lower mass transfer efficiencies (e.g. = 0.3), we find that the merger rate of binary black holes with  BH,1 ≳ 20  ⊙ is not affected at all, unless  = 1.0.
(iii) The merger rate of the rCEE channel is only non-negligible with the convective envelope prescription of Klencki et al. (2020).Furthermore, the GW sources of this channel comprise of relatively lower mass BHs (i.e  BH,1 ≲ 20  ⊙ ) in all of our model variations (see section 3.4).

The impact of angular momentum loss mode on the stable channel
As evident from Fig. 2, the merger rate density of the stable channel ( stable ) depends sensitively on the assumed angular momentum loss mode during the first phase of mass transfer (e.g.compare first with third row).In our  = 2.5 models, this formation channel is particularly efficient and is typcially the dominant out of the three main formation paths considered here.On the other hand, in the  = 1 models, the merger rate of this channel is negligible, unless the mass-loss exponent is increased, i.e.  ad,rad = 7.5.In order to understand the reason for the relation between  and  stable , we need to consider the following three important features of this formation channel: (i) The progenitors of the stable channel have relatively short initial orbital separations.For these systems, the largest values of  ZAMS is a few hundred solar radii (shown in Fig. 4).If  ZAMS is much larger than that, the orbit does not shrink sufficiently by two phases of stable mass transfer, such that a merging binary black hole would be formed.
(ii) There is a minimum initial orbital separation for the binaries of the stable channel.We identify a minimum  ZAMS associated with this channel, below which binaries typically do not form GW sources.This minimum  ZAMS roughly coincides with the separation at which the initial primary would fill its Roche-lobe just when it evolved off MS (see Fig. 4).If  ZAMS is below this minimum value, the first phase of mass transfer (Case A) leads to a stellar merger eventually in the vast majority of cases (shown in Fig. 4).As explained in detail in section C3, even if the initial orbit is wide enough for the system to survive the Case A mass transfer phase, the binary still likely merges in the subsequent mass transfer event.In this case, the first phase of Case A mass transfer leads to a less massive black hole from the initial primary (because the formation of the helium core is halted) and to a more massive secondary (because more mass is transferred to the accretor) with respect to systems experiencing a mass transfer episodes at later stages.This results in a relatively high mass ratio at the onset of the second phase of mass transfer ( MT,2 ) and such systems merge as a result of a dynamically unstable second phase of mass transfer (i.e.typically  MT,2 ≳  crit with  ad,rad = 4).To conclude, systems with a first phase of Case A mass transfer typically merge before forming BH-BH binaries (see also Gallegos-Garcia et al. 2022) and therefore, there exists a minimum initial separation associated with the stable channel. .
(iii) The lower the angular momentum loss is, the wider the orbit becomes during the mass transfer phase.This means that the orbital separation of binaries with  = 1 typically widens more during the mass transfer phase than with  = 2.5.The degree by which the orbit changes due to a first phase of (stable) mass transfer is primarily determined by the initial mass ratio ( ZAMS ),  and  (see right panel of Fig. C1).Systems that form GW sources via the stable channel typically have initially near equal masses (0.7≲  ZAMS ≲1.0, but see e.g.Gallegos-Garcia et al. 2021).This is because such binaries develop sufficiently large mass ratios by the onset of the second phase of mass transfer (see left panel of Fig. C1), and therefore experience efficient orbital shrinking during a stable phase of mass transfer with a black hole accretor (see equation 3).For these initial mass ratio ranges, the orbit of the binary significantly widens in the  = 1 models, while the net change in the orbital separation is very small in the  = 2.5 models (see right panel of of Fig. C1).We note that, if  ZAMS is very small, the orbit shrinks due to the first phase of mass transfer, even with  = 1, but in that case  MT,2 will be small and therefore the orbit will widen due to the second (stable) phase of mass transfer and no GW source will be formed (compare left and right panels of Fig. C1).
Considering these three points, it is possible to understand why the stable channel is inefficient with  = 1.The orbit widens too much due to the first phase of mass transfer, even for those binaries that start out with the minimum  ZAMS associated with this channel.The mass distribution of the primary of the merging binary black hole in the local universe (z = 0) for all of our models with our standard stellar wind models (e.g. wind =  wind,WR = 1, see Table 1).The distribution is shown by stacked histograms, where each colour indicates a different formation channel.In the legend on the upper right corner in each panel, we show the percentage of each channel.In Fig. E1, we provide a comparison between these distributions and the primary mass distribution of BH-BH binaries as inferred from the the third LIGO-Virgo Gravitational-Wave Transient Catalog (GWTC-3) (Abbott et al. 2023) Therefore, the binary black holes that eventually form this way are too wide to merge due to GWs within the Hubble time.
As shown in Fig. 2, stable channel can be efficient with  = 1, if  ad,rad = 7.5.In this case, the significant orbital widening due to the first mass transfer phase can be counteracted by the effect of the second phase of the mass transfer for sources with relatively large mass ratios at the onset of the second mass transfer (3.2 ≲  MT,2 ≲ 5.5, see equation 3).Such sources can form, if the first mass transfer phase is Case B and the mass transfer efficiency is relatively large (e.g. = 0.7) or if the first mass transfer is (very late) Case A, which typically leads to large values of  MT,2 .
The relationship between  stable and  is therefore dependent on the predicted outcome of Case A mass transfer episodes.We emphasise again our point in section 2.1, that the treatement of Case

Secondary mass [M ]
Figure 3. 2D Histograms that show the masses of the merging binary black holes for a few selected models at Z = 0.0007.Each histogram has been normalised to the merger efficiency of their own model variation.The dashed lines show constant mass ratios.We also show the observed GW detections from the third LIGO-Virgo Gravitational-Wave Transient Catalog (GWTC-3) and associated measurement uncertainties (Abbott et al. 2023).We do not show those detections, where any of the inferred merging black hole mass is above  BH ≥ 45  ⊙ .
A mass transfer in stellar evolutionary codes based on Hurley et al. (2000), is extremely simplified and its predictions should be treated with caution.While we should expect that the prediction of large values  MT,2 following a Case A mass transfer phase is qualitatively true (and theferore a minimum  ZAMS could indeed exist for the stable channel), whether indeed the vast majority of them would be above  crit should be further investigated.

The impact of angular momentum loss mode on the CEE channel
Fig. 2 shows that the efficiency of the cCEE channel is also sensitively dependent on the assumed , although in an opposite way as for the stable channel. cCEE is about a factor of 8 larger with  = 1 than with  = 2.5 in our low mass transfer efficiency models.With  = 0.7, this difference is about a factor of two.
Below we explain the reason for this relationship.The binaries of the cCEE channel have  ZAMS of a few thousand solar radii.Only in this case, the binaries are sufficiently wide by the onset of the second phase of mass transfer, such that the donor star (i.e. the initial secondary star) fills its Roche-lobe with a deep convective envelope.
In general, binaries have wider orbital separations at the onset of the second phase of mass transfer ( MT,2 ) with  = 1 than with  = 2.5.Consequently, there are significantly more systems in the latter case for which the second phase of mass transfer is Case Cc (compare the top panels of Fig. C3).This also leads to a higher  cCEE , since in this channel the second phase of mass transfer is by definition Case Cc (see section 1.1).There are two reasons for this: (i) The rate of unstable first phase of mass transfers of the widest interacting binaries sensitively depends on angular momentum loss.The binaries with the longest periods that still exchange mass engage in Case Cc first phase of mass transfer (see Fig. 4).As shown in Fig. 4 Our simulated binaries in the initial mass ratio -initial separation space.We show the orbital separation after tidal circularisation.The systems are simulated with the following parameters;  = 0.3,  ad,rad = 7.5.Furthermore, we have used the assumption of Ivanova & Taam (2004) regarding the development of deep convective envelopes.In the left panel we show systems with  = 2.5, in the right with  = 1.We show the outcome of each simulated binary with difference colours.We distinguish the following types; (i) interacting, detached: these binaries interact via mass transfer phases at least once and form detached BH-BH binaries which do not merge due to gravitational waves within Hubble time, (ii) interacting, disrupted: these binaries interact at least once and get disrupted due to the supernova kick of one of the binary components, (iii) interacting, stellar merger: these systems interact at least once and merge due to one of their mass transfer phases.Their first mass transfer phase is typically stable or unstable with a MS or HG donor, (iv) non-interacting: none of the stars fill their Roche-lobes, they either stay detached or get disrupted, (v) CEE in 1st MT with CHeB donor: these binaries have a core-helium burning donor at the onset of the first mass transfer phase, and they initiate a dynamically unstable mass transfer (vi) GW source, stable channel: gravitational wave progenitors evolving via two subsequent stable mass transfer phases, (vii) GW source, CEE channel: gravitational wave progenitors evolving via a stable and an unstable mass transfer phase.
80 per cent) of Case Cc episodes occur in an dynamically unstable way, if  = 2.5.For these binaries, the orbital separation drastically decreases due to the first phase of unstable mass transfer and these binaries typically do not form GW sources.On the other hand, the majority of the same mass transfer episodes are stable with  = 1, and as a result, the orbit typically widens for these systems.This leads to larger values of  MT,2 and consequently a higher rate of Case Cc second phase of mass transfers compared to the models with  = 2.5.Therefore, the critical mass ratio associated with a first phase of Case Cc mass transfer is sensitively dependent on the assumed  and , (see Fig. 4 and Table B1).This can be understood by considering that with larger angular momentum loss, the orbit shrinks at a faster rate at the beginning of the mass transfer phase and therefore a dynamically unstable mass transfer phase is more easily instigated.
(ii) The periods of binaries engaging in a first phase of Case Cr mass transfer tend to increase more with lower angular momentum loss.As a result, the number of binaries that have large  MT,2 is higher in the models with  = 1 with respect to models with  = 2.5, and consequently, so is the rate of a second phase Case Cc mass transfer, since in wider binaries, the donor star fills its Roche-lobe at a later evolutionary stage and therefore it is more likely that this occurs when the star has already developed a deep convective envelope.

The impact of mass transfer stability parameter
Fig. 2 shows that  stable increases and  rCEE decreases with increasing  ad,rad (compare row 1 and 2, or row 3 and 4 of Fig. 2).This effect is not surprising; larger  ad,rad translates to larger  crit , which implies a larger parameter space for stable mass transfer episodes, in case the donor star is evolved and has a radiative envelope.At the same time, the degree by which the orbital separation shrinks due to a stable phase of mass transfer significantly increases with increasing  MT,2 (see e.g equation 3), which results in an efficient formation of GW sources via the stable channel.For example, with  MT,2 = 3.2 (i.e. crit at  ad,rad = 4), the orbit shrinks typically by ∼ 200  ⊙ due to a stable phase of mass transfer with a BH accretor, while the same with  MT,2 = 5.5 (i.e. crit at  ad,rad = 7.5) is ∼ 1000  ⊙ .This means that the orbital shrinking due to a stable mass transfer with mass ratios  MT,2 ≳ 5 can be as efficient as due to common envelope evolution in our models.
However, Fig. 2 also shows that only the merger rate of relatively lower mass merging BH-BH binaries (i.e. BH,1 ≲ 20  ⊙ ) is significantly affected by an increased  ad,rad (and  = 2.5).In particular, when we increase  ad,rad from 4 to 7.5, the merger rates of systems with  BH,1 ≳ 20  ⊙ remains practically unchanged in our low mass transfer efficiency model, and increases only by a factor of 3 in our high mass transfer efficiency models (while the  stable for the entire mass range increases almost by a factor of 9).The reason for this can be understood by inspecting Fig. C2.At low metallicities, where we expect the vast majority of the GW progenitors to originate from, the maximum  MT,2 that binaries engaging in a first phase of Case B mass transfer develop significantly decreases with increasing  ZAMS,1 .For example, at Z = 0.0007 and for binaries with  ZAMS,1 ≳ 60  ⊙ , this maximum of  MT,2 is about 3 and 3.5 for  = 0.3 and  = 0.7, respectively.Consequently, increasing  ad,rad from 4 to 7.5 (corresponding to increasing  crit from 3.2 to 5) does not have a significant effect on  stable for the most massive GW progenitors.An exception to this can bee seen in the model variations with  = 1.In these models, the most massive GW sources of the stable channel have (a very late) Case A first phase of mass transfer, and develop  MT,2 ≳3.2.In this case, increasing  ad,rad affects the merger rate of the most massive binary black holes too (see last two panels in the 4th row of Fig. 2)

The impact of different convective envelope prescriptions
In our models with the convective envelope prescription of Klencki et al. (2020), the maximum primary mass of the BH-BH binaries of the cCEE channel is about  BH,1 ≈ 20  ⊙ (see e.g.first two columns of Fig. 2).This is about a factor of two lower than with the prescription of Ivanova & Taam (2004).This is because the models of Klencki et al. (2020) predict that most massive stars (i.e. ZAMS ≳ 50  ⊙ ) either never develop deep convective envelopes or they do at a late evolutionary stage, which is not followed by a significant expansion of stellar radius.This implies that the rate of Case Cc mass transfer episodes with donor stars with  ZAMS ≳ 50  ⊙ is negligible and therefore so is The merger rate of rCEE is negligible with the prescription of Ivanova & Taam (2004).In these models, stars develop a deep convective envelope soon after the onset of core-helium burning.Therefore, the occurrence rate of a second phase of unstable Case Cr mass transfer is very low.This is no longer true for the models with the prescription of Klencki et al. (2020), in which the rCEE channel can account for up to 42 per cent of all mergers (see M10 in Fig. 2).However, the primary BH mass is always  BH,1 ≲ 20  ⊙ for this channel.This can be understood again by inspecting Fig. C2.The maximum value of  MT,2 is either below or slightly above  crit for the most massive systems experiencing a Case Cr mass transfer phase (depending on the assumed mass transfer efficiency).Therefore, very few binaries with  1,ZAMS ≳ 50  ⊙ initiate an unstable Case Cr phase.We can also see that the values of  MT,2 are considerably lower for Case C than for Case B mass transfer episodes.This is due to the strong LBV winds that decrease the mass ratios of the binaries over time.Therefore, metallicity independent LBV winds also contribute to the low  rCEE in our models.

Comparison to earlier studies
Relatively early rapid population synthesis studies predicted that merging binary black holes overwhelmingly originate from the CEE channel (see e.g.Dominik et al. 2012;Belczynski et al. 2016;Stevenson et al. 2017), while the contribution from the stable channel is negligible.On the other hand, Neĳssel et al. (2019) found that the stable channel is the dominant source of merging binary black holes.Qualitatively similar results were found by several subsequent studies (i.e.Gallegos-Garcia et al. 2021;Olejak et al. 2021;van Son et al. 2022a).A common interpretation for this difference is that in the latter studies, a significantly higher  ad,rad is assumed (or in case of detailed binary models, computed).For example, Neĳssel et al. (2019) assumes  ad,rad = 6.5 following Ge et al. (2015), which leads to higher rate of stable mass transfer episodes with significant orbital shrinkage when compared to, for example, Stevenson et al. (2017).
Our results confirm that the value of  ad,rad indeed plays an important role for  stable (see discussion in section 3.3).However, our models also suggest that the significance of the stable channel is affected by the assumed angular momentum loss mode as well, and we expect this to be also the reason why the models of Belczynski et al. (2016) predict a negligible  stable .In particular, Belczynski et al. (2016) assumes  = 1 and a critical mass ratio of  crit = 3 for HG donors (for comparison, our assumed  ad,rad = 4 is equivalent to  crit = 3.2, if the accretor is a BH).Therefore, the models of Belczynski et al. (2016) are fairly similar to our M11 and M12 models (see Table 4).Our simulations of M11 and M12 indicate that the stable channel is essentially negligible, in broad agreement with Belczynski et al. (2016).However, in the models with a higher angular momentum loss (i.e. = 2.5), but with the same  ad,rad ,  stable significantly increases and even dominates for  BH,1 ≳ 25  ⊙ .In conclusion, the stable channel can embody a significant formation channel for merging binary black holes not only for  ad,rad ≳ 6.5 but also for strong angular momentum loss modes such as  = 2.5.
In the majority of our model variations, the most massive merging binary black holes originate from the stable channel (i.e. BH,1 ≳ 25  ⊙ ) in broad agreement with Neĳssel et al. (2019); van Son et al. (2022a) and Briel et al. (2022).In particular, our Fig. 2 can be directly compared to Fig. 5 of van Son et al. (2022a) and Fig. 7 of Briel et al. (2023).Our models generally agree more with that of van Son et al. (2022a) than with that of Briel et al. (2023).This is not surpising, as the models of van Son et al. (2022a) were generated by COMPAS (Riley et al. 2022), which also use the fitting formulae of Hurley et al. (2000).Smaller differences between these models most likely can be attributed to different choices in binary physics assumptions of van Son et al. (2022a), such as  ad,rad = 6.5, isotropic angular momentum loss mode, i.e.  =  d / a and a mass transfer efficiency, which is related to the thermal timescale of the accretor star.The differences between our predictions and the results of Briel et al. (2023) are more significant.While they predict that the stable channel dominates at high masses, the overall contribution of this channel is still small (∼6.4 per cent).Furthermore, they find several other channels to be efficient, including ones in which only one phase of mass transfer occurs.These differences are most likely due to the following: i) their model is produced with BPASS (Eldridge et al. 2017;Stanway & Eldridge 2018), which uses detailed binary models, ii) they assume that if the accretor star accretes 5 per cent of its initial mass, it will evolve chemically homogeneously (which allows channels with only one mass transfer episode to be efficient), iii) they assume a super-Eddington accretion for BH accretors.
Furthermore, we find non-negligible differences in the reported final mass ratio distributions.In particular Briel et al. (2023) finds that most systems from the stable channel have  ≲ 4 (most likely can be attributed to their assumption of super-Eddington accretion rate), while van Son et al. (2022a) finds a relatively narrow range of 0.5≲  ≲ 0.8.On the other hand, the predicted mass ratios of the systems from the stable channel in our models typically peak at  ≈1 and drop rapidly beyond  ≈ 0.8 with  ad,rad = 4 and beyond  ≈ 0.5 with  ad,rad = 7.5 (see Fig. 3 and Fig. B1).van Son et al. (2022a) finds a mass ratio distribution of the CEE sources in the range of 0.2≲  ≲ 1.0, which peaks around q≈ 0.3 and gradually decreases from that value with increasing q.While our predicted mass ratios for this channel are typically in the same range, they peak around  ≈ 1 in all of our model variations.

RESULTS: THE IMPACT OF STELLAR WINDS
In this section, we investigate how different assumptions about mass loss rates of line-driven winds affect the evolution of interacting massive binaries (subsection 4.1) and progenitors of GW sources (subsection 4.2).We test three different models, these are Model I with  wind =  wind,WR = 1, Model II with  wind = 1/3,  wind,WR = 1, and finally Model III  wind =  wind,WR = 1/3.The three different stellar wind models are summarised in Table 3 (see also subsection 2.7).Finally, in subsection 4.3, we discuss the importance of LBV winds and the Humphreys-Davidson limit on GW sources.

The effects of stellar winds on binary evolution
In Figure 5, we show the remnant mass as a function of initial mass for single stars and for stars in interacting binaries at Z = 0.01.We see that stars in interacting binaries produce considerably lower mass black holes compared to their single counterparts (see also for similar conclusions: Woosley 2019; Laplace et al. 2021;Vanbeveren et al. 1997).
There are two major reasons for this.Firstly, as the donor star loses its hydrogen envelope as a result of the mass transfer phase, its hydrogen-shell burning is halted and therefore so is the growth of its helium core.As the helium core is expected to grow only slightly during the very short Hertzsprung gap phase, but significantly during the core-helium burning phase (see later Fig. 6), the remnant mass significantly depend on whether the system undergoes Case B or Case C mass transfer.In particular, at metallicities  ≳ 0.005, the radial expansion after the onset of core helium burning of stars with  ZAMS ≳ 40 ⊙ is negligible (Fig. A1).Practically, all interacting massive binaries with  ZAMS ≳ 40 ⊙ initiate Case A or Case B mass transfer episodes at these metallicities (see similar discussion in Vanbeveren et al. 2007).Consequently, the yellow curves in Fig. 5 (corresponding to systems initiating mass transfer at the end of the Hertzsprung gap phase of the donor star) show approximately the maximum remnant mass that stars with  ZAMS ≳ 40 ⊙ from interacting binaries can have at Z = 0.01.Because of the early envelope stripping, the helium cores of the donor stars cannot grow significantly via hydrogen shell burning.As a result, they also form less massive remnants than their non-interacting counterparts.
Secondly, the lifetime of the Wolf-Rayet phase of a star in an interacting binary increases compared to that of the single star.The primary star in the interacting binary spends most, if not all of its The black lines show the remnant mass for single stars (e.g.their hydrogen envelope is not stripped as a result of a mass transfer episode), the blue lines show the remnant mass for stars, which lose their envelopes just at the moment when they cross the Hertzsprung gap (e.g. at zero-age Hertzsprung gap), the yellow lines show stars, which lose their envelopes at the end of the Hertzsprung gap phase (e.g.terminal-age Hertzsprung gap).Here we ignored the accretion by the BH in an eventual second phase of mass transfer, which is a valid assumption, if the accretion is Eddington-limited core helium burning lifetime as a stripped helium star.As a consequence, the star in the binary ends up losing more mass due to Wolf-Rayet winds than its single counterpart.Since the Wolf-Rayet winds directly affect the mass of the helium core, the mass of the black hole sensitively depends on total mass lost during this evolutionary phase.The single star also loses a significant amount of mass via LBV winds.However, this impacts only the hydrogen envelope, and not the mass of the helium core.The difference between interacting and non-interactive binary evolution is illustrated in Fig. 6, where we compare the evolution of a single star with an initial mass of  ZAMS = 100  ⊙ , with the evolution of a binary with the same primary mass as the mass of the single star.The mass of the black hole originating from single stellar evolution is appreciably higher ( BH = 23.3 ⊙ ) than the black hole formed through binary evolution ( BH = 16.2  ⊙ , but also compare the black with the yellow lines in Fig. 5).In the lower panel of Fig. 6, we compare the evolution of the masses of these two systems starting from the Hertzsprung gap phase.During this phase, the helium core of the single star grows only slightly, from  core ∼ 31.2  ⊙ to  core ∼ 31.7  ⊙ (also compare the blue with the yellow lines in Fig. 5), while it increases up to  core ∼ 35.2  ⊙ by the end of the core helium burning.As previously mentioned, at these metallicities, the most massive binaries are only expected to undergo Case C mass transfer phase in negligible numbers.This implies that the masses of stripped helium stars ( HE,binary ), which evolve from an MS star with  ZAMS = 100  ⊙ in an interacting binary is not above  HE,binary ∼ 31.7  ⊙ for the vast majority of the cases.The single star also loses its envelope eventually, mostly due to LBV winds.However, this occurs well after the onset of the  The importance of stellar winds with and without binary interactions are compared by showing the evolution of a star with an initial mass of  ZAMS,1 = 100 ⊙ at metallicity  = 0.01 as a single star, and as a binary with a companion of  ZAMS,2 = 90 ⊙ and with an initial orbital separation of  ZAMS = 400 ⊙ .We show the binaries for three different stellar wind models.Labels 'Binary', 'Binary II' and 'Binary III' correspond to Model I (  wind = 1), Model II (  wind = 1/3) and Model III  wind =  wind,WR = 1/3 ), respectively.We assume circular orbits.The top panel shows the remnant mass and the mass lost due to different mechanisms.The bottom panel shows the mass of the stars as a function of time.The green symbols indicate the starting points of different stellar evolutionary phases; triangle: hydrogenshell burning phase, hexagon: core helium burning phase, diamond: helium star, cross: helium giant.For clarity, we also show the curves starting from the hydrogen-shell burning phase.core helium burning phase.During this phase, the mass of the helium core can grow uninterrupted.As a result, a more massive helium star is formed, e.g. in this case  HE,single = 35.2⊙ .
We also see that while in Model II single stars form appreciably more massive BHs than in Model I, this difference is much smaller for interacting binaries.As shown in Fig. 6, in Model II, the primary star of the binary system develops a more massive helium core before the envelope loss ( HE,binary ≈ 40.3  ⊙ ) than in Model I ( HE,binary ≈ 31.7  ⊙ ).However, this stripped star in the interacting binary loses an enormous amount of mass via Wolf-Rayet winds and ends up with a black hole that is only about 2.4  ⊙ more massive than the black hole evolving from the same system in Model I (compare the solid lines with the dashed lines in Fig. 5).The decrease in the difference in remnant masses for interacting binary systems compared to single systems is the result of the increased amount of mass lost during the Wolf Rayet phase (∼ 20.3  ⊙ as opposed to ∼ 13.7  ⊙ ).In Fig. 6, we also show the evolution of the binary with Model III, in which the Wolf-Rayet winds are also scaled down by a factor of three.In this case, the mass of the final black hole is about  BH ∼ 28.3  ⊙ , which is about 10  ⊙ more than than the black hole forming in the binary with Model I (compare the solid lines with dotted lines in Fig. 5).

The effect of stellar winds on merging binary black holes
In this subsection, we discuss how the maximum mass of merging binary black holes are affected with different stellar wind models.In Fig. 7, we show the primary mass distribution of merging binary black holes for all of our three stellar wind models and with different values of  ad,rad and  at metallicities  = 0.02, 0.01, 0.005.The difference in the maximum mass in Model I and Model II is only ∼ 4  ⊙ and therefore not significant.This is not surprising, as we found similar results for interacting binaries in section 4.1.
Whether the most massive black holes in Model III form GW sources depends on what we assume about the mass transfer stability criteria for giants with radiative donors.With  ad,rad = 4, the most massive systems in Model III form BH-BH binaries that are too wide to merge within Hubble time, or experience stellar merger and never form BH-BH binaries.Therefore, in that case, the maximum masses of the GW sources do not differ significantly between Model II and Model III.However, with  ad,rad = 7.5, the primary masses of merging binary black holes can reach up to  BH ∼ 30  ⊙ for Model III at  = 0.01, which is significantly larger than that of the most massive GW progenitors in Model II ( BH ∼ 20  ⊙ ).In order to understand this in more detail, let us consider how the most massive merging black hole binaries are formed at metallicties  ≳ 0.005: • The most massive GW sources form via the stable channel.This is due to our assumption that envelope ejection during CEE is only possible with core-helium burning donors.At such high metallcities the expansion of radius in such donors is negligible for  ZAMS ≳ 40  ⊙ (see Fig. A1) and therefore so are the rates of Case C mass transfer events for such systems.
• Binaries with the most massive stars only form GW source, if they have  MT,2 ≳ 4-5.The orbit only shrinks efficiently due to a stable phase of mass transfer with a BH accretor, if  MT,2 is large (see equation 3).In relatively high metallicity environments ( ≳ 0.005), the orbital separations of binaries with the most massive initial masses are typically so wide at the onset of the second phase of mass transfer due to stellar winds (see Figure D2), that only systems with  MT,2 ≳ 4-5 form BH-BH binaries that merge within the Hubble time.
Considering these two points, we can understand why binaries with the most massive stars do not form GW sources with  ad,rad = 4.In these model variations  crit = 3.2 at the onset of the second mass transfer phase (see Table B1).Consequently, a binary with  MT,2 ≈ 4-5 experiences an unstable phase of mass transfer and the system merges before it could form a BH-BH binary.If  ≲  crit , then the BH-BH binary that is formed is too wide to merge within the Hubble time.
In Fig. 8, we show a 2D histograms of the masses of the merging binary black holes at Z = 0.01 for each stellar wind model.We distinguish sources based on the type of the first mass transfer episode (Case A is shown by green and Case B is shown by blue).The most massive systems in Model I and Model II experience a Case B first phase of mass transfer.The most massive black holes from these models have mass ratios 0.8 ≲  final ≲ 1.0.On the other hand, the gravitational wave sources with the most massive primaries are predicted to form in a very different way in Model III.These binaries have their first mass transfers with (late) main sequence donors and the mass ratio distribution of the merging binary black holes are in the range of 0.6 ≲  final ≲ 0.8.In Figure D2, we show typically formation histories of the most massive merging binary black holes from each mode and we discuss their evolution in detail in section D.

The effect of LBV winds on the merging binary black hole population
Stars that cross the Humphreys-Davidson limit are predicted to lose a significant amount of mass via LBV stellar winds.However, the .The primary mass distribution of merging binary black holes for different stellar wind models and different assumptions for mass transfer stability criteria for giants with radiative donors ( ad,rad ) and angular momentum loss during mass transfer with non-compact accretor () at metallicities Z = 0.02, 0.01 and 0.005.For all models shown here the convective envelope prescription of Ivanova & Taam (2004) and  = 0.3 is assumed.This histograms are normalised to the merger efficiency.underlying mechanism for the mass loss, the predicted mass loss rates and its metallicity dependence are extremely uncertain (see e.g.Smith 2014).
If the progenitors of merging binary black holes have sufficiently wide initial orbital separations, such that the donor stars cross the Humphreys-Davidson limit, before they initiate a mass transfer episode, then LBV winds will affect the demographics of this merging binary black hole population.This can occur in the following ways.Firstly, the range of the mass ratio distribution is decreased to lower values at the onset of the second mass transfer by the LBV mass loss rates.This affects the number of mergers of the rCEE channel (see also e.g.Fig. C2 and discussion in section 3.4).Secondly, intense mass loss rates widen the orbit.This can increase the number of binary black holes which are too wide to merge within Hubble time.Finally, extremely high LBV mass loss rates can also affect the maximum size that the stars eventually reach.In principle, if the mass loss rates are sufficiently high (  LBV ∼ 10 −3  ⊙ yr −1 ), the red-ward evolution of massive stars could be truncated by LBV winds, as they would lose their hydrogen-rich envelopes soon after the onset of core-helium burning.
The latter is why LBV winds have been associated with the lack of observed red super giants above luminosities of log(L/L ⊙ ≈ 5.8 (see e.g.Lamers & Fitzpatrick 1988, but also see Higgins & Vink 2020;Gilkis et al. 2021;Sabhahit et al. 2021 for different possible scenarios).In the context of gravitational wave progenitors, Mennekens & Vanbeveren 2014 found that, if the LBV mass loss rate is in the order of  LBV = 10 −3  ⊙ yr −1 or higher, the merger rate of binary black holes drastically decreases.However, we show that, if the BH-BH mergers are dominated by the stable channel, then this is not necessarily true.This is because for such binaries both of the mass transfer phases occur in relatively small orbits, typically when the donor is at the beginning of its hydrogen shell burning phase.Therefore, stars in a large fraction of these systems never cross the Humprehys-Davidson limit.
In Fig. 9, we show the donor stars at the onset of the second phase of mass transfer in the Hertzsprung-Russel diagram for the stable and the CEE channel at Z = 0.0007 for two model variations.We chose these particular models because they have the lowest and the largest number of GW progenitors of the stable channel, in which any of the stars cross the Humphreys-Davidson limit during their evolution.
We also show the percentage of the binaries in the stable channel, in which any of the stars evolves beyond the Humphreys-Davidson limit (  Stable,HD ).These are  Stable,HD = 0.9 per cent in the model with  = 2.5 and  Stable,HD = 29.4 per cent in the model with  = 1.This already indicates that in some of our model variations, only a negligible number of GW progenitors of the stable channel is affected by LBV winds.We note that this fraction is dependent on the mass of the stars, In particular, in the model with the  = 1, all stars with donor mass above 80  ⊙ initiate the mass transfer phase after they crossed the Huprheys-Davidson limit.
We also show  StableHD for all model variations at all metallcities in Fig. E3.This clearly demonstrates that especially at lower metallicties (at which most GW sources are predicted to originate from), the majority of the systems in the stable channel never become LBV stars.
As evident from Fig. 9, the most massive GW progenitors in the CEE channel initiate the second phase of mass transfer well beyond the Humphreys-Davidson limit, since the donors of these systems typically have a deep convective envelope before they fill their Rochelobes.
If the LBV mass loss rates were a magnitude higher than our assumed value, all binaries in the CEE channel with stars that cross the Humphreys-Davidson limit would lose their hydrogen rich envelopes before those would become mostly convective.This would decrease the predicted merger rate and the maximum black hole mass of this formation channel, broadly consistent with Mennekens & Vanbeveren (2014).On the other hand, this is not the case for the stable channel, as a significant fraction of these systems never evolve beyond the Humphreys-Davidson limit.Even those systems which do, the stars of these binaries typically spend there less than ∼ 10 per cent of their Hertzprung gap lifetime.This means that even, if the mass loss rate is  LBV = 10 −3  ⊙ yr −1 , the donor would only lose ∼ 1  ⊙ mass, and therefore the predictions for the stable channel are little affected by a steady, intense LBV mass loss.We show this in Fig E4 , where we determined the mass that would be lost by the donor stars of the stable channel due to LBV wind with a mass loss rate of  LBV = 10 −3  ⊙ yr −1 for the  = 2.5,  = 0.3 models at each metallicity.We see that even with such an enormous mass loss rate, the vastr majority of the donor stars would only lose 1-2  ⊙ .We show the positions of the donor stars at the onset of the second mass transfer phase of GW progenitors in the HD diagram at Z = 0.0007 with 2D histograms for two different model variations.The blue 2D histogram shows the systems of the stable channel, while the red shows the systems of the CEE channel.The yellow line shows the Humphreys-Davidson limit.We also show a few stellar tracks with the same masses as in Fig. A1, as well as the their most importnt evolutionary steps with green shapes (see Fig. A1).In the legend  Stable,HD , expresses the number of systems in the stable channel, in which any of the stars cross the Humphreys-Davidson as a fraction of all stable channel GW sources. Stable,LumHD is the same but expressed as a fraction of only those stable channel sources, in which any of the stars in the binary would cross the Humphreys-Davidson limit, if it evolved without any binary interactions

CONCLUSION
We performed a parameter study on the classical isolated binary formation channel of gravitational sources.Our primary aim was to investigate how sensitively the demographics of merging binary black holes depend on current uncertainties related to the first phase of mass transfer (between two hydrogen rich stars) and on stellar winds.For this, we used the rapid population synthesis code SeBa to simulate the evolution of massive binaries over a metallicity range  = 0.0001 − 0.02 with several model variations each with different assumptions about the binary physics.In the first part of the paper (section 3), we varied the (i) angular momentum loss mode and (ii) the mass transfer efficiency during first phase of mass transfer, (iii) the mass transfer stability criteria for giants donors with radiative envelopes and (iv) the effective temperature at which evolved stars develop deep convective envelopes.In the second part of the paper (section 4), we also varied the mass loss rates of line-driven winds (see Table 1 for an overview for all of our model variations).
In our model variations, we identify two dominant evolutionary paths, one involves two stable mass transfer phases (stable channel), and one in which a stable mass transfer episode is followed by an unstable mass transfer phase (CEE channel with two variant: cCEE and rCEE, see see Fig. 1).
We find that current uncertainties related to first phase of (stable) mass transfer have a huge impact on the relative importance of different dominant channels, while the observable properties (i.e.merger rate, mass and mass ratio distribution) of merging binary black holes are not significantly affected.This implies that models with different binary physics assumptions might yield the same predicted demographics of merging binary black holes, the origin of the majority of GW progenitors could be entirely different.This shows why it is very challenging to infer physics of massive binary evolution solely from GW observations in a meaningful way, given the large uncertainties in our current models (see also Belczynski et al. 2022).
Our results are in broad agreement with recent studies that have shown that the stability of mass transfer and the structure of the envelope of the donor star plays a significant role in determining the relative importance of the stable and the CEE channel (see e.g.Neĳssel et al. 2019;Klencki et al. 2020;Bavera et al. 2021;Marchant et al. 2021;Olejak et al. 2021;van Son et al. 2022a;Briel et al. 2023).Furthermore, we have demonstrated that uncertainties regarding the modelling of the first phase of mass transfer, such as the angular momentum loss mode and mass transfer efficiency adds an additional layer of complexity in predicting the dominant formation channel of GW sources from isolated binaries.
In order to break the degeneracy between the uncertainties in modelling mass transfer and the demographics of the population of merging binary black holes, it is clear that we need more constraints than gravitational wave data is currently supplying us, for example, constraints from electromagnetic observations.In particular, regarding the first mass transfer phase, observations of the mass ratios and periods of WR-O/B binary systems could offer invaluable clues about the physics of a mass transfer episode between two hydrogen rich stars.This is because WR-O/B binary systems have not yet experienced core-collapse, and therefore their mass ratios and orbital separations are directly related to the mass transfer efficiency and the angular momentum loss during the first phase of mass transfer.While, only about a few dozens of WR-O/B are known, for which the relevant parameters (such as mass ratio, period, rotational period) are reliable inferred (see e.g.van der Hucht 2001; Crowther 2007), related studies have already lead to very important insights about binary physics (see e.g.Vanbeveren et al. 1997Vanbeveren et al. , 1998b;;Petrovic et al. 2005;Eldridge 2009;Shara et al. 2017;Vanbeveren et al. 2018;Shao & Li 2016;Shenar et al. 2016).Our study demonstrates that a larger, unbiased and more complete catalogue of WR-O/B systems with reliably measured parameters could significantly improve our understanding about the origin of merging stellar mass binary black holes.
Below, we summarise our most important results: • Impact of angular momentum loss mode on the dominant channels (section 3.1 and 3.2): Angular momentum loss during the first phase of mass transfer determines how the orbital separation changes during mass exchange and has an indirect effect on the rate of stable and unstable mass transfer phases.Consequently, it has also a strong impact on the merger rate of different formation channels.We find that typically the stable channel dominates in our models with  = (/ tot )/(/ tot ) = 2.5, while, the merger rate of cCEE channel remains low.On the other hand, in our models with  = 1, the merger rate of the stable channel is typically negligible, while the cCEE channel is efficient.
• Impact of mass transfer stability criteria (section 3.3): The merger rate of the stable channel increases with increasing  ad,rad .This is because the degree of orbital shrinkage due to a stable mass transfer phase with a black hole accretor becomes substantially more efficient with increasing mass ratios (see e.g.van den Heuvel et al. 2017;Pavlovskii et al. 2017;Ge et al. 2020).However, we find that only the merger rate of relatively lower mass BH-BH binaries is significantly affected (i.e. BH ≲ 20  ⊙ ), when  ad,rad is increased from 4 to 7.5 (corresponding to a  crit of 3.2 and 5.5).This is because at lower metallicities (which are the most relevant for GW sources), the maximum mass ratios of the most massive binaries at the onset of the second mass transfer phase is limited to 3-3.5, depending on the mass transfer efficiency.
• Line driven stellar winds (section 4.1 and 4.2): Although, decreasing the mass loss rates of optically thin line driven winds by a factor of three significantly increases the masses of the black holes formed without mass exchange (see Fig. A1), we do not find an appreciable difference in the masses of merging binary black holes (Fig. 7).The primary reason for this is that stars in interacting binaries experience a longer Wolf-Rayet phase due to envelope stripping than their single counterparts.Although, with weaker winds on the main sequence, stars develop a more massive helium core, they also experience significantly higher Wolf-Rayet mass loss rates after the envelope stripping, since the more massive helium core will be more luminous after envelope loss, and this leads to considerably higher mass loss rates (see similar results in Woosley 2019; Laplace et al. 2021).However, if the Wolf-Rayet mass loss rates are simultaneously decreased too, we find a significant increase in the masses of gravitational wave sources, but only in the model variations with  ad,rad = 7.5.
• LBV winds (section 4.3): Intense LBV mass loss rates can have a significant effect on the CEE channel.If the mass loss rates are above  LBV ∼ 10 −3  ⊙ yr −1 , the merger rate of the CEE channel becomes negligible (see also Mennekens & Vanbeveren 2014).On the other hand, the stable channel is not appreciably affected, as these binaries have so short periods, such that their envelopes are stripped as a result of a mass transfer episode before or very soon after they cross the Humpyhreys-Davidson limit.massive as  BH ≳ 30  ⊙ are formed in Model II.Since the initial masses of massive stars follow a distribution of  ∼  −2.3 ZAMS (Kroupa 2001), we therefore expect that only the high mass tail of the black hole mass distribution would be affected significantly when scaling down optically thin line driven stellar winds.Below  ≈ 0.002, the differences between the two stellar wind models become negligible.

APPENDIX B: MERGER RATE DENSITY IN THE LOCAL UNIVERSE
The merger rate density is an important characteristics of a synthetic population of GW sources, as it allows comparison to the catalogue of gravitational wave source detections made by LIGO and Virgo (e.g.Abbott et al. 2021a).We calculate the merger rate density in the local universe (i.e. ≈ 0) for all of our model variations presented section 3.In order to determine this quantity, we follow a formalism similar to the one outlined in Dominik et al. (2015).
First, we define the merger efficiency as the number of mergers originating from a black hole binary with masses  BH,1 and  BH,2 at metallicity  as the fraction of the full parameter space (i.e as the fraction of all the stars born, assuming the initial conditions introduced in subsection 2.8): where  merger is the number of mergers originating from a BH binary with given masses at a given metallicity,  simulated is the number of sampled systems in our simulation and  pm is the simulated parameter space as a fraction of the complete parameter space.We determine the latter as: (B2) here  IMF is the normalised initial mass function of Kroupa (2001), and has been normalised in an interval of 0.08  ⊙ -100  ⊙ ,   is the normalised mass ratio distribution, which has been normalised in the interval of 0-1, and   is the normalised semimajor axis distribution, which has been normalised between 10-10 6  ⊙ ,  bin is the binary fraction, which we assume to be 0.7, following Sana et al. (2012).
The merger efficiency can be converted to a merger rate density by assuming a metallicity specific star formation density model SFRd * (, ): where  d is the so-called delay time, which is the time between the zero-age main sequence of the stars in the binary and the time of the merger due to GWs,  ZAMS is the redshift at zero-age main sequence and it is a function of the time delay, M is the average mass of all stellar systems born, such that SFRd(, ) * / M gives the average number of stars born at redshift  and in the metallicity bin centred around . Delay time can be directly obtained from the population synthesis simulations.In order to determine  ZAMS for a given  and delay time, we use the relationship for lookback time: where The star formation rate density in a metallicity bin centered at  with a width of 2Δ is Here  met (, ) gives the distribution of metallicities of binaries at redshift  as: where  = 0.5 and the mean metallicity varies with redshift as given by Madau & Fragos (2017): where we took  ⊙ = 0.02.Furthermore, the star formation rate density SFRd() is given by Madau & Dickinson (2014): We note again that the metallicity specific cosmic star formation rate is very uncertain and this uncertainty has an important impact on the predicted merging binary black hole population (see e.g.Chruślińska et al. 2020;Tang et al. 2020;Briel et al. 2022;Chruślińska 2022)

APPENDIX C: THE EFFECT OF THE FIRST PHASE OF MASS TRANSFER
In this subsection, we discuss in detail how the first phase of mass transfer affects the mass ratio and orbital separation distributions of binaries.In the following, we restrict the discussion to sources with  ZAMS,1 ≳ 45  ⊙ at  = 0.0007.We do this so that the the impact of the first mass transfer phase can be more easily understood, since at this metallicity the stellar winds are negligible and for this mass range the orbit is not modified by any natal kicks.Therefore, the impact of the first mass transfer episode can be more easily understood..However, we note that results presented here are still qualitatively true for all binaries.
In the next section, we present examples for Case B and Case C mass transfer episodes.For this, we assume that the donor fills its Roche-lobe when it reaches the midpoint of its Hertzsprung gap lifetime.As noted in subsection 2.1, the outcome of Case B mass transfer epsidoe is not sensitivey dependent on when exactly the Roche-lobe overflow occurs (or equivalently, on the initial orbital separation).
For the examples of Case C mass transfer episodes, we assume that the donor fills its Roche-lobe when it reaches an effective temperature of log(T eff ) = 3.73 K.The latter criteria does not strictly imply that the donor is in its core helium burning phase.In particular, at Z = 0.0007, binaries with  ZAMS,1 ≳ 90  ⊙ are still in their Hertzsprung gap phase at this effective temperature.Nevertheless, we still choose this criteria, as we find that it is fairly typical for GW sources of the CEE channel to initiate the first mass transfer phase at this stage.We also note, that outcome of the Case C mass Table B1.Conversion between the mass transfer stability criteria for giants and the critical mass ratio ( crit =   /  ) for different assumptions on the mass transfer efficiency and the angular momentum loss mode during the mass transfer phase.The last column reflects a binary with a black hole accretor whose mass accretion is Eddington-limited and the expelled mass having the specific angular momentum of the accretor.Donors with radiative envelopes are characterised by  ad,rad , whereas donors with convective envelopes are described by  HW .For the latter, we assume a core mass to total mass ratio of 0.45-0.63,which we find is typical for giants at log(T eff ) ≈ 3.73 K with initial masses  ZAMS = 50 − 100  ⊙ with a metallicity independent LBV mass loss rate of  LBV = 1.5 • 10 −4  ⊙ yr −1 . crit with  ad =  HW 0.9 − 0.7 0.9 − 1.1 1.5 − 1.9 1.1 − 1.4 1.4 − 1.7

Mass transfer stability criteria
transfer phases is more sensitively dependent the initial separation (see discussion 2.1).
The distribution of mass ratios and orbital separation of massive binaries at the onset of the second phase of mass transfer are distinct for each formation channel.Therefore, the distribution of these parameters helps us to understand the importance of each channel.The shape of this distribution is primarily determined by the initial conditions and the first phase of mass transfer.The progenitors of each sub-channel introduced in section 1.1 have the following range of mass ratios at the onset of the second phase of mass transfer: • cCEE channel: in this scenario the second phase of mass transfer is Case C, unstable and the donor star has a deep convective envelope, i.e.  MT,2 >  crit with  ad =  HW • rCEE channel: the second phase of mass transfer is Case C, unstable but the donor star has a radiative envelope, i.e.  MT,2 >  crit with  ad =  ad,rad .
• stable channel: The binary has a relatively short period (i.e.Case B first phase of mass transfer).The binary has a relatively large  MT,2 .The orbit only then can shrink efficiently due to a stable mass transfer (see Equation 3).This means there is a minimum value of  MT,2 associated with this channel (for a given  ZAMS,1 ).This value depends on  and  and on the predicted outcome Case A mass transfer episodes in a complicated way, but it is typically  MT,2,min ≳ 2.2 for at Z≲0.005 (see Fig. C3).Therefore:  MT,2,min <  MT,2 <  crit with  ad =  ad,rad .
We note that a stable phase of mass transfer does not typically change the orbit by orders of magnitude.Therefore, the types of the first and the second phases of mass transfers are the same in most cases.For example, most CEE sources have a Case C first phase of mass transfer.

C1 Effect on the mass ratio
In the left panel of Fig. C1, we show the relation between the initial mass ratio  ZAMS ≡  ZAMS,2 / ZAMS,1 and that at the onset of the second mass transfer phase  MT,2 ≡   /  .This figure shows the effect of different mass transfer efficiencies and different initial masses for Case B and Case C mass transfer types.
Fig. C1 tells us what fraction of the binaries can potentially form sources of each channel, and what their initial mass ratio range is.For example, as we will see, GW sources of the stable channel typically require 2.2 ≲  MT,2 ≤  crit .If  = 0.3, we expect the most massive sources of the stable channel to evolve from binaries with  ZAMS ≳ 0.8.On the other hand, the formation via the rCEE channel, which requires  MT,2 >  crit , is not possible at all with mass transfer efficiencies as low as  = 0.3.Even if  = 0.7, the formation becomes only possible if  ZAMS ≳ 0.85 and  ZAMS,1 ≳ 60  ⊙ .
We can see that the curves for Case C mass transfer episodes do not have a simple linear relationship, as is the case for Case B. Instead a dip can be observed in the curves at a given range of initial mass ratio.The position of this dip shifts to lower values with increasing mass.This is due to the effect of LBV mass loss.The dip occurs at the lowest initial mass ratio at which the secondary star is still massive enough after the first mass transfer phase, such that it evolves beyond the Humphreys-Davidson limit and eventually loses a significant fraction of its envelope before the onset of the second phase of the mass transfer.At the same time, the mass of the accretor black hole at the second mass transfer phase (formed by the initial primary star) is not significantly affected by LBV winds.This implies that the range of distribution of  MT,2 of binaries with a Case C first phase of mass transfer are narrowed down by the LBV winds.
The maximum  MT,2 that can be reached by binaries at Z = 0.0007 is shown in Fig. C2.The maximum  MT,2 is given as a function of  ZAMS,1 for a given mass transfer efficiency and mass transfer type.The maximum occurs for near equal mass binaries.We expect the results to depend only weakly on the metallicity as long as  ≲ 0.002, because the  ZAMS −  BH relation does not significantly change at this metallicity range (see discussion later in section 4).
We can use Fig. C2 to asses the importance of  ad,rad .A higher  ad,rad implies a higher critical mass ratio and thereby increases the parameter space for stable mass transfer phase.The degree of orbital shrinkage during stable phase of mass transfer with a black hole accretor increases exponentially with increasing mass ratio (Equation 3).Hence, we may expect that the merger rate of the stable channel could increase significantly with increasing  ad,rad .However, we also have to take into account the limit on the maximum mass ratio at the onset of the second phase of mass transfer.This value is determined by  and the  ZAMS −  BH relation.The maximum mass ratio of binaries with Case B mass transfer phase is below 3.2, if  ZAMS,1 ≳ 60  ⊙ and  = 0.3.Therefore, in our low mass transfer efficiency models, increasing  ad,rad from 4 to 7.5 only affects the evolution of binaries with  ZAMS,1 ≲ 60  ⊙ .Fig. C2 also shows the importance of the mass transfer efficiency on the rCEE channel.We find that with low accretion efficiencies, the formation of GW sources via the rCEE channel is not possible for almost the entire mass range.This can be seen from Fig. C2  The mass ratio distribution of the primary of the merging binary black hole in the local universe (z = 0) for all of our models with our standard stellar wind models (e.g. wind =  wind,WR = 1, see Table 1).The distribution is shown by stacked histograms, where each colour indicates a different formation channel.
up to 7.5, then the formation of such systems become impossible even with larger mass transfer efficiencies.

C2 The effect on the orbital separations
The relationship between the initial orbital separation and the orbital separation at the onset of the second phase of mass transfer ( MT,2 ) can be approximated as: if the effects of supernova are negligible.Here,  wind,post−MT,1 and  wind,pre−MT,1 describes by what factor the orbit widens due to stellar winds from the birth of the binary until the onset of the first mass transfer phase and from the end of the first mass transfer phase until the onset of the second mass transfer phase, respectively.These changes in orbit are described by Equation 4. The terms ) and the mass ratio at the onset of the second mass transfer phase ( MT,2 =   /  ).We have neglected the effects of supernova but not the stellar winds.We show these relations for different mass transfer types, mass transfer efficiencies and initial masses.For Case B, we assume that the mass transfer phase starts at the half of the Hertzsprung gap phase lifetime of the donor.In the other case, we assume that the mass transfer phase occurs just before when the donor reaches an effective temperature of log 10 T eff = 3.73K.This coincides with the stage where donors develop deep convective envelopes according to Ivanova & Taam (2004).For the majority of the binaries, this occurs with a core-helium burning donor ( ZAMS,1 ≳ 90 ⊙ at Z = 0.0007).We chose this criteria, because we find that most of the gravitational wave sources from the CEE channel roughly initiate their first mass transfer phase at this stage.Right panel: The relationship between the the initial mass ratio and the shift in the orbital separation in log space for different different mass transfer types, mass transfer efficiencies, angular momentum losses and initial masses.Here we have always assume stable mass transfer phase regardless of the assumed  ad,rad .We have neglected the effects of supernova but not the stellar winds.The upper horizontal, dotted line shows  crit = 3.2, which corresponds to the critical mass ratio for a mass transfer phase with black hole accretor and giant donors with radiative envelopes with  ad,rad = 4.The lower bottom line show the critical mass ratio when the donor has a deep convective envelope.See Table B1  wind,post−MT,1 and  wind,pre−MT,1 are negligible for binaries evolving through Case B mass transfer at metallicities  ≲ 0.005.On the other hand, if the first mass transfer phase is a Case C, then LBV winds contribute to the widening of the orbit appreciably at all metallicities.The term  MT,1 is the change in the orbit during the first mass transfer phase.In principle, these are due to the com-bined effects of stellar winds and the mass transfer phase itself.However, the effects of the former are typically negligible especially if the mass transfer episode proceeds on thermal timescale.The change in the orbit due to the mass transfer episode is determined by Equation 2. In all cases,  MT,1 is the dominant term, i.In the right panel of Fig. C1, we show log(a MT,2 /a initial ) as a function of  ZAMS for several model variants and for both Case B and Case C mass transfer types.We see that for Case B mass transfer phase, the change in the orbital separation is mostly determined by the mass ratio of the system, and it is only very weakly independent on the masses of the binary.The latter no longer holds if the first mass transfer phase is a Case C.Here the mass of the initially primary star does affect the amount of orbital widening.There is an upward shift in the curves of log(a 2ndMT /a initial ) for donor stars with strong LBV winds.This effect is reduced for more massive donor stars, as the amount of mass lost due to LBV winds before the onset of the second mass transfer phase (at logT eff = 3.73 K) decreases with increasing initial mass.
In general, the orbital separations of systems with low initial mass ratios shrink, whereas the opposite is true for binaries with mass ratios close to one.This is also clear from Equation 2, which shows that as long as   >   the separation of the system decreases during the mass transfer phase.Where the turnover occurs between widening and shrinking, however, is dependent on the assumed angular momentum loss mode  and accretion efficiency .As we increase , the effect of the first mass transfer phase starts to converge to the fully conservative case, and therefore the difference between the angular momentum loss becomes less important.
In the context of gravitational wave sources, there is an important point to make regarding the stable channel.In the initial mass ratio range, where we expect the GW sources of stable channel to originate from (i.e. from relatively similar initial masses), the model with  = 2.5 and  = 0.3 shows only an insignificant change in the orbit.The largest widening occurs in models with  = 1.This shows why the stable channel is the most efficient in the  = 2.5 and  = 0.3 model, and why it is typically inefficient with  = 1.

C3 Case A mass transfer episode
The outcome of a Case A mass transfer phase is much more sensitively dependent on the initial separation and the mass ratio than that of the previously discussed mass transfer types (see also subsection 2.5).This means that it is challenging to show simple relations for these kind of systems as we did in Fig. C1 and C2 for Case B and Case C.
To demonstrate typical outcomes of Case A mass transfer phases, we show the evolution of a few selected binaries until the second mass transfer phase as simulated by SeBa in Fig. D1.These results show two important characteristics of systems evolving through a first phase of Case A mass transfer.Firstly, they develop  MT,2 that are significantly larger than that of their Case B counterparts.These mass ratios are often above  crit .Secondly, if they survive the Case A mass transfer phase, their orbit widens more compared to systems with Case B mass transfers.These two points imply that the majority of these systems initiate an unstable second phase of mass transfer with an HG donor, which results in stellar merger.All of the systems shown in Fig. D1 would indeed merge with  ad,rad = 4, while one of them would merge with  ad,rad = 7.5 as a result of the second mass transfer phase.We, however, stress again that such results from rapid population synthesis codes following the stellar tracks of Hurley et al. (2000) should be treated with caution.These stellar tracks do not track the developing helium core during main sequence and consequently the mass of the stripped donor after the Case A mass transfer might be severely underestimated.

APPENDIX D: EVOLUTION OF MASSIVE BINARIES WITH DIFFERENT STELLAR WIND MODELS
Here, we discuss the evolution of the systems, which form merging binary black holes with the most massive black hole binary in each of our stellar wind models.In order to gain a more detailed understanding of the combined effect of  ad,rad and different stellar wind models, we show the typical formation histories of sources with the most massive primary black holes for each of our stellar wind models in the Appendix in Fig. D2,.In these scenarios, we assume  = 2.5,  = 0.3 and  ad,rad = 7.5.Fig. D2 shows that the formation path of the most massive sources in Model I and Model II are fairly similar.The evolution starts out with an almost equal mass binary, as this leads to the highest mass ratios at the onset of the second mass transfer phase.On the other hand, the minimum initial separation (so that the first mass transfer is Case B) of such sources is different.Decreasing the optically thin stellar winds by a factor of three increases the radii of stars at the end of the main sequence.In particular, at Z = 0.01, for  ZAMS = 100  ⊙ ,  TAMS can increase from 83.2  ⊙ to 110  ⊙ when switching from Model I to Model II.Therefore, binaries of Model II need to start out with a higher initial separation, because of their larger  TAMS .It is also confirmed, that these sources indeed require large mass ratios of  ∼ 4 − 5 at the onset of the second mass transfer.The formation of such systems is not possible in the  ad,rad = 4 model as they would undergo unstable mass transfer and would merge during CEE.We note that primary masses of ∼ 16 − 18  ⊙ for Model I with  ad,rad = 4 can still be seen in Fig. 7, however these systems have secondary black hole masses  2 ∼ 5 − 7  ⊙ .The formation of such systems is possible because they start out with a high initial mass ratio and therefore the orbit significantly shrinks after the first mass transfer (even if  = 1).
Binaries with initial parameters as shown for Model II in Fig. D2 do not form gravitational wave sources in Model III.The evolutionary path of such a binary would be exactly the same as for Model II until the formation of the first naked helium star.However, because of the decreased Wolf-Rayet-like mass loss rates, the initial primary star forms a 28.0  ⊙ black hole.The intially secondary star becomes a 93.5  ⊙ massive HG giant at the start of the second mass transfer.The mass ratio of 3.3 is not sufficiently high enough to shrink the orbital separation (∼ 420  ⊙ ) enough so that the 28.0  ⊙ -30.5  ⊙ BH-BH would merge within Hubble time.
Instead, the sources with the most massive primaries in Model III form with similar initial masses as in Model II but at shorter orbital separations ( ∼ 180 − 200  ⊙ ).At such a separation the first mass transfer occurs with a MS donor.As explained earlier, this halts the growth of the helium core.It is predicted that the mass transfer episode proceeds while the donor eventually evolves off the main sequence and it finally ends with the HG donor losing its hydrogen envelope, leaving behind a 26.8  ⊙ massive naked helium star.With the lowered Wolf-Rayet-like winds, this collapses into a 19.8  ⊙ black hole.At the beginning of the second mass transfer, the donor star is a 96.3  ⊙ HG giant.The mass ratio ∼ 4.9 is high enough to shrink the orbital separation ( ∼ 370.8  ⊙ ) significantly.The same system in Model II would have a larger mass ratio because of the lower mass of the black hole ( 1 ≈ 14.9  ⊙ ) and the mass transfer would become unstable, which results in a stellar merger.As systems with such high masses always require mass ratios  MT,2 > 4, it is clear why there is such a huge difference in the maximum masses of merging binary black holes between  ad,rad = 4 and  ad,rad = 7.5 for Model III.

APPENDIX E: ADDITIONAL FIGURES
In this subsection, we present additional figures that are helpful (but not necessary) to understand the main results of this work.In Fig. E1, we compare the primary BH mass distributions of our models (shown in Table 4) with the inferred distribution of the GWTC-3 catalog (Abbott et al. 2023).We represent the inferred distribution with the "Power Law + Peak" model (see e.g.equation B4 in Abbott et al. 2023).The values of the model parameters are:  = −3.5,  min = 4.6  ⊙ ,  max = 86  ⊙ ,  peak = 0.038,   = 34  ⊙ ,   = 5.7  ⊙ ,   = 4.82  ⊙ .For simplicity, we ignored uncertainties.We normalised the distribution to the inferred merger rate density of BH-BH mergers;  GWTC3 = 28.3Gpc −3 yr −1 .While the upper end of the inferred distribution reaches  BH,1 = 86  ⊙ , the maximum BH mass of our model variations is only  BH,1 = 42  ⊙ , since we do not sample stars more massive than  ZAMS = 100  ⊙ .To make the comparison easier, we only show the figures up to  BH,1 = 60,  ⊙ .None of our models can reproduce the peak at 34  ⊙ .If this feature is indeed related to (pulsational) pair instability supernova ( PPISN Farmer et al. 2019;Marchant et al. 2019), then this is not surprising .Firstly, we do not model PPISN in this study.Secondly, since we do not consider stars with  ZAMS > 100  ⊙ , our PPISN rate would be negligible in any case.All of our models, except M9-M10, yield distributions that are broadly similar to the powerlaw component of the inferred distribution.M9-M10 exhibit much steeper decrease at  BH,1 ≳ 20  ⊙ .Models M5-M16 have a peak around  BH,1 ≈ 10  ⊙ , similarly to the inferred distribution.On Figure D1.A summary of the evolution of a few selected binaries going through Case A mass transfer phases compared with a binary evolving through a Case B mass transfer episode.We assume  = 1 and  = 0.3.We only show the evolution until the initial secondary fills its Roche-lobe The numbers denoting the evolutionary stage are the same as shown in Fig. 1.Therefore Stage 1 refers to the binary at its birth, Stage 2 is at the end of the second mass transfer phase and Stage 3 is just before the second mass transfer episode.We note that for Case A binaries the first mass transfer phase is interrupted at the end of the main sequence of the donor, as its radius starts shrinking.Shortly after, the radius starts to expand rapidly again with the beginning of the Hertzsprung gap phase and therefore refills its Roche-lobe.This second stage ends with the donor losing its hydrogen envelope.For Case A binaries stage 3 refers to stage when the donor becomes a stripped helium star.
the other hand,for models M1-M4, this peak occurs at somewhat higher masses ( BH,1 ≈ 15-20  ⊙ ), which is not supported by observations.In Fig. E2, we show the merging efficiency of each of our model variations (with standard stellar wind models) at each simulated metallicity.In Figure E3, we show the number of binaries of the stable channel in which any of the stars cross the Humphreys-Davidson limit, expressed as a fraction of all systems in the stable channel.In Figure E4, we show the estimated total mass lost due to LBV winds with a mass loss rate of  LBV = 10 3  ⊙ yr −1 of those stars in the systems of the stable channel sources, which cross eventually the Humphreys-Davidson limit in the  = 2.5,  = 0.3,  ad,rad = 4,  eff IT model variation.In Fig. C3, we show density contours, which reflect the mass ratios and orbital separations of binaries at the onset of the second phase of mass transfer.These binaries were simulated with SeBa at a metallicity  = 0.0007.In Fig. C3, we limit the initial primary mass to  ZAMS,1 = 90  ⊙ .We show four model variations with different assumptions on  and .At this stage, the initial primary star has already formed a black hole.If the first phase of mass transfer is Case B, then the mass of the black hole is  BH,1 = 35.44 ⊙ .We distinguish binaries for which the first phase of mass transfer is Case A (shown by the green contour, hereafter 'Case A binaries') and for which it is Case B or Case C (shown by the blue contour, hereafter 'Case B binaries' and 'Case C binaries', respectively).We also show the progenitors of GW sources as 2D histograms over the contours.We distinguish three types: CEE chan-nel, stable channel with Case B binaries ('stable Case B', hereafter) and stable channel with Case A binaries ('stable Case A').

Figure 1 .
Figure1.The most common formation channels of gravitational wave sources from isolated binaries as predicted in this paper.cCEE means that the dynamically unstable mass transfer is initiated by a giant donor, which has a deep convective envelope, whereas for rCEE the donor still has mostly a radiative envelope.
here the period,  is in days and  = −0.55.Although, other studies suggested somewhat flatter distributions (e.g.Kobulnicky et al. 2014or Dunstall et al. 2015 based on observations of B stars in the Tarantula Nebula in the Large Magellanic Cloud).
Figure2.The mass distribution of the primary of the merging binary black hole in the local universe (z = 0) for all of our models with our standard stellar wind models (e.g. wind =  wind,WR = 1, see Table1).The distribution is shown by stacked histograms, where each colour indicates a different formation channel.In the legend on the upper right corner in each panel, we show the percentage of each channel.In Fig.E1, we provide a comparison between these distributions and the primary mass distribution of BH-BH binaries as inferred from the the third LIGO-Virgo Gravitational-Wave Transient Catalog (GWTC-3)(Abbott et al. 2023) Figure 4. Our simulated binaries in the initial mass ratio -initial separation space.We show the orbital separation after tidal circularisation.The systems are simulated with the following parameters;  = 0.3,  ad,rad = 7.5.Furthermore, we have used the assumption ofIvanova & Taam (2004) regarding the development of deep convective envelopes.In the left panel we show systems with  = 2.5, in the right with  = 1.We show the outcome of each simulated binary with difference colours.We distinguish the following types; (i) interacting, detached: these binaries interact via mass transfer phases at least once and form detached BH-BH binaries which do not merge due to gravitational waves within Hubble time, (ii) interacting, disrupted: these binaries interact at least once and get disrupted due to the supernova kick of one of the binary components, (iii) interacting, stellar merger: these systems interact at least once and merge due to one of their mass transfer phases.Their first mass transfer phase is typically stable or unstable with a MS or HG donor, (iv) non-interacting: none of the stars fill their Roche-lobes, they either stay detached or get disrupted, (v) CEE in 1st MT with CHeB donor: these binaries have a core-helium burning donor at the onset of the first mass transfer phase, and they initiate a dynamically unstable mass transfer (vi) GW source, stable channel: gravitational wave progenitors evolving via two subsequent stable mass transfer phases, (vii) GW source, CEE channel: gravitational wave progenitors evolving via a stable and an unstable mass transfer phase.

Figure 5 .
Figure 5.The mass of the remnant as a function of initial mass at Z = 0.01 for stars in interacting binaries and single stars for our three stellar wind models.Model I shown by solid lines, Model II shown by dashed lines and Model III shown by dotted lines.The black lines show the remnant mass for single stars (e.g.their hydrogen envelope is not stripped as a result of a mass transfer episode), the blue lines show the remnant mass for stars, which lose their envelopes just at the moment when they cross the Hertzsprung gap (e.g. at zero-age Hertzsprung gap), the yellow lines show stars, which lose their envelopes at the end of the Hertzsprung gap phase (e.g.terminal-age Hertzsprung gap).Here we ignored the accretion by the BH in an eventual second phase of mass transfer, which is a valid assumption, if the accretion is Eddington-limited Figure6.The importance of stellar winds with and without binary interactions are compared by showing the evolution of a star with an initial mass of  ZAMS,1 = 100 ⊙ at metallicity  = 0.01 as a single star, and as a binary with a companion of  ZAMS,2 = 90 ⊙ and with an initial orbital separation of  ZAMS = 400 ⊙ .We show the binaries for three different stellar wind models.Labels 'Binary', 'Binary II' and 'Binary III' correspond to Model I (  wind = 1), Model II (  wind = 1/3) and Model III  wind =  wind,WR = 1/3 ), respectively.We assume circular orbits.The top panel shows the remnant mass and the mass lost due to different mechanisms.The bottom panel shows the mass of the stars as a function of time.The green symbols indicate the starting points of different stellar evolutionary phases; triangle: hydrogenshell burning phase, hexagon: core helium burning phase, diamond: helium star, cross: helium giant.For clarity, we also show the curves starting from the hydrogen-shell burning phase.
Figure7.The primary mass distribution of merging binary black holes for different stellar wind models and different assumptions for mass transfer stability criteria for giants with radiative donors ( ad,rad ) and angular momentum loss during mass transfer with non-compact accretor () at metallicities Z = 0.02, 0.01 and 0.005.For all models shown here the convective envelope prescription ofIvanova & Taam (2004) and  = 0.3 is assumed.This histograms are normalised to the merger efficiency.

Figure 8 .
Figure8.2D histograms of masses of the merging binary black holes for the three stellar wind models at  = 0.01.We have separated sources that have Case A first mass transfer phase (green) and sources that have Case B first mass transfer phase (blue).The histograms have been normalised to the merger efficiency at  = 0.01.
Figure 9.We show the positions of the donor stars at the onset of the second mass transfer phase of GW progenitors in the HD diagram at Z = 0.0007 with 2D histograms for two different model variations.The blue 2D histogram shows the systems of the stable channel, while the red shows the systems of the CEE channel.The yellow line shows the Humphreys-Davidson limit.We also show a few stellar tracks with the same masses as in Fig.A1, as well as the their most importnt evolutionary steps with green shapes (see Fig.A1).In the legend  Stable,HD , expresses the number of systems in the stable channel, in which any of the stars cross the Humphreys-Davidson as a fraction of all stable channel GW sources. Stable,LumHD is the same but expressed as a fraction of only those stable channel sources, in which any of the stars in the binary would cross the Humphreys-Davidson limit, if it evolved without any binary interactions

Figure A1 .
Figure A1.Hertzsprung-Russel diagrams for massive stars with the mass-loss rate colour-coded in the tracks.The stars have  ZAMS = 20, 30, 50, 60, 80, 100  ⊙ .Different shaded regions denote different stellar wind mechanisms operating on stars.Within the 'Line-Driven' range, the two dashed lines show where the first and the second bi-stability jump occurs.Green figures indicate the starting points of different stellar evolutionary phases; circle: main-sequence, triangle: hydrogen-shell burning phase, hexagon: core helium burning phase, star: AGB, diamond: helium star, cross: helium giant.

Figure A3 .
Figure A3.Zero age main sequence mass and remnant mass functions for different metallicities for single stellar evolution.We assume that black hole formation occurs above  ZAMS = 20 ⊙ .
Fig.C2also shows the importance of the mass transfer efficiency on the rCEE channel.We find that with low accretion efficiencies, the formation of GW sources via the rCEE channel is not possible for almost the entire mass range.This can be seen from Fig.C2as the curve for Case C mass transfer for  ≤ 0.3 and  ad,rad = 4 is below the critical mass ratio for  ZAMS,1 ≳ 40  ⊙ .If  ad,rad is increased

Figure C1 .
Figure C1.Left panel: the relationship between the initial mass ratio ( ZAMS =  ZAMS,2 / ZAMS,1 ) and the mass ratio at the onset of the second mass transfer phase ( MT,2 =   /  ).We have neglected the effects of supernova but not the stellar winds.We show these relations for different mass transfer types, mass transfer efficiencies and initial masses.For Case B, we assume that the mass transfer phase starts at the half of the Hertzsprung gap phase lifetime of the donor.In the other case, we assume that the mass transfer phase occurs just before when the donor reaches an effective temperature of log 10 T eff = 3.73K.This coincides with the stage where donors develop deep convective envelopes according toIvanova & Taam (2004).For the majority of the binaries, this occurs with a core-helium burning donor ( ZAMS,1 ≳ 90 ⊙ at Z = 0.0007).We chose this criteria, because we find that most of the gravitational wave sources from the CEE channel roughly initiate their first mass transfer phase at this stage.Right panel: The relationship between the the initial mass ratio and the shift in the orbital separation in log space for different different mass transfer types, mass transfer efficiencies, angular momentum losses and initial masses.Here we have always assume stable mass transfer phase regardless of the assumed  ad,rad .We have neglected the effects of supernova but not the stellar winds.The upper horizontal, dotted line shows  crit = 3.2, which corresponds to the critical mass ratio for a mass transfer phase with black hole accretor and giant donors with radiative envelopes with  ad,rad = 4.The lower bottom line show the critical mass ratio when the donor has a deep convective envelope.See TableB1

Figure C2 .
Figure C2.The estimated maximum value of mass ratio at the onset of the second mass transfer phase ( MT,2 =   /  ) that a system can reach with a given initial primary mass, mass transfer efficiency and mass transfer phase type.The maximum  MT,2 is achieved by binaries that have initially near equal masses.How exactly we defined 'Case B' and 'Case C' mass transfer episode types are explained in the caption of Fig.C1.

Figure D2 .Figure E2 .Figure E3 .Figure E4 .
Figure D2.Typical formation histories of merging binary black holes with the most massive primary black holes according to our three different stellar wind models.For brevity, we only show a few stages of the most important steps in the evolution of these systems.The stage numbers correspond to the evolutionary stages as shown in Fig.1.For all the sources here, we assume  = 2.5,  = 0.3 and  ad,rad = 7.5.The asterisk for Stage 2 for Model III means that the first mass transfer phase is Case A

Table 2 .
Hjellming & Webbink (1987)sfer stability criteria used in this study for each stellar evolutionary phase.For HG and CHeB stars with radiative envelopes we show two values (4 and 7.5) as we test both values in our study.HW87 stands forHjellming & Webbink (1987)