\section{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. \citealt{Langer_2012}).
Unfortunately, it is challenging to directly observe them as they are rare and short-lived. On the other hand, the population of merging compact binaries, observed as gravitational wave (GW) sources, offer a unique but indirect way to study the evolution of these objects.
Since the first detection of GW, 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 (\citealt{Abbott_2021}).

There are many proposed channels for the formation of merging stellar mass binary black holes. These include the formation via dynamical interactions, such as in globular clusters (e.g. \citealt{Kulkarni1993Natur.364..421K}, \citealt{Sigurdsson1993Natur.364..423S}, \citealt{PortegiesMcMillian2000}, \citealt{Rodriguez_2016}, \citealt{DiCarlo2020}) or in nuclear clusters (e.g. \citealt{AntoniniPerets2012ApJ...757...27A}) or in hierarchical, isolated triples (e.g. \citealt{Antonini2017ApJ...841...77A}, \citealt{Alejandro2021}, \citealt{Silsbee2017ApJ...836...39S}, \citealt{Stegemann2022arXiv220316544S}, \citealt{Martinesz_https://doi.org/10.48550/arxiv.2105.01671}). Formation in AGN discs has also been proposed, where besides dynamical interactions, gaseous drag can also contribute to the hardening of the binaries  (e.g. \citealt{McKerna2020}, \citealt{Bird2017MNRAS.464..946S}, \citealt{Bartos2017ApJ...835..165B}). The possibility of merging binary black holes originating from population III stars (\citealt{Belc2004}, \citealt{Kinugawa10.1093/mnras/stu1022} \citealt{Inayoshi10.1093/mnras/stx757}) or from primordial black holes (\citealt{BirdPhysRevLett.116.201301}, \citealt{Sasaki2018}) has also been studied. Formation from isolated massive binaries could occur via mass transfer episodes (i.e. the so-called classical isolated binary channel) or with chemically homogeneously evolving stars (\citealt{de_Mink_2016}, \citealt{Mandel_2016}, \citealt{Marchant_2016}, \citealt{Buisson2020MNRAS.499.5941D}, \citealt{Riley2021MNRAS.505..663R}).

The classical isolated binary channel is perhaps the most studied formation path (e.g. \citealt{Paczynski1976IAUS...73...75P}, \citealt{Heuvel76}, \citealt{Tutukov1993MNRAS.260..675T}, \citealt{Dominik_2012}, \citealt{Belczynski_2016}, \citealt{Giaccobbo2018MNRAS.480.2011G}, \citealt{Stevenson_2017}, \citealt{Kruckow_2018}, \citealt{KlenckiNelemans2018}, \citealt{Neijssel_2019}, \citealt{marchant2021role}, \citealt{Bavera2021}).
In these studies, typically two main sub-channels are identified. The first one (CEE channel), in which the key step in the formation of close BH-BH binaries is the so-called common envelope evolution (e.g. \citealt{Ivanova_2013}). It is often assumed that envelope ejection is only possible after the donor has reached the core helium burning branch (\citealt{IvanovaTaam2004}, \citealt{Dominik_2012}).
Depending on the assumed common envelope efficiency and the binding energy of the envelope of the donor star, this channel can produce a merger rate that is broadly consistent with the currently observed LIGO rate (\citealt{Abbott_2021}, \citealt{Mandel_2022}).
However, recent studies (\citealt{marchant2021role}, \citealt{Klencki_2020}) pointed out that the binding energy of evolved stars with radiative envelopes can be underestimated by rapid population synthesis studies, which commonly use prescriptions for determining the binding energy. Furthermore, a deep convective envelope could potentially be developed at a significantly cooler effective temperature (\citealt{Klencki_2021}) than previously assumed. These two developments could imply an appreciably lower predicted merger rate for this channel.

The second dominant channel (stable channel) involves two subsequent stable mass transfer episodes. Here, the key point is that the orbit can shrink significantly even due to a stable phase of mass transfer, if the mass ratio is relatively high and angular momentum loss is efficient, such as with a black hole accretor (\citealt{vandenheuvel17}).

Earlier studies did not predict this formation path to be significant, but the detailed simulations of  \citet{Pavlovskii_2016}
 suggest that stable mass transfer is more readily achieved than previously assumed for an evolved donor star with a radiative envelope and a black hole accretor.
 This implies that the parameter space for the formation of GW sources via this channel is also larger (\citealt{KlenckiNelemans2018}) and in fact it can be the dominant channel within the isolated binary scenario (\citealt{olejak2021impact},
 \citealt{Lieke2022ApJ...931...17V},
 \citealt{Bavera2021}, \citealt{marchant2021role}, \citealt{Andrews_2021}).
 This channel is predicted to have considerably longer time delays than the CEE channel. This leads to an increased merger rate density in the near universe (z$\approx$0), compared to that of the CEE channel (\citealt{Neijssel_2019}). This is because the cosmic star formation density is expected to monotonically increase up to $z\sim2$ (see e.g. \citealt{Madau_2014}), .

The properties of the  merging binary black hole population depend on the binary processes that they go through during their lifetime. This, in principle means that observations by LIGO and Virgo (\citealt{ligovirgo}) can be used to constrain these processes. Unfortunately, there are currently too many uncertainties in massive stellar evolution to make meaningful inference, even if it was known which subset of GW observations came from isolated binary evolution (if any). Nevertheless, it is still essential to understand how the uncertainties of the binary 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 the first part of this paper, we concentrate on the uncertainties regarding the first phase of mass transfer. In this case, the accretor is a non-compact object. How much matter is accreted and what is the angular momentum that is lost from the binary is still poorly understood. These two quantities govern the further evolution of mass ratios and orbital separations of the binary systems. This in turn determine which of the previously mentioned sub-channels dominate.
This is because each sub-channel needs to have a specific and distinct range of mass ratios and orbital separations at the time when the initial secondary star is about to fill its Roche-lobe.

For example, the progenitors of the stable channel have relatively low orbital separations (typically in the order of $\sim100\,R_{\odot}$) and relatively high mass ratios, where the latter is required such that the orbit shrinks sufficiently due to the subsequent, stable phase of mass transfer. At the same time, the mass ratios cannot be above the so-called critical mass ratio, above which the mass transfer phase would become dynamically unstable. The latter quanitity depends on the structure of the envelope of the donor star and on how exactly the mass transfer episode proceeds.

The progenitors of the CEE channel, on the other hand, have large separations. These need to be large enough so that the initial secondary star fills its Roche-lobe with a core-helium burning donor. In this case, the mass ratios always have to be above the critical mass ratio.

Therefore, the first mass transfer phase has a significant effect on the predicted merger rates of the two main channels.
To understand this effect better, we test different assumptions for the angular momentum loss and the fraction of mass ejected during a mass transfer phase with a non-compact accretor.
We also vary the mass transfer stability criteria of evolved stars with radiative envelopes to understand the consequences of studies such as \citet{Ge_2020}, \citet{Pavlovskii_2016}.
Furthermore, we also test the implications of recent developments regarding the evolutionary stage at which giant stars develop convective envelopes, following \citet{Klencki_2021}.
We create models using all possible combinations of parameter variations. This allows us to understand the interrelated effects of uncertainties.

In the second part of the paper, we investigate the effects of uncertainties in stellar winds. Both theoretical (\citealt{BJ22https://doi.org/10.48550/arxiv.2203.08218}, \citealt{bjrklund2020new}, \citealt{Krti_ka_2017},  \citealt{Sundqvist_2019}) and observational (e.g. \citealt{Fullerton2006}) studies suggest that mass loss rates of O/B stars could be overestimated by a factor of 2-3 by the prescription of \cite{Vink_2001}, which is commonly used in stellar evolutionary codes. As we will see, the impact of lowered mass loss rates on GW sources significantly depends on our assumption of 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 paper is organised as following. In section \ref{section:isolated_brief}, we briefly review the classical isolated binary formation channel to introduce the terminology used in this paper. In section \ref{section:methodology}, we describe the code used in this study. In section \ref{sub:binzeta}, we show our results that highlight the importance of uncertainties in the first mass transfer phase. In section \ref{sub:stellarwinds}, we show how decreasing the line-driven winds for O/B stars and/or WR stars by a factor of 3 can change the merging binary black hole population. Finally, in \ref{sub:conclusion}, we summarise our main findings.

\section{The classical isolated binary formation channel}
\label{section:isolated_brief}
We now give 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 rather than to give a detailed overview on this subject. For the latter see reviews such as \citet{m2018merging}, \citet{Mapelli2020}, \citet{Postnov_2014}.

Black hole binaries with circularised orbits merge within the Hubble time if their orbital separation is not larger than a few tens of solar radii (\citealt{Peters}, \citealt{m2018merging}). If a massive binary starts out at such a short period, then the donor fills its Roche-lobe on main sequence. Therefore, GW sources of the isolated binary channel must originate from interacting binaries.

We show a schematic drawing of the three most common formation paths of GW sources according to our simulations in Figure \ref{fig:isolatedchannelsl}. At birth, the binaries are not interacting and their orbits widen due to stellar winds (stage 1). Then, the initial primary star fills its Roche-lobe and the first mass transfer phase starts (stage 2). This always occurs in a dynamically stable manner for the channels considered here.
This phase ends with the initial primary star losing its hydrogen envelope. At this stage, the donor star is a stripped helium star, which could experience intense stellar winds and may be observed as a Wolf-Rayet star (stage 3). The mass ratio and the orbit of the system depend on how much matter was accreted by the secondary and how much angular momentum was carried away from the binary by the expelled matter 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 this, the secondary evolves, starts expanding too, 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 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  (\citealt{Ivanova_2013}). As the common envelop 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 in the orders of $\sim 1000 \,R_{\odot}$.
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 (CEE) 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 still shrink sufficiently, given, the mass ratios at the onset of the mass transfer phase are relatively high (stage 5a). For example with mass ratios $q = M_{donor}/M_{accretor} \sim 3$, the orbit for a $M_{accretor} \approx 30\,M_{\odot}$ shrinks by $\sim 100\,R_{\odot}$. Depending on the binary physics assumptions, the stable channel can also contribute to a significant number of GW sources.

All channels discussed above involve two subsequent phases of mass transfer. 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. \citealt{Tout1997MNRAS.291..732T}). For example, systems with an initial primary mass of $M_{ZAMS,1} = 100\,M_{\odot}$ can evolve in such a way, even if their initial mass ratio are as close to unity as $q_{ZAMS} = M_{ZAMS,2}/M_{ZAMS,1} = 0.999$. We also find that only these binaries can form gravitational wave sources in such numbers so that they are compatible with the currently observed merger rate from LIGO and Virgo (\citealt{Abbott_2021}). Therefore, in the rest of the paper, we always assume that in the three dominant channels introduced above, the accretor is a black hole by the onset of the second phase of mass transfer.

\begin{figure}
\hspace*{-0.5cm}
  \includegraphics[scale = 0.40]{isolated_channels.pdf}
  \caption{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.}
  \label{fig:isolatedchannelsl}
\end{figure}

\section{SeBa and model variations}
\label{section:methodology}
For our simulations, we use the rapid population synthesis code SeBa\footnote{ https://github.com/amusecode/seba} (\citealt{Portegies96}, \citealt{Toonen_2012}). An up-to-date description of the code can be found in \cite{Toonen_2012}, therefore we only describe elements which are especially relevant for this study, or which have been changed with respect to \citet{Toonen_2012}.
\subsection{Stellar wind prescriptions}
\label{subsection:swp}
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. We further distinguish line-driven winds of main sequence and evolved stars (i.e. optically thin line-driven winds) and line-driven winds of stripped helium stars (i.e. Wolf-Rayet winds or optically thick line-driven winds). Figure \ref{fig:Wind1} shows approximately where these winds are dominant on the Hertzsprung-Russel diagram. In the following, we briefly summarise the stellar wind prescriptions we use in this study, while Table \ref{tab:mdotapplied} shows at what evolutionary stage these prescriptions are applied. A more detailed overview can be found in the appendix \ref{appendix:subsection:swp}.

\begin{table*}
\caption{Summary of stellar wind prescriptions applied to different evolutionary stages of massive stars in SeBa. HD stands for Humphreys-Davidson limit. The scaling factor $f_{wind}$ is 1 for the default model and 1/3 for the 'scaled-down' model. $\dot{M}_{V01}$ stands for \citet{Vink_2001}, $\dot{M}_{NJ}$ for \citet{Jager}, $\dot{M}_R$ for \citet{Reimers}, $\dot{M}_{WR}$ for \citet{Sander_2020}, $\dot{M}_{VW}$ for \citet{1993ApJ...413..641V} and  $\dot{M}_{LBV}$ is assumed to be $1.5\cdot10^{-4}\,M_{\odot}yr^{-1}$ following \citet{Belczynski_2010}}
\begin{tabular}{|l|l|ll|}
\hline
\textbf{Main sequence} & \textbf{Hertzsprung gap \& Core helium burning} & \textbf{AGB phase}                                     & \textbf{Helium star \& helium giant} \\ \hline
\multicolumn{1}{|c|}{\textit{if $T_{eff} \leq 50 kK$:}} &
  \textit{if $T_{eff} \geq 8 kK$:} &
  \multicolumn{1}{l|}{\textit{By default:}} &
  \textit{Always:} \\
\multicolumn{1}{|c|}{$f_{wind}\cdot\dot{M}_{V01}$} &
  \multicolumn{1}{c|}{$f_{wind}\cdot\dot{M}_{V01}$} &
  \multicolumn{1}{l|}{$\dot{M} = max(\dot{M}_{VW},\dot{M}_{R},\dot{M}_{NJ})$} &
  \multicolumn{1}{c|}{} \\ \cline{1-3}
\textit{otherwise:}    & \textit{if $T_{eff} \leq 8 kK$:}                & \multicolumn{1}{l|}{\textit{if star beyond HD limit:}} & \multicolumn{1}{c|}{}                \\
\multicolumn{1}{|c|}{$f_{wind}\cdot\dot{M}_{NJ}$} &
  \multicolumn{1}{c|}{$\dot{M} = max(\dot{M}_{NJ},\dot{M}_{R})$} &
  \multicolumn{1}{c|}{} &
  \multicolumn{1}{c|}{\textit{$f_{wind,WR}\cdot\dot{M}_{WR}$}} \\ \cline{1-2}
                       & \textit{if star beyond  HD limit:}              & \multicolumn{1}{c|}{$\dot{M}_{LBV}$}                   &                                      \\
                       & \multicolumn{1}{c|}{$\dot{M}_{LBV}$}            & \multicolumn{1}{c|}{}                                  &                                      \\ \hline
\end{tabular}
\label{tab:mdotapplied}
\end{table*}

Line driven winds of O/B stars are modelled in SeBa using the mass-loss rates from \citet{Vink_2001} as long as the star is within the grid defined by \citet{Vink_2001}, otherwise, the empirical formula from \citet{Jager} is adopted. These mass-loss rates are applied throughout the evolution of stars until they reach $T_{eff} \sim 8000K$ or they finish their core-helium burning phase (see Table \ref{tab:mdotapplied}).

 We note that \citet{Vink_2001} estimates the global metallicity dependence to be $\dot{M} \propto Z^{0.69}$ for O stars, by assuming a metallicity dependence of the final wind velocity to be $v_\infty \propto Z^{0.13}$, following \citet{Leitherer92}.
This is consistent with the observations as shown by \citet{Mokiem07} (although their results still depend on the findings of \citet{Leitherer92}). We note that however, \citet{Kriticka06} and \citet{bjrklund2020new} 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 $\dot{M} \propto Z^{0.85}$ (i.e. ignoring the metallicity dependence of $v_\infty$), so that we are consistent with other population synthesis studies  \citep[e.g.][ etc.]{Dominik_2012, Stevenson_2017, Mapelli_2016,Giacobbo_2017}.

If a main-sequence star is outside of the grid defined by \citet{Vink_2001}, the empirical formula of  \citet{Jager} is used.

We note that we assume a metallicity scaling $\sim Z^{0.85}$, which is different from the what was originally suggested by \citet{Kudritzki87}. 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 $ \dot{M}_{LBV}\sim 10^{-5} -10^{-3} M_{\odot}yr^{-1}$, although this is highly uncertain. Even less known about their possible eruptions, in which huge amount of mass is lost in a very short time (\citealt{Humphreys}, \citealt{Smith}, \citealt{lbvvink2012ASSL..384..221V}).
For the mass-loss rates of LBV stars, we follow the assumption of \citet{Belczynski_2010}, i.e. the mass loss rate consant and has a value of $\dot{M}_{LBV} = 1.5\cdot10^{-3}\,M_{\odot}yr^{-1}$.

If the star becomes a cool giant (ie $T_{eff} \leq 8000 K$), we calculate the mass loss rate according to Reimer's empirical formula \citep{Reimers}. We compare it with the mass loss from \citet{Jager} and we take the maximum value.

For the mass loss rates of thermally pulsating AGB stars, we use the prescription of \citet{1993ApJ...413..641V}.

Mass loss rates for helium stars and helium giants (Wolf-Rayet star winds) are calculated according to \citet{Sander_2020}.

\subsection{Treatment of binary interactions}
\label{sec:masstransfers}

In this section, we summarise how binary interactions are treated in SeBa. The physical processes that can significantly influence the orbital evolution of a massive binary are tidal interactions, mass transfer episodes, mass loss due to stellar winds, and gravitational waves.

In SeBa, we assume tidal interactions circularise the orbit by the onset of the mass transfer (\citealt{Portegies96}). Change in the orbital separation and eccentricity due to gravitational waves emission are calculated according to \citealt{Peters}.

As the stars of the binary evolve, their radii grow, which can lead to Roche-lobe overflow. If one of the stars fills its Roche-lobe, the mass transfer episode is initiated. The Roche-lobe is calculated according to Eggleton's formula \citep{EggletonApprox}.

In order to determine the new orbital separation after a mass transfer episode, we use the relation \citep[e.g.][]{BookCompactXray, 1994vandenheuvel}:
\begin{equation}
\label{eq:dotaa}
    \frac{\dot{a}}{a} = -2 \frac{\dot{M}_{d}}{M_{d}}\left[ 1 - \beta \frac{M_{d}}{M_{a}} - \left(\gamma+\frac{1}{2}\right)\frac{(1 - \beta)M_{d}}{M_{d} + M_{a}} \right],
\end{equation}
where $M_{d}$, $M_{a}$ are the mass of the donor and the accretor, respectively, $\beta$ is the mass accretion efficiency, ie. the amount of mass that can be accreted and $\gamma$ is the ratio of specific angular momentum that leaves the system and the total specific angular momentum of the binary, ie. $dJ/dM_{tot} = \gamma J/M_{tot}$, with J being the angular momentum of the binary and $M_{tot}$ being the total mass of the system.
Although, $\beta$ and $\gamma$ is expected to depend on the parameters on the binary, in rapid population synthesis codes, it is commonly assumed that this value is constant.
If $\gamma$ is constant, Equation \ref{eq:dotaa} can be integrated:
\begin{equation}
\label{eq:stable_mt_eq}
\frac{a}{a_i} = \left(\frac{M_{d}}{M_{d,i}}\frac{M_{a}}{M_{a,i}}\right)^{-2}\left(\frac{M_{tot}}{M_{tot,i}}\right)^{2 \gamma +1},
\end{equation}
where the subscript 'i' stands for initial, that is just before the 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 $\gamma = M_{d}/M_{a}$. This mass transfer mode is the so-called isotropical reemission. The change in orbital separation for isotropical reemission in the limit of $\beta \rightarrow0$  \citep[e.g.][]{soberman1997stability, vandenheuvel17}:
\begin{equation}
    \label{eq:iso}
     \frac{a}{a_{i}} = \frac{M_{tot,i}}{M_{tot}}\left(\frac{M_{d,init}}{M_{d}}\right)^2\exp\left(2\frac{M_{d}- M_{d,i}}{M_{a}}\right).
\end{equation}
If the mass ratio is sufficiently high, the shrinking in this mass transfer mode becomes efficient, in which case GW sources can be formed via the stable channel.

When the members of the binary lose mass via stellar winds, we assume that a fraction of it is accreted via Bondi-Hoyle accretion \citep{bondihoyle}, while the rest leaves the system with a specific angular momentum of the donor, which corresponds to $\gamma = M_{a}/M_{d}$.
In this case, the orbit widens as:
\begin{equation}
\label{eq:Jeans}
    \frac{a}{a_{i}} = \frac{M_{tot,i}}{M_{tot}}.
\end{equation}

\subsection{First phase of mass transfer }

There are currently numerous uncertainties regarding mass transfer episodes with non-degenerate accretors (see e.g. \citealt{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. First of all, it is important whether an accretion disc is formed or if the overflowing matter hits the accretor directly (e.g. \citealt{LubowShu1975ApJ...198..383L}). In the former case, the accretor is expected to spin up and reach critical rotation very quickly (\citealt{Packet1981}). Whether accretion really ceases, once this limit is reached is debated (see e.g. \citealt{PophamNarayan1991ApJ...370..604P}). As the centrifugal force cannot remove the unaccreted material entirely from the binary, it is not clear what mechanism (jets near the accretror, radiation pressure, enhanced stellar winds, etc.), if any, would be capable to do so (e.g. \citealt{Langer_2012}).
Another physical process that can hinder the donor star to accrete efficiently is the response of its radius due to mass transfer. If the mass transfer rate occurs on a higher rate than the thermal timescale of the donor star, then it falls out of thermal equilibrium, expands and possibly cannot accommodate more donated mass (see e.g. \citealt{Hurley2002}).

\citet{Langer_2020} used detailed stellar simulations to study O/B star-black hole systems. In this study, it was assumed that stars cannot accrete above critical rotation. It was found that the mass transfer efficiency for systems with donors that already evolved off the main sequence tends to be very low ($\beta \approx 0.1$, on average). On the other hand, if the donor is still a main sequence star, the orbit is so short, that tidal effects can counteract the spinning up of the accretor, which could result in significantly higher mass transfer efficiencies.
This is broadly consistent with the observational studies of \citet{Shao_2016} and \citet{Petrovic}, who inferred a low mass transfer efficiency based on galactic WR-O/B binaries. These systems have likely initiated their mass transfer phase with evolved donors.
Also, \citet{2007deMink} found that it is not adequate to use a single value for the mass transfer efficiency; they discuss indications that mass transfer proceeds more conservatively for shorter period binaries. At the same time, other observations are more consistent with near-conservative mass transfers episodes \citep[e.g][]{Schootemeijer_2018}. \citet{Vinciguerra_2020} found that the observed population of Be X-ray binaries in the SMC can only be reproduced by population synthesis, if a mass transfer efficiency $\beta = 0.5$ or higher is adopted. We note however, that these Be stars have typically somewhat lower masses that the systems discussed in this paper.

If the accretor is not a remnant, we assume that $\gamma = 2.5$, following \citet{Portegies96}. We also test $\gamma = 1$, following \citet{Podsiadlowski92} and \citet{Belczynski2008}. We note that in other rapid population synthesis studies the assumption of $\gamma = M_{d}/M_{a}$ is also commonly made i.e. when the mass expelled from the system is assumed to carry away a specific angular momentum equal that of the accretor  \citep[e.g.][]{Riley_2022, Kruckow_2018,Stevenson_2017}.
In this study, we test two, constant values for mass transfer efficiency when the accretor is a non-compact object; $\beta = 0.3$ and $\beta =0.7$. When the accretor is a neutron star or black hole, we assume that the accretion is Eddington-limited and $\gamma = M_{d}/M_{a}$.

\subsection{Mass transfer episode types based on the evolutionary phase of the donor}
\label{sub:typesofmt}

We distinguish the following mass transfer phase types:

\textit{Case A}: the donor is a main sequence star (see e.g. \citealt{KoushikSen2022A&A...659A..98S}).
If the period of the system is sufficiently short, the secondary might also fill its Roche-lobe, leading to the formation of contact systems (\citealt{Menon2021MNRAS.507.5013M}, \citealt{Wellstein2001A&A...369..939W}, \citealt{Pols1994A&A...290..119P}).
If the binary remains in a semi-detached phase throughout then mass transfer episode, then typically three sub-phases can be distinguished. The first one is the fast Case A phase, during which the mass transfer rate occurs on a thermal timescale. In this case, the donor star is still more massive than the secondary and the orbit shrinks as the result of the mass transfer phase. This is followed by the slow Case A, which roughly starts at the inversion of mass ratio. During this phase, the mass transfer rate occurs on a nuclear timescale. These binaries can be observed as Algol systems. The last phase (Case AB) occurs when the donor star evolves off main sequence and refills its Roche-lobe again as a Hertzsprung gap giant (i.e. hydrogen shell burning giant with an inert helium core). In the rest of the paper when we refer to Case A mass transfer phase, we mean all these three sub-phases together.
 The outcome of Case A mass transfer phase is 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 that of a system with the same initial masses, but which evolved through a mass transfer episode where the donor star had already evolved off the main sequence (e.g. see \citealt{Langer_2020}).
In SeBa, Case A mass transfer episode is the most challenging to model accurately. There are two major reasons for this. Firstly,  the stellar tracks of \cite{Hurley_2000} do not track the mass of the developing helium core of the main sequence star. 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 \ref{subsec:mt_stability}). This means that the radius response of the donor during a mass transfer phase is assumed to be the same, regardless of how much mass had already left the star.
The consequence of these two points is that, in principle, there is no limit to the amount of mass that can be transferred to the accretor as long as the donor keeps refilling its Roche-lobe. Therefore, the mass of the black hole that the donor eventually forms can be severely underestimated by codes like SeBa (for a different apporach in a rapid population synthesis code, see e.g. \citealt{Kruckow_2018}). Therefore, the outcomes of Case A mass transfer episodes predicted by codes based on the stellar tracks of \cite{Hurley_2000} should be treated with caution. Nevertheless, we still show such systems in this work for completeness. However, we always clearly distinguish those systems from the other types.

\textit{Case B}: the donor is on the Hertzsprung gap. During this phase the radius of the donor increases very rapidly, typically by $200-400 \,R_{\odot}$ at  $Z \lesssim 0.001$. In fact, for a large range of metallicities, this is expected to be the most common mass transfer type (assuming a log-uniform distribution of initial periods). The mass transfer phase ends with donor losing its hydrogen envelope, leaving a naked helium star behind \citep[but see][]{Laplace2020A&A...637A...6L}. Since the Hertzsprung gap phase lasts only for a few $\sim 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. Following the reasoning of \citealt{IvanovaTaam2004}, we assume that these donors cannot survive the common envelope phases.

\textit{Case C}: the donor is in its core helium burning phase.
We distinguish two sub-categories. Core helium burning donor with radiative envelope (Case Cr) and with deep convective envelope (Case Cc). We assume that the mass transfer stability criteria is the same for Case B and fore Case Cr. However, the difference is that the core helium burning phase lasts orders of magnitude longer with a slower expansion in radius. Therefore for such donors, the mass of the remnant is dependent on the initial separation. Moreover, the effects of LBV winds are strong; they can widen the orbit significantly and decrease the mass ratio, which in turn changes the outcome of the mass transfer phase. Donors with deep convective envelopes have very different envelope structures and therefore different mass transfer stability criteria. Unstable mass transfer is more readily realised for these systems (see Table \ref{tab:zeta2q}).

\subsection{Mass transfer stability criteria and treatment of mass transfer}
\label{subsec:mt_stability}
The stability of mass transfer is determined by comparing the  so-called mass-radius exponents (\citealt{soberman1997stability}):
\begin{equation}
    \zeta_L = \frac{d \log R_L}{d \log M}
\qquad
    \zeta_{ad} = \left(\frac{d \log R_{d}}{d \log M}\right)_{ad}
\qquad
    \zeta_{th} = \left(\frac{d \log R_{d}}{d \log M}\right)_{th},
\label{eq:zetas}
\end{equation}
where $R_{d}$ is the radius of the donor star, $\zeta_L$ describes how the Roche-lobe radius reacts to mass overflow, while $\zeta_{ad}$ and $\zeta_{th}$ describes how the radius of the donor reacts to mass loss during mass transfer on dynamical, and thermal timescale, respectively.
We then distinguish three different mass transfer modes: stable mass transfer on nuclear time scale ($\zeta_L \leq min(\zeta_{ad},\zeta_{eq})$), stable mass transfer on thermal timescale ($\zeta_{ad}\geq \zeta_L \geq \zeta_{eq}$.) and unstable mass transfer ($\zeta_L > max(\zeta_{ad}, \zeta_{eq}$). In the first two cases, we assume that the mass transfer rate is $\dot{m} = M_{d}/\tau$, where $\tau$ is the nuclear timescale in the first and the thermal timescale in the second case.
When the mass transfer becomes dynamically unstable, we assume that this results in common envelope evolution.

In SeBa, we determine the mode of the mass transfer in the following way. First, $\zeta_L$ is  determined numerically and then it is compared with the adopted $\zeta_{ad}$ and $\zeta_{th}$, which are summarised in Table \ref{tab:zetas}. As a major simplification we assume a constant $\zeta_{ad}$/$\zeta_{th}$ for a given stellar evolutionary phase.
Giants with deep convective envelopes tend to have low $\zeta_{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 an 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 $T_{eff,boundary})$. We test two assumptions. First, that a deep convective envelope develops at an effective temperature $logT_{eff} = 3.73\,K$, following \citet{IvanovaTaam2004} and \citet{Belczynski2008}. In the second model variation, we follow the prescription of \citet{Klencki_2021}. This prescription gives $T_{eff,boundary}$ as a function of luminosity and metallicity. The predicted values of $T_{eff,boundary}$ from \citet{Klencki_2021} are typically considerably cooler than $logT_{eff} = 3.73\,K$.

\begin{table*}
\caption{A summary of the mass transfer 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 for \citet{Hjellming}}
\begin{tabular}{@{}lcccccc@{}}
\toprule
 &
  \textbf{Main sequence} &
  \textbf{Hertzsprung gap} &
  \textbf{Core helium burning phase} &
  \textbf{AGB phase} &
  \textbf{Helium star} &
  \textbf{Helium giant} \\ \midrule
$\zeta_{ad}$ &
  \textit{4} &
  \textit{4 or 7.5} &
  \textit{\begin{tabular}[c]{@{}c@{}}For radiative envelope: 4 or 7.5\\ For fully convective envelope: HW87\end{tabular}} &
  \textit{HW87} &
  \textit{15} &
  \textit{HW87} \\
$\zeta_{th}$ &
  \textit{0.55} &
  \textit{-2} &
  \textit{\begin{tabular}[c]{@{}c@{}}For radiative envelope: -2\\ For fully convective envelope: 0\end{tabular}} &
  \textit{0} &
  \textit{1} &
  \textit{-2} \\ \bottomrule
\end{tabular}
\label{tab:zetas}
\end{table*}

\subsection{Common envelope evolution}

In this study, we model common envelope evolution by adopting the energy formalism \citep[e.g.][]{Webbink84,Heuvel76}. The energy used to unbind the envelope comes from the orbital energy of the binary as it is shrinking. The efficiency by which the energy is transferred to unbind the envelope is described by parameter $\alpha_{CE}$. Then the change in orbital separation can be given as:
\begin{equation}
    \frac{G M_{d} (M_{d}- M_{d,core})}{\lambda R_{d}} = \alpha_{CE} \left( \frac{GM_{a}M_{d,core}}{2 a_f} - \frac{G M_{d} M_{a}}{2 a_i} \right),
\end{equation}
here, the left hand side term is the binding energy of the envelope. $M_{d,core}$ is the mass of the helium core of the donor, $R_{d}$ is the radius of the donor, $a_i$ and $a_f$ are the initial and final orbital separation, respectively The term $\lambda$ describes the structure of the envelope \citep{Kool90}. Several studies published tabulated or fitted data for $\lambda$ for a range of masses and evolutionary stages \citep[e.g.][]{Klencki_2020, Claeys_2014, Xu_2010, Loveridge_2011, dewi}.
In this study, we assume a constant $\lambda = 0.05$, which is in reasonable agreement with the results of \citet{Xu_2010} for core helium burning giants (see also \citet{Dominik_2012}. We adopt $\alpha_{CE} = 1$.

\subsection{Supernova and natal kicks}
The mass of the remnant after black hole formation is computed based on the delayed supernova model from \citet{Fryer_2012}.
This prescription gives the proto-NS mass and the amount of fallback as a function of CO core mass, which in SeBa is obtained from the fits of \citet{Hurley_2000}.
The kick velocity for black holes is calculated as :
\begin{equation}
\label{eq:kick}
    v_{BH} = (1 - f_b)\left(\frac{M_{NS}}{M_{BH}}\right)v_{kick}
\end{equation}
Where $f_b$ is the fallback, $M_{NS}$ is the canonical neutron star mass $M_{NS} = 1.4 M_{\odot}$ and $v_{kick}$ is a random velocity kick drawn from distributions inferred from observations. The standard setting in SeBa for $v_{kick}$ is based on \citet{Verbunt2017} .

\subsection{Initial conditions}
\label{subs:icamv}
Observations suggest that about half of the stars are in binaries (or in higher order, hierarchical systems) and this fraction seems to increase for more massive stars \citep{stellarmult1}.
In particular, \citet{Sana2012} showed
that the binary fraction reaches $f_b \simeq 0.7$ for stars in the mass range $M \simeq 15-60\, M_{\odot}$. These results are based on observations of nearby ($\leq 3 kpc$), young (about 2 Myr old) open clusters and star associations. The same observations showed that the orbital period distribution of these young, massive binaries favour short period systems:
\begin{equation}
     f_p (\log p)  \propto ( \log p )^\pi, \quad \quad \text{for $ \log p\in
[0.15, 5.5]$},
\end{equation}
here $p$ are in days and  $\pi = -0.55$. Although, other studies suggested somewhat flatter distributions (e.g.  \citealt{Kobulnicky_2014} or \citealt{Dunstall} based on observations of B stars in the Tarantula Nebula in the Large Magellanic Cloud).

According to \citet{Sana2012}, the mass ratio distribution of massive stars is nearly flat, $f_q \propto q^{-0.1}$, while the distribution of eccentricities follows $f_e \propto e^{-0.42}$. A similar mass ratio distribution is inferred by \citet{Kobulnicky_2014}, while \citet{Dunstall} finds $f_q \propto q^{-2.9}$ for B stars in the Tarantula Nebula, and \citet{DunstallVIII} finds $f_q \propto q^{-1}$ or O stars in the Tarantula Nebula.

The effects of uncertainties in initial conditions were investigated by \citet{mink2015merger}, and they found that only the change in initial mass function alters significantly the properties of the merging compact binary object population \citep[but see][]{Klencki_2018}.
There are observational and theoretical indications that the initial mass function might not be universal  \citep[see][and references therein]{Chru_li_ska_2020}.
In our simulations, we assume these initial distributions are not correlated, however, in reality this might not be the case \citep[see e.g.][]{Moe_2017, Klencki_2018}.

In this paper we used the following initial conditions for our simuations:
\begin{itemize}
    \item Initial mass function:
    the initial mass of the primary is drawn from a distribution determined by \citet{Kroupa}. In the mass range $M_{ZAMS} = 20 - 100\, M_{\odot}$, this is $N_{IMF}(M_{ZAMS}) \simeq M_{ZAMS}^{-2.3}$.
    \item Initial mass ratio distribution:
    we assume a uniform mass ratio distribution between 0.1 and 1, where the mass ratio is defined as $q = M_{ZAMS,2}/M_{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 \citet{Sana2012}.
    \item Initial separation distribution: we assume a flat distribution in the logarithmic space of binary separation in the interval of $1 R_{\odot}$-$10^4 R_{\odot}$; $N_{a} \simeq log(a)$. This is equivalent to Opik law, ie. a uniform distribution in $\log p$ \citep{OpikLaw1924}.
    \item Initial eccentricity distribution: the initial eccentricity assumed to follow a  thermal distribution \citep{Heggie1975}, i.e. $f_e  \propto e$.
\end{itemize}
\begin{table*}
\centering
\caption{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 \ref{section:methodology}. For the detailed descriptions of stellar wind models see Table \ref{tab:mdotapplied}.
}
\label{tab:models}
\begin{tabular}{@{}lll@{}}
\toprule
\multicolumn{3}{l}{\textbf{Parameters not varied in this paper}} \\ \midrule
\textit{Parameter} &
  \textit{Model/value} &
  \textit{Label} \\
 &
   &
   \\
Single stellar tracks &
  \citet{Hurley_2000} &
  - \\
Tidal interactions &
Orbits are circularised by the time of the\\  mass transfer \citep{Portegies96} &
  - \\
SN prescription &
  Delayed model of \citep{Fryer_2012} &
  - \\
Natal kick velocity distribution &
  \citet{Verbunt2017} &
  - \\
Natal kick velocity scaling for BHs &
  $v_{kick,BH} = (1-f) \cdot v_{kick, NS} \cdot (M_{BH}/M_{NS})^{-1}$ &
  - \\
Common envelope treatment &
  $\alpha$-formalism with $\alpha\lambda = 0.05$ &
  - \\
Angular momentum loss mode with BH or NS accretor &
  $\gamma = M_{d}/M_{a}$ &
  - \\
Accretion efficiency for BH and NS accretors&
  Eddington limited accretion &
  a \\
 \midrule
\multicolumn{3}{l}{\textbf{Model variations}} \\ \midrule
\textit{Parameter} &
  \textit{Model/value} &
  \textit{Label} \\
 &
   &

   \\
   Stellar wind model &
1.) $f_{wind}$ = $f_{wind,WR}$ = 1 and $f_{LBV} = 1.5$ &
  Model I \\
 &
2.) $f_{wind}$ = 1/3, $f_{wind,WR}$ = 1 and $f_{LBV} =1.5$ &
  Model II \\
 &
3.) $f_{wind}$  = $f_{wind,WR}$ = 1/3 and $f_{LBV} = 1.5$  &
  Model III \\
 Angular momentum loss mode with non-compact accretors &
 1.) $\gamma = 2.5$ (\citealt{Portegies96}) &
  $\gamma = 2.5$ \\ &
2.) $\gamma = 1$ (\citealt{Podsiadlowski92}, \citealt{Belczynski2008}) &
  $\gamma = 1$ \\
 Accretion efficiency of non-compact accretors &
 1.) $\beta = 0.3$ &
  $\beta = 0.3$ \\ &
2.) $\beta = 0.7$ &
  $\beta = 0.7$ \\
Mass radius exponent of giants with radiative envelopes \\ &
  1.) $\zeta_{ad,rad} = 4$ &
  $\zeta_{ad,rad} = 4$  \\ &
2.) $\zeta_{ad,rad} = 7.5$ &
$\zeta_{ad,rad} = 7.5$  \\
Boundary of deep convective envelope &
1.) Prescription of \citet{Klencki_2021} &
  $T_{eff,K}$ \\ &
2.) convective above $log_{10}T_{eff} = 3.73\,K $ (\citealt{Belczynski2008}) &
  $T_{eff,IT}$ \\ \bottomrule
\end{tabular}
\end{table*}

\subsection{Merger rate density at the local universe and detection rates}
\label{subs:rates}

The merger rate density and the predicted detection rates are important characteristics of a synthetic population of merging binary black holes as they allow comparison to the catalogue of gravitational wave source detections made by LIGO and Virgo.
In order to determine these two quantities, we follow a formalism  similar to the one outlined in \citet{Dominik_2015}.

First, we define the merger efficiency as the number of mergers originating from a black hole binary with masses $M_{BH,1}$ and $M_{BH,2}$ at metallicity $Z$ 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 \ref{subs:icamv}):
\begin{equation}
\label{eq:mereff}
    \epsilon_{merger} (M_{BH,1}, M_{BH,2}, Z) = f_{pm} \cdot \frac{N_{merger}}{N_{sim, interacting}},
\end{equation}
where $N_{merger}$ is the number of mergers originating from a BH binary with given masses at a given metallicity, $N_{sim, interacting}$ is the number of interacting binaries in our simulation and $f_{pm}$ is the simulated parameter space as a fraction of the complete parameter space. We determine the latter as:
\begin{equation}
    f_{pm} =  f_{interacting}\cdot\int_{20\,M_{\odot}}^{100\,M_{\odot}} N_{IMF}(M) dM,
\end{equation}
here $N_{IMF}$ is the normalised initial mass function of \citet{Kroupa}, $f_{interacting}$ is the number of interacting binaries as a fraction of all massive binaries, which we assume to be 0.7, following \cite{Sana2012}. Therefore, in our rate calculations, we only take interacting binaries into account, and we assume that such systems constitute 70$\%$ of all massive binaries. Furthermore, we also assume that the mass ratios of massive binaries range from $0.1-1$.

The merger efficiency can be converted to a merger rate density by assuming a metallicity specific star formation density model $\textrm{SFRd}^{*}(Z,z)$:
\begin{equation}
    \frac{dN}{dVdzdt}(z) = \iiint\frac{\textrm{SFRd}^{*}(Z,z_{ZAMS})}{\Tilde{M}}\cdot \epsilon_{merger}  dZ dM_{BH,1} dM_{BH,2},
\end{equation}
where $z_{ZAMS}$ is the redshift at the birth of the binary. The time between the birth of a binary and the merger is called time delay. This quantity can be directly obtained from the population synthesis simulations. In order to determine $z_{ZAMS}$ for a given $z$ and time delay, we use the relationship for lookback time:
\begin{equation}
    t_{d} = \frac{1}{H_0}\int_{z}^{z_{ZAMS}} \frac{dz'}{(1 + z') E(z')}dz',
\end{equation}
where $E(z) = \sqrt{\Omega_m(1+z)^3 + \Omega_{\lambda}}$, with $\Omega_M = 0.3$, $\Omega_{\lambda} = 0.7$ and $H_{0} = 70 kms^{-1}Mpc^{-1}$.

The star formation rate density in a metallicity bin centered at $Z$ with a width of $2\Delta Z$ is $\int_{\Delta Z}^{\Delta Z} SFRd(z)\cdot f_{met}(Z,z) dZ$. Here $f_{met}(Z,z)$ gives the distribution of metallicities of binaries at redshift $z$ as:
\begin{equation}
    f_{met}(Z,z) = \frac{1}{\sigma \sqrt{2\pi}}\exp\left(\frac{(\log_{10}(Z) - \mu(z))^2}{2\sigma^2}\right),
\end{equation}
where $\sigma = 0.5$ and the mean metallicity varies with redshift as given by \citet{Madau_2017}:
\begin{equation}
    \mu(z) = \log_{10}(Z_{\odot}\cdot 10 ^{0.153 - 0.074z^{1.34}}) - 0.5ln(10)\sigma^2,
\end{equation}
where we took $Z_{\odot} = 0.02$. Furthermore, the star formation rate density $\textrm{SFRd}(z)$ is given by \citet{Madau_2014}:
\begin{equation}
    \textrm{SFRd}(z) = \frac{0.01\cdot(1 + z)^{2.6}}{1 + ((1+z)/3.2)^{6.2}} \, M_{\odot}yr^{-1}Mpc^{-3}.
\end{equation}

\section{Effect of uncertainties in stable mass transfer}
\label{sub:binzeta}

We now discuss the importance of uncertainties of the stable mass transfer phase and their predicted impact on the merging binary black hole population. We do this by varying four parameters in our simulations, namely:
\begin{itemize}
    \item $\gamma$, 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 $\gamma$ only for non-compact accretors.
    \item $\beta$, which describes the mass transfer accretion efficiency. We vary $\beta$ only for non-compact accretors.
    \item $\zeta_{ad,rad}$, which  is the mass-radius exponent for giant donors with radiative envelopes. It describes the boundary between stable and unstable mass transfer.
    \item $T_{eff,boundary}$, which is the effective temperature at which a deep convective envelope is expected to develop.
\end{itemize}

 The details of these model variations are summarised in Table \ref{tab:models}. Using only our standard stellar wind model, we run simulations with all possible combinations of $\gamma$, $\beta$ $\zeta_{ad,rad}$ and $T_{eff,boundary}$. We simulate $10^6$ binaries at each value of our metallicity grid for each model variations. Our metallicity grid is defined as $Z = [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 $\zeta$ to the critical mass ratio for each model variation in table $\ref{tab:zeta2q}$.

\begin{table*}
\caption{Conversion between the mass transfer stability criteria for giants and the critical mass ratio ($q_{crit} = M_{d}/M_{a}$) for different assumptions on the mass transfer efficiency and the angular momentum loss 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 $\zeta_{ad,rad}$, whereas donors with convective envelopes are described by $\zeta_{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})\approx 3.73\,K$ with initial masses $M_{ZAMS} = 50 - 100\, M_{\odot}$ with a metallicity independent LBV mass loss rate of $\dot{M}_{LBV} =  1.5\cdot10^{-4} M_{\odot}yr^{-1}$.}
\hspace*{1.5cm}
\begin{tabular}{|c|c|c|c|c|c|l|}
\hline
{\textbf{Mass transfer stability criteria}} & \multicolumn{6}{c|}{\textbf{Angular momentum loss}}                                                                                \\ \cline{2-7}
                                                           & \multicolumn{2}{c|}{$\gamma = 2.5$} & \multicolumn{2}{c|}{$\gamma = 1.0$} & \multicolumn{2}{c|}{$\gamma = M_{d}/M_{a}$}    \\ \cline{2-7}
                                                           & \multicolumn{6}{c|}{\textbf{Mass transfer efficiency}}                                                                                \\ \cline{2-7}
                                                           & $\beta = 0.3$    & $\beta = 0.7$    & $\beta = 0.3$    & $\beta = 0.7$    & \multicolumn{2}{c|}{\textit{Eddington limited accretion}} \\ \hline
$q_{crit}$ with $\zeta_{ad,rad} = 4$                         & $3.8$            & $2.9$            & $6.2$            & $3.4$            & \multicolumn{2}{c|}{$3.2$}                                \\ \hline
$q_{crit}$ with $\zeta_{ad,rad} = 7.5$                      & $8.7$            & $5.2$            & $11.6$           & $5.7$            & \multicolumn{2}{c|}{$5.5$}                                \\ \hline
$q_{crit}$ with $\zeta_{ad} = \zeta_{HW}$                  & $0.9 - 0.7$      & $0.9 - 1.1$      & $1.5 - 1.9$      & $1.1 - 1.4$      & \multicolumn{2}{c|}{$1.4 - 1.7$}                          \\ \hline
\end{tabular}
\label{tab:zeta2q}
\end{table*}

Typically, the two most dominant channels within the isolated binary formation scenario are the stable and CEE channel (subsection \ref{subs:icamv}). We find that this is the case for all of our model variations. However, their relative importance drastically changes depending on our assumptions of the binary physics.
In the following, we concentrate only on these two evolutionary channels, as other formation paths typically do not account for more than $\sim3\%$ of the total number of mergers.

\subsection{The effect of the first phase of mass transfer}
\label{sub:effect_of_first_mt}

 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.

In principle, the demographics at the onset of the second phase of mass transfer depend on the initial distributions in a complicated way via tidal effects, stellar winds, the first mass transfer phase (characterised by mass transfer efficiency, the angular momentum loss during mass transfer) and the amount of mass that is lost due to supernova and the supernova kick (see Appendix \ref{sec:app_{accretor}andq}). However, for the most massive stars at low metallicities the picture simplifies considerably for the following reasons. Firstly, the line-driven stellar winds have a negligible effect at metallicities $Z \lesssim 0.001$, as we expect that such mass loss rates scale strongly with metallicity (Sect.\,\ref{subsection:swp}). The supernova models of \citet{Fryer_2012} suggest that at sufficiently large masses ($M_{ZAMS} \gtrsim 40 M_{\odot}$), black hole progenitors form via direct collapse, without any fallback and little kick (Equation \ref{eq:kick}). Furthermore, tides are of less importance in the context of gravitational wave progenitors, as orbits typically  circularise by tidal effects before the onset of mass transfer \citep[e.g.][]{Hurley2002}.
In conclusion, at low metallicities and at sufficiently high masses, the binary demographics  at the onset of the second mass transfer phase are primarily shaped by a single proces: the first mass transfer episode.

\begin{figure*}
  \includegraphics[width=0.5\textwidth]{q2q.pdf}\hfill
  \includegraphics[width=0.5\textwidth]{q2deltaa.pdf}
  \caption{Left panel: the relationship between the initial mass ratio  ($q_{ZAMS} = M_{ZAMS,2}/M_{ZAMS,1}$) and the mass ratio at the onset of the second mass transfer phase ($q_{MT,2} = M_{d}/M_{a}$). 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.73 K$. This coincides with the stage where donors develop deep convective envelopes according to \citet{IvanovaTaam2004}. For the majority of the binaries, this occurs with a core-helium burning donor ($M_{ZAMS,1} \gtrsim 90 M_{\odot}$ 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 $\zeta_{ad,rad}$. We have neglected the effects of supernova but not the stellar winds. The upper horizontal, dotted line shows $q_{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 $\zeta_{ad,rad} = 4$. The lower bottom line show the critical mass ratio when the donor has a deep convective envelope. See Table \ref{tab:zeta2q}}
\label{fig:qsandas}
\end{figure*}

\begin{figure}
  \includegraphics[width=0.55\textwidth]{qmax.pdf}
  \caption{The estimated maximum value of mass ratio at the onset of the second mass transfer phase ($q_{MT,2} = M_{d}/M_{a}$) that a system can reach with a given initial primary mass, mass transfer efficiency and mass transfer phase type. The maximum $q_{MT,2}$ is achieved by binaries that have initially near equal masses. How exactly we defined 'Case B' and '$logT_{eff}$' mass transfer episode types are explained in the caption of Figure \ref{fig:qsandas}.}
\label{fig:qmax}
\end{figure}

In the following, we discuss how the properties of binaries are affected by the first mass transfer phase with different assumptions for $\gamma$ and $\beta$ with an analytically-based approach. We stress that while in this subsection, we make a number of simplifying assumptions (e.g. neglecting the effects of supernovae), in the rest of the paper, the results are directly produced by the rapid population synthesis code SeBa, in which no such simplifications are applied. Our assumptions for this subsection include: 1) We restrict the discussion to sources with $M_{ZAMS,1} \gtrsim 45 \,M_{\odot}$ at $Z = 0.0007$ for reasons discussed above. However, we note that results presented here are still qualitatively true for all simulated binaries. 2) We ignore the effects of supernova on the orbit. 3) We ignore tidal evolution by assuming all binaries are circularised at birth.
4) For the examples of Case B mass transfer \ref{sec:masstransfers}, we assume that the donor fills its Roche-lobe when it reaches the midpoint of its Hertzsprung gap lifetime. 5) Furthermore, 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 $M_{ZAMS,1} \gtrsim 90\,M_{\odot}$ 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 transfer phases is more sensitively dependent the initial separation (see discussion \ref{sec:masstransfers}).

\subsubsection{Effect on the mass ratio}
\label{subsec:effectonthemassratio}
In the left panel of Figure \ref{fig:qsandas}, we show the relation between the initial mass ratio $q_{ZAMS}\equiv M_{ZAMS,2}/M_{ZAMS,1}$ and that at the onset of the second mass transfer phase $q_{MT,2}\equiv  M_{d}/M_{a}$. This figure shows the effect of different mass transfer efficiencies and different initial masses for Case B and Case C mass transfer types.

For Case B mass transfer, the relation between $q_{ZAMS}$ and $q_{MT,2}$ can be approximated as:
\begin{equation}
    q_{MT,2}(q_{ZAMS},M_{ZAMS,1}) \approx \frac{M_{ZAMS,1}\cdot q_{ZAMS} + \beta M_{envelope,1}}{M_{BH,1}},
\label{eq:qmt2}
\end{equation}
where $M_{envelope,1}$ is the mass of the hydrogen envelope of the primary star at the start of the first mass transfer phase, and  $M_{BH,1}$ is the mass of the black hole formed by the primary star. As  $M_{envelope,1}$ and $M_{BH,1}$  are not sensitively dependent on the initial separation for binaries evolving through Case B mass transfer phases (subsection \ref{sec:masstransfers}), we approximate $M_{envelope,1} = 7.03 + 0.478\cdot M_{ZAMS,1}$ and  $M_{BH} = -5.57 + 0.46\cdot M_{ZAMS,1}$ at $Z = 0.0007$ and $M_{ZAMS,1} \gtrsim 40\,M_{\odot}$.
The relation between $q_{MT,2}$ and $q_{ZAMS}$ is more complicated for Case C systems, mostly due to the metallicity independent LBV winds. For such binaries, we use relations presented in Appendix \ref{sec:app_{accretor}andq}.

Figure \ref{fig:qsandas} can be used to gauge the effect of the assumed accretion efficiency on the GW progenitor channels. The progenitors of each sub-channel introduced in section \ref{section:isolated_brief} have the following range of mass ratios at the onset of the second phase of mass transfer:
\begin{itemize}
    \item 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. $q_{MT,2}>q_{crit}$ with $\zeta_{ad} =\zeta_{HW}$

    \item rCEE channel: the second phase of mass transfer is Case C, unstable but the donor star has a radiative envelope, i.e.
    $q_{MT,2}>q_{crit}$ with $\zeta_{ad} =\zeta_{ad,rad}$.
    \item stable channel:
     The binary has a relatively short period (i.e. Case B first phase of mass transfer). This is because the orbit does not shrink as efficiently due to a stable phase of mass transfer as due to CEE.
     The binary has a relatively large $q_{MT,2}$. The orbit only then can shrink efficiently due to a stable mass transfer (see Equation \ref{eq:iso}). This means there is a minimum value of $q_{MT,2}$ associated with this channel (for a given $M_{ZAMS,1}$). As we will see in subsection \ref{subsec:impact_stable}, this value depends on $\gamma$ and $\beta$ and on the predicted outcome Case A mass transfer episodes in a complicated way. Therefore: $q_{MT,2,min}<q_{MT,2}<q_{crit}$ with $\zeta_{ad} =\zeta_{ad,rad}$.
    \end{itemize}
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.

Figure \ref{fig:qsandas} 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\lesssim q_{MT,2} \leq q_{crit}$. If $\beta = 0.3$, we expect the most massive sources of the stable channel to evolve from binaries with $q_{ZAMS} \gtrsim 0.8$. On the other hand, the formation via the rCEE channel, which requires $q_{MT,2} > q_{crit}$, is not possible at all with mass transfer efficiencies as low as $\beta = 0.3$. Even if $\beta = 0.7$, the formation becomes only possible if $q_{ZAMS} \gtrsim 0.85$ and  $M_{ZAMS,1} \gtrsim 60\,M_{\odot}$.

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 $q_{MT,2}$ of binaries with a Case C first phase of mass transfer are narrowed down by the LBV winds.

\subsubsection{The effect on the maximum mass ratio}

The maximum $q_{MT,2}$ that can be reached by binaries at Z = 0.0007 is shown in  Figure \ref{fig:qmax}. The maximum $q_{MT,2}$ is given as a function of $M_{ZAMS,1}$ for a given mass transfer efficiency and mass transfer type. The maximum occurs for near equal mass binaries for Case B binaries, and therefore, can be approximated by calculating $q_{MT,2}$ with $q_{ZAMS} = 1$.
We expect the results to depend only weakly on the metallicity as long as $Z\lesssim 0.002$, because the $M_{ZAMS} - M_{BH}$ relation does not significantly change at this metallicity range (see e.g. Figure \ref{fig:single_bh_spectrum}).

We can use Figure \ref{fig:qmax} to asses the importance of $\zeta_{ad,rad}$. A higher $\zeta_{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 \ref{eq:iso}).
Hence, we may expect that the merger rate of the stable channel could increase significantly with increasing $\zeta_{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 $\beta$ and the $M_{ZAMS} - M_{BH}$ relation.
The maximum mass ratio of binaries with Case B mass transfer phase is below 3.2, if $M_{ZAMS,1} \gtrsim 60 \, M_{\odot}$ and $\beta = 0.3$. Therefore, in our low mass transfer efficiency models, increasing $\zeta_{ad,rad}$ from 4 to 7.5  only affects the evolution of binaries with $M_{ZAMS,1} \lesssim 60\,M_{\odot}$ .

Figure \ref{fig:qmax} also tells us about the effects 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 full mass range. This can be seen from Figure \ref{fig:qmax} as the curve for Case C mass transfer for $\beta \leq 0.3$ and $\zeta_{ad,rad} = 4$ is below the critical mass ratio for $M_{ZAMS,1} \gtrsim 40\, M_{\odot}$.  If $\zeta_{ad,rad}$ is increased up to $7.5$, then the formation of such systems become impossible even with larger mass transfer efficiencies.

\subsubsection{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 ($a_{MT,2}$) can be approximated as:
\begin{equation}
    \frac{a_{MT,2}}{a_{init}} = f_{wind,post-MT,1}\cdot f_{MT,1}\cdot f_{wind,pre-MT,1},
\end{equation}
if the effects of supernova are negligible.
Here, $f_{wind,post-MT,1}$ and $f_{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 are determined according to Equation \ref{eq:Jeans}. These changes are negligible for binaries evolving through Case B mass transfer at metallicities $Z \lesssim 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 $f_{MT,1}$ is the change in the orbit during the first mass transfer phase. In principle, these are due to the combined 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 \ref{eq:stable_mt_eq}.
In all cases, $f_{MT,1}$ is the dominant term, i.e. $|f_{MT,1}| \gg max(|f_{wind,post-MT,1}|, |f_{wind,pre-MT,1}|)$.

In the right panel of Figure \ref{fig:qsandas}, we show $log(a_{MT,2}/a_{initial})$ as a function of $q_{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 $log T_{eff} =10^{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 \ref{eq:stable_mt_eq}, which shows that as long as $M_{d} > M_{a}$ 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 $\gamma$ and accretion efficiency $\beta$. As we increase $\beta$, 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 $\gamma = 2.5$ and $\beta = 0.3$ shows only an insignificant change in the orbit. All other models display orbital widening, up to a factor of a few.
We will see later this has important consequences on the predicted merger rate of the stable channel.

\subsubsection{Case A mass transfer episode}
\label{subsec:CaseA}
 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 \ref{sub:typesofmt}). This means that it is challenging to show simple relations for these kind of systems as we did in Figure \ref{fig:qsandas} and \ref{fig:qmax} for Case B and Case C.  Instead, in Figure \ref{tab:CaseA} we show the evolution of a few selected binaries until the second mass transfer phase as simulated by SeBa. These results show two important characteristics of systems evolving through a first phase of Case A mass transfer. Firstly, they develop $q_{MT,2}$ that are significantly larger than that of their Case B counterparts. These mass ratios are often above $q_{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 Figure \ref{tab:CaseA} would indeed merge with $\zeta_{ad,rad} = 4$, while one of them would merge with $\zeta_{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 \citet{Hurley_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.

\begin{figure*}
  \includegraphics[scale = 0.7]{intial_plot_final.pdf}
  \caption{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; $\beta$ = 0.3, $\zeta_{ad,rad}$ = 7.5. Furthermore, we have used the assumption of \citet{IvanovaTaam2004} regarding the development of deep convective envelopes. In the left panel we show systems with $\gamma = 2.5$, in the right with $\gamma = 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 fills 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.}
  \label{fig:initial}
\end{figure*}

\subsection{The fate of the binary as a function of its initial parameters}
\label{subsec:init}
Figure \ref{fig:initial} show our binary populations as simulated by the population synthesis approach with SeBa for a range of systems with primary masses of $M_{ZAMS,1} = 50\,M_{\odot}$ at metallicity of $Z = 0.0007$. The figure shows that the combination of initial mass ratio and semilatus rectum is a good proxy for the outcome of the evolution for each simulated binary (indicated by the color). For a fraction of systems the outcome is a GW event (shown in reds). The different panels highlight how the distributions change for different $\gamma$.

Most of the Case A binaries  merge. This is either a result of the first or the second mass transfer phase. The latter occurs because at this stage the donor is typically a Hertzsprung gap giant and the mass ratios are above $q_{crit}$ (see discussion in \ref{subsec:CaseA}).

In Figure \ref{fig:initial}, we highlight those binaries that initiate unstable mass transfer with core helium burning donors (shown in yellow).
This is helpful in order to understand how donors with deep convective envelopes interact with their companion under different assumptions of binary physics.
While in model variations with $\gamma = 2.5$, the first phase of Case Cc mass transfer is stable, if $q_{ZAMS} \gtrsim 0.8$, this occurs for a much larger range, when $\gamma = 1$, i.e. $q_{ZAMS} \gtrsim 0.3$. Therefore, we predict that the fraction of binaries with giant donors with deep convective envelopes that undergo unstable mass transfer can be vastly different based on what is assumed for mass transfer efficiency $\beta$ and angular momentum loss $\gamma$.

Figure \ref{fig:initial} shows very clearly that the CEE channel and the stable channel originate from different initial separations. The initial orbital separations of sources of the CEE channel are around $ 1000\,R_{\odot}$, while for the stable channel, they are around $ 50 - 100\,R_{\odot}$.
Although, in this figure, we only show binaries in a restricted primary mass range, we find that the initial separations of gravitational waves progenitors are fairly distinct for these channels regardless of the initial primary mass and assumptions of binary physics.

The maximum of the initial separation of sources of the stable channel is determined by the factor by which an orbit shrinks during a stable mass transfer phase with a black hole accretor. As we will see in subsection \ref{subsec:impact_stable}, the orbit can shrink a few $\sim 100\, R_{\odot}$ at most for the majority of the systems in such a case. This implies that the initial orbital separation must also be in the same order, since the first mass transfer and the natal kick does not typically change the orbit by orders of magnitude.
The minimum value of the initial separation of this channel
coincides with the separation at which the initial primary star would fill its Roche-lobe just when it evolved off main sequence. As we have already seen, if the first mass transfer phase occurs in an earlier stage than late Case A, then the mass ratios at the onset of the second mass transfer are typically above the critical mass ratio.

The range of initial orbital separation distribution of the CEE sources is primarily determined by three factors. Firstly, the adopted  efficiency of the common envelope phase and the assumed binding energies of the donor star. Our $\alpha_{CE}\lambda = 0.005$, typically decreases the orbital separation by a few $\sim 1000\,R_{\odot}$.
Secondly, the assumption that HG donors cannot survive CEE. This puts a constraint on the minimum initial separation of these systems. Thirdly, by the critical mass ratio for a first phase of Case Cc mass transfer (i.e. with a non-compact accretor). This latter introduces a $q_{ZAMS}$ dependent limit on the maximum value of the orbital separation. The role of this critical mass ratio can be understood from Figure \ref{fig:initial}. In the model variations with $\gamma = 1$, binaries initiate a stable, first phase of Case Cc mass transfer, if $q_{ZAMS}\gtrsim 0.3$. This range is $q_{ZAMS}\gtrsim 0.8$ with $\gamma =  2.5$. This implies that for a range of $0.3 \lesssim q_{ZAMS} \lesssim 0.8$, in model variations with  $\gamma = 1$, binaries initiate a stable episode of Case Cc mass transfer, as opposed to unstable episodes in the $\gamma = 2.5$ models. As such binaries of the former model variation retain relatively large orbital separations, many of them have a second phase of unstable, Case Cc mass transfer and therefore can potentially be sources of the cCEE channel. This implies that in the $\gamma = 1$ model, the maximum initial separation of the CEE sources is larger than in the $\gamma = 2.5$ model.

\begin{figure*}
  \hspace*{1cm}
  \includegraphics[width=0.9\textwidth, height=0.8\textheight]{mps_mod2.pdf}
  \caption{We show the density contours of BH-star binaries, at the onset of the second phase of mass transfer in the mass ratio - orbital separation space for four different model variations. Contours with blue colour denote binaries for which the first phase of mass transfer is either Case B or Case C, while contours with green colour show systems for which this is Case A. In this Figure, we limit our population to systems for which the initial primary star had a mass of $M_{ZAMS,1} = 90\,M_{\odot}$. The metallicity is $Z = 0.0007$ and the convective envelope prescription is of \citet{IvanovaTaam2004}. The $a_{GW,final}$ line shows the separation that binaries have by the time they form BH-BH binaries, if they were to merge due to gravitational waves just within Hubble time. For this, we estimate the remnant mass of the donor assuming that the second phase of mass transfer takes place at the midpoint of the hydrogen shell burning phase. The  $a_{GW,MT,2}$ line indicates the separation that the same binaries have at the onset of the second phase of mass transfer. We assume stable mass transfer (i.e. Equation \ref{eq:iso}) and we neglect the effect of supernova kick on the orbit. The $a_{TAMS,2}$ line indicates the orbital separation at which the second phase of mass transfer occurs when the donor just evolved off the main sequence. The $a_{TAHG,2}$ is the orbital separation at which the second mass transfer phase occurs with a donor, which just evolved off the Hertzsprung gap phase, while $a_{conv,2}$ is the orbital separation at which the donor fills its Roche lobe just before it develops a deep convective envelope. Furthermore we also show 2D histograms of the GW sources; CEE channel (dark red), stable channel with Case B second mass transfer phase (orange), stable channel with Case A second mass transfer phase (blue). These histograms have been normalised separately.}
\label{fig:mps}
\end{figure*}

\subsection{Binaries at the onset of second phase of mass transfer}
In Figure \ref{fig:mps}, 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 $Z = 0.0007$.  In Figure \ref{fig:mps}, we limit the initial primary mass to $M_{ZAMS,1} = 90\,M_{\odot}$. We show four model variations with different assumptions on $\gamma$ and $\beta$. 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 $M_{BH,1} = 35.44\, M_{\odot}$.  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 channel, stable channel with Case B binaries (`stable Case B', hereafter) and stable channel with Case A binaries (`stable Case A').
In the following two subsections, we discuss in detail how different assumptions of $\gamma$, $\beta$ and $\zeta_{ad,rad}$ shape the properties of the GW sources. In particular, relating $q_{MT,2}$ to the observable mass ratios of the binary black holes (i.e. $q_{final,obs} =min(M_{BH,1},M_{BH,2})/max(M_{BH,1},M_{BH,2})$), we can understand how the uncertainties in the aforementioned parameters affect the mass ratio distribution of each channel.

\subsection{Impact of uncertainties on the stable channel}
\label{subsec:impact_stable}

In Figure \ref{fig:mps}, we can see that Case A binaries indeed attain significantly larger values of $q_{MT,2}$ than Case B binaries. In fact, most Case A systems initiate an unstable second phase of mass transfer with $\zeta_{ad,rad} = 4$.
Binaries which evolve through Case A mass transfer phases do not typically form gravitational wave sources, except when $\gamma = 1$ and $\zeta_{ad,rad} = 7.5$. These stable Case A systems have typically large values of $q_{MT,2}$ and hence relatively unequal final mass ratios, ranging between $q_{final,obs}\approx 0.5 - 0.8$.

The majority of the progenitors of stable channel originate from Case B systems.
In order to understand the range of $q_{final,obs}$ of the stable Case B channel we need to know the range of $q_{MT,2}$ of their progenitors
The largest value of $q_{MT,2}$ for the stable Case B sub-population for a given $M_{ZAMS,1}$ is determined by $min(q_{crit}, q_{MT,2,max})$. Here $q_{MT,2,max}$ is the maximum value of the mass ratio at the onset of the second phase of mass transfer that can be achieved for a given $\beta$ and $M_{ZAMS,1}$ as shown in Figure \ref{fig:qmax}.
The minimum value of $q_{MT,2}$ of the stable Case B systems depends on the orbital separation distribution at the onset of the second phase of mass transfer and therefore it is primarily set by $\gamma$ and $\beta$.
In order to gain a better understanding of this, let us consider the part of the population shown in Figure \ref{fig:mps} that forms stable Case B binaries. For such systems, the following conditions hold:
\begin{itemize}
    \item \textit{The system forms a black hole binary with sufficiently small orbital separation}: For the systems shown in Figure \ref{fig:mps}, it is possible to determine the orbital separation that a binary has at the onset of the second phase of mass transfer, if it was to form a black hole binary that merges due to GW just within the Hubble time (shown by $a_{GW,MT,2}$ in Figure \ref{fig:mps}). In order to do this, first we estimate the mass of the black hole that the initial secondary star forms. Secondly, we determine the corresponding orbital separation that this binary black hole has if it was to merge just within the Hubble time (shown by $a_{GW,final}$ in Figure \ref{fig:mps}). Finally, we determine $a_{GW,MT,2}$ by assuming that the second phase of mass transfer is stable, Case B and occurs when the donor reaches half of its Hertzsprung gap lifetime.  From this it follows that stable Case B binaries have an orbital separation at the onset of the second mass transfer phase ($a_{MT,2}$) that is smaller than or equal to $a_{GW,MT,2}$.
    \item \textit{The second mass transfer is not Case A}: The second phase of mass transfer is Case A, if $a_{MT,2}<a_{TAMS,2}$, where $a_{TAMS,2}$ is the orbital separation at which the donor star (initial secondary) would fill its Roche lobe at the terminal age MS (shown by red curve in  Figure \ref{fig:mps}). Systems which fulfil the previous condition (i.e. $a_{MT,2} < a_{GW,MT,2}$) but for which the second mass transfer phase is Case A do not necessarily end up with $a_{BHBH} <a_{GWmerger}$. In fact, a binary going through a Case A mass transfer phase often ends up with a wider orbit than the same binary evolving through Case B mass transfer, since relatively more mass is transferred onto the accretor. We find that systems for which the second mass transfer phase is Case A do not typically form GW sources (but see \citealt{marchant2021role}).
    \item \textit{The first mass transfer is not Case A}:
   The first mass transfer phase cannot be Case A, therefore $a_{ZAMS}>a_{TAMS,1}$, where $a_{TAMS,1}$ is the orbital separation at the onset of the first phase of mass transfer, where the initial primary fills its Roche-lobe just at the end of its main sequence phase.
\end{itemize}
As a result of the first and the second conditions, $a_{TAMS,2} \lesssim a_{MT,2} \lesssim a_{GW,MT,2}$ has to hold for the progenitors of the stable channel. When $\gamma = 2.5$, there is a non-negligible fraction of Case B binaries for which the orbit shrinks so much after the first phase of mass transfer that they have a Case A second phase of mass transfer.
In this case, the minimum value of $q_{MT,2}$ is where $a_{TAMS,2} = a_{GW,MT,2}$. This is also shown in Figure \ref{fig:mps}, where for the $\gamma = 2.5$ models, this value is around $q_{MT,2} \approx 2.2$.
When $\gamma = 1$, the orbit widens so much that none of the Case B binaries go through Case A second mass transfer phase. In such cases, the minimum value of $q_{MT,2}$ follows from the first and the third conditions. Since the first phase of mass transfer cannot be Case A, we should expect that the Case B binaries with the shortest orbit (well represented by the lower edge of the blue contour in Figure \ref{fig:mps}) had an orbital separation of $a_{TAMS,1}$ at the onset of the first phase of mass transfer.
Then, the minimum value of $q_{MT,2}$ for the stable Case B binaries occurs at a $q_{MT,2}$ where the edge of the contour and $a_{GW,MT,2}$ cross. For the model variations with $\gamma = 1$, this occurs roughly at $q_{MT,2} \sim 3.1$.

The upper x-axis of Figure \ref{fig:mps} shows the range of final mass ratios of the merging binary black holes for Case B binaries. For example, this is $q_{final,obs}\approx 0.8 - 1$. for the model variation with $\gamma = 2.5$ and $\beta = 0.3$  With our high mass transfer efficiency model ($\beta = 0.7$), this range gets somewhat broader, i.e. $q_{final,obs}\approx 0.7 - 1$, however only if $\zeta_{ad,rad} = 7.5$. This is because the maximum value of $q_{MT,2}$ of the stable Case B population is larger than 3.2 with $\beta = 0.7$. Since in the model variations with $\gamma = 1$, the minimum value of $q_{MT,2}$ of the stable Case B population is around 3.1, we do not expect GW sources from the stable channel with $q_{final,obs} \gtrsim 0.8$.

\subsection{Impact of uncertainties on the CEE channel}
\label{subsec:impact_cee}
We now turn our attention to CEE channel. We assume envelope ejection during common envelope evolution is only possible with core-helium burning donors i.e. with a second phase of Case C mass transfer. Therefore, in Figure \ref{fig:mps}, any GW progenitors of this channel have an orbital separation at the onset of the second phase of mass transfer for which  $a_{MT,2}>a_{TAHG,2}$ holds. Here $a_{TAHG,2}$ represents the minimum orbital separation where the initial secondary fills it Roche-lobe at the end of its Hertzsprung gap phase. Furthermore, to delineate the rCEE and cCEE channel,  we show $a_{conv,2}$  which is the shortest orbital separation of those binaries in which the donor has a convective envelope (i.e. the second phase of mass transfer is Case Cc).

The following conditions hold for the sources of the rCEE channel: $a_{TAHG,2}\leq a_{MT,2} < a_{conv,2}$
and $q_{MT,2}>q_{crit}$. The same for the binaries of the cCEE channel are: $a_{MT,2}> max( a_{TAHG,2}, a_{conv,2})$ and $q_{MT,2}>q_{crit}$.
Note that in the two subchannels $q_{crit}$ has very different values. In the prior case $\zeta_{ad}=\zeta_{ad,rad}$, such that $q_{crit}\approx 3.2$ in the default model, while in the latter case $\zeta_{ad}=\zeta_{HW}$ such that $q_{crit}\approx 1.4-1.7$ (see Table \ref{tab:zeta2q}).

The effective temperature at which a giant develops a deep convective envelope ($T_{eff,boundary}$) has a significant effect on the CEE channel. In Figure \ref{fig:mps}, we only show model variations, where we assume this occurs at $\log(T_{eff}) = 3.73\,K$ (e.g. following assumptions of \citealt{Ivanova_2004}, \citealt{Belczynski2008}).
The models of \cite{Klencki_2021} predict this would occur at considerably cooler temperatures and that the most massive stars ($M\gtrsim 40\,M_{\odot}$) do not even develop deep convective envelopes at low metallicities ($Z \lesssim 0.0007$) at all. The latter naturally implies a lack of massive cCEE sources for these model variations.

Figure \ref{fig:mps} shows that there merger rate of the cCEE channel is considerably larger in model variations with $\gamma = 1$ compared to the ones with $\gamma = 2.5$. In the former case, the binaries have wider orbits and therefore more of them initiate a Case Cc second phase of mass transfer (i.e. $a_{MT,2} \geq a_{conv,2}$). There are two reasons for this.

Firstly, in the $\gamma = 1$ models, more binaries initiate a stable Case Cc first phase of mass transfer compared to the $\gamma = 2.5$ models. This is because the critical mass ratio for a Case Cc mass transfer phase with a MS accretor is much larger with $\gamma = 1$ (see Figure \ref{fig:initial}).
As a result, these binaries do not experience dramatic decrease in their orbits due to CEE and they retain relatively large $a_{MT,2}$.

Secondly, the Case C binaries widen more in model variations with $\gamma = 1$ (see right panel of Figure \ref{fig:qsandas}). Therefore in this case, there are more binaries for which the first phase of mass transfer is stable and Case Cr and the second is unstable and Case Cc. This scenario does not occur, if the orbit shrinks after the first phase of mass transfer because at shorter orbits the initial secondary star fills its Roche-lobe at an earlier evolutionary stage, at which it still has a radiative envelope.

Our assumptions on $\gamma$ and $\beta$ also has an effect on the range of $q_{final,obs}$  of the GW sources of the cCEE channel. As Figure \ref{fig:mps} shows, with $\gamma = 1$, the mass ratios of the merging binary black holes of the cCEE channel can be as low as $q_{final,obs}\approx0.55$. In this case the minimum value of $q_{final,obs}$ is determined by the critical mass ratio of a Case Cc mass transfer phase with a black hole accretor (see Table \ref{tab:zeta2q} or the left vertical dashed blue line in Figure \ref{fig:mps}). On the other hand if $\gamma = 2.5$, there are no binaries with $a_{MT,2} \geq a_{conv,2}$ for $q_{MT,2} \gtrsim 2.0$. This leads to a narrower mass ratio distribution, i.e. $q_{final,obs}\approx 0.7-1.0$.

In the particular examples shown in Figure \ref{fig:mps}, there are no systems that evolve through the rCEE channel. If such systems existed, they would neet to have at least $q_{MT,2} = q_{crit}$. This implies that even if $\zeta_{ad,rad} = 4$, the donor has to be at least $M_{2} = 85\,M_{\odot}$ at the onset of the second phase of mass transfer. At this metallicity, such a massive star ignites the helium in its core at a stage where it already has a convective envelope (i.e $a_{TAHG,2} > a_{conv,2}$, where $q_{MT,2}>q_{crit}$).
The parameter space for Case Cr mass transfer increases with decreasing metallicity or decreasing initial mass, however, only a very few binaries evolve to $q_{MT,2} \gtrsim q_{crit}$ (see Figure \ref{fig:qmax}).
In fact, with the assumption of $\log(T_{eff, boundary}) = 3.73\,K$, we do not find the rCEE channel particularly efficient at any mass range.
With the assumption of $T_{eff,boundary}$ of \citet{Klencki_2021}, we find that the rCEE channel is only efficient for systems with an initial primary mass of $M_{ZAMS,1} \lesssim 40\,M_{\odot}$ and only when  $\zeta_{ad,rad} = 4$.

\begin{figure*}
  \includegraphics[scale = 0.6]{hist2ds_combined_new.pdf}
  \caption{2D Histograms that show the masses of the merging binary black holes for a few selected models at Z = 0.0007. Each histogram have been normalised to the merger efficiency of their own model variation. The dashed lines show constant mass ratios. }
  \label{fig:43_scatter}
\end{figure*}

\begin{figure*}
  \includegraphics[width=1.0\textwidth]{Updated_primaryhist_mod_new.pdf}
  \caption{The mass distribution of the primary of the merging binary black hole in the local universe (z = 0) for all of our models. The distribution is shown by stacked histograms, where each colour indicates a different formation channel.}
  \label{fig:43_hist}
\end{figure*}

\subsection{The properties of the merging binary black hole population }
\label{sub:popprops}

In Figure \ref{fig:43_scatter}, we show 2D histograms of the masses of the merging binary black holes for a few selected model variations at Z = 0.0007 from the CEE and the stable channel. These figures allow us to see the predicted mass ratios of GW sources for the full mass range.
For $M_{BH} \gtrsim 20\,M_{\odot}$, both channels can yield a population of merging binary black holes with a relatively narrow range in mass ratio ($0.8 \lesssim q_{final,obs} \lesssim 1.0$) or a population with a moderately wide ($0.5 \lesssim q_{final,obs} \lesssim 1.0$) mass ratio distribution, depending on the assumptions on the uncertain binary physics. On the other hand, much wider mass ratio distributions are yielded by less massive systems.  We note that for systems with $M_{BH} \gtrsim 20\,M_{\odot}$, the range of mass ratios are in agreement with the discussion in subsection \ref{subsec:impact_stable} and \ref{subsec:impact_cee}.

So far we have only discussed GW formation at a metallicity of Z = 0.0007.
In the appendix in Figure \ref{fig:merger_efficiency} we show the merger efficiency (defined in Equation \ref{eq:mereff}) as a function of metallicity for each model variation.
The rate of mergers of the stable channel monotonically increases with decreasing metallicity regardless of which model variant we consider. The stars at the end of main sequence are more compact at lower metallicities. As previously discussed, a smaller expansion in radius during the main sequence implies that the parameter space for the stable channel increases, as the initial orbits of such sources are bounded by a separation at which the donor fills its Roche-lobe when it just evolved off main sequence. Therefore, the number of mergers of this channel increases with decreasing metallicity.

On the other hand, the merger efficiency of the CEE does not exhibit a monotonic behaviour as a function of metallicity. In some cases, we see a peak around $Z \sim 0.003$. This channel is dominated by the sources of cCEE. The trend shown in Figure \ref{fig:merger_efficiency} is determined by two factors.  Firstly, the onset of core helium burning phase occurs at higher effective temperatures with decreasing metallicity. This implies that at sufficiently high metallicities, where the core helium ignition occurs at an evolutionary stage at which a deep convective envelope is already developed, the merger efficiency of the CEE channel is expected to increase with decreasing metallicities. Secondly, the stellar tracks we use, predict that the minimum effective temperature that the most massive stars eventually reach typically increases with decreasing metallicities. This suggests that at sufficiently low metallicities, at which the core helium burning starts at an evolutionary stage at which the donor has still mostly radiative envelope, the merger rate efficiency shows a decreasing trend with increasing metallicity.

In Figure \ref{fig:43_hist}, we show the primary black hole mass distribution in the local universe (z = 0) for each channel considered here (stable, rCEE, cCEE). Therefore, they are affected by the combined effects of time delay and the cosmic star formation history. In the legend, we also included the estimated cosmic merger rate densities for each channel. These values have been computed as summarised in subsection \ref{subs:rates}. We can see that these values are in a broad agreement with the currently inferred rate from LIGO and Virgo observations (\citealt{Abbott_o3a_2021}).
Figure \ref{fig:43_hist} shows that the relative importance of each sub- channel sensitively depends on the adopted values of $\gamma$, $\beta$, $\zeta_{ad,rad}$, and $T_{eff,boundary}$.

For most models, the cosmic merger rate density is dominated by sources evolving according to the the stable channel. This is true even for those model variations for which the merger efficiency of the CEE channel is larger than that of the stable channel at most metallicities (see e.g. $\zeta_{ad,rad} = 4$ and $\gamma = 2.5$ models in Figure \ref{fig:merger_efficiency}).  This is due to the longer time delays of the stable channel and the fact that the cosmic star formation rate is predicted to increase up to $z \sim 2-3$. This effect was already shown by \citet{Neijssel_2019}.

The contribution of the stable channel becomes negligible with $\gamma = 1$ in model variations with $\zeta_{ad,rad} = 4$ (compare row 1 and 3 of Figure \ref{fig:43_hist}).
However, if $\zeta_{ad,rad}$ is increased up to 7.5, then the stable channel becomes the dominant channel for all model variations, even with $\gamma = 1$. In particular, if $\gamma = 1$ and $\zeta_{ad,rad} = 7.5$, then 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 \lesssim q_{MT,2} \lesssim 5.5$, to be exact). Such sources can form, if the first mass transfer phase is Case B and the mass transfer efficiency is relatively large (e.g. $\beta = 0.7$) or if the first mass transfer is Case A, which typically leads to large values of $q_{MT,2}$. Consequently, the final mass ratios of these merging binary black holes are smaller, i.e. $q_{final} \leq 0.8$ (see Figure \ref{fig:43_scatter}). The lack of GW sources with mass ratios $q_{final} \geq 0.8$ is in broad agreement with \citet{olejak2021impact}, who implemented the mass transfer stability criteria of \citet{Pavlovskii_2016} and assumed $\gamma = 1$.

The effect of $\gamma$ on the merger rate of the cCEE channel is in agreement with the discussion in subsection \ref{subsec:impact_cee}. Decreasing $\gamma$ from 2.5 to 1 increases the merger rate of this channel by a factor of 8 in our low mass transfer efficiency models. With $\beta = 0.3$, this difference is only about a factor of two.

In general, increasing $\zeta_{ad,rad}$ significantly increases the predicted merger rate density of the stable channel (compare row 1 and 2, or row 3 and 4 of Figure \ref{fig:43_hist}). As previously explained, with low mass transfer efficiencies (e.g. $\beta = 0.3$), only the merger rates of the lower mass systems ($M_{BH} \lesssim 20 \, M_{\odot}$) are affected. At the same time the merger rate of the rCEE channel decreases, while the cCEE channel remains practically unaffected.

Lastly, as expected, the merger rate density of the cCEE channel decreases (and that of the rCEE channel increases) when adopting the convective envelope prescription from \citet{Klencki_2021} instead of that of \citet{IvanovaTaam2004} (compare column 1 and 3, or column 2 and 4 of Figure \ref{fig:43_hist}). It also limits the maximum mass of coalescing binary black holes of the cCEE channel to $M_{1,BH} \sim 25\,M_{\odot}$.

\section{Effects of stellar winds}
\label{sub:stellarwinds}
\subsection{The effects of stellar winds on single stellar evolution}
\label{sub:single}

In this section, we describe the effects of uncertainties in stellar winds on the evolution of single, massive stars, with particular attention paid to the final black hole masses. This is achieved by running stellar evolution simulations with SeBa with three different stellar wind models. These three models are the following: 'Model I' with our standard stellar wind model ($f_{wind} = f_{wind,WR} = 1$), 'Model II' with the optically thin line driven winds scaled down by a factor 3 ($f_{wind} = 1/3, f_{wind,WR} = 1$) and 'Model III' where besides the optically thin line driven winds, Wolf-Rayet-like winds are also scaled down ($f_{wind} = f_{wind,WR} = 1/3$). For the exact description of our applied stellar winds prescriptions and scaling factors, see Table \ref{tab:mdotapplied}.
Model II is meant to capture the effect of the findings of studies such as \citet{Sundqvist_2019}, \citet{bjrklund2020new}, \citet{Krti_ka_2017}, which claim that the prescription of \citet{Vink_2001} systematically overpredicts the mass loss by a factor $\sim 2-3$. Modell III attempts to investigate the general uncertainties in the mass loss rates of stripped helium stars (see e.g. \citealt{Sander_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.

In order to understand the importance of different stellar wind prescriptions, we show where different stellar wind mechanisms dominate in the the Hertzsprung-Russel diagram in Figure \ref{fig:Wind1}, along with the predicted mass loss rates throughout the evolution of a few selected massive stars.

\begin{figure*}
  \includegraphics[width=0.45\textwidth]{0.02_HRDWinds.pdf} \hfill
  \includegraphics[width=0.45\textwidth]{0.005_HRDWinds.pdf}
  \caption{Hertzsprung-Russel diagrams for massive stars with the mass-loss rate colour-coded in the tracks. The stars have $M_{ZAMS}$ = 20, 30, 50, 60, 80, 100 $M_{\odot}$. 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.}
  \label{fig:Wind1}
\end{figure*}

\begin{figure*}
  \includegraphics[width=0.45\textwidth]{0.02_MassLostFraction.pdf} \hfill
  \includegraphics[width=0.45\textwidth]{0.005_MassLostFraction.pdf}
  \caption{Fraction of mass lost due to different mass loss mechanisms during the entire life of a non-interacting star at Z = 0.02 and Z = 0.005 as a function of zero-age main sequence mass.}
  \label{fig:Wind2}
\end{figure*}

In Figure \ref{fig:Wind2}, we show the fraction of mass lost due to the different wind mechanisms for single stars at metallicities $Z = 0.02$ and $Z = 0.005$.
Up to $M_{ZAMS} = 20\, M_{\odot}$, the mass of the remnant is almost completely determined by the mass lost due to supernova and dust-driven winds.
Above $M_{ZAMS} = 40-50\, M_{\odot}$, only a negligible amount of mass is lost due to SN as our model \citep[delayed model from][]{Fryer_2012} predicts full fallback for such systems. These models also assume that 10$\%$ of the proto-SN mass is lost due to neutrinos after core-collapse, regardless of the amount of the fallback.

The mass loss due to line-driven winds during the main-sequence phase becomes increasingly significant for stars with initial masses above $M_{ZAMS} = 40-50\, M_{\odot}$.
At this mass range, the stars eventually cross the Humphreys-Davidson limit, which leads to severe mass losses during the LBV phase. The masses of black holes which form via direct collapse (i.e. from progenitors with $M_{ZAMS}\gtrsim 40\, M_{\odot}$) end up losing 20-40$\%$ of their initial mass via line-driven winds during their main sequence phase. This already indicates that lowering the optically thin line-driven winds by a factor of three can have significant effect on the mass spectrum of the most massive black holes.

\begin{figure*}
  \includegraphics[width=0.45\textwidth]{Standard_Model_Mzams_Mremnant.pdf} \hfill
  \includegraphics[width=0.45\textwidth]{ScaledDownModel_Mzams_Mremnant.pdf}
  \caption{Zero age main sequence mass and remnant mass functions for different metallicities for single stellar evolution. We assume that black hole formation occurs above $M_{ZAMS} = 20 M_{\odot}$.}
  \label{fig:single_bh_spectrum}
\end{figure*}

In Figure \ref{fig:single_bh_spectrum}, we show the mass spectrum of black holes from single stellar evolution at a wide range of metallicities for Model I  and Model II.
Model II produces moderately more massive black holes in the mass range $ 40\, M_{\odot} \lesssim M_{ZAMS} \lesssim 70\, M_{\odot}$ compared to the standard wind model at metallicities $Z \gtrsim 0.005$. However, above $M_{ZAMS} \gtrsim 70\, M_{\odot}$, the differences become increasingly more appreciable. In particular, the most massive star we simulated ($M_{ZAMS} = 100 \,M_{\odot}$) using Model I yields  $M_{BH} \approx 14 M_{\odot}$ at solar metallicity. This increases up to $M_{BH} \approx 24 M_{\odot}$ for the Model II.  At $Z = 0.01$ black holes even as massive as $M_{BH}\gtrsim 30\,M_{\odot}$ are formed in Model II.
Since the initial masses of massive stars follow a distribution of $N \sim M_{ZAMS}^{-2.3}$ (\citealt{Kroupa}), 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 $Z\approx0.002$, the differences between the two stellar wind models become negligible.

\begin{figure}
  \includegraphics[width=0.5\textwidth]{section_32_svsb.pdf}
  \caption{The importance of stellar winds  with and without binary interactions are compared by showing the
  evolution of a star with an initial mass of $M_{ZAMS, 1}$ = 100$M_{\odot}$ at metallicity $Z = 0.01$ as a single star, and as a binary with a companion of $M_{ZAMS, 2} = 90 M_{\odot}$ and with an initial orbital separation of $a_{ZAMS} = 400 R_{\odot}$.
  We show the binaries for three different stellar wind models. Labels 'Binary', 'Binary I' and 'Binary II' correspond to Model I ($f_{wind} = 1$),  Model II  ($f_{wind} = 1/3$) and Model III $f_{wind} = f_{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: hydrogen-shell burning phase, hexagon: core helium burning phase, diamond: helium star, cross: helium giant. For clarity, we plotted the curves starting from the hydrogen-shell burning phase.}
  \label{fig:42_single_vs_binary}
\end{figure}

\begin{figure*}
\includegraphics[width=0.9\textwidth]{wind_primary_masses.pdf} \hfill
  \caption{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 ($\zeta_{ad,rad}$ and angular momentum loss during mass transfer with non-compact accretor ($\gamma$) at metallicities Z = 0.02, 0.01 and 0.005. This histograms are normalised to the merger efficiency.}
  \label{fig:wind_primary_first}
\end{figure*}

\subsection{The effects of stellar winds on binary evolution}
\label{subsec:interacting}

In general, non-interactive binaries form more massive black holes than their interacting counterparts.
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, it matters whether the system undergoes Case B or Case C mass transfer.  In particular, at metallicities $Z \gtrsim 0.005$, the radial expansion after the onset of core helium burning of stars with $M_{ZAMS} \gtrsim 40 M_{\odot}$  is negligible  (Figure \ref{fig:Wind1}). Practically, all  interacting massive binaries at these metallicities initiate Case B mass transfer episodes. Therefore 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 binary spends most, if not all of its 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, since this impacts its hydrogen envelope, this does not directly affect the mass of the helium core.

The difference between interactive and non-interactive binary evolution is illustrated in Figure \ref{fig:42_single_vs_binary}, where we compare the evolution of a single star with an initial mass $M_{ZAMS} = 100 \,M_{\odot}$, 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 ($M_{BH} = 23.3 \,M_{\odot}$) than the black hole formed through binary evolution ($M_{BH} = 16.2 \,M_{\odot}$).
In the lower panel of Figure \ref{fig:42_single_vs_binary}, 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 $M_{core} \sim 31.2 \,M_{\odot}$ to $M_{core} \sim 31.7 \,M_{\odot}$, while it increases up to $M_{core} \sim 35.2 \,M_{\odot}$ 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 ($M_{HE, binary}$), which evolve from an MS star with $M_{ZAMS} = 100\, M_{\odot}$  in an interacting binary cannot be above $M_{HE, binary} \sim 31.7 \,M_{\odot}$. The single star also loses its envelope eventually, mostly due to LBV winds. However, this occurs well after the onset of the 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  $M_{HE,single} = 35.2 M_{\odot}$).

The increase in the mass lost due to Wolf-Rayet-like winds for interacting binaries diminishes the importance of weakening the optically thin line-driven winds.
As shown in Figure \ref{fig:42_single_vs_binary}, in Model II, the primary star of the binary system develops a more massive helium core before the envelope loss ($M_{core} \approx 40.3 \,M_{\odot}$) than in Model I. However, it then loses an enormous amount of mass via Wolf-Rayet winds and ends up with a black hole that is only about 2.4$\,M_{\odot}$ more massive than the black hole evolving from the same system in Model I.
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 ($\sim 20.3 \,M_{\odot}$ as opposed to $\sim13.7\, M_{\odot}$).  In Figure \ref{fig:42_single_vs_binary}, 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 $M_{BH} \sim 28.3\, M_{\odot}$, which is about $10\,M_{\odot}$ more than than the black hole forming in the binary with Model I.

\subsection{The effect of stellar winds on merging binary black holes}

\begin{figure*}
\includegraphics[width=1.0\textwidth]{stellarwind_histogram_stable.pdf} \hfill
\caption{2D histograms of masses of the merging binary black holes for the three stellar wind models at $Z = 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 $Z = 0.01$. }
  \label{fig:windscatter}
\end{figure*}

In the last subsection, we discussed how the masses of black holes formed in interacting binaries are affected by different stellar wind models. In this subsection, we turn our attention to the merging binary black holes specifically. In Figure \ref{fig:wind_primary_first}, we plot the primary mass distribution of merging binary black holes for all of our three stellar wind models with different values of $\zeta_{ad,rad}$ and $\gamma$ at metallicities $Z = 0.02,\, 0.01,\,0.005$, but a fixed $\beta = 0.3$. Firstly, we see that decreasing the mass loss rates of optically thin stellar winds by a factor of three only increases the maximum mass of a merging binary black holes by $\sim\,4\,M_{\odot}$ at most. This is not surprising, as we found similar results in subsection \ref{subsec:interacting}.
Secondly, whether the most massive black holes of Model III end up in merging binary systems depends on what we assume about the mass transfer stability criteria for giants with radiative donors. If we assume $\zeta_{ad,rad} = 4$, the maximum masses of black holes do not differ much between Model II and Model III (Figure \ref{fig:wind_primary_first}). However, this is no longer the case for $\zeta_{ad,rad} = 7.5$. In particular, the merging binary black hole masses can reach up to $M_{BH}\sim\,30\,M_{\odot}$ for Model III at $Z = 0.01$.

In order to understand this, let us consider how the most massive merging black hole binaries are formed at metallicties $Z \gtrsim 0.005$:
\begin{itemize}
    \item \textit{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 $M_{ZAMS} \gtrsim 40\, M_{\odot}$ (see Figure \ref{fig:Wind1}) and therefore so are Case C mass transfer events for such systems.
    \item \textit{As opposed to low metallicity environments, the most massive binaries can achieve considerably larger values of $q_{MT,2}$.} This is because a star with a given $M_{ZAMS}$ forms a black hole with significantly lower mass at such metallicities.
    \item \textit{The minimum mass ratio required for the formation of GW sources at the onset of the second mass transfer of the most massive systems is larger compared to the low metallicity case.} Stellar winds are strong enough to widen the orbit significantly and therefore the second mass transfer phase has to shrink the orbit even more efficiently compared to the low metallicity case. This essentially means that the minimum $q_{MT,2}$ for a given $M_{BH,1}$ needs to be larger than for systems at $Z = 0.0007$.
\end{itemize}

We have seen that in the low metallicity case, increasing $\zeta_{ad,rad}$ from 4 to 7.5 does not influences the maximum mass of the GW sources. However at metallicities $Z\gtrsim 0.005$, this is no longer the case, as the formation of these systems typically require $q_{MT,2} > 3.2$ (e.g. see last bullet point). The effect of the lowered stellar winds on the mass spectrum of merging binary black holes winds is only significant, if $\zeta_{ad,rad}$ is large enough, such that the second phase of mass transfer of binaries with the most massive stars and with $3.2\lesssim q_{MT,2}\lesssim q_{crit}$ is still stable.

We also note that the finding - which was discussed in subsection \ref{sub:effect_of_first_mt} - which states that the if first mass transfer is Case A, then the resulting $q_{MT,2}$ typically will be above $q_{crit}$ is still true for $Z \gtrsim 0.005$.

In Figure \ref{fig:windscatter}, we show a  2D histograms of masses of the merging binary black holes for the three stellar wind models at Z = 0.01 for each stellar wind model. We further 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 are indeed from the stable channel and with Case B first mass transfer episode. The most massive black holes from these models have mass ratios $0.8\lesssim q_{final} \lesssim 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 a main sequence donors and the mass ratio distribution of the merging binary black holes are in the range of $0.6 \lesssim q_{final} \lesssim 0.8$.

\subsection{The effects of LBV winds on the merging binary black hole population}
The very uncertain LBV mass loss rates have a strong impact on the merging binary black hole population in various of our model variations. In some cases, it can significantly decrease the merger rate of these systems.
This can occur in the following ways.
Firstly, the range  of the mass ratio distribution is narrowed down to lower values at the onset of the second mass transfer by the LBV mass loss rates. This can decreases the number of mergers of the rCEE channel
 (see also subsection \ref{sub:effect_of_first_mt}).
Secondly, they also affect the maximum size that these stars eventually reach.  According to \citealt{Mennekens_2014} , if the LBV mass loss rate is in the order of $\dot{M}_{LBV} = 10^{-3}\,M_{\odot}yr^{-1}$, the merger rate of binary black holes drops to zero. We find similar results, although only for the CEE channel.
With such high LBV mass loss rates, stars lose their hydrogen envelopes very soon after the onset of core helium burning. Therefore, the number of Case C mass transfer phases drops dramatically, which makes the CEE channel completely inefficient.
Thirdly, intense mass loss rates widen the orbit. This produces binary black holes which are too wide to merge within Hubble time.

While the CEE channel is sensitively dependent on the assumed LBV mass loss rate, we find that this is not the case for the stable channel.
 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.

We find the majority of such sources never evolve beyond the Humphreys-Davidson limit, and even if they do, they would only shed a few solar masses, which does not significantly reduce the mass ratios at the onset of the second mass transfer. In fact, we find that according to our simulations 85-99$\%$ of the gravitational wave sources from the stable channel never evolve beyond the Humphreys-Davidson limit at metallicities $Z\lesssim 0.0007$ depending on the model variation. At higher metallicities, this number decreases. We find the lowest fraction to be 55.5$\%$, which occurs for our  $\gamma = 1$, $\beta= 0.3$, $\zeta_{ad,rad} = 7.5$ models at Z = 0.05. However, even in this case, systems which evolve beyond the Humphreys-Davidson limit typically spend there no more than $\sim10\%$ of their Hertzprung gap lifetime. This means that even if the mass loss rate is $\dot{M}_{LBV} = 10^{-3}\,M_{\odot}yr^{-1}$, the donor would only lose $\sim1\,M_{\odot}$ mass, and therefore the predictions for the stable channel are little affected by a steady, intense LBV mass loss.

Whether or not the LBV mass loss shuts down the formation of merging binary black holes (under the assumption of $\dot{M}_{LBV} = 10^{-3}\,M_{\odot}yr^{-1}$) is therefore dependent on
which sub-channels are dominant.

\section{Conclusion}
\label{sub:conclusion}
We have performed a parameter study of the classical isolated binary formation channel for gravitational sources in order to better understand how the properties of the merging binary black hole population depend on current uncertainties related to stable mass transfer episodes and stellar winds. For this, we have used the rapid population synthesis code SeBa to simulate the evolution of  massive binaries over a metallicity range $Z = 0.0001 - 0.02$ with several model variations each with different assumptions about the binary physics. We have varied the mass transfer efficiency, the specific angular momentum lost from the system during mass transfer episodes and the mass transfer stability criteria for giants donors in order to characterise current uncertainties in our understanding of stable mass transfer episodes. We have also varied the boundary where giants develop a deep convective envelope as this plays an extremely important role for gravitational wave source formation via common envelope evolution.
At higher metallicities, where line-driven winds are assumed to significantly affect the evolution of massive stars, we have also varied the mass loss rates to understand their impact on merging binary black holes. While many of our models (e.g. constant value for mass transfer efficiency, angular momentum loss, mass transfer stability criteria) are very likely not adequate to describe real massive binaries, these results still help our understanding how these parameters shape the properties of the gravitational waves population.

We find that within the isolated binary formation scenario, current uncertainties related to stable mass transfer do not only change the observable properties of merging binary black holes, but they also change which channel dominates the rates. In many of our model variations, the most common formation path involves two stable mass transfer phases (stable channel, \citealt{vandenheuvel17}), whereas in others a stable mass transfer episode followed by an unstable mass transfer phase (CEE channel) is responsible for the formation of most of the GW sources.

We also find that the mass transfer stability criteria strongly influences the merger rate of the stable channel (in broad agreement with recent studies such as \citet{Neijssel_2019}, \citet{olejak2021impact}, \citet{Lieke2022ApJ...931...17V}. However the degree and the quality of the influence are strongly interrelated with other uncertain binary physics parameters, such as angular momentum loss during mass transfer, mass transfer efficiency etc.
We test two assumptions for the mass transfer stability criteria for giants with radiative donors; $\zeta_{ad,rad} = 4$ and $\zeta_{ad,rad}$ = 7.5. These correspond to a critical mass ratios of $q_{crit} = 3.2$ and $q_{crit}$ = 5.5, respectively, if the accretor is a black hole and the accretion is Eddington limited (see subsection \ref{subsec:mt_stability}). The latter assumption is somewhat more consistent with recent studies of \citet{Pavlovskii_2016}, \citet{Ge_2020}, \citet{marchant2021role}. To summarise our main results:
\begin{itemize}
    \item \textit{Effect of mass transfer stability criteria:} The merger rate of GW sources of the stable channel increases with increasing $\zeta_{ad,rad}$. This is because the degree of orbital shrinkage during the stable mass transfer phase with a black hole accretor is crucially dependent on the mass ratio, i.e. it increases exponentially with increasing mass ratios.
    However, we also find that the critical mass ratio for stable mass transfer is not achievable in every model.
    At low mass transfer efficiencies ($\beta = 0.3$) and at metallicites $Z \lesssim 0.005$, the maximum mass ratios at the onset of the second mass transfer are below 3.2 for initial primary masses of $M_{ZAMS,1} \gtrsim 60\, M_{\odot}$.
    As a consequence, in our low mass transfer efficiency models, only the merger rate of black holes with pirimary masses $M_{BH} \lesssim 20\,M_{\odot}$ is affected when $\zeta_{ad,rad}$ is increased.
    \item \textit{Minimum initial separation of progenitors of the stable channel:}
    We find that there is a minimum orbital separation that the progenitors of the stable channel have at the zero-age main-sequence. This minimum initial separation is very close to the separation at which the initially primary star fills its Roche-lobe just when it evolved off main sequence. That is, the formation of GW sources is not typically possible, if the donor is a main sequence star during the first phase of mass transfer\footnote{While this is not strictly true for a few of our model variations, even in those cases the donor has to be very close to the end of its main sequence phase.}. This occurs because these systems develop mass ratios by the onset of the second mass transfer phase which are typically above the critical mass ratio, and consequently they undergo a dynamically unstable mass transfer phase and they merge before they would form a black hole binary.
    \item \textit{Effect of angular momentum loss on the stable channel:} The merger rate of the stable channel  sensitively depends on the assumed angular momentum loss mode during the non-conservative mass transfer phase. The angular momentum loss is most uncertain for the first phase of mass transfer, when the accretor is not a compact star yet.
    We find that the orbits of massive binaries typically widen too much after the first mass transfer phase, if we assume that the expelled matter leaves the system with the specific angular momentum of the binary (i.e. $\gamma  = (dJ/dM_{tot})/(J/M_{tot}) = 1$, following the assumptions of \citealt{Podsiadlowski92}, \citealt{Belczynski2008}).
     On the other hand, if we assume $\gamma  = 2.5$ (\citealt{Portegies96}), the formation of GW systems via the stable channel become possible, since in this case the net change in the orbit of the binary is reduced. We note that while the assumption of constant $\gamma$ is most likely not physical and should only be considered for simplified models.
    \item \textit{Effect of angular momentum loss on the CEE channel:} Both $\gamma$ and $\beta$  significantly affect the merger rate and the mass ratio distribution of the cCEE channel, i.e. the channel which involves a common envelope phase with a donor with convective envelope.
    This occurs for two reasons. Firstly, the critical mass ratio in case of a mass transfer episode between a non-compact accretor and a donor with a deep convective envelope can vary between $q_{crit}\sim 0.9-1.9$ depending on the assumed values of $\gamma$ and $\beta$. Larger values of $q_{crit}$ occur for $\gamma = 1$. In fact, we find that in some of our model variations with $\gamma = 1$, about 70 $\%$ of the giants with deep convective envelopes are predicted to initiate stable mass transfer episodes when the accretor is a non-compact object, as opposed to $\sim 20\%$ in some of the $\gamma = 2.5$ models. On the other hand, during the second mass transfer phase when the accretor is a black hole and the donor has a deep convective envelope, the critical mass ratio is $q_{crit} \approx 1.75$, regardless of the assumed $\gamma$ (since by definition $\gamma$ and $\beta$ is only meant for non-compact accretors and for black holes we assume that the accretion is Eddington-limited). Consequently, in the $\gamma = 1$ models, the merger rate of the cCEE channel is higher than in the $\gamma = 2.5$ model, since in the former case donors with deep convective envelopes more readily instigate a stable mass transfer phases \footnote{The first mass transfer phase in the CEE channel is stable by definition. Even if the binary undergoes an unstable first phase of mass transfer and survives it, we find that such a channel does not contribute to the overall merger rate in appreciable numbers.}. Secondly, if the first mass transfer phase typically widens the orbit (as is the case for $\gamma =1$), then a system in which the initial primary fills its Roche-lobe with a radiative envelope is more likely to have the second mass transfer phase with a giant donor with a deep convective envelope. This is because a widened orbit implies that the initially secondary star fills its Roche-lobe at a later evolutionary stage, at which a deep convective envelope might have already been developed.
    \item \textit{Line driven stellar winds:} While decreasing the mass loss rates of  optically thin line driven winds by a factor of three increase the black hole masses in single stellar evolution, we do not find an appreciable difference in the masses of merging binary black holes. This is because at metallicities where line-driven stellar winds are deemed important, the stars typically undergo mass transfer episodes on the Hertzsprung gap (that is mass transfer episodes with core-helium burning donors are almost non-existent due to negligible expansion of radius of such stars). As the donor loses its hydrogen envelope during hydrogen shell burning, the growth of the helium core is halted very early on. This is followed by a longer Wolf-Rayet phase, which results in larger mass loss with respect to their single counterparts. The halted growth of the helium core and increased Wolf-Rayet mass losses diminish the importance of lowered mass loss rates of optically thin line driven stellar winds.
    However, if we simultaneously decrease the Wolf-Rayet winds too, we find a significant increase in the masses of gravitational wave sources for the $\zeta_{ad,rad}=7.5$ model alone.
    \item \textit{LBV winds:}
    Intense LBV mass loss rates can have significant effect on the CEE channel. If the mass loss rates are above $\dot{M}_{LBV} \sim 10^{-3}\,M_{\odot}yr^{-1}$, the merger rate of the CEE channel becomes negligible. On the other 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 they would evolve beyond the Humpyhreys-Davidson limit.
\end{itemize}

Our results suggest that the current uncertainties
in the first mass transfer phase have a huge impact on the predicted merging binary black hole population.
We have also shown that seemingly unrelated parameters can have a complex, interrelated effects on the population properties of GW sources.
This implies that along with the poorly understood common envelope evolution, core collapse and natal kick of black hole progenitors, it is currently extremely difficult to infer parameters related to binary or stellar physics directly from LIGO and Virgo observations (see also \citealt{Startrack2022ApJ...925...69B}). In order to break the degeneracy between the uncertain parameters and the predicted population of merging binary black holes, it is clear that we have to analyse the gravitational wave data along with electromagnetic observations.
In particular, regarding the first mass transfer phase, observations of WR-O/B binaries and high mass X-ray binaries could offer invaluable clues about the how the first mass transfer episode proceeds. Especially an extensive catalogue of the former type would be useful, as these systems did not go through core-collapse, and their mass ratios and orbital separations are directly related to the mass transfer efficiency and the angular momentum carried away by the expelled matter.

Finally, we note that with recent developements in stellar physics, it has been shown that the stable channel could be important for the formation of merging stellar mass binary black holes. As these sources originate from relative short period binaries, a better modelling of Case A mass transfer phase in rapid population synthesis codes based on stellar tracks of \citet{Hurley_2000} is required.

\section*{Acknowledgements}

 Computational work was performed University of Birmingham's BlueBEAR HPC service.
 ST acknowledges support from the Netherlands Research Council NWO (VENI 639.041.645 and VIDI 203.061 grants

\bibliographystyle{mnras}