© D. Pauli et al. 2022

1. Introduction

Wolf-Rayet (WR) stars are helium-rich massive stars with strong, optically thick winds which lead to emission-line dominated spectra (Crowther 2007). Their high surface temperatures place many of them to the left of the zero-age main sequence in the HR diagram (Abbott & Conti 1987; Schmutz et al. 1989; Koesterke & Hamann 1995; Hainich et al. 2014), allowing them to emit copious amounts of ionizing radiation (Walborn et al. 2004), which may strongly influence star-forming regions in galaxies (Ramachandran et al. 2018; Crowther 2019), and the whole appearance of star-forming so-called WR galaxies (Conti 1991; Brinchmann et al. 2008).

The strong winds and ionizing photons of WR stars drive shocks through the interstellar medium and sweep up interstellar bubbles (Garcia-Segura et al. 1996a,b; Weaver et al. 1977). Their ejection of hydrogen- or helium-burning products (Maeder 1983; Dray et al. 2003) and dust (Cherchneff et al. 2000) drives the chemical evolution of galaxies. Finally, WR stars may also give rise to hydrogen-free (Dessart et al. 2011; Groh et al. 2013) and superluminous supernovae (Inserra et al. 2013; Dessart et al. 2020; Aguilera-Dena et al. 2020, 2022b), to long-duration gamma-ray bursts (Woosley 1993; Yoon & Langer 2005), and black hole binaries (Remillard & McClintock 2006; Langer et al. 2020).

Despite their importance, the formation process of WR stars is still poorly constrained. There are two main ways for massive stars to have their envelope stripped. In the single star scenario, the star develops strong winds that remove its hydrogen-rich envelope on a timescale comparable to its life time (De Loore et al. 1978; Maeder & Meynet 1987; Langer 1987). Sufficiently strong winds are only achieved by the most massive stars. Due to the metallicity dependence of hot star winds (Mokiem et al. 2007), it is unknown whether self-stripping of single stars is possible at low metallicity. Alternatively, stars of any mass and metallicity may lose most of their hydrogen envelope due to interaction with a close binary companion (Paczyński 1967; Vanbeveren & Conti 1980; Wellstein & Langer 1999). A third, but not dominant channel of forming WR stars would be through internal mixing, which requires quasi-chemically homogeneous evolution (Maeder 1987; Yoon et al. 2006; Woosley & Heger 2006; Marchant et al. 2016; Aguilera-Dena et al. 2018; Hastings et al. 2020). It is suggested only to be important at low metallicity, high mass, and extreme rotation.

It is unclear as of yet which of the two envelope stripping channels, the single star or the binary channel, is more important for forming WR stars. The main reason for this is that the mass-loss rate of cool massive stars are not yet well understood. In particular, the role of the Luminous Blue Variables stage (Weis & Bomans 2020) in removing the H-rich envelope in stars near the Eddington-limit (Langer et al. 1994; Smith & Owocki 2006; Grassitelli et al. 2021) and the metallicity dependence of LBV mass loss (Kalari et al. 2018), are uncertain. Similarly, the mass-loss rates of red supergiants are weakly constrained (Mauron & Josselin 2011) as well as, again, the impact of metallicity (Kee et al. 2021).

The single star path for forming WR stars has been explored extensively (Maeder & Meynet 2012; Langer 2012, and references therein). Meynet & Maeder (2003, 2005) found that single star models which include strong rotationally induced internal mixing can explain the major properties of the WR populations in the Galaxy and the Magellanic Clouds. However, the magnitude of rotational mixing has been questioned since then, based on the spectroscopic analysis of large samples of O and early B stars (Hunter et al. 2008; Grin et al. 2017; Markova et al. 2018), and the realization that the majority of stars massive enough to potentially form WR stars rotate relatively slowly (Bestenlehner et al. 2014; Sabín-Sanjulián et al. 2017; Ramírez-Agudelo et al. 2017).

Over the last decade, several works have argued that binary interactions are crucial for our understanding of stellar evolution and that their importance cannot be neglected, as its stellar products – compared to single star models – contribute to different stellar populations (e.g., de Mink et al. 2013; Schneider et al. 2015; Wang et al. 2020). Nonetheless, it is still unclear to which extend binaries contribute to these populations (e.g., Shenar et al. 2020; Massey et al. 2021). Recent observations of massive stars show that their evolution is dominated by binary interaction (Sana et al. 2012, 2014; Moe & Di Stefano 2017), meaning that at least half of all massive stars are born with a close companion such that mass transfer is unavoidable. The immediate implication is that binarity plays a crucial role for WR star formation. Fortunately, the uncertainties in the process of the formation of WR stars through binary systems are much smaller than those for single stars. The result of binary induced mass stripping is well defined; the donor in mass transferring binary systems contract and end the mass stripping once most of their hydrogen envelope is gone. Even though binary evolution in general needs to consider a larger number of initial parameters and contains additional physics uncertainties (Eldridge et al. 2008), the predictions of the properties of WR star models formed through binary interaction may be more credible than comparable single star predictions, since the uncertain wind mass-loss rates for cool stars, or their metallicity dependence, hardly matters.

Here, we explore WR star formation in single stars and close binaries in order to bring clarity to the question on how the individual channels contribute to an entire WR population (e.g., Götberg et al. 2018; Shenar et al. 2020). For this purpose, we use a dense grid of detailed binary evolution models. The computation of a complementary grid of single star models using the same input physics allows us to address the relative importance of both WR star formation channels. The ideal testing ground for our models turns out to be the Large Magellanic Cloud, because with more than 100 WR stars, it provides a statistically significant WR sample (Breysacher et al. 1999; Bartzakos et al. 2001a; Hainich et al. 2014). At the same time, the LMC WR sample is thought to be essentially complete (Neugent et al. 2018) and likely originates from a more chemically homogeneous O star population than the Galactic WR stars.

Our paper is organized as follows. In the next section, we describe the set-up and choice of physics parameters for our single star and binary evolution models and provide examples of evolutionary tracks in the HR diagram for both. In Sect. 3, we outline our procedure to obtain a synthetic LMC WR star population from our stellar evolution models. Sect. 5 provides a discussion of the obtained synthetic populations of WR stars, separating between WN stars with and without hydrogen, as well as WC stars and compares the result with the observed WR populations. In Sect. 5.3, we discuss the properties of our our computed WR binary systems and compare them with the properties of observed WR binaries. Section 6 gives a discussion of the uncertainties involved in our analysis, before we provide our final conclusions in Sect. 7.

2. Methods

In this section, we describe the methods employed to produce and analyze a grid of detailed stellar models. In Sect. 2.1 we provide a brief overview of the chosen input parameters of our single and binary evolutionary models. In Sect. 2.2 we introduce a criterion based on the optical depth to differentiate stripped-envelope (He-stars) and WR stars. Finally, in Sect. 2.3 we describe the assumptions made to create a synthetic WR stellar population.

2.1. Input parameters for our stellar evolution models

We used version 10 398 of the Modules for Experiments in Stellar Astrophysics (MESA) code (e.g., Paxton et al. 2011, 2013, 2015, 2018, 2019) to calculate a grid of detailed binary stellar evolution models at LMC metallicity. The initial chemical composition of our models by mass fraction is X_H = 0.7383, X_He = 0.2569, and Z = 0.0048, where Z includes all elements heavier than helium. Solar-scaled abundances for the heavy elements in metal poor dwarf galaxies, like the LMC, are commonly used in stellar evolutionary models (i.e., Eggenberger et al. 2021; Eldridge et al. 2017). However, dwarf galaxies have different metallicity distributions than galaxies like the Milky Way (i.e., Chruslinska & Nelemans 2019). Therefore, we employed for our models a tailored initial chemical compositions, similar as in the work Brott et al. (2011). For the initial C, N, and O mixtures we used the abundances as they are observed from H II regions (Kurt & Dufour 1998), namely log(C/H)+12 = 7.75, log(N/H)+12 = 6.90, and log(O/H)+12 = 7.35. For Mg and Si we used log(Mg/H)+12 = 7.05, and log(Si/H)+12 = 7.20 based on observations from B-type stars in the LMC (Trundle et al. 2007; Hunter et al. 2007, 2009). As initial iron abundance we used log(Fe/H)+12 = 7.05 (see Brott et al. 2011) which is in agreement with the observed abundances in the LMC (e.g., Ferraro et al. 2006). For the remaining metals we adopted the initial chemical abundances of Asplund et al. (2009) reduced by 0.4 dex, accounting for the lower metallicity in the LMC.

The resulting model grid covers initial primary masses in the range of M_1, i ≃ 28 M_⊙ − 89 M_⊙ in steps of Δlog(M_1, i/M_⊙) = 0.05. Stellar models with lower initial mass are not expected to contribute to the WR population (cf. Fig. A.2) and thus are not included in this work. The initial orbital periods cover the range of P_i ≃ 1.4 d − 10 000 d in steps of Δlog(P_i/d) = 0.05 and the initial mass ratios cover the range of q_i = 0.25 − 0.95 in steps of Δq_i = 0.05. Our grid contains a total of 10 780 models. We assumed that our binary models are initially tidally synchronized at the zero-age main-sequence (ZAMS). As the models with the highest initial orbital periods (i.e., P_i ≃ 10 000 d) never interact during their evolution, they serve as our grid of quasi nonrotating single star models. While initial tidal synchronization is not realistic for wide binary systems, it ensures that both components evolve essentially as our single star models before the onset of mass transfer via Roche-lobe overflow (RLOF).

Stellar wind mass-loss is included as in Brott et al. (2011) with slight modifications, described below. The mass-loss rates of Vink et al. (2001) are used for hydrogen-rich main-sequence stars. For temperatures below the bi-stability jump (Vink et al. 2001, their Eqs. (14) and (15)), the maximum value of either Vink et al. (2001) or Nieuwenhuijzen & de Jager (1990) wind prescriptions is adopted. To model mass loss during the WR phase, we use the mass-loss rates of Nugis & Lamers (2000) with a clumping factor of D = 2.9 for the hydrogen-rich WN phase where the surface hydrogen abundance is below X_H < 0.4, and for the hydrogen-free WN and WC phase we used the mass-loss recipe of Yoon (2017) with a clumping factor of D = 4, meaning that this leads to a variation of the original mass-loss rate as Ṁ ∝ D^−1/2. We note, that some of the mass-loss recipes are originally calibrated for D = 10 and thus when using a clumping factor of D = 4 or D = 2.9 the mass-loss rates are increased by a factor of 1.58 or 1.86, respectively. For surface hydrogen abundances between 0.7 and 0.4, the wind mass-loss rate is linearly interpolated between those of Vink et al. (2001) and Nugis & Lamers (2000). For all used mass-loss rates a solar baseline of Z_⊙ = 0.017 is used.

Rotational mixing is modeled as a diffusive process including the effects of dynamical and secular shear instabilities, the Goldreich-Schubert-Fricke instability, and Eddington-Sweet circulations (Heger et al. 2000). The efficiency of rotational mixing is calibrated following Brott et al. (2011), with efficiency factors f_c = 1/30 and f_μ = 0.1. In addition, internal angular momentum transport via magnetic fields is included according to Spruit (2002).

Convection is modeled using the Ledoux criterion and the standard mixing length theory (Böhm-Vitense 1958) with a mixing length parameter of α_mlt = l/H_P = 1.5. As such, the effect of envelope inflation near the Eddington limit (Sanyal et al. 2015) is not suppressed. However, to avoid convergence failure during the late evolutionary stages of the primaries and secondaries in our binary models, the option MLT++ (see Sect. 7.2 of Paxton et al. 2013) is turned on for stellar models at core helium ignition if their helium core mass exceeds > 14 M_⊙. For core hydrogen burning models, step overshooting is used so that the convective core is extended by 0.335 H_P (Brott et al. 2011) where H_P is the pressure scale height at the boundary of the convective core. Thermohaline mixing is included with an efficiency parameter of α_th = 1 (Kippenhahn et al. 1980) and semiconvective mixing is included with an efficiency parameter of α_sc = 1 (Langer et al. 1983).

Mass transfer is modeled by using the ‘contact’ scheme from MESA (Marchant et al. 2016). This is an implicit method in which the mass-transfer rate is adjusted in such a way that the radius of the donor star stays inside the Roche lobe. Moreover, as the name already suggests, it is capable of modeling a contact phase where both stars overfill their Roche volumes, as long as mass outflow via the 2nd Lagrange point is avoided (Marchant et al. 2016; Menon et al. 2021).

The mass donors in our binary models are calculated until core helium depletion, after which we assume that they form a compact object. We assume that the system becomes disrupted and thus set the orbital separation to infinity, implying that from then on the mass gainer is modeled without further binary interaction, until core carbon depletion.

2.2. The optical depth of WR winds

Classical WR stars are hydrogen-poor stars with optically thick winds that show strong emission lines. Shenar et al. (2020) found that WR stars have a minimum luminosity that depends on metallicity. Stripped-envelope stars below this luminosity are believed to have optically thin atmospheres and do not posses the emission features that characterizes WR stars (i.e., Götberg et al. 2018). Stars below the minimum luminosity are also believed to have mass-loss rates that are lower than those that would be inferred from extrapolating empirical mass-loss rate prescriptions to lower luminosities (Vink 2017; Sander & Vink 2020).

Therefore, we employed a criterion to select WR stars from our models and discard transparent-wind stripped-envelope stars. Following Aguilera-Dena et al. (2022a) we estimated the surface optical depth of our stellar models and only included those with large enough optical depths in our synthetic WR population. This criterion assumes that the velocity structure in a stellar wind can be described by a β-law with β = 1, and that the opacity κ can be approximated by the electron scattering opacity κ_es = 0.20 (1 + X_H) cm² g⁻¹. We then compute the optical depth according to (Langer 1989) as

$$\begin{array}{c}\hfill \tau (R)={\displaystyle \frac{-\kappa \dot{M}}{4\pi R(v_{\infty }-v_0)}}ln\left({\displaystyle \frac{v_{\infty }}{v_0}}\right).\end{array}$$(1)

Here, Ṁ is the wind mass-loss rate, R the radius of the star, v_∞ is the terminal wind velocity, and v₀ the expansion velocity at the base of the wind, which is on the order of the sound speed v₀ = 20 km s⁻¹.

Following Gräfener et al. (2017), the terminal wind velocity v_∞ of H-free WR stars is related to their escape velocity v_esc and can be expressed as

$$\begin{array}{c}\hfill v_{\infty }=1.3v_{\mathrm{esc}}=1.3\sqrt{\frac{2GM}{R}(1-\mathrm{\Gamma }}),\end{array}$$(2)

where Γ is the Eddington factor. For WC stars a steeper relation of v_∞ = 1.6 v_esc is used (Gräfener et al. 2013). For hot OB stars (27 500 K ≤ T_eff ≤ 50 000 K) we used a steep relation of v_∞ = 2.6 v_esc and a more shallow relation for cool OB stars (12 500 K ≤ T_eff ≤ 22 500 K), given by v_∞ = 1.3 v_esc (Vink et al. 2001). Regarding H-rich and H-poor WN stars, no such relation has been derived yet. We assume that v_∞ = 1.3 v_esc is a valid choice for this phase, as it is also used for the later H-free WN phase. For the supergiant stage the maximum value of the mass-loss rates of Vink et al. (2001) or Nieuwenhuijzen & de Jager (1990) are used. As in Nieuwenhuijzen & de Jager (1990) no such correlation is given, the correlation of v_∞ = 0.7 v_esc for temperatures below T_eff ≤ 10 000 K is used (Vink et al. 2001).

Aguilera-Dena et al. (2022a) calibrated this method for hydrogen-free WR stars using the minimum luminosities inferred from the SMC, LMC, and Milky Way by Shenar et al. (2020) and find that a wind optical depth of τ = 1.45 describes the borderline between WR stars and helium stars with transparent winds best (see their Fig. 1). For simplicity in this work τ_thresh = 1.5 is used. For WR stars with hydrogen, no such calibration is available. However, using the above scheme to estimate the optical depth of observed LMC WR stars with hydrogen indicates optical depth values between 0.3 and 1. Therefore, this criterion is employed in our models to distinguish between classical WRs, hydrogen burning WRs, and transparent-wind stripped-envelope stars (see also Appendix A).

2.3. Population synthesis

To create a synthetic WR population we employed the intrinsic mass, period, and mass-ratio distributions obtained by Sana et al. (2013) for stars in the LMC. To summarize, we assumed that the initial masses of single and donor stars are given by an initial mass function of the form f(M/M_⊙)∝(M/M_⊙)^−α with an exponent of α = 2.3 (Salpeter 1955), that the intrinsic orbital period distribution is given by f(log(P/d)) ∝ log(P/d)^π with an exponent of π = −0.45 (Sana et al. 2013), and that the intrinsic mass-ratio distribution is given by f(q)∝q^κ with an exponent of κ = −1.0 (Sana et al. 2013). The intrinsic binary fraction derived by Sana et al. (2013) is f_bin = 0.51 ± 0.04, but for simplicity we used f_bin = 0.5 in the following analysis.

To minimize statistical fluctuations, we drew 15 000 000 binary and 15 000 000 single star models according to the initial distributions. To link the drawn initial parameters to a model in our grid, a nearest neighbor approach is used.

It is assumed that there was a constant star formation rate in the LMC for the last 10 × 10⁶ yr thus the age of all stars follows a uniform distribution. The time resolution of our population model is limited by the time resolution of the stellar evolution models. With typical time steps of ≲1000 yr in our models, it is ensured that also short lived evolutionary phases are covered. A model only contributes to the synthetic WR population if its calculated optical depth is τ ≥ 1.5 (see Sect. 2.2).

To make the synthetic WR population comparable to observations it is divided into four subgroups: H-rich WN, H-poor WN, H-free WN, and WC. The synthetic H-free WN population is used to obtain a normalization factor such that the number of synthetic H-free WN stars matches the observed number of H-free WN stars in the LMC. The resulting normalization factor is then applied to the other WR subgroups. Harris & Zaritsky (2009) found that within the last 50 × 10⁶ yr the average star formation rate in the LMC was SFR_obs ≈ 0.4 M_⊙ yr⁻¹. With the normalization required to reproduce the WR population we obtain a similar result of SFR_req ≈ 0.4 M_⊙ yr⁻¹.

We note, that our synthetic population is not complete, in a sense, as not all possible evolutionary scenarios can be modeled with the MESA code yet. The supernova explosion of the primary and its effect on the binary orbit is not included. Afterwards the secondary is evolved as a single star. Common envelope evolution and its descendant stellar products are not modeled, and thus not included in the final population. We note, that common envelope evolution, especially for high mass stars, is not well understood, and the final contribution of systems that have undergone a common envelope evolution to the WR population remains highly uncertain (Ivanova et al. 2020). We discuss possible impacts of post common envelope evolution systems on our synthetic WR population in more detail in Sect. 6.3.

3. Single star and binary evolution models

In the following section, we discuss the most important aspects of the evolutionary tracks of our single star and binary models, as function of the initial stellar masses and the initial orbital periods. For a comparison with other stellar evolutionary models and their predictions on the WR population, we refer to Appendix B.

3.1. Single star models

We present five selected evolutionary tracks of our single star models in Fig. 1 that are representative for the possible outcomes of our single star grid. All models start their life on the zero-age main-sequence and they expand on the nuclear timescale during core hydrogen burning. Core hydrogen exhaustion can be identified by a hook in the evolutionary track, which corresponds to a short contraction phase before hydrogen shell ignition. For the two tracks with the highest initial masses this contraction occurs at effective temperatures below 10 kK (see also Brott et al. 2011; Köhler et al. 2015).

Fig. 1. Evolutionary tracks of selected single star models in the Hertzsprung-Russell diagram. The tracks are labeled by the initial stellar mass and color-coded by the surface helium abundance XHe (see color bar to the right). A star symbol (also color-coded), marks the point of core helium depletion. Lines of constant radius are drawn as thin straight lines.

All stellar models evolve into cool supergiants, with large envelope masses and effective temperatures below 10 kK. Here, models below ≲30 M_⊙ become classical red supergiants, while the more massive models develop He-enriched surfaces and become yellow supergiants. With the mass-loss rates adopted as described above, only models with masses above M_i ≳ 50 M_⊙ manage to remove most of their hydrogen-rich envelope and develop into a WR star. Since the mass-loss rates of cool supergiants are uncertain, so is the mass limit above which our single star models can become WR stars. We discuss the consequence of this for our synthetic WR population in detail in Sect. 6.4.

In order for our models to reach a WR phase, they need to achieve a surface helium mass fraction of Y_s ≃ 0.4 or more (see Schootemeijer & Langer 2018, for a model-independent view on this). Only then, their surface temperatures can exceed 30 kK. During their evolution as WR stars, the models go first though a phase during which the surface still contains hydrogen (H-poor WN phase). In case they lose their hydrogen envelope completely, as the two most massive models shown in Fig. 1, they settle on the helium main sequence (we do not correct the surface temperature of our WR models for optical depth effects of their winds) as H-free WN stars. In case the models also lose their helium-envelope, they may appear as WC stars. We do not consider a possible WO phase here (see Aguilera-Dena et al. 2022a).

3.2. Binary star models

The evolution of a star in a binary system differs drastically from the single star scenario as it is possible that both components interact during their lifetime. In a binary system mass transfer is initiated when one of the stellar components grows and exceeds the so-called Roche-lobe radius R_RL as the overflowing material is gravitational attracted by the other component (Tauris & van den Heuvel 2002). In Fig. 2 the evolutionary tracks of the mass-donors and mass-gainers for different initial orbital periods and mass ratios are shown.

Fig. 2. Evolutionary tracks of mass donors (panels a–d) and of mass gainers (panels e and f) for selected binary models in the Hertzsprung-Russell diagram. The tracks are labeled with the selected initial stellar masses. The mass gainers shown in panels e and f are the counterparts to the mass donors shown in panels a and b, respectively. The panels a–c depict the donor evolution for three different initial orbital periods, for binaries with a fixed initial mass ratio of qi = 0.7. Panels b and d: donors for the same orbital period (Pi = 40 d) but for two different initial mass ratios. In the panel a diagram, arrows identify different parts of the same tracks. The tracks and star symbols are color-coded as in Fig. 1. Double rings marks a predicted merger situation; our calculations end there.

Panel a of Fig. 2 shows the evolutionary tracks of the mass donors for five initial donor masses of systems with an initial orbital period of P_i = 4 d and an initial mass ratio q_i = 0.7. All five donors start mass transfer while they are still undergoing core hydrogen burning. From the indicated helium surface abundances it becomes evident that most of the hydrogen rich envelope is removed from the donor star model during the mass transfer. While the donor models with initial masses M_i = 28.2 M_⊙, 39.2 M_⊙, 50.1 M_⊙, and 63.1 M_⊙ are calculated up to core helium depletion, the binary model with an initial donor mass of 79.4 M_⊙ enters a common envelope phase which leads to a merger of both components.

Panel e shows the evolution of the mass-gainers corresponding to the mass-donors of panel a. These models accrete a large amount of the transferred matter, as indicated by the large luminosity increase compared to single star models. Their surface helium abundance is enriched during the mass transfer, but later reduced by thermohaline mixing.

The evolutionary tracks of mass-donors of systems with initial orbital period of P_i = 40 d are shown in panel b of Fig. 2. Here, the two lowest mass models start mass transfer after core hydrogen exhaustion. Models above ≳50.1 M_⊙ inflate during their main-sequence evolution and initiate mass transfer during core hydrogen burning. However, the effect of the mass-transfer phase is similar for all models, as almost all of the hydrogen-rich material is removed from the donor star model, followed by a blueward evolution. A comparison to the models shown in panel a shows that the resulting helium star models cover a larger luminosity range. The more limited luminosity increase of the corresponding mass gainers (panel f) compared to the case of the P_i = 4 d indicates that the mass transfer in the wider binaries is less efficient and most of the mass is lost from the systems.

The tracks of the mass-donors in binaries with an initial orbital period of P_i = 400 d are depicted in panel c of Fig. 2. As expected, the mass-transfer phase is initiated at a later evolutionary phase. The models have a rapid mass-transfer phase near or shortly after the end of their main-sequence. Due to their convective envelopes, the three more massive models are expected merge during the mass transfer. In contrast to the short period binary models which are shown in panel a the intermediate and long period models can produce more luminous WR stars if they are able to avoid a common envelope phase.

Panel d shows the evolutionary tracks of mass-donors with an initial orbital period of P_i = 40 d, similar as panel b, but for an initial mass ratio of q_i = 0.4. By comparing the stellar evolution tracks shown in panels b and d it becomes evident that the mass-transfer phase starts approximately at the same evolutionary stage. Even though the tracks differ slightly during the mass transfer phase, the resulting helium star models are very similar. This shows that the mass ratio is of secondary importance, as long as a binary merger can be avoided, which is inevitable for too small initial mass ratios (see Appendix C).

In summary, from the comparison between panels a and b, as well as e and f, we have learned two things: First, the initial orbital period indeed has a direct impact on the final products of the mass-donors as well as mass-gainers and second, the accretion efficiency decreases with increasing initial orbital period. Panels b and d show that the initial mass ratio is expected to have only a small effect on the outcomes of a binary evolution.

4. The observed WR population in the LMC

Based on the LMC WR cataloge of Breysacher et al. (1999), the WN population in the LMC as analyzed by Hainich et al. (2014) contains 107 stars. Among those are 17 confirmed and 22 suspected binary system. Hainich et al. (2014) analyzed the WN spectra with the Potsdam Wolf-Rayet (PoWR) atmospheric models and determined their stellar properties. The most important parameters for our analysis are listed in Table D.1.

Even though the so-called stellar temperatures of the WN stars are estimated by Hainich et al. (2014) we do not include them in our analysis, for two reasons. First, the temperature estimation of WR stars is a non trivial process as they are covered by optically thick material thus the star cannot be observed directly. The differences between the observed stellar temperatures and the ones predicted by stellar evolutionary models is still debated in the literature (e.g., Gräfener et al. 2012). Secondly we do not have temperature estimates for all WR stars, especially for the WC stars.

The luminosities are more reliable than the temperatures and are therefore used for comparison with our models. The uncertainties on luminosity are not discussed in Hainich et al. (2014). Therefore, a typical uncertainties of Δlog(L/L_⊙) = 0.2 is adopted which is used in similar works like the studies of galactic WN stars (Hamann et al. 2019).

Binary parameters, namely the orbital period and mass ratio, of known WN binaries are adopted from Shenar et al. (2019) and are listed in Table D.1 as well. The uncertainties of the orbital period are below 1% and the typical uncertainties in the mass ratio is in the order of 10%.

In addition to the WN stars, the LMC contains 24 known WC type stars. Only 8 WC have been analyzed with stellar atmosphere codes by Crowther et al. (2002), Ramachandran et al. (2018), and Hillier et al. (2021). Bartzakos et al. (2001a) determined the absolute visual magnitudes and when applying a typical bolometric correction of M_bol = 4.5 (Smith et al. 1994) it is possible to obtain the luminosity of all 24 WC stars. The calculated luminosities, as well as the luminosities obtained by stellar atmosphere codes which are used for comparison, are listed in Table D.2. The uncertainties of the calculated luminosities are about Δlog(L/L_⊙) = 0.3, which is similar to the typical uncertainties from stellar atmospheric models (Sander et al. 2012). One can see that the discrepancy for the single stars is rather small and within the given uncertainties. Only the value of the WC binary BAT99 53 differs from the calculated luminosity by Δlog(L /L_⊙) = 0.5.

The orbital periods of known WC binaries are adopted from Bartzakos et al. (2001b) and are listed in Table D.2. The uncertainties of the orbital period are below 1%.

Complementary to the WN and WC stars, the LMC hosts 3 known WO type stars, namely LH41-1042, LMC195-1, and Sanduleak 2. Typically one assumes that these stars have enriched surface oxygen abundances, but in the recent work of Aadland et al. (2022) it is shown that WO type stars only have a higher C content than WC type stars while their surface O abundance is comparable to those of WC type stars. It is still believed that these stars reflect the final stages of core He- and C-burning, and hence cannot be not fully covered by our model grids (Tramper et al. 2015; Sander et al. 2019; Aguilera-Dena et al. 2022a). Nonetheless, we want to give a brief comparison of the fundamental stellar parameters of these WO type stars and the predictions of our model grid.

LH41-1042 was studied by Tramper et al. (2015) and has a luminosity of $logL/L_{\odot }=5.{26}_{-0.14}^{+0.12}$ and a surface carbon abundance of X_C = 0.60. Aadland et al. (2022) recently reanalyzed LMC195-1 and Sanduleak 2, finding that both have similar luminosities of log L/L_⊙ = 5.41 and surface carbon abundances of X_C = 0.62. Our models predict a maximum surface carbon abundance of X_C = 0.53. Aadland et al. (2022) have shown in their work that also other stellar evolutionary models have problems in predicting this high carbon surface abundances. They find that it is linked to a too efficient ¹²C+⁴He→¹⁶O nuclear reaction rate which already becomes important during core He-burning. We believe that this might also be the case in our models. However, a more in-depth study of this problem would be beyond the scope of our paper. Despite the disagreement of the carbon surface abundance, we find that our models with X_C ≳ 0.50 have luminosities in the range of log L/L_⊙ = 5.3 − 5.5 which is close to those of the observed WO stars.

Our models only cover a luminosity range until log L/L_⊙ ≤ 6.35, thus observed stars with higher luminosity are excluded from the sample. Furthermore, we exclude stars with spectral types Of as they are on the verge of becoming WR stars but in most of the cases they are quite cool and have a high amounts of H in their envelopes, which is not predicted by our stellar models when applying the optical depth criterion.

5. The synthetic WR population

5.1. Luminosity distributions

In this section, we evaluate the number of WR stars of different types, obtained from our models via the single and the binary channels, and derive the luminosity distributions of the different WR types from our synthetic WR populations. We then compare these findings with the observed WR population in the LMC. In Sect. 5.3, we show the main properties of our WR binary models and compare them with the – more limited – observed properties of WR binaries in the LMC.

As discussed above, we restrict ourselves to compare the luminosity distributions of the synthetic and observed WR population in the LMC, temperatures and radii are not always available and are subject to larger uncertainties. In the work of Hainich et al. (2014) the WN population of the LMC was studied spectroscopically by using a model atmosphere grid with surface hydrogen abundances of X_H = 0.0, 0.2, and 0.4. Here, we follow their differentiation and distinguish “H-free WN stars” as having X_H < 0.05 (undetectable), “H-poor WN stars” with X_H ≲ 0.4, and “H-rich WN stars” with X_H > 0.4. Notably, the observed WN population in the LMC shows a strong drop in the number of objects with hydrogen abundances above 0.4 (see below).

According to Hainich et al. (2014), four genuine WN stars show hydrogen abundances above 0.4, with a hydrogen mass fraction of either 0.6 or 0.7. Two of them have luminosities of log L/L_⊙ = 6.7 and log L/L_⊙ = 6.8, well above the luminosity limit of our model grid (log L/L_⊙ = 6.35), and thus not considered here. The other two have a luminosity close to this limit (BAT99-49 with log L/L_⊙ = 6.34 and BAT99-111 with log L/L_⊙ = 6.25). BAT99-49 is member of an eccentric (e = 0.35) binary system (Foellmi et al. 2003), which excludes a previous Roche-lobe overflow phase. BAT99-111 is a bright X-ray source and as such likely a colliding wind binary. In agreement with Hainich et al. (2014) and others, we conclude that these four stars are not stripped envelope stars, but rather core H-burning stars with extreme mass and luminosity such that their stellar wind becomes optically thick.

Of the stars with luminosities below log L/L_⊙ = 6.35 in the list of Hainich et al. (2014), 37 have a WN spectral type and a H-abundance in the range 0.1–0.4, 38 have a WN spectral type and no hydrogen detected, and 19 stars have a WC spectral type. The resulting luminosity distributions of the observed and predicted WR population are shown in Fig. 3. The predicted WR population contains 24.3 H-poor WN star models, 38 H-free WN star models and 19.3 WC star models, where the matching number of H-free WN stars is due to our applied normalization.

Fig. 3. Luminosity distribution of observed H-poor WN (orange), H-free WN (green), and WC (blue) star populations (colored shading) with our synthetic WR population (black histogram, τ ≥ 1.5). The synthetic WR population distinguishes WR models obtained from single star (dotted bins) and binary models (dashed bins). The observed H-poor WN population excludes stars with spectral types Of and Of/WN. In the upper right corner of each plot, we give the number of observed stars, the numbers of predicted single stars and binary produced WR stars of the respective WR subtype. For the synthetic H-poor WN populations an additional histogram (brown) for lower optical depths (τ ≥ 0.5) is shown.

The synthetic H-free WN and the (also H-free) WC star populations show an overall agreement with the observed luminosity distribution. All of these models are predicted to be core helium burning which is in agreement with the overall picture of the evolution of WR stars. Our synthetic populations predict that the majority of the less luminous stars (log L/L_⊙ ≲ 6.0) are expected to be in binary systems and that the more luminous WR stars originate from stars evolved in isolation.

The synthetic luminosity distribution of the H-poor WN stars with an optical depth above τ ≥ 1.5 contains only 24.3 stars and clearly underestimates the observed population. Our synthetic population shows again, that the majority of these stars is expected to be in binary systems and that the single star models dominate at luminosities above log L/L_⊙ ≳ 6.0. The observed distribution of H-poor WN stars also shows a double peaked distribution which might be linked to the binary nature of the systems. However, the observed distribution contains many more stars at luminosities around log L/L_⊙ ∼ 5.7 which are not reproduced by the synthetic population.

The threshold of τ ≥ 1.5 was determined by using theoretical prescription on the stellar properties and mass-loss rates of H-free WN stars. However, the values used as input in Eq. (1) may differ for H-poor WN stars, especially the core hydrogen burning ones. Therefore, the wind optical depths for such stars may be different (i.e., lower) as well as their threshold value. We explore this possibility by modifying the optical depth threshold for the H-poor WN distribution. This is roughly equivalent to performing a new calibration for the optical depth threshold in this population, which is otherwise impossible to perform due to the uncertainties in their luminosity distribution. The best agreement between the number of H-poor WNs and their shape can be found when using τ ≥ 0.5. The resulting synthetic population predicts that almost 50% of the H-poor WN stars are currently core hydrogen burning.

In all predicted luminosity functions, the single stars produce a peak above 10⁶ L_⊙. Except perhaps for the WN stars with a bit of hydrogen, this peak does not find a correspondence in the observed luminosity distributions. An initial binary fraction of the most massive stars above 50% could remedy this situation (cf. Sect. 6.5).

A further discrepancy between the observed and predicted distributions is the following. Our synthetic WR population predicts that about 80% of all WR stars are in binary system, while only ∼40% of the WR stars in the LMC are known to be in a binary system (Hainich et al. 2014). According to Hainich et al. (2014), it is possible that a significant fraction of the apparently single WR stars are in fact in a binary system with a yet undetected companion. We return to this point in Sect. 6.6.

5.2. Surface hydrogen abundances

Figure 4 shows the predicted distribution of surface hydrogen abundances and luminosities of our WN models with hydrogen. We find that the most likely hydrogen abundances are in the range X_H = 0.15 − 0.25 throughout most of the considered luminosity range. Furthermore, our models predict the possibility of higher hydrogen abundances for more luminous stars. Focussing on the distribution obtained for a limiting optical depth of τ = 0.5, we predict luminous WN stars with hydrogen abundances of up tp X_H = 0.6, of which, according to Fig. 3 the majority still undergoes core hydrogen burning. We note that for luminosities below 10⁶ L_⊙, where our single star models do not contribute, H-rich WR stars may be produced during the slow Case A mass-transfer phase, as detailed in Sen et al. (2022).

Fig. 4. Surface hydrogen abundances as a function of luminosity of the observed H-poor WN stars (orange dots) and of our synthetic population (histogram). Black lines indicate borderlines of the theoretical distributions obtained for different optical depth thresholds (τ ≥ 0.5, 1.0, and 1.5).

As for the lumnosity distributions, when comparing to the observed distribution, the agreement is fair when τ = 0.5 is used to defined the WR star models. Except for two very luminous ones, all observed stars fall within the predicted range. From the trend for the borderline of the predicted distributions as function of threshold optical depth, one can anticipate that these two objects would also be recovered if a smaller optical depth threshold would be adopted (cf. Sect. 2.2).

The area where we expect the largest number of objects in Fig. 4 is well populated with observed WR stars. However, we see perhaps more stars at X_H = 0.4 than expected and we predict more stars with low hydrogen abundances X_H ≃ 0.1 than observed. On the other hand, we note that the observed hydrogen abundances are derived by Hainich et al. (2014) based on three grids of atmosphere models, using X_H = 0, 0.2, and 0.4, which are the values preferentially found in the observed stars (X_H = 0 is not shown in Fig. 4). We conclude that the hydrogen abundances in the models of our synthetic WR population agrees rather well with the observations.

5.3. WR binary properties

In the previous sections we focused on the stellar properties of the WR stars in the LMC. In this section, we investigate their binary properties of our WR population. We derive the predicted distributions of orbital periods and mass ratios for WR binaries with different WR subtypes and compare the with the available observations. In contrast to luminosity and hydrogen abundance, where the available information is rather complete, orbital periods are known for only 18 LMC WR binaries and mass ratios only for 8 of them.

5.3.1. H-poor WN stars

The orbital period distribution of the observed and predicted H-poor WN population in the LMC can be seen in the upper right panel of Fig. 5. Most of the observed periods of the H-poor WN binaries are below P ≲ 10 d. This is linked to an observational bias toward smaller orbital periods as the companion in these systems has a stronger impact on the appearance of the spectrum, namely by a stronger Doppler-shift of the emission and absorption lines of each stellar component. Due to this observational bias, it is reasonable to assume that the sample is nearly complete for the shortest periods and nearly incomplete for long periods. In the sample of observed periods there are two binary systems with orbital periods above P ≳ 100 d which are brighter than our most luminous models (log L/L_⊙ ≥ 6.35). Both binaries show X-ray emission, most probably originating from a colliding wind, which is an additional indicator that these stars have a binary companion.

Fig. 5. Predicted mass-ratio and period distributions of the H-poor WN stars. Upper left: predicted mass-ratio distribution of the H-poor WN stars in our synthetic population with optical depths τ ≥ 0.5 (brown dashed histogram) and τ ≥ 1.5 (black dashed histogram), where we distinguish WR stars during core hydrogen and core helium burning. Overlayed is the histogram of the observed mass ratios in LMC WN stars with hydrogen (orange). Upper right: same as in the upper left panel, but for the orbital periods. Lower left: 2D histogram of mass ratios and orbital periods of our synthetic population with optical depths τ ≥ 0.5, represented by colored squares, where the color indicates the expected number of objects (see color bar to the right). Individual observed binaries are shown as orange dots. Lower right: same as the lower left panel but now for optical depths τ ≥ 1.5. The observational data is take from Shenar et al. (2019).

By comparing the small sample of observed orbital periods to the predicted distribution of the H-poor WN binaries several similarities become apparent as well as some discrepancy. Our models with optical depths τ ≥ 1.5 are under-predicting the number of binaries with orbital periods of about P ≈ 2 d, hinting toward a unconsidered formation channels. Possible formation channels are discussed in Sect. 6. However, our models with τ ≥ 0.5 can reproduce the short period binaries and predict that they are still core hydrogen burning stars. Another discrepancy is that our models predict a large amount of systems with long period (P ≳ 100 d) that are barely observed. This is most likely linked to the aforementioned observational bias.

The mass ratio distributions, found in the upper left panel of Fig. 5, of the observed H-poor WN population is as small as the sample of the orbital periods. The sample is spread over a wide range of mass ratios and with respect to the small sample size is in good agreement with the predicted distribution, as there are no outliers as in the case of the orbital period distribution.

In the lower panels of Fig. 5 2D histograms of the combination of the period and mass ratio distributions with the two different optical depth criteria are shown. The predicted distributions show a wide variety of the mass ratio q = 0.4 − 2.4 at small orbital periods P ≲ 1.0 d. When going to larger orbital periods the mass ratio distribution gets shallower and closer to q = 1, which is linked to the efficiency of mass transfer. It appears that the predicted distribution of mass ratios and orbital periods is in rough agreement with the observed distribution with respect to the small available sample size and their respective uncertainties.

5.3.2. H-free WN stars

Figure 6 contains the observed and predicted period and mass ratio distributions of the H-free WN binaries. Similar to the orbital period distribution of the H-poor WN binaries, the predicted distribution under-estimates the amount of short period WR stars and over-estimates the binaries with long orbital periods. This can also be a consequence of observational biases. The under-prediction of the short period binaries hints at an unconsidered formation channel which seems to be independent of the WR subtype (see Sect. 6.2).

Fig. 6. Same as Fig. 5, but for H-free WN stars. In the right panel, the 2D histogram, systems with unknown mass ratio but known orbital period are indicated by arrows. The observational data is take from Shenar et al. (2019).

Compared to the predicted orbital period distribution of the H-poor WN binaries, the orbital period distribution of the H-free WN binaries is shifted toward longer orbital periods. This shift can be explained by the mass that is lost from the system during the WR phase, leading to a change in the orbital angular momentum and, therefore, to an longer orbital period.

The observed mass ratio distribution is too small to get any useful information, but at least none of the observed mass ratios disagrees with the predicted mass ratio distribution. Compared to the predicted mass ratio distribution of the H-poor WN binaries the mass ratio distribution of the H-free WN binaries is shifted toward more extreme mass ratios q > 1. This shift is also connected to the strong mass-loss rates during the WR phase, which drastically alters the mass of the WR star on a timescale which is short compared to the main-sequence lifetime of the companion.

Inspecting the 2D histogram of the predicted mass ratio and orbital period distribution of the H-free WN binaries in the right panel of Fig. 6 one can see that this distribution now spreads over a large range of mass ratios and orbital periods with a small tendency toward shorter orbital periods P ≈ 20 d. Comparing this with the small sample of observed H-free WN binaries (except of the binaries with the shortest periods) no discrepancies can be found between predictions and observations.

5.3.3. WC stars

The orbital period and mass ratio distributions of stars in the WC phase are shown in Fig. 7. The sample of observed WC binaries contains three short period binaries, of which two (those with the shortest periods) cannot be explained by our models. However, there is a lack of intermediate orbital periods around ∼10 d which is also found by Dsilva et al. (2020) for the galactic WC binaries.

Fig. 7. Same as Fig. 5, but for WC stars. The observational data is take from Bartzakos et al. (2001b).

By comparing the predicted orbital period distribution with the ones of the H-free and H-poor WN binaries, a shift toward even longer orbits can be seen. As in the case of the H-free WN, we explain this shift by the strong mass-loss rates that lead to a widening of the orbit. Furthermore, it becomes evident that the predicted number of WC binaries is much smaller than the number of WN binaries. This is linked to the evolution of the stars themselves, as not every WN star will evolve into a WC star because the helium rich layers will not be removed for all WN stars and they spend less time as WCs than as WNs in this mass and metallicity range (see Aguilera-Dena et al. 2022a).

Unfortunately, the mass ratios of the WC binaries have not been measured yet. Therefore, we can only show the predicted distribution and compare it with the ones of the other WR subtypes. As in the case of the H-free WN binaries, the distribution of the mass ratios is again shifted toward more extreme mass ratios q > 1 as the donor star looses several solar masses via its strong WR wind.

The 2D histogram in the right panel of Fig. 7 shows that the expected distribution of mass ratios and orbital periods is rather uniform for orbital periods in the range of P = 10 d − 1000 d and for mass ratios in the range of q = 0.8 − 2.0. The distribution clearly lacks of short orbital periods P ≲ 10 d and mass ratios below q ≲ 0.8.

6. Discussion

In the previous section we found that the observed WR population of the LMC with luminosities below 10⁶ L_⊙ cannot be explained by our singe star models. Therefore, we compare our findings with the most commonly used single and binary evolution models to see whether they lead to similar conclusions or not. In this Section, we briefly summarize the predictions of the WR population of the different stellar evolutionary models. A more detailed comparison can be found in Appendix B. Moreover, we discuss possible formation channels which are not included in our analysis that can explain the discrepancies between the predicted and observed WR population.

6.1. Comparison with previous single and binary star evolution models

We first focus on the similarities and differences of nonrotating stellar evolution models. Commonly used are the so-called Geneva models (e.g., Ekström et al. 2012; Georgy et al. 2012, 2013; Eggenberger et al. 2021). These models are available for Galactic (MW), LMC and SMC (Z = 0.014, 0.006, and 0.002) metallicity which makes it possible to compare these models with our models and to use them as a proxy for other model grids that do not cover LMC metallicity. Their lower WR luminosity limits are log(L /L_⊙)≈5.4, log(L /L_⊙)≈5.75, and log(L /L_⊙)≳6.4 (see Fig. B.1) for MW, LMC, and SMC models, respectively. Their 40 M_⊙ LMC model spends a very short time in a H-free WR phase, which is not predicted by our grid. Their 60 M_⊙ LMC model predicts similar WR phases, however, their model is less luminous. The 85 M_⊙ model is predicted to evolve similarly like our 80 M_⊙ model. Our models seem to fit the general trend, although the WR luminosity limit of the Geneva models for the LMC is lower than our prediction which is linked to the more efficient red supergiant (RSG, T_eff ≤ 4800 K) winds in their models.

The MIST (Choi et al. 2016) stellar evolutionary models are also calculated with MESA but use less efficient overshooting and semiconvective mixing. It is a known issue, that the choice of the mixing parameters can have a deep effect on the resulting stellar populations, particularly after the main-sequence (e.g Schootemeijer et al. 2019; Gilkis et al. 2021; Klencki et al. 2022). Their LMC models ([Fe/H]= − 0.5) only produce WR stars with luminosities above log(L /L_⊙)≳6.4 (see Fig. B.1) as their models spend most of their time in the RSG phase.

Finally, the FRANEC models (Limongi & Chieffi 2018) are only calculated for MW and SMC metallicity ([Fe/H] = 0 and [Fe/H]= − 1 ). Their lower WR luminosity limits are log(L /L_⊙)≈5.4 and log(L /L_⊙)≈6.2 for MW and SMC, respectively. As for the Geneva models, our models seem to fit the general trend. Even though we are unable to directly compare the models, their predicted H-poor, H-free WN, and WC phases appear to agree better with our models (cf. Fig. B.1).

Comparing the different slow and the fast rotating models of the different codes, presented in Fig. B.1, one can see two major effects: First, rotation may extend the time spent in the WR phase significantly as seen from the Geneva models, but this needs not to be the case, as the FRANEC models show. Second, rotation leads to more efficient mixing in the stars, which then have a shortened H-free WN phase.

From the MIST models one can conclude that the choice of mixing and the stronger RSG mass-loss rates can affect the appearance of the WR population drastically. Even though our models spend less than a quarter of their core helium burning time in the cold supergiant phase (T_eff ≤ 12 500 K) and are in relative good agreement with the other single star models, we are still overestimating the RSG population. Gilkis et al. (2021) studied this phenomenon in more detail using MESA models and show that more efficient mixing can shorten the RSG phase and avoid an overprediction of cool supergiants.

We conclude that single star models with enhanced mixing and/or with enhanced mass loss in the cool supergiant regime, may predict WR stars down to luminosities of about log(L /L_⊙)∼5.6, at LMC metallicity. However, they experience difficulties reproducing the observed fainter population of WR stars in the LMC (see also Shenar et al. 2020).

For the binary models we compare our models with those from Eldridge et al. (2017). Their binary models are calculated for LMC metallicity (Z = 0.004) and are based on detailed stellar evolution models. We note, that typically a metallicity of Z = 0.006 is used for the LMC, but as already argued in Sect. 2.1, observations support a lower iron abundance which is taken into account in our model grid. Thus for a comparison to our stellar evolutionary models we prefer to use the BPASS models calculated at the lower metallicity of Z = 0.004. We find, that most of their models spend similar times in the individual WR phases and have similar luminosities compared to our models (cf. Fig. B.2). We find that some of their models that undergo a Case A mass transfer are predicted to produce WR stars that are by far more luminous than our models predict. This is because their models enter a common envelope phase, during which they merge and, therefore, leading to the formation of more massive and luminous WR stars. Besides that we find that their models spend to some extend different amount of time in the individual WR phases. The binary models of Eldridge et al. (2017) differ to our models by four assumptions: First, in BPASS the WR mass-loss recipe of Nugis & Lamers (2000) is used only when X_H < 0.4 and T_eff > 10 kK. Second, RLOF is modeled using a simplified formula to ensure numerical stability (Eldridge et al. 2017, their Eqs. (2) and (3)). Third, the secondary is approximated by the single star stellar evolution equations of Hurley et al. (2000) and only retrospectively modeled in detail. Fourth, their models do not include rotation. We find that the different physical assumptions can well explain the differences in the time a stellar model spends in the individual WR phases.

6.2. Short period WR binaries

In Sect. 5.3 it was shown that the observed shortest-period WR binaries, with orbital periods below P_orb < 3 d are not explained by our binary models. Below, we discuss three possible implications of this finding.

Common envelope evolution. The initially widest interacting binaries in our model grid are expected to undergo a common envelope evolution, which we do not model. In this scenario, the binary goes through an unstable mass-transfer phase, during which the primary’s envelope engulfs the secondary star. This is followed by an spiral-in process during which the orbital energy of the secondary is deposited in the envelope of the primary. This might lead to the ejection of the envelope and the emergence of a compact binary, consisting of the core of the primary (now a WR star) and the mostly unaffected secondary (Ivanova et al. 2013). While we can not exclude that the observed shortest-period WR binaries evolved through this path, simplified one-dimensional estimates which compare the binding energy of the primary’s envelope with the released orbital energy indicate that for massive stars, a successful common envelope ejection appears difficult as long as the companion star is a main sequence star (e.g., Kruckow et al. 2016). Pauli (2020) investigated the binary models presented in this paper and found that none of them with an initial primary mass above M_1,  i ≳ 28 M_⊙ is expected to eject the envelope. However, it is possible that there are unconsidered phenomena that contribute to a successful ejection, for example pulsations during the supergiant phase or LBV eruptions (e.g., Langer et al. 2020), such that we can not exclude this formation channel.

Chemically homogeneous evolution. Assuming that close binaries at the ZAMS are tidally locked leads to rapidly rotating stars in the binaries with the shortest orbital periods. This scenario has been suggested by de Mink et al. (2009) to lead to chemically homogeneously evolving binary components and investigated with MESA models using a similar input physics as our models for a large initial parameter space by Marchant et al. (2016) and Hastings et al. (2020). Their conclusion was that this scenario applies only to very metal poor binaries and only to the most massive systems (e.g., possible to the extremely massive SMC binary HD 5980; Koenigsberger et al. 2014). Our binary models with initial orbital periods below ≲2 d all result in a merger during their contact phase and binary systems with larger initial orbital period cannot explain the observed WR binaries, which renders this scenario unlikely.

Angular momentum loss by winds. When wind particles leave a massive star and the binary system, they may still interact gravitationally with both stellar components and thereby enhance or reduce the orbital angular momentum of the binary. This mechanism, which is not included in our models, may be of relevance when the wind speed is comparable to the orbital velocities of the two stars, which may be particularly applicable to the tightest binaries (Brookshaw & Tavani 1993; de Mink & Mandel 2016). MacLeod & Loeb (2020) show that in binaries with identical components a strong orbital shrinkage may occur. Whether this scenario can explain the shortest period WR binaries remains unclear. However, since the terminal wind velocities of WR stars are typically larger than v_∞ ≳ 2000 km s⁻¹, while the orbital velocities in our short-period models rarely exceed 400 km s⁻¹ we assume this effect to be small.

In the end, we can not rule out any of the three channels and conclude that all three deserve deeper study. An understanding of the shortest period WR binaries my be of significant relevance for a reliable prediction of the population of black hole mergers in the universe.

6.3. Disregarded binary products

Due to the limitations of our method, we need to disregard various binary evolution products in our synthetic WR population. Here, we briefly discuss the various neglected channels and give an estimate of their relevance.

As shown by our example plots giving an overview of the fates of individual binary models for a given primary mass (cf. Appendix C), we expect our binaries to merge in various corners of the parameter space spanned by the initial orbital period and the initial mass ratio. In particular, very short and very large initial periods, as well as extreme mass ratios are found to lead to merger conditions. Overall, we see from these plots that perhaps one quarter of all computed binaries are thought to have this fate, which, for flat initial distributions in log period and mass ratio might represent the corresponding fraction in a binary population.

As the merger product is a single star, the mergers will enhance the single star population. Thus, to first approximation, neglecting the merger products may be simply compensated by adopting a smaller initial binary fraction – which is still uncertain anyways. In case one of the merging components already undergoes core helium burning, the merger product may strongly differ from a single star. However, the fraction of these cases in a given population is estimated be small, due to the small lifetime of these merger products (Justham et al. 2014).

In our binary evolutionary models the secondary star is evolved as a single star after the primary has depleted helium in its core. This might be in general a good approach when assuming that the BH formation is associated with a birth kick as that of neutron stars (Janka 2013) and be representative for what happens in a sizable fraction of the systems. However, this might not be true for all systems in which the primary forms a BH. Some systems might not get disrupted, even though the orbit could widen significantly. For those systems, we neglect in our final population the potential WR stars that would be formed when the secondary interacts with compact companion via RLOF.

The number of WR stars formed in this way will be smaller than the number of WR stars formed from the first mass transfer, because binaries can break up when the compact object forms, or merge upon mass transfer from the initial secondary to the compact object. In fact, in case the formed compact object is a neutron star, both processes are very likely, such that WR stars are not expected. However, when the compact object is a BH – which is more likely in the mass range we consider here – the first process depends on the unknown magnitude of BH formation kicks (Janka 2013; Mirabel 2017; Chan et al. 2018). The second process depends on the weakly constrained rate at which the BH can expel the transferred matter from the binary (King & Begelman 1999; King et al. 2000). It is therefore difficult to quantify the number of WR stars formed from the inverse mass-transfer process. Currently, with Cyg X-3, there is only one candidate WR+BH binary in the Milky Way (van Kerkwijk et al. 1992).

While we are not able to quantify the fraction of WR stars originating from secondaries via stable RLOF, the above arguments substantiate that the WR stars originating from primaries are dominating in any given population. We note that, while this implies that the comparison of our model results with the observed LMC WR population is still meaningful, it does not preclude that a fraction of the LMC WR stars hosts a compact object, as speculated by van den Heuvel et al. (2017).

6.4. WR mass-loss rates

Stellar wind mass loss of massive stars is still poorly understood, especially those of the later evolutionary stages. In the case of WR stars, empirical mass-loss rates have been proposed, but the true mass-loss rates strongly depend on the assumed physics, like wind clumping. In our models we use a clumping factor of D = 4 as suggested by Yoon (2017) in order to be able to reproduce the luminosity range of the observed WC stars in the LMC. In the literature, for a long time a clumping factor of D = 10 was used for spectral analyses of WR winds, but this value is not accurately constrained (Nugis et al. 1998; Hainich et al. 2014). However, increasing the clumping factor from D = 4 to D = 10 would imply that the mass-loss rate drops by a factor of ≈1.58. This would imply that many of the WC stars predicted by our synthetic WR population could not have been formed. If one would use an unclumped wind with D = 1, the mass-loss rates would increase by a factor of 2 compared to our case, which would lead to the formation of too faint WC stars. Evidently, in particular number and properties of WC stars depends sensitively on the adopted WR mass-loss recipe. We note that while using the one proposed by Yoon (2017) allows us to recover the lower luminosity limit of LMC WC stars, their absolute number and luminosity function are not prescribed by selecting this recipe. While a better understanding of the winds from WR stars is highly desirable, we note that the recipe used here is in good agreement with the recent findings of Sander & Vink (2020).

6.5. Single star contribution

While our predicted WR luminosity functions represent the corresponding observed luminosities rather well, the high-luminosity peaks produced by the single star contributions appear to be too large (see Fig. 3). The largest discrepancies occur for the H-free WN stars and for the WC stars, for which the single star models provide large contributions above ≳10⁶ L_⊙ which are not reflected in the observed WR star population. A simple solution is to reduce the single star fraction would be to increase the binary fraction. As Fig. 3 shows, reducing the single star contribution by a factor of two, that is assuming a binary contribution of 75%, would yield a much better agreement.

A larger binary fraction would also help to alleviate another problem of our single star models. In our synthetic population, we predict 3 red supergiants (T_eff ≤ 4800 K) and 5 yellow supergiants (T_eff ≲ 7500 K) with luminosities above log(L /L_⊙) = 5.6. Davies et al. (2018) and Gilkis et al. (2021) find that the number of observed stars of this types is smaller by roughly 50%. While Gilkis et al. (2021) conclude that stronger mixing in stars can reduce the number of luminous cool supergiants, they do not investigate which effect this has on the WR luminosity functions. From their Fig. 6, we expect that they would indeed overproduce very luminous WR stars (log(L /L_⊙)≥6). The same is true for models with strong rotationally induced mixing (cf. Fig. B.1).

Our results indicate that for the most luminous stars, with initial masses above ∼40 M_⊙, a binary fraction larger than >50% might resolve the two issues discussed above. An increased initial binary fraction would be in agreement with the results from Moe & Di Stefano (2017) who find for stars in our Galaxy with initial masses above > 16 M_⊙ an intrinsic close binary (log(P/d)≤3.7 and q > 0.1) frequency of f_bin = 1.0 ± 0.2. On the other hand, we can not exclude the possibility that in the large parameter space of mass loss and mixing, there may be single star models which also overcome the mentioned problems. In any case, our results show that the contribution of binary stars to our understanding of the populations of the various types of massive stars is essential and that any calibration of the uncertain parameters in massive star evolution based on comparing only single star models with observations may lead into the wrong direction.

6.6. The WR binary fraction

Several efforts to determine the fraction of WR stars which reside in binary systems are documented in the literature. Foellmi et al. (2003) and Schnurr et al. (2008) investigated the WN stars in the LMC and found a binary fraction on the order of 30%. More recently, Hainich et al. (2014) and Shenar et al. (2019) derive a WN binary fraction between 16 and 36%. In the following, we refrain from discussing the binary fraction of WN stars with hydrogen, because in the corresponding observational samples we find many objects which may not be genuine WN stars, but rather fall into the Of category. From the 38 H-free LMC WR stars listed in Hainich et al. we find 5 certain and 6 possible binaries, leading to a binary fraction in the range 13%−29%. Bartzakos et al. (2001a) investigated the WC stars in the LMC and found binary induced radial velocity variations in 13% of their sample, but signatures of an O star in the spectrum for 65%.

The numbers appear difficult to reconcile with to our predicted WR binary fractions. In our models, assuming a single star fraction of 50%, we find an overall WR binary fraction of ∼80%, while for WR stars with luminosities below ≲10⁶ L_⊙ we expect nearly all WR stars to reside in binary systems. However, the quoted numbers are detection fractions and bias corrections are difficult. When considering the measured orbital periods of WR binaries (cf. Figs. 6 and 7), we find an average of 13 d for the H-free WN stars and 6.6 d for the WC stars. In contrast, average predicted orbital periods are above 50 d for both groups. Since radial velocity variations are larger for short-period binaries, conceivably, many long period WR binaries may have escaped detection so far. Similarly, though based on even less data (Fig. 6), the two mass ratios measured in H-free WN binaries indicate O star companions with nearly twice the mass of the WR star, while our models predict many lighter companions, with spectral types well into the B star regime (see also Langer et al. 2020).

We conclude that our binary models can not be ruled out based on the mismatch between predicted and detected WR binary fractions. Notably, in a recent study of northern Galactic WC stars, Dsilva et al. (2020) find a bias corrected minimum binary fraction of 72%.

7. Conclusions

In this work, a large grid of detailed single star and binary evolution models at LMC metallicity is used to create a synthetic WR star population. These models are calculated with the MESA code and include the physics of mass-loss, differential rotation, angular momentum transport by magnetic fields, inflation, tides, and binary mass and angular momentum transfer. The grid covers initial primary masses in the range of M_1, i ≃ 28 M_⊙ − 89 M_⊙. We compare our results with the observed LMC WR population, which is thought to be complete.

The WR stars in our the synthetic population reproduce the number of observed WR stars of different subtype (WN with hydrogen, WN without hydrogen, WC) and their luminosity distributions well (see Fig. 3). Similarly, we find a good overall agreement of the observed and predicted hydrogen mass fractions in WN stars (see Fig. 4). This is remarkable, since in contrast to single star models, the properties of WR star models produced by mass transfer are much less affected by the uncertain physics of stellar wind mass loss. Based on the observed high fraction of main sequence binaries, a large impact of binary mass transfer on the predicted WR population is expected and confirmed. This renders predictions for massive star populations which are solely based on single star models questionable.

Our results also raise several new questions. When comparing the distributions of our WR model binaries, in particular orbital periods and mass ratios, with the (sparse!) available observations (Figs. 6 and 7), we find that our models do not reproduce the few shortest-period WR binaries (P_orb < 3 d). We discuss several possible reasons for this (Sect. 6.2), and point out that discriminating these may be important for the predictions of BH mergers.

When adopting an initial binary fraction of 50%, we find that the large impact of single stars on the predicted number of very luminous WR stars (> 10⁶ L_⊙) is not seen in the observational data (Fig. 3), while halving the initial single star fraction for stars above ∼40 M_⊙) can remedy this and at the same time prevents an overproduction of luminous red and yellow supergiants (Sect. 6.5).

A large initial binary fraction leads to a large expected binary fraction of the WR stars. Our models predict an overall binary fraction of ∼80%. While this may seem to be at odds with rather modest binary detection fractions in observed WR populations, our models imply that the bias correction may be large. This view is supported by a recent study of Galactic WC stars, where a bias corrected minimum binary fraction of 72% is derived, while based on actual detections, a binary fraction of 13% has been suggested for the LMV WC stars (cf. Sect. 6.6).

Clearly, our work is just a small step to assess the impact of binary evolution on the formation of WR stars and on massive star evolution as a whole. The WR stars, due to their outstanding spectra, appear suitable to provide well defined populations for robust tests of our models. However, in the end, we need to compare the same models with populations of many other types of massive stars (contact systems, massive Algols, X-ray binaries, BH-merger, etc.) and at different metallicities, before robust conclusions can be drawn.

Acknowledgments

The authors thank Koushik Sen for his useful discussions and advise when analyzing the binary model grid. D.P. acknowledges financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50 OR 2005. D.R.A.-D. is supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of “Science and Society” Action Always strive for excellence – “Theodoros Papazoglou” (Project Number: 01431). P.M. acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N

References

Appendix A: The stellar wind optical depth

In Sect. 2.2 we calculated the optical depth of the winds with the formula given in Eq. 1 for the main-sequence and post main-sequence evolution of our single and binary models. In the following we want to show exemplary phase diagrams of the theoretical optical depth of the wind throughout the different evolutionary stages of each model would have.

Figure A.1 shows the calculated optical depth of the winds of our single star models, and of the surface hydrogen abundance, as function of mass and normalized time during the main-sequence and post main-sequence evolution. During their post main-sequence evolution, models with initial mass M_i ≳ 50 M_⊙ are able to efficiently strip off their hydrogen-rich layers during the supergiant phase and go through a WR phase. During the transition from the supergiant stage to the WR phase, the optical depth of the wind rises quickly. The wind transitions from an optically thin to an optically thick wind occurs synchronously with the drop of the surface hydrogen abundance close to X_H ≈ 0.3. Therefore, the produced helium stars are suggested to be observed as WR stars according to our criterion. Due to the quick transition from an optically thin to an optically thick wind a lower threshold of τ ≥ 0.5 does not lead to a significant change on the contribution of the single stars to the H-poor WN population (see Fig. 3).

Fig. A.1. Optical depth (left, color-coded) and surface hydrogen abundance (right) of our single star models as function of initial mass and normalized time, during their main-sequence (lower) and core helium burning (upper) evolution. The transition line from an optically thin to an optically thick wind is marked by a solid black line labeled with “1.500”. The unlabeled solid black lines divide the different considered WR phases as labeled in each region. In the right plot, threshold values of the surface hydrogen abundances in steps of ΔXH = 0.1 are indicated by thin solid lines, which have dark red color for the range XH = 0.8 − 0.5 and brown color for XH = 0.4 − 0.1. The borderlines of XH = 0.4 and XH = 0.0 are marked by solid black lines, as well as the borderline beyond which stellar models have a surface carbon mass fraction in excess of XC = 0.01, which is our threshold for defining WC stars.

None of our single star models develops an optically thick wind while the surface hydrogen abundance is above X_H ≳ 0.3. This means, that our models predict that all observed WR stars with surface hydrogen abundances above X_H ≳ 0.3 are formed in binary systems. The quick transition from an optically thin to an optically thick wind in the models is linked to the switch from the Vink et al. (2001) OB wind scheme to the Nugis & Lamers (2000) WR wind scheme when the surface hydrogen abundance drops from X_H = 0.7 to X_H = 0.4. Of course as discussed in Sects. 6.4 and 6.5 this picture can change when using different mass-loss rates and mixing efficiencies.

Figure A.2 shows the same diagrams but now for the mass donors of our binary models with an initial orbital period P_i = 10 d and an initial mass ratio q_i = 0.5. For these values of initial orbital period and mass ratio, all systems undergo stable mass transfer and are calculated until core helium depletion. While the amount of mass that is transferred depends on both parameters, the dependence is weak in the sense that in all cases most parts of the hydrogen-rich envelope are removed during the mass-transfer phases.

Fig. A.2. Same as in Fig. A.1, but for our binary models with an initial orbital period of P i = 10 d and an initial mass ratio of q i = 0.5. In contrast to Fig. A.1, the optical depths here mark fixed values of the wind optical depths in the range from τ = 1.5 − 3.5 and are labeled accordingly.

Similar to the single star models, the binary models are expected to have an optically thin wind during most of their main-sequence lifetime. However, during a mass-transfer phase more than 70% of the H-rich envelope is removed. The models shown have a mass-transfer phase during the main-sequence (Case A mass transfer). After the onset of mass transfer the surface hydrogen abundance drops rapidly from X_H = 0.7 to X_H = 0.4, and the wind of the models is predicted to become optically thick as soon as the surface hydrogen abundance drops below X_H ≈ 0.4.

The most notable difference to the single star models is that for most massive primaries we expect them to develop an optically thick wind, and thus a WR-like spectrum, already during their late phases of core hydrogen burning. We would like to note here, that here it can be seen that a lower threshold on the optical depth of τ ≥ 0.5 is expected to strongly impact the H-poor WN population, as it leads to a sustainable longer H-poor WN phase which starts during core hydrogen burning and is able to explain surface hydrogen abundances above X_H ≥ 0.4 (see also Fig. 4).

Appendix B: Comparison with previous works

Evolutionary models of massive stars contain several weakly constrained physics parameters that can strongly affect their evolution and the fate. In the literature various stellar evolutionary models computed with different stellar evolution codes as well as different input parameters, in particular mass-loss recipes and mixing efficiencies, are known and used to explain observed population of stars at different metallicities. Here, we compare our single and binary star models with the predictions of other works.

B.1. Single star models

At first we compare our models to the Geneva (Ekström et al. 2012; Georgy et al. 2012, 2013) stellar evolutionary models which cover a large range of masses and several metallicities. Unfortunately, the recent Geneva models are calculated for solar (Z_⊙ = 0.014), LMC (Z = 0.006), and SMC (Z = 0.002) metallicity. In Fig. B.1 we show the predicted times a model spends as WR star as a function of luminosity, while the tracks are color-coded by the different WR phases predicted by the models.

Fig. B.1. Predicted time a model spends in the different WR phases τWR as a function of luminosity for different initial parameters. The tracks taken from literature cover initial masses M i = 30 M⊙, 40 M⊙, 60 M⊙, and 80 M⊙ (or 85 M⊙ for the Geneva models). Due to the logarithmic spacing of our initial masses our models cover initial masses M i = 32 M⊙, 39.2 M⊙, 63.1 M⊙, and 79.4 M⊙. The tail of each track is labeled with its initial mass. The tracks are colored in orange, green, and blue according to their surface abundances which are used in the later analysis to respectively differentiate the H-poor WN, H-free WN, and WC phase. In panel a and d the nonrotating and rotating Geneva tracks at galactic (dashed) and SMC (dotted) metallicity are compared to our model tracks (solid). Panels b and e show the nonrotating and rotating stellar evolution MIST tracks at galactic (dashed) and LMC (dotted) metallicity, in comparison to our tracks (solid). The last two panels c and f emphasize the differences of the galactic (dashed) and SMC (dotted) of the FRANEC models to our LMC models (solid).

The tracks of the nonrotating Geneva models are shown in panel a of Fig. B.1. In our considered luminosity range their nonrotating SMC models avoid a WR phase at all, while their Galactic models can explain WR stars with luminosities above log(L/L_⊙)≳5.4. Their LMC models are a somewhat fainter than our models with luminosities above log(L/L_⊙)≳5.7. This is mainly because their adopted mass-loss rate during the RSG phase is larger. The Geneva models use the mass-loss rates of de Jager et al. (1988), which is similar to the one of Nieuwenhuijzen & de Jager (1990) used by us, except for luminosities above log(L/L_⊙) > 5.5, where the mass-loss rates of de Jager et al. (1988) become larger by a factor of 2 or more (Mauron & Josselin 2011). Additionally, in the Geneva models the mass-loss rates during the RSG are increased by a factor of 3 if the luminosity in the envelope is five times higher than the Eddington luminosity in the convective envelope. This has a strong impact on the formation of WR stars at different luminosities. The Geneva models use the mass-loss recipe of Nugis & Lamers (2000) with a clumping factor of D = 10 for the WR phase, which makes their WR mass-loss less efficient than the one we adopt. This impacts the predicted times a stellar model can spend in a WR phase.

Results from the Geneva models with rotation can be found in panel d of Fig. B.1. These models drastically differ from the nonrotating case. First, they spend three times more time in a WR phase, most notably as H-poor WN. Second, all Galactic models can evolve now into WC stars with luminosities as low as log L/L_⊙ ≈ 5.4. Third, their most massive SMC model with initial mass 85 M_⊙ now evolves into a H-poor WN type star.

In contrast to the nonrotating models, these models differ strongly from ours. As we argue in Sect. 6.5, their WR lifetimes are exceptionally high because of their efficient rotational mixing, which results in larger convective cores, which leads to an overestimate of the number of the most luminous WR stars.

The second set of stellar evolutionary models we compare to are the MESA Isochrones and Stellar Tracks (MIST) models from Choi et al. (2016). As these models are also calculated with MESA like our models, one might expect relative similar results. However, their nonrotating Galactic and LMC models (panel b) of Fig. B.1) produce a WR phase only for luminosities above log(L/L_⊙) = 6.3. Their models are calculated with a smaller overshooting parameter of 0.2H_P which leads to smaller He-cores, and their treatment of semiconvection is also less efficient than that in our models with α_sc = 0.1. This inefficient mixing leave their models spend four times longer in the cool supergiant regime (T_eff ≤ 12500 K) compared to our models, which explains the scarcity of predicted WR stars.

Their rotating LMC models show similar features, as they still can only explain WR stars with luminosities larger than log(L/L_⊙)≳6.3. It appears, that their Galactic models are in better agreement with our LMC models, as they match roughly the luminosity range, the times spend in the different WR phases.

The third set of single star models with which we want to compare our models to are the models calculated with the FRANEC code from Limongi & Chieffi (2018). As in the case of the Geneva models the FRANEC models only cover Galactic ([Fe/H] = 0) and SMC ([Fe/H]= − 1) metallicity. In contrast to the Geneva models, the nonrotating FRANEC models at SMC metallicity, shown in panel c of Fig. B.1, produce WR stars with luminosities above log(L/L_⊙)≳6.2. Their Galactic models have longer H-poor and H-free WN phases than those from the Geneva code. Overall these models seem to be in better agreement with our models, as they fit into the general trend of their models, while covering the expected luminosity range and WR phases. Contrary to the rotating Geneva models, the fast rotating FRANEC models spend only slightly more time in a WR phase, which appears to be in better agreement with the observations.

We emphasize that in the case of the single stars, the predictions for the WR population depend sensitively on the adopted physics parameters. However, all of the presented single star models struggle explaining the observed WR stars in the LMC with luminosities below log(L/L_⊙)≲5.7.

B.2. Binary star models

In the literature, only few comprehensive sets of detailed binary evolution models are available and the most commonly used ones are the Binary Population and Spectral Synthesis (BPASS) models from Eldridge et al. (2017), covering wide ranges of masses, periods, mass-ratios, and metallicities, including LMC metallicity. A comparison of selected BPASS tracks with our tracks for different initial orbital periods is shown in Fig. B.2.

Fig. B.2. Predicted time a model spends in the different WR phases τWR as a function of luminosity for different initial masses and orbital periods of binary models at LMC metallicity. The BPASS tracks, depicted as dashed lines, cover initial donor masses M1,  i = 30 M⊙, 40 M⊙, 60 M⊙, and 80 M⊙. The solid lines correspond to our binary models which cover initial donor masses M1,  i = 32 M⊙, 39.2 M⊙, 63.1 M⊙, and 79.4 M⊙. Each track is labeled at its end with its initial mass. The tracks are colored in orange, green, and blue according to their surface abundances which are used in the later analysis to respectively differentiate the H-poor WN, H-free WN, and WC phase. All binary models have an initial mass-ratio of q i = 0.7 and the tracks of initial orbital periods P i = 4 d, 40 d, and 400 d in panel a, b, and c, respectively.

By comparing the binary models with initial periods of P_i = 4 d shown in panel a of Fig. B.2, one can see similarities and at the same time strong differences. The donor models with initial masses 30 M_⊙ and 40 M_⊙ evolve similar. They show an agreement with luminosities and, therefore, masses, as well as in the time they spend in the individual WR phases. However, the predictions of the donor models with initial masses 60 M_⊙ and 80 M_⊙ differ strongly, with respect to their luminosity, lifetime, and surface abundances. This can easily be explained, the BPASS models enter during the RLOF a common envelope phase which results in a merger. As a consequence, their models evolve into luminous RSGs (log(L/L_⊙)≥6.5), where a second mass-transfer phase is initiated and the remaining H-rich envelope is stripped, leading to the formation of the appearently “overluminous” WR stars.

The BPASS models with initial periods of P_i = 40 d (panel b) of Fig. B.2) spend similar times as WR stars as for our models. However, the models predict different surface abundances and thus different WR phases at different luminosities. For instance, their model with initial mass 30 M_⊙ spends it entire WR phase as H-poor WN star, while our model with similar initial mass is able to full strip off its H-rich envelope and spends more than >50% of its time as H-free WN star. This is likely linked to the effects of the different WR mass-loss recipes used.

A discrepancy for the different predicted WR subtypes can also be found for the longest period binaries, with initial periods P_i = 400 d (panel c) of Fig. B.2). Here, the lack of very luminous WR stars in our calculations occurs due to mass-transfer rates exceeding the threshold beyond which a common envelope evolution is initiated. In BPASS there is a prescription to model systems undergoing unstable mass transfer and common envelope evolution. Therefore, it is likely that these differences occur due to their different treatment of RLOF.

We see that for close binaries, where BPASS assumes complete mixing during mass transfer and for wide binaries, where BPASS appears to adopt a different merger criterion, the models which are not affected by these two issues agree well with our binary models.

Appendix C: Model overview

Here, we provide three figures which allow to obtain an overview of the different final states of our binary models. For three exemplary initial primary masses (M_1, i = 31.6 M_⊙, 63.1 M_⊙, and 79.4 M_⊙), Figs. C.1, C.2, and C.3 show the covered initial parameter space in orbital period and mass ratio. Each pixel in these plots corresponds to one computed detailed binary evolution sequence and its color, and hatching gives information about its evolution (see below). They further indicate the borderline between Case A and Case B mass transfer (blue dashed line) as well as the borderline between interacting and noninteracting binaries (blue dotted line). The information coded in each pixel is as follows.

• Both dep. C:

  Both stars are calculated until carbon depletion in their core. Due to convergence problems stellar models with helium core masses above M_He > 14 M_⊙ are calculated to core helium depletion instead. When the primary depleted carbon the companion is modeled as a single star, ignoring the possibility of the formation of a binary with compact object due to the uncertainties of a supernova kick.

• L2 overflow:

  During a contact phase both stars overfill their Roche lobe and mass is overflowing the second Lagrangian point L₂. The system is suspected to merge.

• Post MS + Inv. MT:

  Inverse mass transfer from the secondary to a post main-sequence primary. These systems are suspected to have an unstable mass transfer and to merge, leading to the formation of a post main-sequence star.

• Upper Ṁ limit:

  The transferred material can neither be expelled nor accreted anymore (equation 5.3 with R_RL, 2 instead of R₂ in Marchant (2016) has been reached). The transferred material can neither be expelled nor accreted anymore. Therefore, the formation of a common envelope is assumed.

• Lower Ṁ limit:

  The transferred material might no-longer be able to be radiated away from the system (equation 5.3 in Marchant (2016) has been reached). We note that the models reaching this condition continue the calculation. This label has to be understood as a warning that radiation may not be able to drive a wind to expel the transferred material from the system.

• MT with q_i < 0.25:

  It is assumed that a mass-transfer phase for systems with extreme mass ratios below q_i < 0.25 will fill the Roche lobe of the secondary quickly leading to a common envelope phase and a potential merger.

• max MT limit:

  This is the additional label mentioned above. Binary systems that trigger this condition have a mass-transfer rate that exceeds 0.1 M_⊙ yr⁻¹. These systems are suggested to form a common envelope.

• convergence error:

  The calculation has stopped due to numerical issues and did not converge. This likely happens in the late phases of the evolution or late stages of mass transfer.

• Had contact:

  Both stars had a contact phase. We note that this does not imply that this system results in a merger.

• Not-int. boundary:

  Systems above this line are not interacting during their entire life. They evolve as single stars.

• Case A/B boundary:

  Boundary line between systems undergoing Case A and Case B mass transfer.

Fig. C.1. Phase diagram of binary models with a fixed initial donor mass of M1, i = 31.6 M⊙. The initial mass ratio q i is plotted versus the initial orbital period P i. The different phases are indicated in the legend and a more detailed description can be found in the text.

Fig. C.2. Same as Fig. C.1, but for initial donor mass of M1, i = 63.1 M⊙.

Fig. C.3. Same as Fig. C.1, but for initial donor mass of M1, i = 79.4 M⊙.

To get a better understanding of the implications of the phase diagram the most characteristic regions of Fig. C.1 will be described in more detail. Binary systems with initial orbital periods beneath the “Case A/B” boundary already have a mass transfer phase during their time on the main-sequence. One can see that most of the systems with very tight orbits experience a contact phase. During this phase a large fraction of the donors envelope is transferred and accreted, leading to a mass ratio close to q = 1 after the mass-transfer phase. Therefore, it is possible for most systems that either the secondary fills its Roche lobe and inverse mass transfer sets in, or mass is overflowing through the second Lagrangian point during the contact phase. Both scenarios are suspected to result in a merger. For systems with initially more extreme mass ratios q_i → 0, the accretor is not luminous enough so that infalling material cannot be transferred to infinity by a radiation driven wind. These systems are also suspected to form a common envelope and merge. Only a certain fraction of the binary systems with wide enough orbits and large enough mass ratios can have a successful mass-transfer phase where both stars can be modeled until carbon depletion. It is worth mentioning that the accretion efficiency is larger for case A systems and decreases with increasing initial orbital period whereas the mass-transfer rate grows with increasing initial orbital period.

Because of the effect of envelope inflation near the Eddington limit the most massive stars of our grid expand during their time on the main-sequence. This shifts the Case A/B boundary to higher initial orbital periods for increasing initial masses. The Case A/B boundary reaches its maximum at M_1, i ≈ 63.1 M_⊙, where almost all systems undergo Case A mass transfer. For higher initial masses we predict the Case A/B boundary to go down to lower initial orbital periods again, as the winds become efficient enough to remove large parts of the envelope, allowing the donors to remain more compact during their main-sequence evolution (cf. Fig. 1).

Binary systems that have a mass-transfer phase while they evolve off the main-sequence are located between the Case A/B boundary and the “not-interacting” boundary. These systems are too wide to have a contact phase and only binary models with initial mass ratios close to q_i ≃ 1 are expected to have an inverse mass-transfer phase. Models with extreme mass ratios q_i → 0 are still unable to drive a wind that transfers the infalling material to infinity as they reach the upper limit on mass transfer. For the systems with the widest orbits (P_i ≳ 300 d), the maximum value for mass transfer is triggered. These systems are expected to have donors that have evolved into a red supergiant phase and have large convective envelopes which drastically increase the mass transfer. The stability of mass transfer of these systems is questionable and in this work it is assumed that these systems form a common envelope. Only a certain fraction of binary systems is suspected to be able to evolve both stars until core carbon depletion.

Binary systems with an initial orbital period larger than the not-interacting boundary (P_i ≳ 3000 d) evolve as single stars. As the stars are tidally synchronized their initial rotation is rather small and they can be considered as quasi nonrotating single stars.

In total less than < 5% of our systems encounter “convergence errors” of different kinds. By inspecting Figs. C.1, C.2, and C.3 one can see that the fraction of systems encountering these convergence errors increases with increasing initial donor mass. At the highest masses (M_ini,  1 = 89) the fraction goes up to 30%.

At higher initial donor masses some of our primary models have numerical issues, typically during mass transfer. For instance, when the primary already evolved to a RSG (log(P/d)≳2.5 it has a big convective envelope which rapidly expands, leading to too small time steps and thus to convergence errors. Other models in this region of the phase diagram end their evolution under the “max MT limit”. We believe that most of the systems, which have numerical issues during mass transfer, would have triggered one of the criteria associated with a following common envelope evolution (e.g., convergence errors or “Post MS + Inv. MT”) anyways, meaning that these systems are not expected to contribute to our synthetic WR population even if they would not have encountered the convergence errors.

With increasing initial donor mass also the masses of our secondaries increase. The secondary, which is evolved as a single star after the primary has depleted helium in its core, evolves in most cases into a RSG star. During this late evolutionary stage it can be in some cases that these models encounter numerical instabilities and trigger the convergence errors flag. Implying that for these systems the primary who predominantly contributes to the WR population is fully modeled and included in our synthetic population. Consequently, the true fraction of systems that encounter numerical issues and are neglected in our final synthetic WR population is thought to be below < 1%.

Appendix D: Tables with the observed data of WN and WC stars in the LMC

All Tables

All Figures

	Fig. 1. Evolutionary tracks of selected single star models in the Hertzsprung-Russell diagram. The tracks are labeled by the initial stellar mass and color-coded by the surface helium abundance XHe (see color bar to the right). A star symbol (also color-coded), marks the point of core helium depletion. Lines of constant radius are drawn as thin straight lines.
In the text

Fig. 2. Evolutionary tracks of mass donors (panels a–d) and of mass gainers (panels e and f) for selected binary models in the Hertzsprung-Russell diagram. The tracks are labeled with the selected initial stellar masses. The mass gainers shown in panels e and f are the counterparts to the mass donors shown in panels a and b, respectively. The panels a–c depict the donor evolution for three different initial orbital periods, for binaries with a fixed initial mass ratio of qi = 0.7. Panels b and d: donors for the same orbital period (Pi = 40 d) but for two different initial mass ratios. In the panel a diagram, arrows identify different parts of the same tracks. The tracks and star symbols are color-coded as in Fig. 1. Double rings marks a predicted merger situation; our calculations end there.

In the text

Fig. 3. Luminosity distribution of observed H-poor WN (orange), H-free WN (green), and WC (blue) star populations (colored shading) with our synthetic WR population (black histogram, τ ≥ 1.5). The synthetic WR population distinguishes WR models obtained from single star (dotted bins) and binary models (dashed bins). The observed H-poor WN population excludes stars with spectral types Of and Of/WN. In the upper right corner of each plot, we give the number of observed stars, the numbers of predicted single stars and binary produced WR stars of the respective WR subtype. For the synthetic H-poor WN populations an additional histogram (brown) for lower optical depths (τ ≥ 0.5) is shown.

In the text

	Fig. 4. Surface hydrogen abundances as a function of luminosity of the observed H-poor WN stars (orange dots) and of our synthetic population (histogram). Black lines indicate borderlines of the theoretical distributions obtained for different optical depth thresholds (τ ≥ 0.5, 1.0, and 1.5).
In the text

Fig. 5. Predicted mass-ratio and period distributions of the H-poor WN stars. Upper left: predicted mass-ratio distribution of the H-poor WN stars in our synthetic population with optical depths τ ≥ 0.5 (brown dashed histogram) and τ ≥ 1.5 (black dashed histogram), where we distinguish WR stars during core hydrogen and core helium burning. Overlayed is the histogram of the observed mass ratios in LMC WN stars with hydrogen (orange). Upper right: same as in the upper left panel, but for the orbital periods. Lower left: 2D histogram of mass ratios and orbital periods of our synthetic population with optical depths τ ≥ 0.5, represented by colored squares, where the color indicates the expected number of objects (see color bar to the right). Individual observed binaries are shown as orange dots. Lower right: same as the lower left panel but now for optical depths τ ≥ 1.5. The observational data is take from Shenar et al. (2019).

In the text

	Fig. 6. Same as Fig. 5, but for H-free WN stars. In the right panel, the 2D histogram, systems with unknown mass ratio but known orbital period are indicated by arrows. The observational data is take from Shenar et al. (2019).
In the text

	Fig. 7. Same as Fig. 5, but for WC stars. The observational data is take from Bartzakos et al. (2001b).
In the text

Fig. A.1. Optical depth (left, color-coded) and surface hydrogen abundance (right) of our single star models as function of initial mass and normalized time, during their main-sequence (lower) and core helium burning (upper) evolution. The transition line from an optically thin to an optically thick wind is marked by a solid black line labeled with “1.500”. The unlabeled solid black lines divide the different considered WR phases as labeled in each region. In the right plot, threshold values of the surface hydrogen abundances in steps of ΔXH = 0.1 are indicated by thin solid lines, which have dark red color for the range XH = 0.8 − 0.5 and brown color for XH = 0.4 − 0.1. The borderlines of XH = 0.4 and XH = 0.0 are marked by solid black lines, as well as the borderline beyond which stellar models have a surface carbon mass fraction in excess of XC = 0.01, which is our threshold for defining WC stars.

In the text

	Fig. A.2. Same as in Fig. A.1, but for our binary models with an initial orbital period of P i = 10 d and an initial mass ratio of q i = 0.5. In contrast to Fig. A.1, the optical depths here mark fixed values of the wind optical depths in the range from τ = 1.5 − 3.5 and are labeled accordingly.
In the text

Fig. B.1. Predicted time a model spends in the different WR phases τWR as a function of luminosity for different initial parameters. The tracks taken from literature cover initial masses M i = 30 M⊙, 40 M⊙, 60 M⊙, and 80 M⊙ (or 85 M⊙ for the Geneva models). Due to the logarithmic spacing of our initial masses our models cover initial masses M i = 32 M⊙, 39.2 M⊙, 63.1 M⊙, and 79.4 M⊙. The tail of each track is labeled with its initial mass. The tracks are colored in orange, green, and blue according to their surface abundances which are used in the later analysis to respectively differentiate the H-poor WN, H-free WN, and WC phase. In panel a and d the nonrotating and rotating Geneva tracks at galactic (dashed) and SMC (dotted) metallicity are compared to our model tracks (solid). Panels b and e show the nonrotating and rotating stellar evolution MIST tracks at galactic (dashed) and LMC (dotted) metallicity, in comparison to our tracks (solid). The last two panels c and f emphasize the differences of the galactic (dashed) and SMC (dotted) of the FRANEC models to our LMC models (solid).

In the text

Fig. B.2. Predicted time a model spends in the different WR phases τWR as a function of luminosity for different initial masses and orbital periods of binary models at LMC metallicity. The BPASS tracks, depicted as dashed lines, cover initial donor masses M1,  i = 30 M⊙, 40 M⊙, 60 M⊙, and 80 M⊙. The solid lines correspond to our binary models which cover initial donor masses M1,  i = 32 M⊙, 39.2 M⊙, 63.1 M⊙, and 79.4 M⊙. Each track is labeled at its end with its initial mass. The tracks are colored in orange, green, and blue according to their surface abundances which are used in the later analysis to respectively differentiate the H-poor WN, H-free WN, and WC phase. All binary models have an initial mass-ratio of q i = 0.7 and the tracks of initial orbital periods P i = 4 d, 40 d, and 400 d in panel a, b, and c, respectively.

In the text

	Fig. C.1. Phase diagram of binary models with a fixed initial donor mass of M1, i = 31.6 M⊙. The initial mass ratio q i is plotted versus the initial orbital period P i. The different phases are indicated in the legend and a more detailed description can be found in the text.
In the text

	Fig. C.2. Same as Fig. C.1, but for initial donor mass of M1, i = 63.1 M⊙.
In the text

	Fig. C.3. Same as Fig. C.1, but for initial donor mass of M1, i = 79.4 M⊙.
In the text