# Wolf-Rayet – compact object binaries as progenitors of binary compact objects

Erika Korb, Michela Mapelli, Giuliano Iorio, Guglielmo Costa

► To cite this version:

Erika Korb, Michela Mapelli, Giuliano Iorio, Guglielmo Costa. Wolf-Rayet – compact object binaries as progenitors of binary compact objects. Astronomy & Astrophysics - A&A, 2025, 695, pp.A199. 10.1051/0004-6361/202452585 . hal-04757918

HAL Id: hal-04757918

https://hal.science/hal-04757918v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Wolf–Rayet–compact object binaries as progenitors of binary compact objects

Erika Korb<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*, Michela Mapelli<fnref n="3" />,<fnref n="1" />,<fnref n="2" />,<fnref n="4" />,*, Giuliano Iorio<fnref n="1" />,<fnref n="2" />,<fnref n="4" />,†,
Guglielmo Costa<fnref n="5" />,<fnref n="4" />,‡, and Marco Dall'Amico<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,§

Received 11 October 2024 / Accepted 4 February 2025

## ABSTRACT

Binaries with a Wolf–Rayet star and a compact object (WR–COs), either a black hole (BH) or a neutron star (NS), have been proposed as possible progenitors for the binary compact object mergers (BCOs) observed with gravitational wave (GW) detectors. In this work, we use the open source population synthesis code SEVN to investigate the role of WR–COs as BCO progenitors. We consider an initial population of $5 \times 10^6$ binaries, and we evolve it across 96 combinations of metallicities, common envelope efficiencies, core-collapse supernova models, and natal kick distributions. We find that WR–COs are the progenitors of most BCOs, especially at high and intermediate metallicity. At $Z = 0.02, 0.014$, and $0.0014$, more than $\ge 99\%$ of all the BCOs in our simulations evolved as WR–COs. At $Z = 0.00014$, inefficient binary stripping lowers the fraction of BCOs with WR–CO progenitors to $\approx 83–95\%$. Despite their key role in BCO production, only $\approx 5–30\%$ of WR–COs end their life as BCOs. We find that Cyg X-3, the only WR–CO candidate observed in the Milky Way, is a promising BCO progenitor, especially if it hosts a BH. In our simulations, about 70–100% of the Cyg X-3-like systems in the WR–BH configuration (BH mass $\le 10 M_\odot$) are BCO progenitors, in agreement with the literature. Future observations of WR–COs similar to Cyg X-3 may be the Rosetta stone to interpret the formation of BCOs.

**Key words.** gravitational waves – methods: numerical – binaries: general – stars: black holes – stars: Wolf–Rayet

## 1. Introduction

The LIGO–Virgo–KAGRA (LVK) collaboration has detected almost one hundred gravitational wave (GW) event candidates as a result of its first three observing runs (Abbott et al. 2023a; Wadekar et al. 2024). All the confirmed events are associated with binary compact object mergers (BCOs); most of them are binary black holes (BBHs), but there are also binary neutron stars (BNSs; e.g., GW170817 and GW190425; Abbott et al. 2019, 2020), and black hole-neutron stars (BHNSs; e.g., GW200105 and GW200115; Abbott et al. 2021). The fourth observing run of LVK is ongoing and has already led to more than a hundred significant GW event candidates<fnref n="1" />. Among them, the BHNS system GW230529 comprises a black hole (BH) in the mass range $3–5 M_\odot$, confirming that there is no empty gap between the maximum neutron star (NS) and the minimum BH mass (Abac et al. 2024).

Binaries hosting a Wolf–Rayet star (WR) and a compact object (CO), either a BH or an NS, are regarded as promising BCO progenitor candidates (Belczynski et al. 2012). Tight WR–CO binaries already host a CO, the WR can end its life as an NS or BH (Limongi & Chieffi 2010; Chen et al. 2015), and these systems likely interacted via at least one mass transfer episode – as expected for BNSs, BHNSs, and BBHs progenitors (see,

e.g., Bethe & Brown 1998; Belczynski et al. 2008; Dominik et al. 2012; Giacobbo & Mapelli 2018; Vigna-Gómez et al. 2020; Belczynski et al. 2020; Bavera et al. 2021; Santoliquido et al. 2021; Broekgaarden et al. 2022a; Iorio et al. 2023; Boesky et al. 2024).

Mass transfer episodes can partially or totally remove the external layers of a star, producing a binary stripped star with a modified internal structure and evolution (e.g., Laplace et al. 2020; Schneider et al. 2021, 2023). Stripped stars can form via binary stripping but also via self-stripping, if they have strong stellar winds (typical mass loss rates of WRs are up to $\dot{M} \approx 10^{-4} M_\odot$; see, e.g., Sen et al. 2021; Sabhahit et al. 2022; Sander et al. 2023; Götberg et al. 2023). Regardless of the stripping process, stripped stars can be visible as WRs if they are sufficiently hot (temperature at the sonic point $\ge 40 \times 10^3$ K; Nugis & Lamers 2002; Gräfener & Vink 2013; Gräfener et al. 2017), luminous ($L \ge 10^5 L_\odot$ at solar metallicity; Shenar et al. 2020), and exhibit emission line-dominated spectra (e.g., helium, carbon, nitrogen lines; see the review by Crowther 2007). Observations in the Milky Way (MW), Small Magellanic Cloud, and Large Magellanic Cloud (Van Der Hucht 2001; Bartzakos et al. 2001; Foellmi et al. 2003a,b; Schnurr et al. 2008) indicate that binary stripping is a common formation channel for WRs (binary fractions $\ge 40\%$), but it does not become dominant at decreasing metallicity, where stellar winds are weaker. Shenar et al. (2020) argue that the WR binary formation channel has a non-trivial dependence on metallicity and is heavily affected by WR spectral definitions and wind models.

<footnote n="1">¹ Physics and Astronomy Department Galileo Galilei University of Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy</footnote>

<footnote n="2">² INFN – Padova, Via Marzolo 8, 35131 Padova, Italy</footnote>

<footnote n="3">³ Institut für Theoretische Astrophysik, ZAH, Universität Heidelberg, Albert-Ueberle-Straße 2, 69120 Heidelberg, Germany</footnote>

<footnote n="4">⁴ INAF – Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="5">⁵ Univ Lyon, Univ Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France</footnote>

<footnote>* Corresponding authors; erika.korb@studenti.unipd.it; mapelli@uni-heidelberg.de</footnote>

<footnote n="1">¹ Updated list of significant event candidates at https://gracedb.ligo.org/</footnote>
Wolf-Rayet-compact object systems are one of the possible electromagnetic progenitors of BCOs. So far, seven WR-CO candidates have been observed (see Esposito et al. 2015 and references therein), but their characterization is challenging. The strong WR winds prevent a precise dynamical CO mass measurement, leaving uncertainties on the BH or NS nature. Uncertainties in the CO nature also affect the only WR-CO candidate known in the MW: Cyg X-3 (e.g., Zdziarski et al. 2013; Laycock et al. 2015; Koljonen & Maccarone 2017; Binder et al. 2021). Belczynski et al. (2013) showed that observational and theoretical uncertainties, such as models of core-collapse supernovae (CCSNe) and wind mass loss, significantly affect the possible fates of WR-COs similar to Cyg X-3 and should be taken into account in evolutionary studies.

In this work, we investigate the role of binaries hosting a WR and a CO as BCO progenitors in the isolated evolution scenario. We used the population-synthesis code SEVN (Iorio et al. 2023) to simulate the evolution of a binary population similar to the MW one under different conditions. We probed a large parameter space to characterize BCO formation and bracket uncertainties in theoretical models. Key parameters for modeling evolutionary processes, such as mass transfer and CCSNe, are still poorly constrained, and their impact on the demography of BCOs is large (see the reviews by Mapelli 2021 and Mandel & Broekgaarden 2022, and references therein). Here, we systematically explore the role of stellar metallicity, common envelope (CE) ejection efficiency, CCSN physics, and natal kicks. We also discuss the impact of uncertainties in mass transfer physics on the WR-CO and BCO formation channels. Eventually, we consider the implications for BH spins in the case of tidal spin-up.

This paper is organized as follows. In Section 2, we describe the features of SEVN relevant to this work and the assumptions on the parameter space. In Section 3, we present the results of our simulations and a comparison with Cyg X-3. In Section 4, we discuss the impact of mass transfer assumptions, BH tidal spin-up, and other aspects of our parameter-space exploration. In Section 5, we summarize our results.

## 2. Methods

### 2.1. The sevn code

We performed our simulations using the SEVN population-synthesis code (Spera et al. 2015; Spera & Mapelli 2017; Spera et al. 2019; Mapelli et al. 2020; Iorio et al. 2023). SEVN interpolates stellar properties from a grid of pre-computed stellar tracks and implements binary evolution processes through analytic and semi-analytic models. We report hereafter only the assumptions most relevant for this work; we refer to Iorio et al. (2023) for a more detailed discussion about SEVN.

### 2.2. WR definition

Wolf-Rayet stars are commonly defined based on their spectral characteristics (Crowther 2007), but SEVN does not model the spectral emission of stars directly. Thus, we do not use a spectrum-based definition for WRs but we approximate it with a mass-based definition. We assume that a WR star is a helium star: a core-He or shell-He burning star with helium core mass constituting most of its stellar mass ($M_{\text{He-core}} \ge 97.9\% M_{\text{star}}$).

From an observational point of view, helium stars are stripped-stars that may or may not already exhibit the spectral features of WRs (Shenar et al. 2020). As the stripping proceeds, the hotter internal layers of the He star are exposed and the

surface abundance of advanced burning products is expected to increase. Eventually, the spectrum of a He star will exhibit strong helium, carbon or oxygen emission lines, resembling the one of a WR star of the WN, WC or WO sub-type, respectively (Crowther 2007). Our definition of WR is based on the properties of the WN sub-type at Galactic metallicity and on the WR spectral models elaborated with the PoWR code (Hamann & Gräfener 2004; Hainich et al. 2014).

### 2.3. Stellar tracks

We interpolate stellar properties from two set of tracks: H-rich stars and pure-He stars. These tracks were calculated with the PARSEC stellar evolution code (Bressan et al. 2012; Chen et al. 2015; Costa et al. 2019; Nguyen et al. 2022; Costa et al. 2025), with the version described in Costa et al. (2025).

H-rich stars evolve from the pre-main sequence (pre-MS) to the end of the carbon burning phase. In such models, the mass-loss rate of hot massive stars depends on both stellar metallicity and Eddington ratio (Vink et al. 2000, 2001; Gräfener & Hamann 2008; Vink et al. 2011). Specifically, winds of WRs follow the prescriptions presented by Sander et al. (2019), with the additional metallicity dependence adopted by Costa et al. (2021).

Pure-He stars evolve from a He zero-age main sequence (ZAMS) to the end of the carbon burning phase. He-ZAMS stars are initialized accreting mass on a $0.36 M_{\odot}$ He core, obtained removing the H-rich envelope from the progenitor star when it started the He-core burning phase. The initial metallicity (Z) determines the initial helium mass fraction ($Y = 1 - Z$) of stars in the He-ZAMS. Winds of pure-He stars follow Nugis & Lamers (2000). In this work, we used this set of pure-He tracks to interpolate the properties of the He-stars that we used as proxy for WRs.

Since the definition of WR star depends on the helium core mass, we recall here the most important assumptions for its evolution adopted in PARSEC. We define the radius of the He core as the most external shell with hydrogen mass fraction $X < 10^{-3}$. For the convection treatment, we adopt the Schwarzschild criterion (Schwarzschild 1958) and the mixing-length theory (Böhm-Vitense 1958) with a solar-calibrated value for the parameter $\alpha_{\text{MLT}} = 1.74$ (Bressan et al. 2012). We consider a penetrative overshooting scheme above the core boundary limit and we parametrize it with the mean free path $\lambda_{\text{ov}}$ of the unstable elements. In the tracks adopted in this work, $\lambda_{\text{ov}} = 0.5 H_p$, where $H_p$ is the pressure scale height; according to the ballistic approach (Bressan et al. 1981). This assumption corresponds to a value of about $f_{\text{ov}} = 0.025$ in the exponential decay overshooting formalism (diffusive scheme; Herwig et al. 1997). For more details about the PARSEC tracks adopted in this work we refer to Costa et al. (2025) and Iorio et al. (2023).

### 2.4. CCSN models

In this work, we tested three CCSN models for CO formation: the rapid and delayed ones proposed by Fryer et al. (2012) and a compactness-based one that uses a compactness parameter (e.g., O'Connor & Ott 2011), as done in Mapelli et al. (2020). The three models calculate the remnant mass as a function of the carbon-oxygen core mass at the end of the carbon burning phase ($M_{\text{CO-core}}$).

Regardless on the CCSN model, in this work, we used the fiducial SEVN classification for COs (Iorio et al. 2023): we considered a remnant a BH if its mass is $M_{\text{rem}} \ge 3 M_{\odot}$, an NS otherwise.
### 2.4.1. Rapid and delayed CCSN models

The rapid and delayed models assume that the shock is revived by the energy accumulated by neutrinos in a convective region between the infalling external layers and the surface of the proto-NS. Fryer et al. (2012) distinguish rapid or delayed explosions considering shocks revived <250 ms or >500 ms after core bounce, respectively. In both models, the remnant mass $M_{rem}$ is expressed as

$$M_{rem} = M_{proto} + f_{fb} (M_{star} - M_{proto}), \quad (1)$$

where $f_{fb}$ (fallback parameter) is the fraction of stellar mass $M_{star}$ that falls back to the proto-compact object ($M_{proto}$). Both $f_{fb}$ and $M_{proto}$ are functions of $M_{CO-core}$. The heaviest stars ($M_{CO-core} \ge 11 M_{\odot}$) always undergo a direct collapse, accreting all the ejected material ($f_{fb} = 1$), and forming a BH.

### 2.4.2. Compactness-based CCSN model

The compactness-based CCSN model tested in this work uses the compactness parameter $\xi_{2.5}$ (O'Connor & Ott 2011) as a proxy for the explodability of stellar cores. The parameter $\xi_{2.5}$ is defined as the ratio between a reference mass (2.5 $M_{\odot}$ here) and the radius that encloses such mass at the onset of the core-collapse. In SEVN, we calculate $\xi_{2.5}$ following the fit derived by Mapelli et al. (2020) on non-rotating stellar tracks generated with the FRANEC code (Limongi & Chieffi 2018):

$$\xi_{2.5} = \frac{M = 2.5 M_{\odot}/M_{\odot}}{R(M = 2.5 M_{\odot})/1000 \text{ km}} = 0.55 - 1.1 \left( \frac{M_{\text{CO-core}}}{1 M_{\odot}} \right)^{-1.0}. \quad (2)$$

Several studies suggest that there is a complex relation between the compactness parameter and the carbon-oxygen core mass (e.g., Sukhbold et al. 2018; Patton & Sukhbold 2020; Chieffi & Limongi 2020; Schneider et al. 2021, 2023). In this work, we followed the conservative approach by Mapelli et al. (2020), based on the results of Limongi & Chieffi (2018), and considered a monotonic relation.

Here, we assumed that if $\xi_{2.5} \le 0.35$ ($M_{CO-core} \le 5.5 M_{\odot}$), the CCSN is successful and produces an NS. The NS mass was drawn from a Gaussian distribution centered at 1.33 $M_{\odot}$, with standard deviation of 0.09 $M_{\odot}$ and minimum NS mass set to 1.1 $M_{\odot}$, in agreement with NSs observations in binaries (Özel & Freire 2016). If $\xi_{2.5} > 0.35$, the star collapses directly to a BH, with mass $M_{BH}$,

$$M_{BH} = M_{\text{He}} + f_{\text{H}} (M_{\text{star}} - M_{\text{He}}), \quad (3)$$

where $M_{\text{He}}$ is the He-core mass at the ignition of the carbon burning and $f_{\text{H}}$ is the fraction of the residual hydrogen-rich envelope that collapses to a BH. The amount of $f_{\text{H}}$ is substantially unconstrained. The work by Fernández et al. (2018) shows that it depends on the binding energy of the envelope at the onset of core collapse. Hence, our choice of $f_{\text{H}} = 0.9$ should be regarded as a crude approximation.

### 2.5. Natal kicks

Observations of CO proper motions indicate that compact remnants receive a kick at their birth (Hobbs et al. 2005; Verbunt et al. 2017; Atri et al. 2019). We tested four models for the magnitude of the kick following the work by Hobbs et al. (2005),

Atri et al. (2019), Fryer et al. (2012), and Giacobbo & Mapelli (2020). We refer to the corresponding models with the acronyms M265, M70, Mfb, and G20.

Hobbs et al. (2005) fitted the proper motions of young pulsars in the MW with a Maxwellian function with root-mean-square (rms) velocity $\sigma = 265 \text{ km s}^{-1}$. Similarly, Atri et al. (2019) measured the proper motions of BHs in X-ray binaries and fitted their velocity distribution with a Maxwellian with a mean of $107 \pm 16 \text{ km s}^{-1}$, corresponding to a rms of $\sigma = 70 \text{ km s}^{-1}$. Therefore, in the M265 (M70) model we assign kick magnitudes drawing random numbers from a Maxwellian distribution with rms of $\sigma = 265 \text{ km s}^{-1}$ ($\sigma = 70 \text{ km s}^{-1}$).

Fryer et al. (2012) assumed that the kick is distributed according to a Maxwellian function with $\sigma = 265 \text{ km s}^{-1}$, modulated by the fraction $f_{fb,CCSN}$ of mass that falls back onto the proto-NS:

$$v_{\text{kick}} = f_{\sigma=265} (1 - f_{\text{fb,CCSN}}), \quad (4)$$

where $f_{\sigma}$ is a random number extracted from the Maxwellian curve. In our simulations with the Mfb model, we choose $f_{fb,CCSN}$ consistently with the CCSN model (see Appendix A.1).

In the G20 model, we followed Giacobbo & Mapelli (2020) and re-scaled the kicks drawn from the M265 model:

$$v_{\text{kick}} = f_{\sigma=265} \frac{\langle M_{\text{NS}} \rangle M_{\text{NS}}}{M_{\text{rem}} \langle M_{\text{ej}} \rangle}, \quad (5)$$

where $\langle M_{\text{NS}} \rangle$ and $\langle M_{\text{ej}} \rangle$ are, respectively, the average NS mass and ejecta mass (their values are sensible to the single stellar evolution and CCSN model assumed).

### 2.6. Binary evolution and mass transfer

Binary evolution processes – such as mass transfer, collisions, and tides – are modeled in SEVN following Hurley et al. (2002), with a few differences thoroughly described by Iorio et al. (2023). Hereafter, we summarize the assumptions for mass transfer processes.

A star expanding beyond its Roche lobe donates its mass to the companion in a Roche lobe overflow (RLO) process. As Hurley et al. (2002), we determine RLO stability by computing the donor-to-accretor mass ratio $q = M_{d}/M_{a}$ and comparing it to a set of critical thresholds $q_{c}$. If $q > q_{c}$, mass transfer becomes unstable and we either assume that the two stars merge directly (e.g., if they both fill their Roche lobe), or we start a phase of CE. For a donor star in the main sequence or the Hertzsprung gap phase, we assume $q_{c} = \infty$, since such stars have radiative envelopes, which are less likely to develop unstable mass transfer (Ge et al. 2010, 2015). We adopt the same $q_{c}$ values as Hurley et al. (2002) for stars in more advanced burning stages. We follow their fit on models of condensed polytropes (Hjellming & Webbink 1987) for giant stars with deep convective envelopes; we set $q_{c} = 3.0$ for core-Helium burning stars (both H-rich and pure-Helium ones), and we use $q_{c} = 0.784$ for pure-Helium stars with active Helium-shell burning.

We describe the orbital changes due to CE with the usual two parameters $\alpha_{CE}$ and $\lambda_{CE}$ (Webbink 1984). Here, the parameter $\alpha_{CE}$ is the fraction of the orbital energy that is transferred to the envelope, while $\lambda_{CE}$ describes the concentration of the donor’s profile. For H-rich stars, we use the fits to $\lambda_{CE}$ proposed by Claeys et al. (2014), while for pure-He stars we assume $\lambda_{CE} = 0.5$ (Iorio et al. 2023). According to the original formalism (Webbink 1984), the envelope is ejected when the orbital energy lost by the binary equals the envelope binding energy. Thus, $\alpha_{CE}$
<table><caption>**Table 1.** Explored combinations of CCSN and natal kick models.</caption><thead><tr><th>CCSN</th><th>Kick</th><th>M265</th><th>M70</th><th>Mfb</th><th>G20</th></tr></thead><tbody><tr><td>Delayed</td><td></td><td>dM265</td><td>dM70</td><td>dMfb</td><td>dG20</td></tr><tr><td>Rapid</td><td></td><td>rM265</td><td>rM70</td><td>rMfb</td><td>rG20</td></tr><tr><td>Compactness</td><td></td><td>cM265</td><td>cM70</td><td>cMfb</td><td>cG20</td></tr></tbody></table>

cannot be larger than one. However, the orbital energy is not the only source of energy that can be used to eject the CE. Efficient recombination of H and He inside the envelope lowers its binding energy, facilitating its removal. This and other additional sources of energy should be taken into account when modeling CE evolution (e.g., Paczyński & Ziółkowski 1968; Han et al. 1994, 1995; Ivanova et al. 2013; Klencki et al. 2021; Röpke & De Marco 2023). Thus, in this work, we also explored values of $\alpha_{CE} > 1$.

## 2.7. Simulations and parameter space

We wanted to characterize the impact of different physical assumptions on the formation and fate of WR–COs. Thus, we fixed our initial binary population and we evolved it across 96 combinations of parameters and models. We generated a population of $5 \times 10^6$ massive binary stars with initial mass, mass ratio, orbital period, and eccentricity drawn from observationally motivated distributions (Kroupa 2001; Sana et al. 2012; Moe & Di Stefano 2017), as described in Appendix A.2.

We chose four possible values for the metallicity ($Z = 0.02, 0.014, 0.0014, 0.00014$) and two possible values for the fraction $\alpha_{CE}$ of the orbital energy required to eject the CE ($\alpha_{CE} = 1, 3$). For each ($Z, \alpha_{CE}$) pair, we explored twelve combinations of CCSN and natal kick models: the rapid, delayed, and compactness-based CCSN models, and the M265, M70, Mfb, and G20 natal kick models. In Table 1 we summarize these combinations and we report the correspondent acronyms.

# 3. Results

## 3.1. Most BCOs form from WR–CO systems

For every combination of natal kick, CCSN model, metallicity, and CE efficiency explored in this work (Table 1), we find that most BNSs ($\ge 99\%$), BHNSs ($\ge 93\%$), and BBHs ($\ge 79\%$) evolved as WR–BHs or WR–NSs before becoming BCOs (Figure 1). In each set, more than $\ge 83\%$ of the total BCO population had a WR–CO progenitor.

Metallicity ($Z$) and CE efficiency ($\alpha_{CE}$) are the quantities that influence more the fraction $f_{from}$ of BCOs with a WR–CO progenitor. Here, we calculated $f_{from}$ as the ratio between the number of BCOs with a WR–CO progenitor (WR–BH or WR–NS) and the total number of BCOs (BBHs, BHNSs or BNSs) in each set. If we fix a pair of ($Z, \alpha_{CE}$) values and consider the distribution of $f_{from}$ across the explored combinations of CCSN and natal kick models, we find that nearly every ($\ge 99\%$) BCO at $Z \ge 0.0014$ evolved through the WR–CO configuration. At $Z = 0.00014$, only $\approx 83–95\%$ of BCOs followed this formation scenario (the range indicates values within the 68% credible intervals, Figure 2). The decrease of the fraction of BCOs with a WR–CO progenitor at low $Z$ is determined by inefficient stripping of the WR and by the combined effect of CCSN and natal kicks models.

## 3.1.1. BHNS progenitors

For all the metallicities explored in this work, we find that BHNSs generally form a BH as first CO (Figure 1). Evolution through the WR–BH configuration is increasingly favoured for BHNSs at lower metallicities, becoming the dominant ($\ge 95\%$) formation channel at $Z = 0.00014$. For decreasing metallicities, stellar winds are weaker and stars are able to retain more mass (Vink et al. 2011). At $Z = 0.00014$, BH production is favoured while NSs require one or more CE events to be produced. Assuming a low CE efficiency ($\alpha_{CE} = 1$) completely suppresses the WR–NS channel for BHNSs at $Z = 0.00014$.

Remnant masses, thus their BH or NS nature, are sensible to the model of CCSN adopted. Assuming different natal kick models influences the natal kick magnitudes and the possibility that a binary is able to survive or merge via GW emission (Figure 3). Thus, the abundance of the WR–BH channel over the WR–NS one is sensible to different combinations of CCSN and natal kick models. For each pair of ($Z, \alpha_{CE}$), we find variations up to 40% in the fraction of BHNSs that evolved as WR–BHs in place of WR–NSs. These variations are larger at solar metallicity ($Z = 0.02, 0.014$).

## 3.1.2. “Missing” WR–BHs

For each pair of ($Z, \alpha_{CE}$), binary evolution is determined by the initial orbital properties only until the formation of the first CO. CCSNe impart a natal kick on the newly born CO, modifying the post-CCSN orbital properties, the following evolution of the system, and the possibility that it becomes a BCO and merges within a Hubble time (Peters 1964). As shown in Figure 1, the influence of CCSN and natal kicks is usually negligible ($\Delta f_{from} \le 0.01$), except for BHNSs and BBHs at $Z = 0.00014$, and BBHs at solar metallicity $Z = 0.02, 0.014$ (variations can be up to $\Delta f_{from} \le 0.15$). These systems can avoid the WR–CO configuration but still produce BCOs.

Weak stellar winds disfavour the production of WRs via self-stripping, especially at the lowest metallicity explored in this work: $Z = 0.00014$. In our simulations, most ($\approx 60–80\%$) BCOs are BBHs at sub-solar metallicity ($Z = 0.0014, 0.00014$). Thus, the lower fraction of BCOs with WR–CO progenitors at $Z = 0.00014$ is determined by the lower fraction of BBHs with WR–BH progenitors.

Binary black holes can form without passing through a WR–BH phase if they experience the second CCSN event while the non-degenerate star is still transferring mass to its BH companion. At low metallicity, if the mass transfer stops before the donor’s envelope is completely stripped, the donor star cannot become a WR and the system becomes a BBH avoiding the BH-WR stage. We find that these “missed WRs” still retain a hydrogen-envelope at the end of their life. Residual envelope masses range from $\approx 50\%$ to almost $\approx 2\%$ of the total stellar mass, indicating that the most stripped stars were almost WRs according to the definition adopted in this work (He-core mass $M_{He-core} \ge 97.9\% M_{star}$, Section 2.2).

## 3.2. Formation channels of BCOs from WR–CO systems

When present, the WR–CO configuration constitutes always the last evolutionary stage before BCO formation. Figure 4 shows the formation channels of BCOs with WR–CO progenitors for the $\alpha_{CE} = 3$ case (the $\alpha_{CE} = 1$ case is shown in the Appendix B.1). We distinguish four channels according to the
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Fraction $f_{\text{from}}$ of BCOs with WR-CO progenitors for all the metallicities (columns), $\alpha_{\text{CE}}$ efficiencies (rows), and combinations of CCSN and natal kick models tested in this work (horizontal entries, see Table 1). We distinguish BCOs into BBHs, BNSs and BHNSs; WR-CO progenitors into systems with a WR star and a BH or NS.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Median fractions $f_{\text{from}}$ of BCOs with WR-CO progenitors as a function of metallicity $Z$ and common envelope efficiency $\alpha_{\text{CE}}$. Left (right) plot shows the $\alpha_{\text{CE}} = 1$ ($\alpha_{\text{CE}} = 3$) case. We distinguish fractions $f_{\text{from}}$ of BBHs evolved as WR-BHs (purple dotted lines, plus sign markers), BNSs evolved as WR-NSs (orange dashed lines, cross sign markers), BHNSs evolved as either WR-NSs or WR-BHs (blue dash-dotted lines, hexagon markers). We also show the fraction of all BCOs with a WR-CO progenitor (black lines, round markers). For each $(Z, \alpha_{\text{CE}})$ pair and BCO type, we consider the distribution of the 12 $f_{\text{from}}$ fractions of the correspondent simulated sets (combinations of CCSN and natal kick models, see Table 1), and we show their median values (scatter points) and 68% credible intervals (error bars).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Upper (lower) panels: Fraction $f_{\text{fate}}$ of WR–NSs (WR–BHs) producing a BCO (darkest shades), a bound but not merging BCO (intermediate shades), or resulting in binary ionization (lightest shades). The remaining fraction of WR–NSs (WR–BHs) merges after a collision or a CE episode (not plotted). Each set of panels shows the same fraction $f_{\text{fate}}$ for all the metallicities (columns), $\alpha_{\text{CE}}$ (rows), and combinations of CCSN and natal kick models (horizontal entries; see Table 1).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Fraction of BCOs with WR-CO progenitors evolved through different formation channels for different metallicities (rows) and combinations of CCSN and natal kick models explored in this work (Table 1). Here, we show the $\alpha_{CE} = 3$ case. We distinguish four main evolutionary paths based on mass transfer events that occurred before the WR-CO formation: stable mass transfer only, CEs only, at least one stable and one unstable mass transfer, no mass transfer. Hatched bars indicate WR-CO systems that become BCOs experiencing at least one additional mass transfer event (stable or unstable RLO) after the WR-CO phase.</figcaption></figure>
type of mass transfer events occurring prior to the WR-CO formation: evolution through stable mass transfer only, CE phases only, at least one stable mass transfer and one CE event, no mass transfer at all.

Every family of WR-COs to BCOs (WR-BHs to BBHs, WR-BHs to BHNSs, WR-NSs to BHNSs, and WR-NSs to BNSs) exhibits a preferential formation pathway for a fixed metallicity (Fig. 4). CCSNe, natal kick assumptions, and CE efficiency $\alpha_{CE}$ do not impact the relative fraction of possible formation channels as much as metallicity does. Thus, in Fig. 5 we considered the fiducial sets (simulated with delayed CCSN, G20 natal kick models, and $\alpha_{CE} = 3$) to illustrate the main features of the formation channels for these WR-COs at solar ($Z = 0.02$) and sub-solar ($Z = 0.00014$) metallicity.

WR-COs that are BCO progenitors generally form after one or more mass transfer episodes. Stable mass transfer episodes or CE events favour the stripping of the WR star and can reduce the orbital separation, facilitating the production of both a WR-CO binary and a system that merges via GW emission within a Hubble time. Only a few ($\le 8\%$) BBH progenitors at solar metallicity form without exchanging mass. This sub-population of non-interacting binaries has wide orbits ($a \gtrsim 10^2 R_\odot$) and is able to merge only if the natal kick induced by the second CCSN boosts the eccentricity significantly ($e \gtrsim 0.9$). Such lucky kicks are more frequent for the M70 and G20 natal kicks models at $Z = 0.02$.

Stability and duration of the mass transfer phases depend on several factors (e.g., donor mass, mass ratio, stellar radius, stellar phase, orbital separation) and are influenced by metallicity (e.g., Ge et al. 2024). Thus, as shown in Fig. 5, there are many possible formation pathways for WR-COs, with large variability in the number and type of stable mass transfer or CE episodes required. However, we can distinguish two main scenarios, corresponding to WR-COs that are produced in a “direct” or “indirect” way.

In the “direct” scenario, both stars have a similar ZAMS mass ($M_{2,ZAMS} \gtrsim 0.9 M_{1,ZAMS}$). The primary star initiates a CE event (eventually after other mass transfer episodes) when it is already burning helium in its core, while the secondary is either close to ignite it (end of the red giant branch phase) or is also already burning it. If the CE is successfully ejected, the post-CE configuration consists in two helium cores orbiting about each other: a WR-WR binary has formed. Since the two progenitor stars were already close to the end of their lives, shortly after the end of the CE ($\le 3 \times 10^5$ yr) the initially more massive star (primary at ZAMS) becomes a CO, forming a WR-CO. If the post-CE orbit is tight enough ($a \le 3-4 R_\odot$), the most massive WR can overfill its Roche lobe and initiate additional short ($\le 5 \times 10^4$ yr) stable mass transfer events before the first CO formation.

In the “indirect” scenario, the WR-CO does not evolve through the WR-WR configuration. The primary star can still experience one or more mass transfer phases as donor star, but at the time of the first CO formation, the secondary companion is not yet a WR. Thus, in this scenario the secondary star can donate mass to the CO before entering the WR-CO configuration. This scenario is the most common one to produce WR-COs and generally requires the WR progenitor to start a CE episode, in particular at sub-solar metallicity (Figure 5).

Once the WR-CO phase is reached, binaries might undergo wind-fed accretion but they will not experience other mass transfer events to produce BCOs. RLO and CE episodes after the WR-CO phase are common only among WR-COs with light COs, namely WR-NSs or WR-BHs with light BHs, as are the ones produced with the delayed CCSN model ($\sim 3-6 M_\odot$).

Light compact remnants imply long timescales to merge via GW emission (Peters 1964). Thus, to produce a BCO able to merge within a Hubble time, the orbit often has to shrink via mass transfer events: a stable mass transfer phase is generally sufficient for WR-BHs hosting heavier BHs ($\gtrsim 6 M_\odot$), especially at lower metallicities (Fig. 4), while a CE is preferred if an NS is involved.

### 3.3. Only a small fraction of WR-CO systems become BCOs

We showed that the vast majority of BCOs evolved via a WR-CO stage. However, only few WR-COs end their life as BCOs. At sub-solar metallicity ($Z = 0.0014$, $Z = 0.00014$), about $\approx 30\%$ of WR-CO systems produce BCOs; at solar metallicity ($Z = 0.02$, $Z = 0.014$) only $\approx 5\%$ (Figure 6).

WR-COs are not common: in each simulated set, we find that only $\approx 1-5\%$ of the initial binary population evolves as a WR-CO system. In our simulations, the relative abundance of WR-BHs (WR-NSs) in the WR-CO population is not directly correlated with the abundance of BBHs (BNSs) in the BCO population. WR-BHs constitute always $\approx 50-80\%$ of the WR-CO population, with the relative abundance to WR-NSs that is sensible to the CCSN and natal kick models. In contrast, the relative abundance of BBHs in the BCO population is a strong function of metallicity. For instance, BBHs constitute the dominant BCO family at sub-solar metallicity ($\approx 60-80\%$) but are rare at solar metallicity ($\le 10\%$). Similar arguments apply to the relative abundance of WR-NSs and BNSs, even if BNS production is also sensible to $\alpha_{CE}$.

To investigate what determines the fate of WR-COs, we considered the fraction $f_{fate}$ of WR-CO systems with a specific fate. For instance, to calculate the fraction $f_{fate}$ of WR-COs becoming BCOs, we consider the number of BCOs with a WR-CO progenitor (WR-BHs, WR-NSs, or WR-BH plus WR-NSs; for BBHs, BNSs or BHNSs, respectively) and we divide it for the correspondent total WR-CO population (all WR-BHs, WR-NSs or WR-COs, respectively). WR-COs could have various possible fates: the system could prematurely coalesce (as a result of a collision or CE event); the WR could form the second CO after the binary split; the binary could be loose and not merging via GWs; the binary could be tight enough to merge via GW emission within a Hubble time. In Figure 3, we show the fractions $f_{fate}$ calculated for WR-NSs and WR-BHs for all the possible fates and across the parameter space explored in this work (Table 1). In Figure 6, we show median values and 68% credible intervals for the distributions of $f_{fate}$ leading only to BCO production for a fixed pair of ($Z, \alpha_{CE}$) values.

#### 3.3.1. WR-COs and lucky kicks

As shown in Figure 3, most of the WR-NSs that are able to survive in bound orbits have also a final orbital configuration that is optimized to merge via GW emission within a Hubble time. These systems were selected by lucky kicks: natal kicks that are strong enough to result in highly eccentric orbits but not too strong to ionize the binary.

Figure 7 shows that the G20 and M70 models have natal kick magnitude distributions that peak in the optimal velocity regime of lucky kicks for NSs ($v \approx 100$ km s$^{-1}$). Since BHs are heavier than NSs, the optimal regime for lucky kicks of BHs is shifted to higher kick velocities ($v \approx 200$ km s$^{-1}$), where none of the natal kick distributions considered in this work peaks. Therefore, in our simulations, lucky kicks boost only the production efficiency
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Sankey diagram illustrating formation channels of WR–COs that are also BCO progenitors. We show our fiducial sets, simulated assuming a delayed CCSN model, G20 natal kick model, and $\alpha_{CE}=3$; for solar ($Z=0.02$, left) and sub-solar ($Z=0.00014$, right) metallicity. We distinguish “direct” WR–COs from “indirect” WR–COs. “Direct” WR–COs are produced by a WR-WR post-CE binary; “indirect” WR–COs avoid the WR-WR evolution. Unlike “direct” WR–COs, WR progenitors in “indirect” WR–COs have the possibility to start a stable mass transfer (indicated with SMT) or a CE phase between the first CO formation and the WR–CO phase. Few “indirect” systems avoid further mass transfer (indicated with MT) episodes.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Median fraction $f_{\text{fate}}$ of WR-COs that were able to produce BCOs as a function of metallicity. The left-hand (right-hand) plot refers to $\alpha_{\text{CE}} = 1$ ($\alpha_{\text{CE}} = 3$). In each set, we consider all WR-CO systems and we count how many of them form any type of BCOs (gray lines, circle markers), only BBHs (purple lines, lower triangle markers), only BHNSs (blue lines, diamond markers), and only BNSs (yellow lines, upper triangle markers). For each (Z, $\alpha_{\text{CE}}$) pair and BCO type, we consider the distribution of the 12 $f_{\text{fate}}$ fractions of the corresponding simulated sets (combinations of CCSN and natal kick models, see Table 1), and show their median values (scatter points) and 68% credible intervals (error bars).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Distribution of natal kick magnitudes received by all BHs (left) and NSs (right) produced by the initial binary population at Z = 0.02 ($\alpha_{\text{CE}} = 3$, delayed CCSN) for the four natal kick models explored in this work (Table 1).</figcaption></figure>

of the lightest second-born COs, explaining the higher production efficiency of BNSs and BHNSs with respect to BBHs at Z = 0.02 (Figure 6).

We highlight that selection effects due to *lucky kicks* are responsible also for the higher fraction of BBHs without a WR-BH progenitor at Z = 0.02 (Figure 1). In this case, the selection effect is relevant only for the population of lightest BHs, especially for BHs with mass $\le 6$ M⊙ produced by the delayed CCSN model. For decreasing metallicity, the BH mass increases and the effect of *lucky kicks* is negligible.

### 3.3.2. WR-NSs and CE evolution

CE evolution is the most common channel to form binaries with one or two NSs (WR-NSs, BHNSs, BNSs) and is necessary to

produce BNSs at every metallicity (Figure 4). Unstable mass transfer events may completely remove the stellar envelope and facilitate the creation of NSs. However, entering a CE phase may cause a reduction of the semi-major by $\Delta a \approx 2000$ R⊙ and a premature collision. This possibility is enhanced for low $\alpha_{\text{CE}}$ values. The sets evolved with $\alpha_{\text{CE}}=1$ need to extract more orbital energy to expel the CE with respect to sets evolved with $\alpha_{\text{CE}}=3$, often leaving no surviving system.

Collisions driven by CE events reduce the fraction of WR-NSs that can successfully form BHNSs or BNSs. At the lowest metallicity (Z = 0.00014), NSs are produced via binary stripping and the effect of inefficient CE ejection is more evident. At $\alpha_{\text{CE}}=1$, $\gtrsim 90\%$ WR-NSs coalesce prematurely, limiting the production of BNSs (Figure 6) and suppressing the formation of BHNSs from WR-NSs (Figure 1).
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Orbital properties of WR-COs and Cyg X-3 at the beginning (left) or end (right) of the WR-CO configuration. We highlight the properties of WR-COs (triangular markers) and Cyg X-3-like binaries (non-triangular markers) that are BCO progenitors, comparing them to the properties of the whole sample of WR-CO (gray areas in the background). We distinguish binaries ending up as BBHs (purple upward triangles; red circles), as BNSs (gold downward triangles; yellow thin diamonds), or as BHNSs, either forming the BH first (green rightward triangles; thick turquoise diamonds) or the NS (blue leftward triangles; turquoise pentagons). The black rectangular box indicates the region of orbital parameters that we use to select Cyg X-3 candidates (Section 3.4). Black diagonal lines indicate orbital periods of 1 hour, 1 day or 1 year. Here we show the set evolved with $\alpha_{CE} = 3$, $Z = 0.02$, G20 natal kick model and delayed CCSN model. In Appendix B.2 we show other sets at $\alpha_{CE} = 3$ and $Z = 0.02$.</figcaption></figure>

### 3.4. Cyg X-3: A possible BCO progenitor

Cyg X-3 is the only WR-CO candidate known in the MW (Zdziarski et al. 2013; Esposito et al. 2015; Veledina et al. 2024) and has been indicated as a possible BCO progenitor by Belczynski et al. (2013). Here, we investigate this possibility, looking for Cyg X-3-like binaries in our simulations and analyzing their fates. In this work, we define a Cyg X-3-like system as a WR-CO system if its orbital properties satisfy three criteria: orbital period $P_{orb} \in [4.5, 5.1]$ h, WR mass $M_{WR} \in [8-14]$ M$_{\odot}$, and CO mass $M_{rem} \le 10$ M$_{\odot}$ (Singh et al. 2002; Koljonen & Maccarone 2017). We take a conservative approach on the mass and nature of the CO hosted in Cyg X-3 and we allow for both the BH ($M_{rem} = 3-10$ M$_{\odot}$) and NS hypothesis ($M_{rem} < 3$ M$_{\odot}$), given the typical uncertainties in the CO mass determination of COs with WR companions (Laycock et al. 2015; Binder et al. 2021).

Cyg X-3 lies in the Galactic plane at a distance of $d_{CygX-3} = 7.41 \pm 1.13$ kpc from the Sun (McCollough et al. 2016). According to the metallicity gradient of the MW derived from Gaia DR2 (Lemasle et al. 2018), Cyg X-3 is expected to have solar-like metallicity ([Fe/H] $\approx -0.088$). Given the uncertainties in the metallicity determination, we investigated both solar values considered in this work: $Z = 0.02$ and $0.014$.

#### 3.4.1. Cyg X-3: Optimal properties for BCO formation

We find that only $\approx 0.1-0.01\%$ of WR-CO binaries at solar metallicity have a Cyg X-3-like configuration. Typical WR-CO separations at solar metallicity in our simulations are of $a \gtrsim 10^2$ R$_{\odot}$ (Figure 8). Such large orbital separations usually require lucky kicks in the second-CO formation to increase the eccentricity and allow for BCO production. In contrast, Cyg X-3 has a semi-major axis that is at least two order of magnitudes

smaller ($a \approx 3-5$ R$_{\odot}$), increasing the possibility that the system will merge via GW emission within a Hubble time. Figure 8 shows that Cyg X-3 occupies one of the regions of the space of orbital properties that strongly favors BCO production.

About $\approx 20-80\%$ of the Cyg X-3 candidates in our simulations are BCO progenitors. More than $\ge 70\%$ ($\approx 100\%$ for some combinations of the parameter space, Figure 9) of the Cyg X-3-like binaries hosting a BH end their life as BCOs. Even if we consider evolution with the M265 natal kick model, that generates high-velocity kicks typical of Galactic pulsars (Hobbs et al. 2005), we find that more than $\ge 30\%$ of Cyg X-3 candidates hosting a BH become BCOs. In contrast, only $\le 60\%$ of our Cyg X-3-like systems composed of a WR and an NS become BCOs, because most of the systems break during the second CCSN. Thus, we find that Cyg X-3 has a higher probability to become a BCO progenitor if it hosts a BH, in agreement with Belczynski et al. (2013).

#### 3.4.2. Formation efficiency and evolution

Figure 9 shows the formation efficiency $\eta_{CygX-3}$ of Cyg X-3 for all the sets evolved at solar metallicity. For each set, we calculated $\eta_{CygX-3}$ as the ratio between the number of Cyg X-3 candidates and the total simulated mass, correcting for the incomplete initial mass sampling of the primary (Appendix A.2) and assuming an initial binary fraction of 1 for simplicity.

The delayed CCSN model is the only CCSN model explored in this work that produces COs in the $\approx 2-6$ M$_{\odot}$ range. Therefore, the delayed model produces more Cyg X-3-like binaries ($M_{rem} \le 10$ M$_{\odot}$) with respect to the rapid and compactness-based models. Assuming the delayed CCSN model, the formation efficiency of Cyg X-3 can be up to one order of magnitude higher with respect to values obtained in other sets, reaching values of $\eta_{CygX-3} \approx 5.6 \times 10^{-7}$ M$_{\odot}^{-1}$.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Formation efficiency of Cyg X-3 candidates (Singh et al. 2002; Koljonen & Maccarone 2017) for each set of binaries simulated at the two solar metallicity values considered in this work: Z = 0.02 (left) and Z = 0.014 (right). We distinguish Cyg X-3 candidates hosting a BH (upper panels) or an NS (lower panels). We highlight with darker shades the Cyg X-3 candidates that are BCO progenitors.</figcaption></figure>

We find that Cyg X-3-like binaries evolve similarly to WR–COs at the same metallicity (Section 3.3). Cyg X-3 progenitors experience one or two mass-transfer events (with preference for two, one stable and one unstable), but there are no further mass transfer episodes after the system becomes a WR–CO, except for some WR–NSs. At the beginning of the Cyg X-3-like phase, the WR star does not fill more than ≤80% of its Roche lobe: the binary can be wind-fed and the orbital separation can increase, as it is visible in Figure 8. The second compact object forms after less than ≤0.15 Myrs of evolution in the Cyg X-3 configuration.

The G20 and M70 natal kick models select a larger number of Cyg X-3 candidates, especially if they host an NS (Figure 9). These binaries entered in the Cyg X-3 region of the parameter space as a consequence of the natal kick. Most of these systems will break after the second CCSN or coalesce prematurely, because of the additional mass transfer events that WR–NSs may experience.

# 4. Discussion

## 4.1. WR formation

In this work, WR stars can form either via self-stripping, as a result of wind mass loss, or via binary stripping, as a result of mass transfer events. WR formation via self-stripping requires powerful stellar winds. Since wind mass loss rates depend on the mass of the star and its metallicity (Vink et al. 2001, 2011; Nugis & Lamers 2000; Sander et al. 2023), in our simulations self-stripping can produce WRs at solar metallicity (Z = 0.02, 0.014) but not at sub-solar metallicity (Z = 0.0014, 0.00014). For instance, according to the PARSEC tracks adopted in this work, a star evolved in isolation Z = 0.02 becomes a WR ($M_{\text{He-core}} \ge 97.9\% M_{\text{star}}$) if it has an initial mass $M_{\text{ZAMS}} \ge 30 M_{\odot}$ (WR mass of $M_{\text{WR}} \approx 13 M_{\odot}$). The same

$M_{\text{ZAMS}} \approx 30 M_{\odot}$ star, if evolved at Z = 0.0014, will still retain ≈50% of its mass in its hydrogen-envelope at the end of the carbon burning phase.

We expect that changes in the input physics, in the internal mixing properties, and in the choice of the stellar evolution code adopted would modify mass losses, the internal structure, and the core growth of the star, influencing the thresholds for WR formation via single stellar evolution (e.g., Agrawal et al. 2022 and references therein). For instance, Dall’Amico et al. (2025) showed that accretion-induced chemically quasi-homogeneous evolution can enhance the formation of WR stars over red supergiants in secondary stars, increasing the fraction of binaries that are able to form WR–CO systems.

## 4.2. Mass transfer uncertainties

In SEVN, we use critical donor-to-accretor mass ratios $q_c = M_{\text{donor}}/M_{\text{accretor}}$ to determine mass transfer stability, according to the formalism adopted in BSE (Hurley et al. 2002) and other population-synthesis codes. Several works in the literature highlight that this is a simplistic approach to a complicated problem, and has to be revised. For instance, the works by Ge et al. (2010), Ge et al. (2015), Ge et al. (2020), and Ge et al. (2024) indicate that there is no universal $q_c$ threshold for a fixed evolutionary phase of the donor star; rather, their values are strongly influenced by assumptions on the stellar physics, metallicity, and conservation of mass and angular momentum. Willcox et al. (2023) find that prescriptions calculated assuming conservative mass transfer could underestimate the number of mergers produced by unstable mass transfer. Gallegos-Garcia et al. (2021) showed that setting a $q_c$ threshold favors CE evolution as the dominant formation channel for BBHs, while more detailed stellar evolution calculations favour stable mass transfer evolution as the dominant channel. The work of Gallegos-Garcia et al. (2021)
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Left: probability distribution function (PDF) of the orbital periods $P_{WR-BH}$ at the end of the WR-BH configuration for progenitors of BBHs. The dashed line at $P_{WR-BH} = 1$ day highlights the systems that are tight enough to induce efficient tidal spin-up ($P_{WR-BH} \le 1$ day) according to the model by Bavera et al. (2021b). Center: probability distribution function of the dimensionless spin magnitude $\chi_{BH,2}$ of second-born BHs in WR-BHs evolving into BBHs. Spins are calculated with the model of Bavera et al. (2021b) that assumes tidal spin-up in close binaries and efficient angular momentum transport. All the 96 simulated sets (Table 1) exhibit a bimodal distribution. Right: probability distribution function of the effective spin parameter $\chi_{eff}$ obtained in all 96 sets, color-coded according to the assumed metallicity $Z$.</figcaption></figure>

indicates that this approach to mass transfer processes in population synthesis codes could overestimate the merger rate density of BBHs (e.g., Neijssel et al. 2019; van Son et al. 2022; Olejak et al. 2024; Picco et al. 2024).

In this work, we did not put a $q_c$ threshold for MS and HG donor stars: we assumed that they experience only stable mass transfer events, avoiding CEs. This setup corresponds to the fiducial model presented by Iorio et al. (2023). The resulting BCO formation efficiency and channels are compared with the ones obtained assuming, for instance, the $q_c$ setup of Hurley et al. (2002). From Figs. 14 to 20 of Iorio et al. (2023) we can see that forcing mass transfer to be stable for MS and HG donor stars has minor effects on BBH properties, while it increases BNS merger efficiency at solar metallicity (by about one order of magnitude).

### 4.3. Spin distributions in BBHs

#### 4.3.1. Spin magnitude of second-born BH

According to our models, most BBHs emerge from a WR-BH configuration. The first-born BH can be non-spinning ($\chi_{BH,1} \approx 0$) if the progenitor star loses angular momentum efficiently. In close orbits, as are the ones leading to BBHs, the BH can tidal spin up the WR star. After the CCSN, the second-born BH could be highly spinning and become the dominant contribution to the effective and precessing spin parameters ($\chi_{eff}$ and $\chi_p$) that we can measure from GW observations (Abbott et al. 2023b).

Bavera et al. (2021) used MESA (Paxton et al. 2015) to follow the tidal spin up in WR-BH binaries. They assumed efficient angular transport, resulting in non-spinning first-born BH. Bavera et al. (2021) derived a fitting formula to predict the dimensionless spin magnitude of the second-born BH $\chi_{BH,2}$ as a function $f$ of the orbital period of the WR-BH binary $P_{WR-BH}$ and the WR mass, $M_{WR}$, at He or C depletion (Bavera et al. 2021b):

$$ \chi_{BH,2} = \begin{cases} f(P_{WR-BH}, M_{WR}) & P_{WR-BH} \le 1 \text{ day} \\ 0 & \text{otherwise} \end{cases}. \qquad (6) $$

If we apply this relation to the population of WR-BHs in our simulations, we find a bimodal distribution for the spins

of second-born BHs (Figure 10): they are either highly spinning ($\chi_2 \approx 1$) or no spinning ($\chi_2 \approx 0$). In the model by Bavera et al. (2021b), the spin magnitude is determined by tidal spin-up. Thus, it is strongly influenced by the orbital period of the WR-BH system.

Figure 10 (left panel) shows that most WR-BHs at solar metallicity ($Z = 0.02, 0.014$) have orbital periods too large ($P_{WR-BH} > 1$ day) to be spun up by tides according to Equation (6). Only BBHs produced with the G20 natal kick model (Giacobbo & Mapelli 2020), the model that generates the lowest kick magnitudes for BHs (Figure 7), exhibit a significant fraction of binaries with $P_{WR-BH} \le 1$ day.

At sub-solar metallicity ($Z = 0.0014, 0.00014$), the lower wind efficiency allows WR production only via mass transfer interactions and binary stripping. Hence, WR-BHs at these metallicities have shorter periods resulting in an efficient spin up of the WR that will produce the second-born BH.

#### 4.3.2. Comparison with spins from GW observations

If one of the two BHs is non-spinning ($\chi_1 \approx 0$), the effective spin parameter $\chi_{eff}$ is extremely sensible to the mass ratio $q = M_{BH,L}/M_{BH,H}$ of the lighter-to-heavier BH ($M_{BH,H} \ge M_{BH,L}$). By construction, this<fnref n="2" /> mass ratio is $q < 1$ and can be used to calculate $\chi_{eff}$ with an explicit dependence on the possibility that the second-born BH $M_{BH,2}$ is the heaviest or the lightest one:

$$ \chi_{eff} \approx \frac{m_2 \chi_2 \cos \theta_2}{m_1 + m_2} \approx \begin{cases} \frac{q}{1+q} \chi_2 \cos \theta_2 & M_{BH,2} \le M_{BH,1} \\ \frac{1}{1+q} \chi_2 \cos \theta_2 & M_{BH,2} \ge M_{BH,1} \end{cases}. \qquad (7) $$

GW observations suggest that most BHs have relatively low spins ($\chi \approx 0.2$) and that it might exist a $q - \chi_{eff}$ correlation that favors higher $\chi_{eff}$ values for more asymmetric binaries (Callister et al. 2021; Abbott et al. 2023b). In the case of efficient angular momentum transport and Eddington-limited accretion, the first-born BH could be non-spinning and $\chi_{eff}$ would be determined mainly by the properties of the second-born BH. Assuming also

<footnote n="2">² Other mass ratio definitions could be larger than 1, for instance the donor-to-accretor mass ratio $q = M_d/M_a$ used to determine mass transfer stability.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Mass ratio $q = M_{BH,L}/M_{BH,H}$ of the lighter-to-heavier BH and effective spin parameter $\chi_{\text{eff}}$ for all BBHs produced in the 96 sets simulated in this work (Section 2). Colored thick lines indicate the analytic relations (Equation (7)) that delimit the allowed parameter space for maximally spinning ($\chi_2 = 1$) second-born BH, for spins aligned ($\cos\theta_2 = 1$) or anti-aligned ($\cos\theta_2 = -1$) with the binary orbital angular momentum. Black thin lines show the analytic relation expected for low spinning ($\chi_2 \sim 0.1-0.5$) and aligned ($\cos\theta_2 = 1$) second-born BH. The analytic relations assume that the second-born BH is the heaviest one ($M_{BH,2} = M_{BH,H}$), unless different specifications.</figcaption></figure>

spin aligned to the orbital momentum ($\cos\theta_2 = 1$), the analytical approximation of $\chi_{\text{eff}}$ (Equation (7)) indicates an intrinsic $q-\chi_{\text{eff}}$ (anti-) correlation if the second-born BH the (most) least massive one. The magnitude of the second-born BH spin influences the shape of the (anti-) correlation in the $q-\chi_{\text{eff}}$, as we show in Figure 11.

In the case of maximally spinning second-born BH ($\chi_2 = 1$), the approximated formula for $\chi_{\text{eff}}$ divides the $q-\chi_{\text{eff}}$ plane in three regions: an oval-shaped central region, an outer wing, and a more external one. The regions are symmetric with respect to the $\chi_{\text{eff}} = 0$ center, indicating (anti-) aligned spins for (negative) positive $\chi_{\text{eff}}$ values. In Figure 11 we superimpose also the properties of all the BBHs simulated in this work (96 sets, see Table 1) and spun up via Eq. (6). We find no BBH in the more external region, thus, no systems with high $\chi_{\text{eff}} \approx 1$ and similar masses ($q \approx 1$). Most of our BHs are non-spinning ($\chi_{\text{eff}} = 0$) or have spin aligned with the orbital momentum ($\cos\theta_2 \approx 1$)$^3$. BBHs with second-born BH that is lighter than the first-born one ($M_{BH,2} = M_{BH,L}$) are confined in the central oval-shaped region, with $|q| < 0.5$. BBHs with heavier second-born BHs ($M_{BH,2} = M_{BH,H}$) distribute both in the central oval-shaped region and in the outer wings. In this scenario, the wings centered at $|q| = 0.5$ identify a region that can be populated only by BBHs evolved

through a mass-ratio-reversal evolution. The distribution of our BBHs in the $q-\chi_{\text{eff}}$ plane is in agreement with other works in the literature that studied the mass-ratio-reversal scenario and considered tidal spin up (e.g., Olejak & Belczynski 2021; Broekgaarden et al. 2022b; Olejak et al. 2024; Banerjee & Olejak 2024).

## 5. Summary

We investigated the role of WR-COs as the possible progenitors of BCOs. We used the population synthesis code SEVN to account for uncertainties in theoretical models, evolving the same initial binary population through 96 combinations of metallicity, CE efficiency, CCSN, and natal kick models.

We found that most ($\ge 83\%$) of the simulated BCOs evolved as WR-CO systems in the considered parameter space. More than $\ge 99\%$ ($\sim 83-95\%$) of BCOs at $Z \ge 0.0014$ ($Z = 0.00014$) have a WR-CO progenitor. The fraction of BCOs that derive from WR-CO systems is mostly affected by mass transfer (when it stops before complete envelope stripping we do not form a WR at low $Z$), metallicity (which influences stellar winds, thus the role of self and binary stripping in WR production), and the combined effect of natal kicks and CCSN models.

Only $\approx 1-5\%$ of the simulated massive binaries passes through the WR-CO phase. Despite their key role as BCO progenitors, only a small fraction of WR-COs ($\approx 5-30\%$, depending on the metallicity) becomes a BCO. When formed, WR-COs constitute the last evolutionary configuration prior to BCO formation. These progenitors experience all the RLOs and CE events before the WR-CO phase. Only WR-NSs require at least one additional CE event in order to produce BCOs.

We characterized the possible fate of Cyg X-3, the only Galactic WR-CO candidate, finding that it will likely evolve into a BCO if it is a WR-BH system ($\approx 70-100\%$). If Cyg X-3 is a WR-NS binary, the BCO fate is disfavoured ($<40-60\%$), in agreement with Belczynski et al. (2013). Cyg X-3 possible formation channels are representative of WR-COs at $Z = 0.02, 0.014$, suggesting that the observation and characterization of similar WR-CO systems could help us to constrain the models to interpret the BCO formation history.

*Acknowledgements.* EK and MM acknowledge support from the PRIN grant METE under contract No. 2020KB33TP. MM, GC, GI, and MDA acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. EK and MM also acknowledge financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) STRUCTURES. GI acknowledges financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.4, - Call for tender No. 3138 of 18/12/2021 of Italian Ministry of University and Research funded by the European Union - NextGenerationEU. GC acknowledges financial support from the Agence Nationale de la Recherche grant POPSYCLE number ANR-19-CE31-0022. Numerical calculations have been made possible through a CINECA-INFN (TEONGRAV) agreement, providing access to resources on Leonardo at CINECA. The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grants INST 35/1597-1 FUGG and INST 35/1503-1 FUGG.

## References

Abac, A. G., Abbott, R., Abouelfettouh, I., et al. 2024, ApJ, 970, L34

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 011001

Abbott, B., Abbott, R., Abbott, T., et al. 2020, ApJ, 892, L3

Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 915, L5

Abbott, R., Abbott, T., Acernese, F., et al. 2023a, Phys. Rev. X, 13, 041039

Abbott, R., Abbott, T., Acernese, F., et al. 2023b, Phys. Rev. X, 13, 011048

Agrawal, P., Szécsi, D., Stevenson, S., Eldridge, J. J., & Hurley, J. 2022, MNRAS, 512, 5717

<footnote>$^3$ The spin orientation $\cos(\theta_2)$ was calculated following Equation (1.16) of Mapelli (2021), assuming that only the natal kicks can change the spin orientation.</footnote>
Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2019, MNRAS, 489, 3116

Banerjee, S., & Olejak, A. 2024, A&A, submitted [arXiv:2411.15112]

Bartzakos, P., Moffat, A., & Niemela, V. 2001, MNRAS, 324, 18

Bavera, S. S., Fragos, T., Zevin, M., et al. 2021, A&A, 647, A153

Bavera, S. S., Zevin, M., & Fragos, T. 2021b, RNAAS, 5, 127

Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223

Belczynski, K., Bulik, T., & Fryer, C. L. 2012, ApJ, submitted, [arXiv:1208.2422]

Belczynski, K., Bulik, T., Mandel, I., et al. 2013, ApJ, 764, 96

Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104

Bethe, H. A., & Brown, G. E. 1998, ApJ, 506, 780

Binder, B. A., Sy, J. M., Eracleous, M., et al. 2021, ApJ, 910, 74

Boesky, A. P., Broekgaarden, F. S., & Berger, E. 2024, ApJ, 976, 23

Böhm-Vitense, E. 1958, ZAp, 46, 108

Bressan, A., Chiosi, C., & Bertelli, G. 1981, A&A, 102, 25

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Broekgaarden, F. S., Berger, E., Stevenson, S., et al. 2022a, MNRAS, 516, 5737

Broekgaarden, F. S., Stevenson, S., & Thrane, E. 2022b, ApJ, 938, 45

Callister, T. A., Haster, C.-J., Ng, K. K. Y., Vitale, S., & Farr, W. M. 2021, ApJ, 922, L5

Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068

Chieffi, A., & Limongi, M. 2020, ApJ, 890, 43

Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J., & Verbunt, F. W. M. 2014, A&A, 563, A83

Costa, G., Girardi, L., Bressan, A., et al. 2019, MNRAS, 485, 4641

Costa, G., Bressan, A., Mapelli, M., et al. 2021, MNRAS, 501, 4514

Costa, G., Shepherd, K. G., Bressan, A., et al. 2025, A&A, 694, A193

Crowther, P. A. 2007, Annu. Rev. Astron. Astrophys., 45, 177

Dall'Amico, M., Mapelli, M., Iorio, G., et al. 2025, A&A, 695, A221

Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52

Esposito, P., Israel, G. L., Milisavljevic, D., et al. 2015, MNRAS, 452, 1112

Fernández, R., Quataert, E., Kashiyama, K., & Coughlin, E. R. 2018, MNRAS, 476, 2366

Foellmi, C., Moffat, A., & Guerrero, M. 2003a, MNRAS, 338, 360

Foellmi, C., Moffat, A., & Guerrero, M. 2003b, MNRAS, 338, 1025

Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91

Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, ApJ, 922, 110

Ge, H., Hjellming, M. S., Webbink, R. F., Chen, X., & Han, Z. 2010, ApJ, 717, 724

Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2015, ApJ, 812, 40

Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2020, ApJ, 899, 132

Ge, H., Tout, C. A., Chen, X., et al. 2024, ApJ, 975, 254

Giacobbo, N., & Mapelli, M. 2020, ApJ, 891, 141

Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011

Götberg, Y., Drout, M. R., Ji, A. P., et al. 2023, ApJ, 959, 125

Gräfener, G., & Hamann, W. R. 2008, A&A, 482, 945

Gräfener, G., & Vink, J. S. 2013, A&A, 560, A6

Gräfener, G., Owocki, S. P., Grassitelli, L., & Langer, N. 2017, A&A, 608, A34

Hainich, R., Rühling, U., Todt, H., et al. 2014, A&A, 565, A27

Hamann, W. R., & Gräfener, G. 2004, A&A, 427, 697

Han, Z., Podsiadlowski, P., & Eggleton, P. P. 1994, MNRAS, 270, 121

Han, Z., Podsiadlowski, P., & Eggleton, P. P. 1995, MNRAS, 272, 800

Herwig, F., Bloecker, T., Schoenberner, D., & El Eid, M. 1997, A&A, 324, L81

Hjellming, M. S., & Webbink, R. F. 1987, ApJ, 318, 794

Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974

Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897

Iorio, G., Mapelli, M., Costa, G., et al. 2023, MNRAS, 524, 426

Ivanova, N., Justham, S., Chen, X., et al. 2013, A&A Rev., 21, 59

Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, A&A, 645, A54

Koljonen, K. I. I., & Maccarone, T. J. 2017, MNRAS, 472, 2181

Kroupa, P. 2001, MNRAS, 322, 231

Laplace, E., Göteborg, Y., De Mink, S., Justham, S., & Farmer, R. 2020, A&A, 637, A6

Laycock, S. G. T., Maccarone, T. J., & Christodoulou, D. M. 2015, MNRAS, 452, L31

Lemasle, B., Hajdu, G., Kovtyukh, V., et al. 2018, A&A, 618, A160

Limongi, M., & Chieffi, A. 2010, in Journal of Physics Conference Series, 202, 012002

Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13

Mandel, I., & Broekgaarden, F. S. 2022, Liv. Rev. Relativ., 25, 1

Mapelli, M. 2021, in Handbook of Gravitational Wave Astronomy, 16

Mapelli, M., Spera, M., Montanari, E., et al. 2020, ApJ, 888, 76

McCollough, M., Corrales, L., & Dunham, M. 2016, ApJ, 830, L36

Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15

Neijssel, C. J., Vigna-Gómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740

Nguyen, C. T., Costa, G., Girardi, L., et al. 2022, A&A, 665, A126

Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227

Nugis, T., & Lamers, H. J. G. L. M. 2002, A&A, 389, 162

O'Connor, E., & Ott, C. D. 2011, ApJ, 730, 70

Olejak, A., & Belczynski, K. 2021, ApJ, 921, L2

Olejak, A., Klencki, J., Xu, X.-T., et al. 2024, A&A, 689, A305

Özel, F., & Freire, P. 2016, ARA&A, 54, 401

Paczyński, B., & Ziółkowski, J. 1968, Acta Astron., 18, 255

Patton, R. A., & Sukhbold, T. 2020, MNRAS, 499, 2803

Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15

Peters, P. C. 1964, Phys. Rev., 136, B1224

Picco, A., Marchant, P., Sana, H., & Nelemans, G. 2024, A&A, 681, A31

Röpke, F. K., & De Marco, O. 2023, Liv. Rev. Computat. Astrophys., 9, 2

Sabrahit, G. N., Vink, J. S., Higgins, E. R., & Sander, A. A. C. 2022, MNRAS, 514, 3736

Sana, H., De Mink, S., de Koter, A., et al. 2012, Science, 337, 444

Sander, A. A. C., Hamann, W.-R., Todt, H., et al. 2019, A&A, 621, A92

Sander, A. A. C., Lefever, R. R., Poniatowski, L. G., et al. 2023, A&A, 670, A83

Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021, MNRAS, 502, 4877

Schneider, F. R. N., Podsiadlowski, P., & Müller, B. 2021, A&A, 645, A5

Schneider, F. R. N., Podsiadlowski, P., & Laplace, E. 2023, ApJ, 950, L9

Schnurr, O., Casoli, J., Chené, A.-N., Moffat, A., & St-Louis, N. 2008, MNRAS, 389, L38

Schwarzschild, M. 1958, Structure and evolution of the stars (Princeton University Press)

Sen, K., Xu, X. T., Langer, N., et al. 2021, A&A, 652, A138

Shenar, T., Gilkis, A., Vink, J. S., Sana, H., & Sander, A. A. C. 2020, A&A, 634, A79

Singh, N., Naik, S., Paul, B., et al. 2002, A&A, 392, 161

Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739

Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086

Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889

Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93

Van Der Hucht, K. A. 2001, New Astron. Rev., 45, 135

van Son, L. A. C., de Mink, S. E., Renzo, M., et al. 2022, ApJ, 940, 184

Veledina, A., Muleri, F., Poutanen, J., et al. 2024, Nat. Astron., 8, 1031

Verbunt, F., Igoshev, A., & Cator, E. 2017, A&A, 608, A57

Vigna-Gómez, A., MacLeod, M., Neijssel, C. J., et al. 2020, PASA, 37, e038

Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295

Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574

Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132

Wadekar, D., Venumadhav, T., Mehta, A. K., et al. 2024, Phys. Rev. D, 110, 084035

Webbink, R. F. 1984, ApJ, 277, 355

Willcox, R., MacLeod, M., Mandel, I., & Hirai, R. 2023, ApJ, 958, 138

Zdziarski, A. A., Mikolajewska, J., & Belczynski, K. 2013, MNRAS, 429, L104
# Appendix A: Models and assumptions

## A.1. Fallback in the Mfb natal kick model

The Mfb natal kick model assigns the kick magnitudes drawing a random number $f_σ$ from a Maxwellian distribution ($σ$ = 265 km s⁻¹ rms) and modulating it with the fraction $f_{fb,CCSN}$ of mass falling back onto the proto-NS:

$$v_{\text{kick}} = f_{\sigma=265} (1 - f_{\text{fb,CCSN}}). \qquad (\text{A.1})$$

Here, we choose the fallback fraction $f_{fb,CCSN}$ consistently with the CCSN models (Section 2.4).

In the rapid and delayed case, we set the fallback fraction equal to the fallback parameter $f_{fb,CCSN} = f_{fb}$ of (Fryer et al. 2012). In the compactness-based case, we distinguish successful and failed supernovae. For BHs formed via direct collapse, no material is ejected ($f_{fb,CCSN} = 1$) and no kick is imparted. For NSs, we calculate a fallback fraction in analogy with Fryer et al. (2012): we assume that the remnant mass ($M_{\text{rem}}$) is determined by the fraction $f_{fb,CCSN}$ of the pre-SN mass ($M_{\text{pre-SN}}$) that a proto-CO of $M_{\text{proto}} = 1.1 M_{\odot}$ is able to accrete, where

$$f_{\text{fb,CCSN}} = \frac{\max(M_{\text{rem}} - M_{\text{proto}}, 0)}{\max(M_{\text{pre-SN}} - M_{\text{proto}}, 0)}. \qquad (\text{A.2})$$

## A.2. Initial conditions for binary population

The initial primary masses $M_1$ were drawn from the Kroupa (2001) initial mass function (IMF) and the initial mass ratios $q$ from Sana et al. (2012):

$$\xi(M_1) \propto M_1^{-2.3} \qquad M_1 \in [5, 150] M_{\odot}, \qquad (\text{A.3})$$

$$\xi(q) \propto q^{-0.1} \qquad q \in [0.1, 1], \qquad (\text{A.4})$$

where $q = M_2/M_1$ and $M_1 \ge M_2$. Combining $\xi(q)$ and $\xi(M_1)$, we obtained the distribution for the initial secondary masses $\xi(M_2) = \xi(M_1)\xi(q)$. We choose the lower cuts at $M_{1,\text{min}} = 5 M_{\odot}$ and $q_{\text{min}} = 0.1$ because we focus on BH and NS progenitors.

We generated a total mass of $M_{\text{IC}} \approx 1.08 \times 10^7 M_{\odot}$. Accounting for incomplete mass sampling of the primary from the Kroupa IMF requires a correction factor of $f_{\text{IMF}} = 0.289$, indicating a parent population of $M_{\text{tot}} \approx 3.7 \times 10^8 M_{\odot}$.

We generated the initial orbital periods $P$ from the $P = \log(P/\text{days})$ distribution by Sana et al. (2012):

$$\xi(P) \propto P^{-0.55} \qquad P \in [0.30, 5.5], \qquad (\text{A.5})$$

setting $P = 2$ days as lower limit (we assumed that binaries with shorter period already circularized following Moe & Di Stefano 2017). Moe & Di Stefano (2017) showed that there is a complex correlation between the mass ratio $q$, period $P$ and eccentricity $e$ of early-type binaries. For simplicity, we assumed independent distributions for $q$ and $P$ and we adopted the Moe & Di Stefano (2017) prescription for $e$:

$$\xi(e(P)) \propto 1 - (P/\text{days})^{-2/3} \qquad P \ge 2 \text{ d.} \qquad (\text{A.6})$$

# Appendix B: Additional plots

## B.1. Formation channels and CE efficiency

In Figure B.1 we show the formation channels of BCOs for $\alpha_{\text{CE}} = 1$. The most notable difference with respect to the $\alpha_{\text{CE}} = 3$ case (Figure 4, Section 3.2) occurs for WR-NS formation and final fate. WR-NSs require one or more CE events to form, and at least an additional one to produce a BCO. At $Z = 0.00014$, the role of CE evolution in BNS creation is enhanced at $\alpha_{\text{CE}} = 1$, with $\approx 60-90\%$ of BNSs exchanging mass only via CE events. At these metallicity, NSs may form only after efficient envelope stripping, as it is the one induced by CE evolution. Thus, evolution through stable RLOs is disfavoured. However, a low CE efficiency determines inefficient CE ejection and can induce premature coalescence of the binary members. For instance, at $Z = 0.00014$ and $\alpha_{\text{CE}} = 1$ no WR-NS was able to survive the CE evolution and produce BHNSs, unlike the $\alpha_{\text{CE}} = 3$ case.

## B.2. Orbital properties of Cyg X-3

In Section 3.4, we showed that Cyg X-3 is a possible BCO progenitor, with formation pathways representative of WR-COs evolved at $Z = 0.02, 0.014$.

Figure B.2 shows the orbital properties of WR-COs and Cyg X-3 in a set evolved with $\alpha_{\text{CE}} = 3$, $Z = 0.02$, the G20 natal kick model and the compactness-based CCSN model. This binary population was evolved under the same assumptions of Figure 8, except for the CCSN model. Despite the same $\alpha_{\text{CE}}$ and $Z$ assumptions, the elevated probability of receiving a lucky kick through the G20 model couples with the different CCSN model (compactness-based in Figure B.2, delayed in Figure 8) and determines a different distribution of WR-COs in the space of orbital parameters.

For further comparison, we show in Figure B.3 how the occupation of space of orbital parameters changes if we assumed the kick model to be the Mfb one. In this plane, changing the initial models modifies the number of WR-COs with Cyg X-3-like properties. Nevertheless, binaries similar to Cyg X-3 remain favoured as BCO progenitors, suggesting that that further observations of WR-COs with orbital properties similar to Cyg X-3 may constrain the properties of BCO progenitors.
<figure><img src="image_13.png" /><figcaption>**Fig. B.1.** Fraction of BCOs with WR-CO progenitors evolved through different formation channels for different metallicities (rows) and combinations of CCSN and natal kick models explored in this work (Table 1). Here we show the $\alpha_{CE} = 1$ case. We distinguish four main evolutionary paths on the basis of mass transfer events occurred before the WR-CO formation: stable mass transfer only, CEs only, at least one stable and one unstable mass transfer, no mass transfer. Hatched bars indicate WR-CO systems that become BCOs experiencing at least one additional mass transfer event (stable or unstable RLO) after the WR-CO phase.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. B.2.** Orbital properties of WR-COs and Cyg X-3 at the beginning (left) and end (right) of the WR-CO phase. We highlight the properties of WR-COs (triangular markers) and Cyg X-3-like binaries (non-triangular markers) that are BCO progenitors, comparing them to the properties of the whole sample of WR-COs (gray areas in the background). We distinguish binaries ending up as BBHs (purple upward triangles; red circles), as BNSs (gold downward triangles; yellow thin diamonds), or as BHNSs, either forming the BH (green rightward triangles; thick turquoise diamonds) or the NS first (blue leftward triangles; turquoise pentagons). The black rectangular box indicates the region of orbital parameters that we use to select Cyg X-3 candidates (Section 3.4). Black diagonal lines indicate orbital periods of 1 hour, 1 day or 1 year. Here we show the set evolved with $\alpha_{CE} = 3$, $Z = 0.02$, G20 natal kick model and compactness CCSN model. In Figure 8 we show the same set, but with the delayed CCSN model.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. B.3.** Same as Figure B.2, but here we show the set evolved with $\alpha_{CE} = 3$, $Z = 0.02$, Mfb natal kick model, and delayed CCSN model. In Figure 8 we show the same set, but with the G20 natal kick model.</figcaption></figure>