© D. R. Aguilera-Dena et al. 2022

1. Introduction

Many open questions in contemporary astrophysics are intertwined with the evolution and fate of massive stars. Massive stars shape their environments through the many feedback processes they experience, and they are the main source of chemical enrichment throughout cosmic history. The first generations of massive stars are believed to have contributed significantly to the re-ionisation of the Universe (e.g. Götberg et al. 2020). The populations of massive stars in the more chemically enriched local Universe have a complex morphology (e.g. Evans et al. 2011; Pantaleoni González et al. 2021), and they are responsible for the formation of many astrophysical transients we observe (e.g. Modjaz et al. 2019).

At the end of their evolution, massive stars can explode as core-collapse supernovae (SNe), leaving behind a neutron star (NS) or black hole (BH) remnant. Massive stars that do not produce successful explosions collapse, leaving behind BHs. Very massive stars can undergo pair instability, and above a certain mass they explode without leaving a remnant behind (Langer 2012). They input a large amount of mechanical and radiative energy in their environments, both through stellar winds on relatively long timescales, and through SNe on shorter timescales. As successive generations of stars are formed, their initial composition becomes more enriched in heavy elements produced by their predecessors. Understanding how their composition influences their evolution, energy output, and final fate is therefore important for understanding galaxy evolution, and for understanding the properties of stellar populations, astrophysical transients such as SNe (e.g. Modjaz et al. 2019), and compact object mergers (e.g. Vigna-Gómez et al. 2018) throughout the history of the Universe.

Most massive stars in the Universe are thought to be in stellar binaries (e.g. Sana & Evans 2011; Sana et al. 2012) or multiples (e.g. Sana et al. 2014; Moe & Di Stefano 2017) that will interact at some point during their lifetime. Our understanding of the evolution of massive stars, and the role that multiplicity plays in their evolution in different environments, has taken grand leaps in the past decades (e.g. Podsiadlowski et al. 1992; Yoon et al. 2010; Langer 2012; Schneider et al. 2021; Laplace et al. 2021). However, many gaps remain in our understanding of the effect of metallicity and multiplicity on their evolution and fate.

Simulations of interacting binary systems indicate that at least one component star is likely to lose most or all of its hydrogen-rich envelope (e.g. Podsiadlowski et al. 1992; Wellstein et al. 2001; Yoon et al. 2017; Gilkis et al. 2019; Laplace et al. 2020; Klencki et al. 2020; Vigna-Gómez et al. 2022). This results in the formation of stripped-envelope stars, characterised by a surface that is hydrogen deficient and helium enriched.

Envelope stripping can occur due to mass transfer, for example during the main sequence (MS) evolution (Case A mass transfer; see Sen et al. 2022) or following the depletion of hydrogen or helium in the core (Case B and C mass transfer, respectively; see Tauris & van den Heuvel 2006, for a review). Envelope ejection can also occur on dynamical timescales when a binary system undergoes a common envelope (CE) event that is often accompanied by a significant tightening of the binary orbit (Tauris & Dewi 2001; Ivanova 2011; Kruckow et al. 2016). Furthermore, the most massive and luminous stars may lose their hydrogen envelopes without interacting with a companion star, through winds during the MS evolution (Conti 1975), winds during the red supergiant phase (Maeder 1981), or massive explosive outbursts, when their outer layers are near the Eddington limit (e.g. Langer et al. 1994).

Several formation mechanisms of stripped-envelope stars may contribute to their numbers in the Universe, but their efficiency is uncertain (Shenar et al. 2020). Regardless of the formation mechanism, stripped-envelope stars will have different surface properties and evolve differently than their non-stripped counterparts (Woosley 2019; Schneider et al. 2021; Laplace et al. 2021).

Some stripped-envelope stars will become classical Wolf-Rayet (WR) stars (Groh et al. 2014), a class of massive stars that is distinguished by spectra with strong emission lines, a product of their high mass loss rates that lead to an optically thick wind above the stellar photosphere, and that show little or no hydrogen (Crowther 2007). WR stars are observed in several sub-types that are indicative of their surface composition, temperature, and evolutionary stage. The broadest classification divides them into WN, WC, and WO types. WN-type stars are characterised by the presence of nitrogen in their spectra, while WC stars have strong carbon and oxygen lines. Stars with these spectral types are believed to be undergoing helium burning in their core. WO stars have spectra dominated by oxygen lines, and are believed to have helium-depleted cores, which places them closer to the end of their evolution (Tramper et al. 2015). The mass loss rates of WR stars have been found to be several orders of magnitude higher than those of main sequence stars, and they are strongly metallicity dependent (e.g. Gräfener & Hamann 2005, 2008; Vink & de Koter 2005; Hamann et al. 2006; Hainich et al. 2014; Tramper et al. 2015).

The luminosities of WR stars have been found to have a lower limit, below which no WR stars are observed. This lower luminosity has been observed to be metallicity dependent. Understanding the formation channel of stars above and below this limit is not trivial, as they may originate both from single wind-stripped stars, and from stars that have been stripped via an interaction with a companion (Shenar et al. 2020). Theoretically, optically thick WR winds are expected to form above a given minimum luminosity, which depends on their metallicity (Gräfener et al. 2017; Sander & Vink 2020). Below this limit, the mass loss rates of stripped-envelope stars are expected to drop substantially (Vink 2017), leading to a strong reduction or disappearance of their WR line emission. While this drop in wind mass loss rate still needs to be confirmed observationally, a similar effect is expected for continuous mass loss relations, where WR winds gradually become optically thin for lower masses. Götberg et al. (2018) have shown that, due to this effect, low mass stripped-envelope stars emit most of their radiation in the UV, and are therefore difficult to detect. The study of WR stars through direct observation has been a very active field since their discovery (Neugent & Massey 2019), but the study of stripped-envelope stars below the WR luminosity limit is a field that is expected to grow with upcoming observational campaigns. However, detailed studies are limited to nearby galaxies where systems can be observed individually.

In a recent study, Yoon (2017) found that an increase in the mass loss rate observed in WC stars could account for the formation of the faintest WC- and WO-type stars in the LMC and in our Galaxy. Wolf-Rayet winds are known to have a strong metallicity dependence (Hainich et al. 2014; Tramper et al. 2015), which strongly impacts their populations in different galaxies.

With this in mind, we study the evolution of helium stars as a proxy for stripped-envelope stars, including wind mass loss, in a similar fashion to Yoon (2017) and Woosley (2019). We extend the considered range of initial conditions to cover different metallicity environments, focusing on high metallicities, which are not often addressed in the literature, and can explain several of the observed population properties of WR stars and Type I core-collapse SN (discussed in Aguilera-Dena et al. 2022, hereafter Paper II).

We pay particular attention to the empirical metallicity-dependent lower luminosity limit of WR stars (Shenar et al. 2020), apparent when comparing the WR populations in the Small and Large Magellanic Clouds (e.g. Massey et al. 2003, 2014; Neugent et al. 2018) and the Galaxy (e.g. Hamann et al. 2019). This limit is likely related to the predicted drop in mass loss rate for lower luminosity helium stars (Gräfener et al. 2017; Vink 2017; Sander & Vink 2020).

We have organised our paper as follows. In Sect. 2 we describe the physical and numerical treatment we use in our simulations, and propose a method for finding the minimum luminosity of WN and WC stars as a function of metallicity. In Sect. 3 we describe how metallicity affects the time-dependent properties of our stripped-envelope star models. In Sect. 4 we contrast the results of our numerical models to the inferred minimum luminosities of WR stars. In Sect. 5 we show the implications of our results for the morphology of WR populations. We conclude our paper with a discussion in Sect. 6 and conclusions in Sect. 7. In a companion paper (Paper II) we discuss the properties of our stripped-envelope stellar models at core collapse, and present an analysis of the transients and the remnants we expect from them. Additionally, some of the ramifications for compact-object mergers are also discussed in Antoniadis et al. (2022).

2. Method

In this paper we followed several approaches to study the effect of metallicity on the evolution of stripped-envelope stars, all of which are detailed below. In Sect. 2.1 we describe the method employed to construct a grid of stellar evolution models of helium stars, meant to approximate the evolution of stripped-envelope stars. In Sect. 2.2 we present a semi-analytic method we devised to estimate the wind optical depth of stripped-envelope stars. We then use this criterion (in Sect. 4) to classify our models as either WRs (with optically thick winds) or transparent-wind helium stars. Finally, we combine the models and the semi-analytical method to produce synthetic populations of stripped-envelope stars in different metallicity environments through the method detailed in Sect. 2.3.

To simplify the presentation of our results, we introduce the notation used in this paper for various useful quantities in Table 1.

2.1. Stellar evolution models

We computed grids of evolutionary sequences of non-rotating helium stars using version 10398 of the Modules for Experiments in Stellar Astrophysics (MESA) code (Paxton et al. 2011, 2013, 2015, 2018), from the beginning of core helium burning to the onset of core collapse (when possible), defined as the point where the infall velocity at the edge of the iron core reaches 1000 km s⁻¹.

To test the effect of initial chemical composition in hydrogen depleted stripped-envelope stars, we performed calculations for seven different metallicities (Z_init), from 0.01 to 0.04 in steps of 0.005, scaled from the solar abundances found by Grevesse et al. (1996), and using initial masses between 1.5 and 70 M_⊙ in steps of 0.5 M_⊙, which correspond to the helium core masses of stars with zero age main sequence (ZAMS) masses between roughly 10 and 150 M_⊙ (Woosley 2019, see Eq. (9)).

The initial models were generated from hydrogen-rich pre-main sequence models at the aforementioned metallicities. They were evolved imposing chemical homogeneity by means of an artificially imposed, very high mixing coefficient throughout the star, and without mass loss, until the end of hydrogen burning, but before the ignition of the ¹⁴N(α,γ)¹⁸F(β⁺,ν)¹⁸O(α,γ)²²Ne reaction in the core, a reaction that takes place on a short timescale and at a lower temperature than helium burning, and depletes the core of nitrogen. This was done to ensure that nitrogen was not depleted in the stellar surface at the beginning of the calculation. The models were evolved further until they settled into thermal equilibrium from helium burning, relaxing the condition of homogeneity. This provides the abundances of CNO isotopes in the envelope that correspond to those of the core of a massive star after hydrogen burning (enhanced N, reduced C and O), and provides a helium content given by Y_init = 1 − Z_init.

Although binary evolution models typically retain residual hydrogen in their envelope after having been stripped due to Roche lobe overflow (Yoon et al. 2017; Gilkis et al. 2019; Laplace et al. 2020), we assume that the remaining amount of hydrogen is small, and will be removed on a short timescale by winds or case B mass transfer, and therefore does not significantly affect the evolution of helium stars (see Sect. 6.1.2 for a discussion).

Convection was modelled using the standard mixing length theory (Böhm-Vitense 1958), using the default MESA value of α_MLT = 2.0, and adopting the Ledoux criterion for instability. We employed efficient semi-convection, using α_SC = 1.0, following the results of Schootemeijer et al. (2019). Convective overshooting was not included during any of the phases of evolution. We included thermohaline mixing, as prescribed by Cantiello & Langer (2010). To improve the numerical stability of our calculations, we employed MESA’s predictive mixing in helium-burning shells (Paxton et al. 2015). We calculated the energy generation rates and chemical composition changes using MESA’s approx21 nuclear network, from the initial model to core collapse.

Since helium stars often experience instabilities in their envelopes due to the proximity to the Eddington limit (e.g. Sanyal et al. 2015), our models were calculated using efficient energy transport through MESA’s mlt++ (Paxton et al. 2013), effectively increasing the efficiency of convection to transport energy in convective layers in the envelope. In some cases, in order to produce models that converge up to core collapse, we exclude radiative acceleration in the envelope by setting the radial velocity to 0 in layers with T < 10⁸ K during the late evolution when time steps become smaller than 0.1 yr. Additional numerical details are available online¹.

To compute the mass loss that our models experience during their evolution, we use the mass loss recipes suggested by Yoon (2017), which prescribe different mass loss rates for WN- and WC-type stars. These recipes are based on the work of Hainich et al. (2014) for WN stars, using the metallicity dependence obtained by Hamann et al. (2006) and Tramper et al. (2016) for WC and WO stars. The mass loss rate recipe employed has the form

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{WN}}=f_{\mathrm{WR}}{\left(\frac{L}{L_{\odot }}\right)}^{1.18}{\left(\frac{Z_{\mathrm{init}}}{Z_{\odot }}\right)}^{0.6}{10}^{-11.32}\frac{M_{\odot }}{\mathrm{yr}}\end{array}$$(1)

for Y = 1 − Z_init, while for Y < 0.9 it is given by

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{WC}}=f_{\mathrm{WR}}{\left(\frac{L}{L_{\odot }}\right)}^{0.85}{\left(\frac{Z_{\mathrm{init}}}{Z_{\odot }}\right)}^{0.25}{Y}^{0.85}{10}^{-9.2}\frac{M_{\odot }}{\mathrm{yr}}·\end{array}$$(2)

As suggested by Yoon (2017), we interpolate between them linearly in the regime between Y = 1 − Z_init and Y < 0.9 using

$$\begin{array}{c}\hfill \dot{M}=(1-x){\dot{M}}_{\mathrm{WN}}+x{\dot{M}}_{\mathrm{WC}},\end{array}$$(3)

with x = (1 − Z_init − Y)/(1 − Z_init − 0.9). We take the recommended value of f_WR = 1.58, corresponding to a clumping factor of D = 4, and scale these equations assuming Z_⊙ = 0.02. These wind mass loss rates for WR stars were determined empirically, and are likely to be improved in the future. Some of the helium stars we model are outside of the observed luminosity range where WR stars have been observed. We are therefore effectively extrapolating WR mass loss rates to lower luminosities, which may have lower mass loss rates (Vink 2017). We take this into account a posteriori during the analysis of our results.

2.2. Stellar wind optical depth

Shenar et al. (2020) found the minimum luminosity of WR stars in the LMC, the SMC, and the Galaxy, and propose an analytical expression for this quantity as a function of metallicity. They derive it by means of the ‘transformed radius’ R_t (Schmutz et al. 1989), which depends on the stellar radius, terminal wind velocity, mass loss rate, and wind clumping factor. They assume that the mass loss rate depends on luminosity and metallicity with a functional form given by Ṁ ∝ L^αZ^β, and that all other quantities are independent of either L or Z. Finally, they find the minimum luminosity of WR stars by assuming that there is a metallicity-independent value of R_t above which stripped-envelope stars will have WR spectra. They determine that this minimum luminosity is proportional to Z^{4α/(3−4β)}. Using values for α and β derived by Vink (2017) in the context of low mass helium stars, with luminosities below the minimum observed WN luminosity, they find a power-law relation between the luminosity threshold and Z with an exponent roughly equal to −1, which provides a reasonable fit to the observed values.

Here, we take a similar approach, but instead calculate the optical depth at the base of the winds of WR stars, τ(R), using the analytical approach of Langer (1989). The minimum luminosity of WR stars should then be determined by the minimum luminosity at which WR winds are optically thick (i.e. have τ(R)≳1). In this formalism, for a wind velocity law with β = 1 (Eq. (8) in Langer 1989) the optical depth at the base of the wind is given by

$$\begin{array}{c}\hfill \tau (R)=\frac{\kappa |\dot{M}|}{4\pi R(v_{\infty }-v_0)}ln\frac{v_{\infty }}{v_0},\end{array}$$(4)

where κ is the opacity, R the stellar radius, and v_∞ and v₀ are respectively the terminal wind velocity and the velocity at the base of the wind. We then combine this with the metallicity- and luminosity-dependent wind mass loss rates for both WN- and WC-type stars with Eqs. (1) and (2), thereby extending the scope of Shenar et al. (2020) to also include WC-type stars.

We assume that the terminal wind velocity is given by the escape velocity at the surface, multiplied by a factor K (Gräfener et al. 2017):

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

Here Γ is the Eddington factor. Substituting these values in Eq. (4) yields

$$\begin{array}{cc}\hfill \tau (L,Z)=& \frac{\kappa |\dot{M}(L,Z)|}{4\pi R(L)}{(K\sqrt{\frac{2GM(L)}{R(L)}(1-\frac{\kappa L}{4\pi GcM(L)})}-v_0)}^{-1}\hfill \\ \hfill & \times ln\left(\frac{1.3\sqrt{\frac{2GM(L)}{R(L)}(1-\frac{\kappa L}{4\pi GcM(L)})}}{v_0}\right).\hfill \end{array}$$(6)

We employ the mass-luminosity and radius-luminosity relations from Langer (1989), assume that the opacity at the base of the wind is given by the electron scattering opacity, κ = 0.2 cm² g⁻¹, and that v₀ = 20 km s⁻¹, representative of the sound speed in the envelopes of WR stars. After making these assumptions, the optical depth at the base of the wind of a WR star is only a function of luminosity and metallicity, and it is different for WN- and WC-type stars.

The result of this calculation assuming the mass loss rate of WN stars (Eq. (1)), and using K = 1.3 (Gräfener et al. 2011) is presented in Fig. 1. The dash-dotted line corresponds to the region where τ = 1. This line separates the optically thick region (green), where we expect stripped-envelope stars to have WR spectra, and the optically thin region (purple). Since this model relies on several simplifying assumptions, and the uncertainties they introduce cannot be accounted for in detail, we find the values of the minimum luminosities of each type of WR by fitting a line of constant optical depth to the observed lower luminosity values of WR stars as a function of metallicity. For WN-type stars, we calibrate our calculation to the observed values of the minimum luminosity of WN stars in the SMC, LMC, and the Galaxy from Shenar et al. (2020).

Fig. 1. Optical depth at the base of the wind of a stripped-envelope star as a function of its metallicity and luminosity, obtained from Eq. (6) (based on the optical depth derivation, mass-luminosity, and radius-luminosity relations of Langer 1989). The mass loss rate of WN stars is used, which depends on metallicity and luminosity, as in Eq. (1). The star symbols represent the lowest luminosity observed for WN stars in the SMC, LMC, and Milky Way, according to Shenar et al. (2020). Error bars with a width of 0.05 dex are included as these values are rounded. The lines represent the solutions of Eq. (6) for constant optical depth, including those that best fit the observed values of the observed minimum luminosities of WN stars.

We find that the observed minimum luminosities of WN stars lie at optical depths greater than 1 but of order unity, and conclude from this that our approach can be used to estimate the threshold WN luminosity as a function of metallicity. We find that the curve described by setting τ = 1.45 provides the best fit to the data. This is equivalent to changing the value of any of the constants whose value we (somewhat arbitrarily) set in Eq. (6) to make the curve with τ = 1 coincide with the data. We highlight that neither approach is preferable to the other, but we argue that by setting reasonable values for the unknown quantities and finding relative good agreement with our expectations justifies the validity of this method.

We then find the minimum luminosity for WC stars as a function of metallicity by assuming they are described by a curve with the same optical depth, but with the mass loss rate given by Eq. (2), and a value of K = 1.6 (Gräfener et al. 2011). The values of the minimum luminosities of WN and WC stars found using this method, labelled $L_{\min ,\mathrm{WN}}^{\mathrm{tau}}$ and $L_{\min ,\mathrm{WC}}^{\mathrm{tau}}$, respectively, are shown in Fig. 2.

Fig. 2. Minimum luminosity at which stripped-envelope stars can be observed as WN-type (green) and WC type WC-type (red), as a function of metallicity, according to Eqs. (7) and (8), respectively.

Figure 1 shows that the curves of constant optical depth are well described by power laws between the luminosity and the metallicity. For WN-type stars we find that the best fit to the observations is given by

$$\begin{array}{c}\hfill L_{\min ,\mathrm{WN}}^{\mathrm{tau}}=6.85\times {10}^4{\left(\frac{Z}{0.02}\right)}^{-0.71}{L}_{\odot },\end{array}$$(7)

and for WC-type stars the curve with τ = 1.45 is described by

$$\begin{array}{c}\hfill L_{\min ,\mathrm{WC}}^{\mathrm{tau}}=7.51\times {10}^3{\left(\frac{Z}{0.02}\right)}^{-1.55}{L}_{\odot }.\end{array}$$(8)

We find that curves of constant optical depth provide a better fit to the observed minimum WN-type star luminosities than a power law with slope of −1. Furthermore, by extending this analysis to WC stars, we are also able to compare the minimum luminosity of WC stars predicted by our optical depth model to the minimum luminosity at which WC stars are observed, as well as the minimum luminosity of the carbon-rich models that are produced in our grids.

The proposed method used to find the minimum luminosities of WR stars has been employed by Pauli et al. (in prep.), who find a similar calibration, and extend it to classify H-rich WR stars, obtained in binary evolution models. Their binary models are made at LMC metallicity, and they find that the number and the luminosity distribution of hydrogen-rich WN stars, hydrogen-poor WN stars, and WC stars are well reproduced by their models. They also predict core hydrogen burning stars that develop optically thick winds. Models in their grid that expose the products of helium burning are found to always be luminous enough to be classified as WC stars. Our method is in agreement with their findings as no carbon-rich stars with transparent winds are predicted by our models (see Sect. 4). We note that this method forv distinguishing WR stars from low mass stripped-envelope stars with transparent winds is agnostic to whether the mass loss rates change abruptly between one type and the other or not.

2.3. Population synthesis models

Using the models described in Sect. 2.1 and the semi-analytical method described in Sect. 2.2, we constructed simple synthetic populations of stripped-envelope stars in different metallicities. We distinguish between models with surface optical depth larger or smaller than τ = 1.45, classifying the former as WR stars and the latter as transparent-wind stripped-envelope stars. We classify a model as a WN star if it has τ > 1.45 and shows nitrogen on the surface. If surface nitrogen is depleted and surface carbon enhanced, it is classified as WC. When a model has a surface oxygen abundance X_O > 0.05 and the model has completed core helium burning, it is classified as a WO star (Tramper et al. 2015).

We use these model stellar populations to estimate how the number of stars (WR stars of the different types we consider and transparent-wind stripped-envelope stars) changes as a function of metallicity. We also compute how the average number ratios of WC to WN stars, ⟨N_WC/N_WN⟩, WO to WN stars, ⟨N_WO/N_WN⟩, and WO to WC stars, ⟨N_WO/N_WC⟩ vary as a function of metallicity.

We construct the stellar populations assuming a constant star formation rate and an initial mass function (IMF) given by $\xi (M_{\mathrm{ZAMS}})\propto M_{\mathrm{ZAMS}}^{-2.35}$ (Salpeter 1955). For each metallicity we draw 5000 stars with masses given by the IMF every 1000 yr, for a total of 60 Myr, totalling 3 × 10⁷ stars. We assume that the ZAMS mass of a star that produces a helium core of a given mass (which we later assume becomes a helium star) is not a strong function of metallicity, but only of the initial helium star mass. We employ the formula of Woosley (2019), given by

$$\begin{array}{c}\hfill \frac{M_{\mathrm{He},\mathrm{ini}}}{M_{\odot }}\approx \{\begin{array}{cc}0.0385M_{\mathrm{ZAMS}}^{1.603},\hfill & \mathrm{if}{M}_{\mathrm{ZAMS}}<30M_{\odot }\hfill \\ 0.5M_{\mathrm{ZAMS}}-5.87,\hfill & \mathrm{if}{M}_{\mathrm{ZAMS}}\ge 30M_{\odot }.\hfill \end{array}\end{array}$$(9)

We then assume that the ZAMS lifetime of these stars is given by their nuclear timescale, and that all stars are immediately stripped at the beginning of core helium burning, and evolve as single helium stars. We use the values of t_He − N, t_He − C, and t_He − O found from our models to calculate the number of each type of stripped-envelope star as a function of time.

The maximum ZAMS mass in our populations is taken as 151.74 M_⊙, corresponding to the ZAMS mass of our most massive helium star model, with 70 M_⊙. Due to the slope of the IMF, these results are not sensitive to the exact value of the maximum mass as long as it is high. The minimum mass for the transparent-wind stripped-envelope stars that are tracked in our populations is 3 M_⊙, which corresponds to a ZAMS mass of approximately 15.14 M_⊙. This is chosen to reflect the minimum mass of stripped-envelope stars that undergo core collapse (Woosley 2019; Chanlaridis et al. 2022).

3. Impact of metallicity on the evolution of stripped-envelope stars

In this section we study the effect of metallicity on the evolution of stripped-envelope stars by analysing the grids of evolutionary calculations described in Sect. 2.1. We have subdivided this section as follows. In Sect. 3.1 we review the evolution of helium star models of different masses and metallicities with a few characteristic examples. In Sect. 3.2 we characterise their lifetimes and the duration of their different evolutionary phases. Finally, in Sect. 3.3 we present the luminosity ranges at which carbon- and oxygen-rich stellar models are produced as a function of metallicity.

3.1. Representative examples

Helium stars that are massive enough to burn helium in their core through the triple alpha process begin their core helium burning lifetime with a convective core and a predominantly radiative envelope. The latter is composed mostly of helium and nitrogen at this stage. Its carbon, nitrogen, and oxygen abundances are set during hydrogen burning by CNO equilibrium. The core is mostly composed of helium, but is enriched with neon and depleted of nitrogen since nitrogen burns at a lower temperature than helium. As helium burning progresses the convective helium-burning core becomes increasingly rich in carbon, but also becomes enriched in oxygen through alpha captures on carbon. The total amount of oxygen produced from this reaction remains uncertain, however, because this reaction rate in particular is poorly constrained (e.g. Farmer et al. 2020).

In our models the transition period between the phases where a nitrogen-rich envelope is present and the carbon-rich layers are exposed is short. Nitrogen is thereafter absent from their surface, replaced by the products of helium burning. As mass loss continues to uncover deeper layers, the amount of carbon in the surface decreases and the amount of oxygen increases since deeper layers of the stars spent a longer time exposed to the ash of the ¹²C(α,γ)¹⁶O reaction.

Depending on their mass and metallicity, helium stars can reach different outcomes at the end of their evolution. Possible outcomes have been characterised and studied in detail at solar metallicity (Woosley 2019; Ertl et al. 2020; Chanlaridis et al. 2022). Helium stars that are massive enough to burn helium but with initial masses below ∼1.9 M_⊙ become CO or ONeMg white dwarfs (WDs). A small mass window may exist where helium stars go through more advanced nuclear burning processes that result in the formation of thermonuclear and electron-capture SNe (Antoniadis et al. 2020; Chanlaridis et al. 2022). Helium stars more massive than about ∼2.8 M_⊙ form an iron core that will eventually collapse. Helium stars with final masses higher than ∼35 M_⊙ may form either pulsational pair-instability SNe, leading to core collapse and the formation of a BH, or a pair-instability SN, leading to an explosion that leaves no remnant behind (e.g. Woosley 2017).

With the physical and numerical choices we have taken, mass loss is the process that most drastically affects the evolution of helium stars of different metallicities. The effect on envelope inflation is not captured by our treatment using MESA’s mlt++.

Some models lose enough mass to expose the layers where nitrogen is depleted and carbon and oxygen are enhanced. At this point, they experience an increase in mass loss rate, corresponding to the WC stage. This creates a dichotomy in surface chemical composition of helium star models at core collapse, as well as a distinct distribution of final masses (Yoon 2017; Paper II).

Since helium stars lose mass at a high rate, above a certain initial mass their convective cores tend to decrease in mass. Therefore, no semi-convective layers are formed in regions with chemical composition gradients. This, along with the fact that our models do not include rotational mixing or mixing by tides, implies that stars that could be observable as WN- or WC-type stars, with surface layers rich in both carbon and nitrogen, are not produced in our grids (Langer 1991; Hillier et al. 2021).

The evolution of helium star models of different mass and metallicity is illustrated by the Kippenhahn diagrams of four example evolutionary calculations in Fig. 3. The upper left panel shows our 8 M_⊙ model with a metallicity of 0.01 (corresponding to a star with ZAMS mass of about 28 M_⊙). It has a convective core that initially spans ∼4.8 M_⊙, and which gradually decreases in mass as the total mass of the star decreases due to wind mass loss. This trend continues until the end of core helium burning, when core convection stops. When helium is depleted in the core, helium burning is ignited in a shell above the core. Afterwards, this particular evolutionary sequence goes through carbon burning in a series of convective flames that appear at subsequently higher mass coordinates before igniting neon in the core. The layers enriched by helium-burning ash are never exposed by mass loss, and thus this model has Y_surf = 1 − Z_init and a relatively high surface nitrogen abundance throughout its entire evolution.

Fig. 3. Kippenhahn diagrams following the energy generation and loss rates and the structure of convective regions as a function of time remaining until core collapse, from core helium burning until a few days before core collapse. Represented are the evolutionary calculations with initial helium star masses of 8 (top) and 20 (bottom) M⊙, with metallicities of Z = 0.01 (left) and Z = 0.04 (right). The colour denotes the intensity of the energy generation (red) and loss (blue) rate. The hatched regions denote convective regions. The helium and carbon burning regions are indicated with arrows. The point where the products of helium burning are exposed at the surface, and the model transitions from WN to WC mass loss rates are indicated with horizontal arrows.

In comparison, an evolutionary sequence with the same initial mass and metallicity of 0.04 (upper right panel in Fig. 3) loses mass at a higher rate, and thus its convective helium-burning core also decreases in mass at a faster pace. During core helium burning, the layers depleted of nitrogen and enriched in carbon and oxygen are exposed. At this point the wind transitions from the WN regime to the WC regime, resulting in an increase in mass loss rate, which also accelerates the retreat of the convective core. Subsequent episodes of nuclear burning in the core occur similarly, but the final mass is smaller, and the final surface composition is different, having a lower helium content, and with the ash of helium burning also present.

To illustrate the effect of initial helium mass, we also show an evolutionary sequence with initial mass of 20 M_⊙ and metallicity of 0.01 (corresponding to a star with ZAMS mass of about 52 M_⊙; see lower left panel in Fig. 3). This model experiences the transition between mass loss regimes early in its evolution, similar to the 8 M_⊙ evolutionary sequence at a metallicity of 0.04. However, it has a different core evolution due to its higher mass. Carbon burning in the core occurs radiatively, and as it progresses the burning region moves to higher layers of the star until it settles at the base of a convective zone that encompasses almost the entire carbon-oxygen core. This is common for stars of high mass, and affects the core structure, likely leading to the formation of a BH instead of a SN explosion (Brown et al. 2001; Sukhbold & Woosley 2014).

At a metallicity of 0.04, the evolutionary sequence with an initial mass of 20 M_⊙ (see lower right panel in Fig. 3) has a stronger mass loss and experiences a wind-regime transition at an earlier evolutionary phase than its lower mass counterpart. Due to its reduced mass at the end of core helium burning, it experiences carbon burning in a series of convective flames like the lower mass examples, linked to a higher probability of ending in a successful SN explosion at the end of its evolution (cf. Paper II).

As illustrated by these examples, populations of stripped-envelope stars in different metallicity environments will appear different. Stripped-envelope stars in high metallicity environments will spend a longer fraction of their lives with their carbon- and oxygen-rich layers exposed; they will have lower final masses compared to populations in lower metallicity environments; and the population of SNe they produce will be different, both in terms of total SN rates, of relative rate of different SN sub-types, and of observable SN properties (cf. Paper II).

3.2. Lifetimes and duration of evolutionary stages

Hydrogen-free classical WR stars can be classified as either WN, WC, or WO type. Each of these types corresponds to an evolutionary stage, and they are characterised by different surface abundance patterns and luminosity-to-mass ratios. Their occurrence, properties, and relative rate is known to be different in environments of different metallicities, and has been used as a testbed of massive star evolution (e.g. Massey & Johnson 1998; Massey & Holmes 2002; Meynet & Maeder 2005; Vanbeveren et al. 2007; Georgy et al. 2012; Eldridge et al. 2017; Neugent & Massey 2019; Pauli et al., in prep.). The probability of observing a WR star either as a WN, WC, or WO depends strongly on the duration of each of these stages. In this section we analyse the duration of each stage in our helium star models according to the surface abundances and evolutionary stage, regardless of the surface optical depth. Whether the surface optical depth is large enough to be observed as WR-type stars is discussed in Sect. 4.

The total lifetimes of our helium star models are shown in Fig. 4, and are available online². Similar to main sequence lifetimes, helium burning lifetimes of stripped-envelope stars are a monotonically decreasing function of initial mass. The lifetimes of models that do not expose the products of helium burning in their surfaces are well described by their helium burning nuclear timescale at the moment of stripping. The lifetime between core helium depletion and core collapse is only 2−6% of the total lifetime (see Fig. 5, bottom panel). Models at lower metallicities are shown to have shorter lifetimes because they retain a greater amount of mass at any given evolutionary stage, but the difference is limited to a few percent in this regime. Metallicity has a larger effect on helium star models that spend time with a carbon- and oxygen-rich surface (to the right of the circles in Fig. 4) since they lose mass at a faster rate, and therefore decrease their convective helium-burning core mass early in their evolution. Helium star models that at some point have carbon- and oxygen-rich surfaces at the highest metallicity in our grid have lifetimes that are up to ∼14% shorter than their low metallicity counterparts, particularly in the regime of initial helium star masses of 10−20 M_⊙.

Fig. 4. Total lifetime of helium stars, from the start of helium core burning to core collapse, as a function of initial helium star mass. The different coloured lines represent different metallicities. The downward pointing triangles at the top indicate the lowest initial helium star mass needed for a model to transition into the WC stage before the end of its evolution. The dots on each line indicate the minimum mass at which stripped-envelope stars may be classified as WN stars, obtained via Eq. (7), and using the mass-luminosity relation of Langer (1989). Helium stars below this limit do not have optically thick WR winds, and have overestimated mass loss rates.

In our simplified picture, where stellar models are assumed to be fully stripped of their hydrogen envelope at helium ignition, the duration of the different evolutionary phases depends on their total lifetime and the strength of their winds. The initial surface helium mass fraction of our models is given by Y_surf = 1 − Z_init. The surface abundances of carbon, nitrogen, and oxygen are set by CNO equilibrium, and are consequently initially nitrogen-rich, and carbon- and oxygen-poor. Therefore, all stellar models spend at least a fraction of their lifetime with nitrogen- and helium-rich envelopes (defined as t_He − N), and are thus candidate WN stars. Models that have helium-burning ash in their surface are candidate WC- or WO-type stars. We distinguish them by defining WC star candidates as those that are carbon-rich in the surface, and candidate WO stars as those that have a surface oxygen abundance of at least 0.05 and have finished core helium burning, after which their luminosity-to-mass ratio sharply increases. The duration of the carbon-rich phase (defined as t_He − C) is set by how quickly stars become stripped of their helium-nitrogen envelope. The duration of the oxygen-rich phase after core helium depletion (defined as t_He − O) is set by how quickly stellar models reach the threshold surface oxygen abundance in the lower mass regime, and by the remaining lifetime after core helium depletion in the higher mass regime. The fraction of the total lifetime that our models spend in each of the phases is summarised in Fig. 5.

Fig. 5. Ratios of lifetimes of different evolutionary phases to total lifetime as a function of initial mass, for helium stars of different metallicities. Top panel: ratio of the time spent with a nitrogen-rich envelope to total lifetime. Middle panel: ratio of time spent with a carbon-rich envelope to total lifetime. Bottom panel: ratio of lifetime spent with an envelope with surface oxygen abundance higher than 0.05 after core helium burning to the total lifetime. The dashed lines in the bottom panel indicate the ratios of the lifetime spent with an envelope with surface oxygen abundance higher than 0.05 to the total lifetime, and the dot-dashed lines indicate the ratios of the lifetime after core helium depletion to the total lifetime.

The minimum mass at which models spend a fraction of their total lifetime as both WC- and WO-type star candidates decreases sharply with metallicity. The fraction of the lifetime that stellar models of the same initial mass spend in these more advanced stages is also an increasing function of metallicity. Combined with the increasing total lifetime, the relative number of WC and WO stars is expected to be strongly dependent on metallicity.

3.3. Minimum luminosities of carbon- and oxygen-rich stripped-envelope stars

The minimum luminosity at which stripped-envelope stars are produced depends exclusively on the efficiency of their formation channel, and its dependence on ZAMS mass, binary fraction, and metallicity. Whether such stars are WN or transparent wind helium stars depends on their surface optical depth.

Similarly, the emergence of carbon- and oxygen-rich stripped-envelope stars will depend on the efficiency of mass loss in stripped-envelope stars to lose their nitrogen-rich layer. Therefore, the minimum luminosity of carbon-rich and oxygen-rich stripped-envelope stars will depend on the metallicity of the environment.

Figure 5 shows that stars with an enhanced carbon and oxygen surface abundance will typically have higher initial masses than their nitrogen-rich counterparts, and will therefore be found at preferentially higher luminosities. Carbon- and oxygen-rich stellar models are only produced above a threshold luminosity, $L_{\min ,\mathrm{WC}}^{\mathrm{evo}}$. Stellar models with oxygen surface abundance above 0.05 after core helium burning are produced above a threshold luminosity $L_{\min ,\mathrm{WO}}^{\mathrm{evo}}$. The value of these luminosities is determined from our models, and is indicated in Fig. 6. Both limits decrease with increasing metallicity due to the stronger winds being more efficient at exposing the formerly convective cores of our helium star models.

Fig. 6. Evolution of surface luminosity as a function of age, between the He-ZAMS and core carbon depletion, for helium star models with different initial masses, coloured by surface helium abundance. Each panel represents a set of models with different metallicity, from 0.01 (top left) to 0.04 (bottom). Arrows on the right side of the panel with Z = 0.01 represent the luminosities of observed WC stars in the LMC, inferred from the sample of Bartzakos et al. (2001), assuming a bolometric correction of Mbol = 4.5 (e.g. Moffat 1989; Smith et al. 1994), and have typical uncertainties of 0.3 dex. Arrows on the right side of the panel with Z = 0.02 represent the luminosities of observed WC stars in the Galaxy, taken from Sander et al. (2019). The thin horizontal dashed lines represent the minimum luminosity at which carbon-rich models are formed at each metallicity, $L_{\min ,\mathrm{WC}}^{\mathrm{evo}}$. The solid black lines indicate the moment at which stellar models end core helium burning, and have surface oxygen abundance higher than 0.05. The horizontal thick dotted lines indicate the minimum luminosity at which the optical depth of nitrogen-rich stripped-envelope stars is expected to be higher than 1.45, according to Eq. (7). The horizontal thick dash-dotted lines indicate the minimum luminosity at which the optical depth of carbon-rich stripped-envelope stars is expected to be higher than 1.45, according to Eq. (8).

As shown in Fig. 3, massive models lose mass at rates high enough that their convective helium-burning cores will decrease in size. This leads to a gradual decrease in luminosity during helium burning, as shown in Fig. 6. The drop in luminosity is greater at higher metallicities, and faster for models with carbon-rich surfaces that transition to WC mass loss rates. Models with lower masses have a more steady luminosity during their evolution as their winds are not strong enough to decrease their mass significantly, and they do not transition to WC-type mass loss.

Much like in models of stars with hydrogen rich envelopes, helium star models will experience an increase in luminosity a few thousand years before the end of their evolution, caused by the ignition of helium shell burning when helium is depleted in their core. The mass of the least massive model that has oxygen surface abundance above 0.05 determines the value of $L_{\min ,\mathrm{WO}}^{\mathrm{evo}}$.

For models with luminosities below the minimum luminosity at which helium-burning products are exposed at the surface, the final surface helium abundance corresponds to the initial value. Above this limit helium stars expose their formerly convective layers, and the helium surface abundance decreases steadily, following the smooth composition gradient that is left above the convective helium-burning core. In helium star models with intense enough mass loss, a saturation value of around Y_surf ∼ 0.3 is reached after helium burning.

Figure 7 shows the minimum initial helium star mass above which our models are found to expose the products of helium burning in their surface before core collapse. This minimum mass decreases with increasing metallicity, due to the increase in mass loss rates in the WN stage. Since the mass loss rates also increase with increasing luminosity, which is a proxy of mass (e.g. Langer 1989; Gräfener et al. 2011), helium stars that are massive enough to expose layers of their formerly convective helium-burning core spend a fraction of their lifetime as WC candidates. This fraction increases the more massive they are. Therefore, as illustrated in Fig. 5, for a fixed initial mass, the fractions of stellar lifetime spent as a WC and a WO candidate increase as a function of metallicity. This implies that WC and WO stars ought to be more common in environments of higher metallicity. The value of the transition mass can be fitted by a power law. This fit is given by

$$\begin{array}{c}\hfill M_{\min ,\mathrm{WC}}^{\mathrm{evo}}=8.93{\left(\frac{Z}{0.02}\right)}^{-0.362}{M}_{\odot },\end{array}$$(10)

Fig. 7. Minimum initial helium star mass $M_{\min ,\mathrm{WC}}^{\mathrm{evo}}$ at which a stellar model exposes its carbon- and oxygen-enriched layers (potential WC-type stars) as a function of metallicity (dots). The minimum initial helium star mass can be described as a power law, shown as a solid line (see Eq. (10)).

and is shown in Fig. 7.

4. Wolf-Rayet stars versus transparent wind stripped-envelope stars

In this section we contrast the minimum luminosity of WR stars of different types, calculated employing the criterion based on the optical depth of the stellar wind, described in Sect. 2.2, to the outcome of the models presented in Sect. 3. In Sect. 4.1 we discuss what type of WR star can be expected to be the least luminous as a function of metallicity. In Sect. 4.2 we discuss the minimum luminosities of WN, WC, and WO stars in the context of the results of our models. Finally, in Sect. 4.3 we discuss how our results compare to observed samples of WR stars in the local Universe.

4.1. The least luminous Wolf-Rayet stars as a function of metallicity

The luminosity below which WR stars transition into transparent wind stripped-envelope stars can be estimated as a function of metallicity via the optical depth at their surface, as outlined in Sect. 2.2. Helium- and nitrogen-rich stripped-envelope stars with luminosities below the minimum WN luminosity may evolve with lower mass loss rates than those employed in our models (see Sect. 1). If they are stripped by a binary companion and have no further interactions after being stripped, their final masses are expected to be closer to their mass at the moment of stripping, their luminosities are expected to decrease less during their evolution, and their lifetimes are expected to be shorter.

Stellar models with helium- and nitrogen-rich surfaces in our grids with luminosities higher than $L_{\min ,\mathrm{WN}}^{\mathrm{tau}}$ are therefore representative of WN stars. As shown in Fig. 2, in the metallicity range covered by our models, we find that the minimum luminosity above which stripped-envelope stars with nitrogen-rich helium envelopes are likely classified as WN stars is higher than the minimum luminosity above which stripped-envelope stars with carbon- and oxygen-rich envelopes are likely classified as WC stars (i.e. $L_{\min ,\mathrm{WN}}^{\mathrm{tau}}>L_{\min ,\mathrm{WC}}^{\mathrm{tau}}$). This implies that the lowest luminosity WR that can exist in these environments is of WC type, regardless of whether they are formed by wind mass loss or not.

By extrapolating Eqs. (7) and (8) to higher metallicities, we find that this trend persists since the minimum luminosity of WC stars has a steeper metallicity dependence than the minimum luminosity of WN stars. However, the roles become inverted below metallicities of 0.0014, where according to our best fit models $L_{\min ,\mathrm{WN}}^{\mathrm{tau}}=L_{\min ,\mathrm{WC}}^{\mathrm{tau}}$. We predict that in environments with metallicities below this value, the least luminous WR star that can exist will be of WN type instead.

The considerations made to derive the minimum luminosities of WN and WC stars carry no information about how they are formed, but rather just give the luminosity ranges where we could potentially observe stripped-envelope stars as WR stars. Whether they are formed or not depends first of all on the probability of stripping hydrogen-rich stars, and subsequently on the intensity of their winds.

If the metallicity dependence of WN, WC, and WO winds can be extrapolated to environments of low metallicity, our findings imply that both WN and WC stars become very rare in such environments, as the minimum luminosity increases quite quickly with decreasing metallicity, and they will only originate in very massive stars with high luminosities. Furthermore, the least luminous WR star that can exist in environments with metallicities below about 0.0014 (according to Eqs. (7) and (8)) will be a WN star, as opposed to more metal-rich environments where the least luminous WR stars that can possibly exist are of WC type.

As an example, in a galaxy like I Zwicky 18, with a metallicity of 0.0002 (Szécsi et al. 2015), we predict that the minimum luminosity at which a nitrogen-rich stripped-envelope star has τ ≤ 1.45 is about 1.8 × 10⁶ L_⊙, and for a carbon-rich stripped-envelope star it is about 9.5 × 10⁶ L_⊙, corresponding to minimum masses of about 47 M_⊙ and 147 M_⊙, respectively.

As a final note, we emphasise that the metallicity dependence of the minimum luminosities of WN and WC stars found through the optical depth model presented in Sect. 2.2 have the empirical mass loss rates of WR stars as an input, so we caution that the behaviour of these quantities might change as our understanding of the mass loss rates of WR stars improves.

4.2. Luminosity range of WC- and WO-type stars as a function of metallicity

The least luminous carbon-rich stripped-envelope models produced in our grids of evolutionary calculations are dictated only by the effect of mass loss. As shown in Fig. 8, these models have luminosities well above the luminosity at which carbon-rich stripped-envelope stars would transition into having optically thin winds. The minimum luminosity carbon-rich models in our grids have optical depths τ ≳ 3, and more luminous models have even larger optical depths. Therefore, they are representative of WC stars. Since the mass loss rate of the WC and WO stars we employ has the same functional form, we find that all of our oxygen-rich core helium depleted models are representative of WO stars after core helium depletion as well.

Fig. 8. Analytical and numerical predictions of the minimum luminosity of WC- and WO-type stars as a function of metallicity, compared with observed populations. The solid line in the left panel shows the minimum luminosity at which WC stars are predicted to be observable as a function of metallicity, according to Eq. (8). The black pentagons indicate the minimum luminosity at which our stellar models expose their carbon- and oxygen-enriched layers. The panels to the right show the observed luminosity distribution of WC and WO stars in the LMC (Hainich et al. 2014) and in the Galaxy (Sander et al. 2019). The horizontal lines compare the minimum luminosity carbon- and oxygen-rich stars produced in our models at Z = 0.01 and Z = 0.02 with the observed distributions in the LMC and the Galaxy, respectively.

Furthermore, we find that the luminosity ranges at which WC and WO stars are formed are similar. However, we expect that the latter are more common at low luminosities since the remaining lifetime between helium core depletion and core collapse decreases as helium star mass increases (see Figs. 5 and 6). Models representative of WC stars, on the other hand, spend a greater fraction of their lifetimes as WC stars with increasing mass, so the luminosity distribution of WC stars is expected to result from a competition between the formation probability of a stripped-envelope stars of a given mass (which is partly determined by the IMF, and decreases as mass increases), and the amount of time they spend in the WC phase (which increases with increasing helium star mass).

WC and WO stars with luminosities lower than those found in our models can potentially be formed if there are additional channels that can strip helium stars of their nitrogen-rich helium envelopes, making their surfaces carbon-rich and nitrogen-depleted, such as subsequent episodes of binary mass transfer. Alternatively, they can be a consequence of mixing processes that either make the convective helium-burning cores of stripped-envelope stars larger than what we find in our models, or efficiently exchange material between the envelope and the convective helium-burning core.

The least luminous carbon-rich models produced in our grids reach luminosities below the minimum luminosity of WN at metallicities of about 0.03 and above. This implies that WC stars that are less luminous than the least luminous WN star are formed via wind mass loss only in these environments.

In the lower metallicity regime, since we expect that the minimum luminosity at which carbon- and oxygen-rich stripped-envelope stars can be observed as WC stars will increase, WN stars may be capable of forming carbon-rich and nitrogen-poor helium stars with transparent winds via wind mass loss, but further modelling of stripped-envelope stars in these environments is required to confirm this conjecture.

We note that evolutionary tracks in our grids that expose the products of helium burning during their evolution always have luminosities above the minimum luminosity of WN stars inferred with the wind optical depth method at all metallicities. This means that our numerical modelling employing the Yoon (2017) mass loss rates is correctly applicable to all evolutionary tracks that spend time in the WC regime, as they do not transition into the regime where we expect their mass loss rates to be lower. We note that the minimum luminosity of WN stars is independent of the evolutionary state of our models, whereas the minimum luminosity of WC in our grids of evolutionary calculations is set by the evolutionary history of our models, which decrease in luminosity during their evolution as a consequence of mass loss (see Fig. 6), and does not correspond to the initial luminosity of our models.

4.3. Comparison with observations

Figure 8 shows the luminosity distribution of WC stars in the LMC and in the Galaxy, compared to the predicted minimum luminosities from our optical depth model, and those achieved in our evolutionary calculations. As shown in this Fig. 8 (and in Fig. 6), all of the WC and WO stars observed in the LMC and in the Galaxy are above the minimum luminosity we predict using our optical depth model, and the bulk of their populations are within the luminosity range in which they are formed by our models. This implies that most of them are consistent with having been formed from a WN star that lost its nitrogen-rich helium envelope due to winds, and do not require an additional source of mass loss to be formed. In the LMC, the fact that the only WO star known is close to the lowest luminosity at which such stars are produced in our models is also consistent with the fact that the least luminous WO stars also have the longest lifetimes.

In our Galaxy, however, a tail in the luminosity distribution of WC stars populates the region where we predict WC and WO stars could exist according to our optical depth model, but where no candidates are found in our evolutionary calculations. The tail of the luminosity distribution of WC stars in the Galaxy contains few stars, and the luminosities determined for these stars have significant uncertainties. A possible explanation is that stars that populate this tail may be found in environments within our Galaxy with higher metallicity, but if the tail is significant for the commonly adopted Galactic metallicity value, it may imply that WC- and WN-type winds alone are unable to produce WC stars in the luminosity regime of the least luminous observed WC and WO stars. This might either imply that stars below this limit are formed by a different mechanism, such as case BB mass transfer, that their luminosities are underestimated, or that either WN or WC mass loss rates are higher than what we employed in our models. The least luminous star in the Galactic sample is WR 119, with $logL/L_{\odot }=4.7_{-0.2}^{+0.25}$. Sander et al. (2019) noted that it is the least luminous object in their sample, but the uncertainty on its parallax is large, implying that the real uncertainty on its luminosity could be up to 0.3 dex, still well below $L_{\min ,\mathrm{WC}}^{\mathrm{evo}}$. However, this object is located towards the centre of the Galaxy, where metallicity is potentially higher than in the solar neighbourhood.

A larger sample of observed WR stars, in a larger metallicity range, and a more stringent determination of their distances and parameters will help settle the question of whether or not transparent-wind carbon-rich stars can be formed via wind mass loss (as opposed to case BB mass transfer, see Tauris et al. 2013), and whether or not additional channels for forming WC and WO stars are necessary to populate the gap between the least luminous stars of these classes produced by WN mass loss, and the least luminous carbon-rich stripped-envelope stars.

Recent work has been made to predict the observational characteristics of helium stars with optically thin winds through spectral modelling (Götberg et al. 2018), and this theoretical work has inspired observational campaigns that have produced candidate transparent wind stripped-envelope stars. These pursuits will provide a complementary perspective on the nature of stripped-envelope stars, by giving us a glance at the population of stripped-envelope stars of a lower mass. Information about them will constrain the stripping efficiency of stars as a function of mass, hopefully adding the missing pieces of the puzzle.

As a final note, we highlight that the existence of a minimum luminosity limit for WR stars has been established observationally (Shenar et al. 2020), and therefore the results and discussion we put forward are independent of the model we use to calculate the minimum luminosity of WR stars.

5. Wolf-Rayet populations in the Local Group

Using the method described in Sect. 2.3, we created simplified synthetic populations of stripped-envelope stars. We computed the number of transparent wind stripped-envelope stars with initial masses above 3 M_⊙, WN stars, and WC stars as a function of metallicity, presented in Fig. 9. We normalised these so that the number of WC stars in the population with metallicity of 0.01 matches the 25 known WC stars in the LMC. We then computed the average number ratios of WC to WN stars, ⟨N_WC/N_WN⟩; WO to WN stars, ⟨N_WO/N_WN⟩; and WO to WC stars, ⟨N_WO/N_WC⟩, as a function of metallicity. These ratios are shown as thick blue lines in Fig. 10. In these calculations it is assumed that the star formation rate, binary fraction, and stripping probability are the same at any given mass, and thus their uncertainties cancel out when calculating number ratios of different types of WR stars.

Fig. 9. Numbers of transparent-wind stripped-envelope stars, WN stars, and WC stars obtained from our population models as a function of metallicity. The numbers are normalised so that the population with metallicity of 0.01 has 25 WC stars, matching the population of the LMC (Neugent et al. 2018).

As shown in Fig. 9, the total number of stripped-envelope stars increases slightly as a function of metallicity. This is mainly due to the fact that as metallicity increases, the lifetime of a stripped-envelope star of a given initial mass increases as well (see Fig. 4). The populations are dominated by transparent-wind stripped-envelope stars at every metallicity, followed by WN-type stars, then WC stars, and finally WO-type stars. Qualitatively, this is in agreement with the known populations of WR stars, where WN-type stars represent the majority of WR stars, whereas WC-type stars are progressively more rare, and WO-type stars are very uncommon.

The decrease in the number of transparent-wind stripped-envelope stars, and the corresponding increase in the number of WN-type stars as a function of metallicity, is caused by the decrease in the minimum luminosity of WN stars. As shown in Fig. 9, the number of transparent-wind stripped-envelope stars exceeds the number of WN stars by a factor of 8 in the lowest metallicity regime in our simulation, and by a factor of about 3 at the highest metallicity.

As shown in Figs. 9 and 10, the number of WC-type stars relative to WN-type stars increases as a function of metallicity. This is a consequence of two effects. The first is that, as metallicity increases, models with the same initial helium star mass spend more time as WC stars, as shown in Fig. 5. The second is the decrease in the minimum initial helium star mass above which WC stars are produced in our evolutionary calculations. These two phenomena are both a consequence of the higher mass loss rates during the WN phase with increasing metallicity. However, that this ratio increases as a function of metallicity is not trivial since the minimum mass at which stellar models are classified as WN-type stars also decreases with increasing metallicity, implying that lower mass helium stars, favoured by the IMF in our simple picture, are also produced in greater abundance. This is the reason why this ratio stalls above a metallicity of about 0.03 at a value of about 0.6.

Fig. 10. Average number ratio of WC to WN stars (top panel), WO to WN stars (middle panel), and WO to WC stars (bottom panel) as a function of metallicity, as predicted from our models (solid lines). This quantity strongly depends on the minimum mass for a star to be observable as a WR. The different lines in the top and middle panels correspond to employing different values of the minimum mass of WN stars, found adopting different values of wind optical depth for the transition between WN stars and transparent wind stripped-envelope stars. The solid line corresponds to the fit to observations in Eq. (7), adopting τ = 1.45 (blue). The dashed lines are calculated adopting τ = 1 (purple) and τ = 2 (green). Data shown in triangles correspond to the observed ratios of Neugent & Massey (2019), and are lower limits due to the inclusion of hydrogen-rich WN stars. Data shown in circles for the LMC and SMC correspond to the same ratio, with WN stars with XH > 0.3 removed (according to the analyses of Hainich et al. 2014 for the LMC, and Hainich et al. 2015; Shenar et al. 2016 for the SMC).

Similarly, the increase in the number of WO-type stars relative to both WN-type and WC-type stars follows from the fact that models with higher metallicity have higher mass loss rates, so that the oxygen-rich layers become exposed earlier in their evolution. Another factor that contributes to the increase in these two ratios in a population is that the time between core helium depletion and core collapse increases as a function of metallicity as it is inversely proportional to the mass at the end of helium burning. Finally, this ratio also increases because the minimum mass at which a stellar model reaches the WO stage decreases as a function of metallicity, therefore becoming more favoured by the IMF.

As shown in Fig. 10, we find that the number of WO-type stars relative to both WN- and WC-type stars increases monotonically, without stalling at a fixed value. This occurs because the minimum mass at which stars enter the WO stage decreases faster than the minimum mass at which they become WC stars. Furthermore, as opposed to WC-type stars, the least massive WO stars tend to also be the ones that spend the greatest fraction of their lifetime in this phase. Therefore, on average, a larger proportion of WO stars in a population of WRs is a better indicator of metallicity than the number of WC stars relative to WN stars alone.

The number of WC- and WO-type stars relative to WN-type stars is sensitive to the value of the minimum mass above which stripped-envelope stellar models would be classified as WN-type stars. To show this, we repeated the calculations using a range of minimum WN masses that correspond to different assumptions about the optical depth at which the transition between WN and transparent wind stripped-envelope stars occurs. We adopt τ = 1 and τ = 2, shown respectively with a dot-dashed and a dashed line in Fig. 1. The resulting ratios are shown in the purple and green dashed lines of Fig. 10.

For all three choices of minimum WN mass, the numbers of both WC- and WO-type stars relative to WN-type stars are found to increase as a function of metallicity with roughly the same slope. The increase is less than a factor of 2 in all cases. This behaviour qualitatively reproduces the behaviour of the WC/WN ratio in observed populations (Neugent & Massey 2019). The dependence on the choice of minimum mass is also steep due to the steepness of the IMF. Increasing the value of the minimum mass of WN stars excludes lower mass stars, favoured by the IMF, and more likely to spend time as WN stars than WC stars. Therefore, this choice will increase the number of WC stars at any given metallicity.

On one hand, because of the steepness of the IMF, if the stripping probability is constant, most WN stars formed will be of relatively low mass, and cluster around the minimum WN luminosity. On the other hand, because the mass-luminosity relation of WN stars is also steep, a small variation in the minimum luminosity of WN stars results in a large variation in the number of WN stars. The combination of these two effects, as shown in Fig. 10, has a very strong impact on the morphology of WR populations. We find that after a small variation in the adopted minimum mass of WN stars, this ratio varies by almost a factor of 2 regardless of the metallicity. The variation that results from changing the lowest WN mass is larger than the variation that results from changes in metallicity by a factor of 4, even though more stars evolve to the WC phase at higher metallicity and spend a longer fraction of their remaining lifetime in this phase.

This implies that the minimum luminosity at which stripped-envelope stars are observable as WR stars has a larger impact on WR populations than the metallicity dependence of WR mass loss rates. Since all carbon- and oxygen-rich stars in the metallicity range of our calculations are expected to be observable as WR-type stars, we expect that these ratios are also indicative of the minimum luminosity of WN-type stars. Careful consideration of these findings is paramount to understanding the evolution of WR stars and transparent-wind stripped-envelope stars in the Universe. Previous theoretical studies of WR stellar populations where this quantity has been computed have found that the stripping probability of main sequence stars, either single (e.g. Meynet & Maeder 2005) or binaries (e.g. Eldridge et al. 2017), is key to understanding this ratio, but we argue that including the behaviour of the minimum luminosity of WR stars as a function of metallicity can be just as important.

The picture becomes more complicated in low metallicity environments, since other poorly constrained quantities, such as the minimum luminosity of WC stars, and the efficiency at which they are formed by WR winds become important, and complicate the matter further. However, we hope to see these considerations accounted for in future population synthesis models, to build a more accurate picture of the WR populations in different galaxies, as they have a very large impact in their environments.

The ratio of WO to WC stars does not depend on the wind optical depth model because all WC and WO candidates in our models are also found to have an optical depth that is significantly higher than 1.45. However, it is sensitive to the rest of the uncertainties that our models are subject to, such as the mass loss rates assumed, and to the definition of the boundary between WC and WO stars.

To assess the validity of our results, we compare our population with metallicity of 0.01 to the WR population of the LMC, which is known to be complete (Neugent et al. 2018), and we compare the ratios of WC- to WN-type stars to the those inferred by Neugent & Massey (2019). The observed values of log (O/H) + 12 from Neugent & Massey (2019) were converted to metallicities assuming that the ratio of oxygen mass fraction to total metallicity is constant, and corresponds to the solar value (Asplund et al. 2009), and that helium mass fraction increases linearly with metallicity from its primordial value to the value observed in the Sun (Peimbert et al. 2007). Observed values of N_WC/N_WN are made according to spectral class, which means that hydrogen-rich WN stars are included in the WN sample. Therefore, these measurements represent a lower limit to the results in our calculations. For comparison, the values of WN/WC are shown with hydrogen-rich WR stars removed from the sample for the SMC (Hainich et al. 2014) and the LMC (Hainich et al. 2015; Shenar et al. 2016).

We find that by imposing that the population with metallicity of 0.01 has 25 WC-type stars, matching the number of WC stars in the LMC, the number of WN stars is underestimated. Our population model has 48.4 WN stars, while the LMC has a total of 124. However, of the 97 WN stars analysed by Hainich et al. (2015), 74 have a surface hydrogen abundance of less than 0.3. The discrepancy might be due to the fact that the metallicity in the LMC is lower than 0.01, because our model excludes WN stars that originate from single stars, and because of the simplicity in our assumptions. Using this population model we predict that the number of transparent-wind stripped-envelope stars in the LMC should be around 350.

We find that the number ratios in Fig. 10 agree qualitatively with our findings, but not many environments have observed populations of WR stars in the metallicity regime of our calculations. At lower metallicities, the effect of remaining hydrogen in the envelope after binary interactions will likely become important in determining this quantity, so the comparison cannot be made directly with models similar to those presented here, performed at lower metallicity.

The number ratios involving WO stars are more difficult to assess in these populations due to an intrinsic rarity of WO stars. In the environments where observations have been performed, only very few WO stars are found. In the SMC no WC-type stars are known, and there is only one WO-type star known. Using our models with metallicity of 0.01 as a proxy of the LMC, we predict that there should be between 2 and 3 WO-type stars accompanying the 124 WN-type stars known. This is in line with the two WO-type stars known in this environment.

There are no known WO stars in M 31 (Neugent et al. 2012) and M 33 (Neugent & Massey 2011). M 31 contains about 92 WN stars, which would imply that, according to our models, there should be between 2 and 3 WO-type stars. Similarly, we would expect that the 45 WN stars in the inner part of M 31 are accompanied by about 3 WO-type stars.

In the Milky Way the WR type ratios are difficult to determine since the observed WR population is not complete. Future observations of WR stars in our own Galaxy and others can contribute to measuring these ratios in environments of different metallicities, and this pursuit can potentially be used to put constraints on several aspects of stellar evolution.

6. Discussion

The results presented in the above sections were derived using 1D stellar evolution models of helium stars with metallicities in the range 0.01−0.04. Taking this approach, we were able to study the effects of WR mass loss on stripped-envelope stars, without the additional complications and uncertainties of the previous evolutionary stages. This range includes the metallicity of the solar neighbourhood, but extends to environments of high metallicity, which have often been neglected in the literature. The metallicities included in this study are meant to explore the range where the amount of hydrogen remaining in the envelope after envelope stripping is more likely to be negligible. Our aim was to characterise the effect of metallicity-dependent mass loss rates, as well as the metallicity dependence of the luminosity at which stripped-envelope stars are separated between WR stars and transparent-wind helium stars in the populations of WR stars, particularly through the ratio of WC to WN stars.

In the following subsections we address the limitations and uncertainties of the approach we took to obtain them, and compare them to previous models of stripped-envelope stars.

6.1. Uncertainties in helium star evolution

6.1.1. Uncertainties in Wolf-Rayet mass loss rates

The main effect of metallicity on the evolution of helium stars of equal initial mass, given the physical treatment that we employ, is due to the metallicity-dependent winds. The wind mass loss rates used in this work were adapted from Yoon (2017). These were obtained from separate empirical studies of WN (Hainich et al. 2014) and WC stars (Tramper et al. 2016), and calibrated to have matching clumping factors (namely D = 4), which Yoon (2017) found would help reconcile the lowest luminosity WC stars observed to our understanding of WR mass loss.

The empirically obtained mass loss rates of WR stars used in our simulations were extrapolated beyond the luminosity and mass regimes at which the WR mass loss rate observations were carried out. Therefore, we caution that the results of our simulations, particularly those with masses near and below the obtained values of the minimum mass of WN stars at a given metallicity have mass loss rates that may be overestimated.

Currently there is no consensus in the community as to what are the true mass loss rates of helium stars with different luminosities and compositions. Recently, 1D hydrodynamically consistent simulations of wind acceleration in WR stars have been carried out (Sander et al. 2020; Sander & Vink 2020). These studies suggest, contrary to what has been found empirically, that the mass loss rates of WR stars of similar parameters have a shallower dependence on metallicity than that determined from the observed sample of WN and WC stars. Furthermore, these studies suggest that the mass loss rates of WC-type stars are in fact lower than those of WN-type stars of similar luminosity and mass. If this is indeed the case, then the final masses, ejecta masses, and BH masses predicted for the stars in our sample that spend a fraction of their lifetimes as WC stars are also underestimated. Self-consistent theoretical calculations of WR mass loss rates are uncertain inasmuch as the temperatures (and therefore the properties of the wind-accelerating regions) are poorly constrained.

6.1.2. Structure of massive stars after stripping

The initial condition of our simulations is an approximate representation of the state of stars that have lost their hydrogen-rich envelopes through winds, binary stable mass transfer, or CE ejection. This approach has been employed in many studies of the late evolution of stripped-envelope stars (Arnett 1974; Nomoto & Hashimoto 1988; Woosley et al. 1995, 2020; Tauris et al. 2015; McClelland & Eldridge 2016; Yoon 2017; Kruckow et al. 2018; Woosley 2019; Ertl et al. 2020; Higgins et al. 2021), and has been found to be a powerful tool to study a complex stage in massive star evolution, which is very relevant in the study of Type I SNe, WR stars, and gravitational wave sources.

Detailed simulations of binary systems that become stripped through stable mass transfer show that this interaction does not remove the entire envelope of the donor star in the binary system, and that a remaining hydrogen layer, even of a very low mass, can have significant consequences in the following evolution of these stars (Gilkis et al. 2019; Laplace et al. 2020; Pauli et al., in prep.). Furthermore, stars that lose their hydrogen envelopes before the beginning of core helium burning cannot reproduce phenomena such as WN/WC-type stars, which are observed to be rich in both nitrogen and carbon, and require a mixing mechanism that likely arises from helium gradients that only occurs in convective helium-burning cores that grow in mass (Langer 1991).

Recent simulations of binary evolution at LMC metallicity (Pauli et al., in prep.) show that interacting binary systems will retain a hydrogen-rich layer after stripping. At this metallicity, and using the mass loss rates from Yoon (2017), these layers are only expected to be present for ∼10−20% of the helium burning lifetime of the donor stars, depending on their mass, in short period systems. For the systems with the longest periods and for stars evolved in isolation, their models enter the hydrogen-free WN stage after 50% of the helium core has burned. At higher metallicities the winds remove the hydrogen-rich layers more efficiently, therefore reaching the hydrogen-free WN stage earlier in their evolution.

Binary evolutionary calculations with different initial metallicities (Klencki et al. 2022) suggest that this layer is efficiently removed in systems with solar metallicity and above, suggesting that binary stripping is more efficient as metallicity increases. Therefore, we are confident that our models reproduce the properties of stars that have been stripped by Case A or Case B mass transfer, as well as successful CE ejections that occur soon after the ignition of hydrogen shell burning at the metallicities employed in our grids. Detailed binary stellar models are required to study stripped-envelope stars at metallicities lower than those in our grids to properly account for the effect of the remaining hydrogen layer in the internal structure of stripped-envelope stars, and on their orbital evolution. Furthermore, as noted by Schneider et al. (2021), stars that interact earlier or later have different structures at the beginning of helium burning, and conclusions about such systems can therefore not be derived from our models.

However, the simulations of Gilkis et al. (2019) show that if stars that retain a small amount of hydrogen in their envelopes do not lose winds as WR stars (but instead follow the Vink 2017 prescription made for stripped-envelope stars below the WR luminosity limit), then the remaining hydrogen is less likely to be removed, leading to the preferential formation of Type IIb SN progenitors instead of Type Ib progenitors. If mass loss rates of stripped-envelope stars are indeed not high enough to efficiently remove the remainder of the hydrogen envelope after a binary interaction, then approximating them as pure helium stars from the beginning of their evolution may lead to overestimating their final mass, and to overestimating the number of WC stars in a given population.

The structure of stars that have lost their hydrogen envelope through CE interactions is even more uncertain. CE evolution is typically modelled using the so-called energy formalism (Webbink 1984), which provides a method to calculate whether a binary system will merge or will survive a CE interaction and successfully eject the CE. In the latter case, it allows for the calculation of the final separation of the system. However, this formalism depends on an arbitrary choice of boundary between the stellar core, which will be the surviving part of the star after the interaction, and the envelope that will be ejected. Therefore, the surface composition of stars after the interaction is unknown. The range of possible choices and the consequences of variations in this choice have been discussed by several authors (e.g. Tauris & Dewi 2001; Ivanova 2011; Kruckow et al. 2016; Vigna-Gómez et al. 2022), but no consensus has been found on the outcome of CE events. Recent studies of CE interactions (e.g. Fragos et al. 2019; Law-Smith et al. 2020) find that the final structure and separation of the binary components after the interaction is uncertain even in detailed simulations, and that the binary components might be subject to subsequent interactions during the stripped star’s core helium burning phase. However, it is likely that most or all of the hydrogen envelope will be ejected if the event does not lead to a merger.

Single stars can also become stripped of their hydrogen envelope via their stellar winds. Helium stars below a certain mass are probably only formed due to binary interactions as their winds during the main sequence are unlikely to be strong enough to strip them of their hydrogen envelope. However, this depends heavily on uncertain parameters such as main sequence and post-main sequence mass loss, rotation, and potentially unresolved physical phenomena (e.g. Leung et al. 2021); different sets of stellar evolution models predict different values for the threshold mass above which stripped-envelope stars can be formed from single-star evolution, particularly as a function of metallicity (for a detailed discussion, see Shenar et al. 2020).

6.1.3. Physical and numerical uncertainties

Uncertainties in the evolution of massive stars are known to appear in the late stages of evolution due to the numerical resolution and the choice of nuclear network. We model energy generation rates and chemical abundances with MESA’s approx21 nuclear network, which is the smallest nuclear network that can be used to model massive stars to core collapse. This choice may affect the very late phases, but the lifetimes, compositions, and other global evolutionary properties discussed in this paper are not be affected (Farmer et al. 2016).

The efficiency of convection and the presence of overshooting in helium stars are also uncertain. Varying these parameters may lead to more massive convective helium-burning cores. This results in different lifetimes, and leads to significant differences in evolution, particularly in the lower mass regime. Overshooting in low mass helium stars most significantly affects the locations of the boundaries in initial mass that lead to different final outcomes (Chanlaridis et al. 2022), but since low mass helium stars (with luminosities below the WN limit) likely have very low mass loss rates regardless of their metallicity, we argue that the minimum mass at which stripped-envelope stars experience core collapse instead of becoming thermonuclear or electron-capture SNe is not strongly dependent on metallicity, and make a conservative choice for this limit, setting it at 3 M_⊙.

Another uncertainty is introduced in our results from the implementation of MESA’s mlt++ scheme (Paxton et al. 2013). It artificially increases the efficiency of convective energy transport in the envelopes of stars, and is commonly implemented in calculations of massive star evolution to increase numerical stability. However, this tool also suppresses the effect of envelope inflation, which can occur in WR stars, as they are close to the Eddington limit (e.g. Sanyal et al. 2015; Petrovic et al. 2006). Since the mass-luminosity relation for WR stars is well defined (Langer 1989), and the empirically obtained mass loss rates do not depend on any other parameters defined at the surface of these stars, we argue that using mlt++ does not have an impact on the core evolution of helium stars, or on the amount of mass that they lose due to winds. However, this results in computed radii in our models that are likely underestimated. Furthermore, the proximity to the Eddington limit has been associated with violent mass eruptions, as in the case of luminous blue variable stars (e.g. Gräfener et al. 2012). Therefore we argue that the final masses of our models, particularly those above the minimum WN luminosity, are an upper limit, and that the real distribution of SN ejecta masses is likely to be shifted to lower masses. Additionally, the ejecta masses of stripped-envelope SNe can be further displaced to lower values if helium stars expand significantly due to inflation, and experience a second episode of interaction with their companion (Tauris et al. 2013), or if they experience mass eruptions due to their high luminosities.

6.1.4. Rotation in stripped-envelope stars

The models presented in this paper are computed without accounting for the effects of rotation. Rotation can induce mixing of chemical elements in the stellar interior, thereby modifying its structure, evolution, surface composition, and final outcome (see e.g. Heger et al. 2000). However, Yoon et al. (2010) found that angular momentum loss during mass transfer is significant, and most helium stars produced in binary systems by mass transfer are slow rotators. Helium stars produced in single stars, stripped by winds alone, lose most of their angular momentum in the winds, and therefore are also expected to be slow rotators (e.g. Meynet & Maeder 2003). Therefore, we argue that the evolution of stripped-envelope stars is nicely approximated by non-rotating models of helium stars.

Some helium stars are thought to end their lives with significant rotation, particularly when stripping leads to short-period (< 2 days) binaries where tidal spin-up is particularly efficient (Qin et al. 2018). However, this mainly affects low metallicity systems, and the binary configuration that leads to their formation is rare.

Very fast rotation during early hydrogen burning can lead to chemically homogeneously evolving (CHE) stars, an evolutionary channel that diverges from canonical massive star evolution (Maeder 1987). CHE stars transport most of the hydrogen in their surfaces into their cores, where it is burnt into helium, and form critically rotating helium stars. CHE stars might be formed in mass accretors in binary systems at very low metallicity, but it is not expected to be common in the metallicity range considered in our models (Cantiello et al. 2007). They may form in very close binaries (de Mink et al. 2009), but are also expected to be rare and more prominent at metallicities lower than those of our models.

Woosley (2019) argues that helium stars approximately evolve like CHE stars. However, Aguilera-Dena et al. (2018, 2020) find that rotational mixing is still efficient after core hydrogen depletion. This allows the convective helium-burning core to grow, instead of decreasing in mass as in the case of helium stars (see Fig. 3), and prevents the formation of strong chemical composition gradients in the interior of CHE stars during and after core helium burning. This results in different structures and surface chemical compositions between helium stars and CHE stars throughout their evolution. Since they are more prevalent in low metallicity environments, CHE stars lose very little mass due to winds. However, they lose large fractions of mass (both steadily and in bursts) due the combined effect of contraction and rotation. This mass likely remains near the star as it is lost via centrifugal acceleration, and likely contributes to the photometric and spectral properties. Therefore, we argue that their observable properties will also be significantly different. Furthermore, their core structures differ from those of non-rotating helium stars, and they lead to a different distribution of parameters such as core compactness and carbon abundance as a function of initial mass, that have an effect on their final outcome.

6.1.5. Minimum luminosity of Wolf-Rayet stars

The formalism presented in this work to calculate the optical depth at the surface of WR stars, adapted from Langer (1989), makes several simplifying assumptions on the structure and opacity of WR winds. The choices of the values of parameters in the model, such as setting κ to the electron scattering opacity, the density distribution of wind following a β law with β = 1, and the surface and terminal velocities are a simplification, since these may not represent WR stars accurately. For instance, the opacity at the base of the wind includes lines and may be a factor of few greater than the electron scattering opacity. Furthermore, its value is likely metallicity dependent. This contribution would also increase at higher metallicity. However, we justify the use of these choices because of the agreement they produce with the observed values of the minimum luminosities of WN-type stars in the Galactic, LMC, and SMC populations.

The mass-radius relation employed in this calculation (Langer 1989) is calibrated for models of WR stars that do not experience inflation. Inflation can occur in luminous WR stars, or in WR stars with high envelope opacities (e.g. Petrovic et al. 2006; Gräfener et al. 2012). However, it becomes important only for evolved WR stars, and WR stars with luminosities that are considerably higher than the minimum luminosity WR stars predicted by our optical depth model, even at the highest metallicities in our grid, at least during the helium burning lifetime. The mass luminosity relation employed (Langer 1989) was obtained through fits of hydrostatic, helium burning models. It was found to depend on composition, such that pure helium stars of a given luminosity were predicted to be less massive if they were enriched in carbon and oxygen. In our calculations, we employed only the mass-luminosity relation obtained for pure helium stars, but we found that employing a mass-luminosity relation for carbon and oxygen rich stars (Eq. (19) in Langer 1989) only changes the fit parameters of Eq. (8) by about 2%.

The minimum luminosities of WN and WC stars are not trivial to determine from the observed samples of WR stars. WR stars are rare objects since massive stars only spend a short fraction of their total lifetime in this phase. The populations in the Magellanic Clouds are thought to be complete; there are only 12 WR stars in the SMC (Massey et al. 2003) and 154 in the LMC (Neugent et al. 2018), compared with the hundreds of massive stars in the SMC (e.g. Castro et al. 2018) and thousands in the LMC (e.g. Doran et al. 2013). The value of the minimum luminosity of WC stars is more complicated to determine because WC stars are much more scarce than WN stars, particularly at low metallicity. Of the 12 WR stars in the SMC, 11 are WN type, and only 7 are hydrogen-poor (with X_H < 0.3, but not completely hydrogen free). The value of the minimum luminosity of WC stars in the SMC is very uncertain since there are no WC stars in the SMC, and only one WO star. However, according to our estimate, this star appears to be close to the minimum luminosity limit. In the LMC 28 of the 154 WR type stars are of types WC and WO, making this limit easier to determine.

More WR stars are found in the Galaxy (e.g. Hamann et al. 2019; Sander et al. 2019), but the population is still very incomplete, and parameters of many WR stars, such as mass loss rate, surface composition, temperature, and bolometric luminosity, have not been determined due to the lack of high resolution spectra. Furthermore, there is a significant metallicity gradient in the Galaxy (Anders et al. 2017), which increases the uncertainty in determining the values of the minimum luminosities of WN and WC stars within the Galaxy.

Despite the uncertainties that go into the determination of the minimum luminosities of WN and WC stars, both theoretically and observationally, we argue that our results highlight the importance of taking the luminosity distribution of WR stars in different metallicities into account, as well as the mass loss rates that these stars experience during the WN and WC evolutionary stages. The fits provided in Eqs. (7) and (8) can be employed in stellar evolution models and population synthesis studies that study WR stars to properly account for WR mass loss rates, and to distinguish between helium stars with optically thin winds and WR stars of different types.

6.1.6. Uncertainties in our stellar population models

The calculations obtained from the population model described in Sect. 2.3 are a crude approximation, since not every star that is formed becomes stripped, and this probability likely depends on ZAMS mass. Furthermore, many WRs will likely have evolved from single stars that perhaps retain a hydrogen envelope during a substantial fraction of their helium burning lifetime, and WRs that evolve from binaries are not always stripped exactly at the beginning of helium burning, nor necessarily completely stripped. According to Pauli et al. (in prep.), WR stars in the LMC are most commonly formed by binary stripping. Therefore, our models are a good proxy for most of the formed WR stars, at least at this metallicity. Additional stripped-envelope stars are likely formed by reverse stripping, when the least massive component of a binary system is stripped by its more massive companion, which is not accounted for by Pauli et al. (in prep.). The number of systems that undergo this process, however, is uncertain as it depends on natal kicks (Vigna-Gómez et al. 2020). The approximations made in our population model can help us guide our understanding of how varying strengths of WR winds can influence the components of a stellar population.

Since we assume that stars are fully stripped by the beginning of core helium burning, the fractional lifetime that stars spend as a WC star represents an upper limit (unless other effects such as a second mass transfer occur, as suggested by Laplace et al. 2020 for sub-solar metallicity). However, this general trend occurs regardless of the extent of remaining hydrogen shells because even if the WR lifetime is only a fraction of the helium burning lifetime, the effect of metallicity will be the same (Meynet & Maeder 2005).

As discussed in Sect. 5, in the metallicity range that our simulations cover, every time a star has its carbon and oxygen rich layers exposed it will also be luminous enough to appear as a WC star. Therefore, the ratio of WC to WN stars is set by the metallicity dependence of the minimum luminosity of WN stars, the efficiency of WN winds in forming WC stars, the lifetimes of WR stars in each phase, and the stripping probability of stars as a function of mass. In our population model we assume that this efficiency is constant as a function of mass, effectively isolating the influence of winds and the minimum WN luminosity in shaping a population of WR stars.

6.2. Comparison to previous stellar evolution models

The study of the late phases of evolution of stripped-envelope stars is becoming a very active field of stellar astrophysics. Stripped-envelope stars have been modelled as single helium stars (Woosley et al. 1995; Yoon 2017; Kruckow et al. 2018; Woosley 2019; Higgins et al. 2021), and as the late phases of evolution of massive single stars (Massey & Holmes 2002; Meynet & Maeder 2005), although models of full binary evolution are quickly becoming available as well (Eldridge et al. 2017; Schneider et al. 2021; Laplace et al. 2021). The results of these simulations have evolved significantly in recent years, partly because our understanding of WR mass loss has increased considerably in recent years.

We include the WR mass loss rates of Yoon (2017), which differ from what most works on stripped-envelope stars employ, particularly in that we consider that winds of WN and WC stars are different. However, in agreement with Woosley (2019), Schneider et al. (2021) and Laplace et al. (2021), we find that the retreat in convective helium core mass has a large effect on the evolution of stripped-envelope stars.

The final masses of our models with Z = 0.02 are in agreement with Yoon (2017) and Dessart et al. (2020). The final masses and lifetimes of our Z = 0.015 models are also consistent with the Z = 0.0145 models from Woosley (2019) at masses where the WC wind is unimportant. However, we observe systematically lower final masses for helium stars that evolve as WC than what is found by Woosley (2019), at similar metallicity. The core evolution, carbon-oxygen core masses, and luminosities of our models agree well with the results of Woosley (2019) in every mass regime, even though he employed a much more complete nuclear network in his calculations.

The theoretically computed mass loss rates of stripped-envelope stars of Sander et al. (2020) naturally account for the transition between WR stars with strong mass loss rates and transparent wind helium stars, which they predicted to have mass loss rates that are several orders of magnitude lower than those of WR stars. They find that the transition between them is not abrupt, but the mass loss rates of stripped-envelope stars are predicted to quickly decrease below a certain luminosity. They also find that the metallicity dependence of mass loss rates is less steep than that used in this work. Therefore, the helium star models of Higgins et al. (2021), calculated using the mass loss rates of Sander et al. (2020), yield higher final masses than those of Yoon (2017) and Woosley (2019), and those of the models in this work. However, this becomes significant only in the range of masses of the most massive BH progenitors, and the discrepancies in the mass range where the grids overlap are small.

Although the mass loss rates of WR stars employed by Higgins et al. (2021) are different from those used in this work, the consensus is that erroneous results can be obtained in the analysis of stripped-envelope stars due to potential caveats, which result from extrapolating the empirically derived mass loss rates of WR stars beyond their applicability range, as is the case for progenitors of gravitational wave sources.

A detailed analysis of the core properties of our models, and a comparison with other works that have addressed this in the context of progenitors of Type I SNe is deferred to Paper II. Some of the ramifications for compact-object mergers are also discussed in Antoniadis et al. (2022).

7. Conclusions

We have modelled grids of evolutionary sequences of helium stars to mimic the late evolution of massive stars that have been stripped of their hydrogen envelopes. Envelope stripping in massive stars can occur if they are in binary or multiple systems, where a companion is close enough to have an interaction (including stable mass transfer and CE ejections). Alternatively, stars can lose their envelopes via stellar winds or massive ejections during and after the main sequence, but this is only effective for stars with the highest initial masses, and those born in high metallicity environments. We have covered the range of initial helium masses between 1.5 and 70 M_⊙, and metallicities between 0.01, similar to that observed in the LMC, and 0.04, similar to the environments of highest metallicity known, such as the region near the centre of our Galaxy or in nearby massive galaxies.

Evolved helium stars contain a strong chemical composition gradient. Layers initially outside the convective helium-burning core contain abundant nitrogen, enhanced during hydrogen burning due to CNO processing. The convective helium-burning core, on the other hand, is nitrogen-free and has a progressively decreasing amount of helium, and will instead be enriched in carbon and oxygen. This gives rise to the observational dichotomy between WN and WC stars. We follow Tramper et al. (2015) in assuming that WO stars are formed from stripped-envelope stars that have finished core helium burning, and expose their oxygen-rich layers.

Through our models, we find the evolution of massive stars is deeply affected by stripping, and by the amount of mass loss they experience after stripping, which is greater at high metallicity. We characterised the lifetime of helium stars, finding that it is weakly dependent on metallicity. However, the duration of the phases where our models spend time as WN-, WC-, and WO-type candidates is strongly dependent on metallicity, and will determine the composition of WR populations across the Universe. We found that the minimum initial helium star mass at which stars develop a WC phase and a WO phase varies strongly with metallicity.

We find that the duration of the WO phase is mediated by the remaining lifetime after core helium depletion, and the time at which oxygen-rich layers become exposed, both of which increase with metallicity. Therefore, the minimum initial helium star mass at which models are found to spend some time in the WO phase also increases with increasing metallicity.

We have proposed a method to find the minimum luminosity of WR stars, which is different for WN- and WC-type stars, as a function of metallicity. This is accomplished by combining the empirical mass loss rates of Yoon (2017) and the wind optical depth, as derived by Langer (1989). We have provided a quantitative way to determine this limit in different environments.

Using this method, we find that many stars in our grid do not experience a WN-type phase, but rather evolve as transparent wind helium stars. We also find that helium stars that are enriched with carbon and oxygen in their surface are always observable as WR stars (i.e. have optically thick winds) at the metallicities that we cover in our models. However, extrapolating to low metallicities, we find that stars with transparent winds and with surface layers composed of helium-burning products, should exist (like the models of Szécsi et al. 2015).

Our models reproduce the WC star luminosity distributions in the LMC and our Galaxy well, except for very few low luminosity outliers. Even these have luminosities considerably above the predicted transition between optically thin and optically thick winds. It is possible that these stars are a product of isolated helium star evolution, but could also be formed through rare binary evolution channels, for example white dwarf mergers (Gvaramadze et al. 2019) or instabilities in the envelopes of some helium stars.

Through simplified population synthesis, we find that the number of stripped-envelope stars increases slightly with metallicity. We find that transparent-wind stripped-envelope stars outnumber WR stars by a considerable fraction. Using the metallicity in our grids that most closely resembles the metallicity of the LMC, we predict that it will contain about 350 transparent-wind stripped-envelope stars, which may eventually produce Type Ib SNe. Furthermore, the expected number ratio of WC- to WN-type stars is a function of metallicity, and we find that it increases as a function of metallicity, since stars at higher metallicities spend a longer fraction of their life as WC stars, even after taking into account that the number of WN stars might increase due to the decrease in the WN star luminosity limit. However, this increase is not dramatic and stalls at Z ∼ 0.03 due to the shorter lifetime of WC stars compared to WN stars.

We also compute the expected number ratio of WO to WN stars, and of WO to WC stars. Both quantities increase with increasing metallicity because the lifetimes of WO stars increase as a function of metallicity, while the minimum mass at which they are formed decreases. Our estimates are in agreement with the observed number of WO stars in the LMC. For M 31 and M 33, our models predict two to three WO-type stars, while none have been found so far. We note that these number ratios are sensitive to the minimum WN star luminosity, which has not been considered in previous population synthesis works.

With this study we find that populations of stripped-envelope stars are not only determined by the efficiency of their formation channels (i.e. main sequence star winds and binary interactions), they are also influenced by the metallicity-dependent minimum luminosity at which stripped-envelope stars develop optically thick winds and WR characteristics. Furthermore, the metallicity-dependent wind mass loss of stripped-envelope stars is key in determining the WR sub-type distribution. In combination with our findings in Paper II and those in Antoniadis et al. (2022), we emphasise the need for a holistic approach to the study of the evolution of stripped-envelope stars, both through observations of individual systems and stellar populations, and through the study of stripped-envelope SNe across cosmic time.

Acknowledgments

D.R.A.-D. and J.A. are 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). D.P. acknowledges financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50 OR 2005. S.C.Y. has been supported by the National Research Foundation of Korea (NRF) grant (NRF-2019R1A2C2010885). D.R.A.-D. would like to acknowledge valuable discussions with Eva Laplace and Jakub Klencki. This research made extensive use of NASA’s ADS and MESA (http://mesastar.org) (Paxton et al. 2011, 2013, 2015, 2018).

References

All Tables

All Figures

Fig. 1. Optical depth at the base of the wind of a stripped-envelope star as a function of its metallicity and luminosity, obtained from Eq. (6) (based on the optical depth derivation, mass-luminosity, and radius-luminosity relations of Langer 1989). The mass loss rate of WN stars is used, which depends on metallicity and luminosity, as in Eq. (1). The star symbols represent the lowest luminosity observed for WN stars in the SMC, LMC, and Milky Way, according to Shenar et al. (2020). Error bars with a width of 0.05 dex are included as these values are rounded. The lines represent the solutions of Eq. (6) for constant optical depth, including those that best fit the observed values of the observed minimum luminosities of WN stars.

In the text

	Fig. 2. Minimum luminosity at which stripped-envelope stars can be observed as WN-type (green) and WC type WC-type (red), as a function of metallicity, according to Eqs. (7) and (8), respectively.
In the text

Fig. 3. Kippenhahn diagrams following the energy generation and loss rates and the structure of convective regions as a function of time remaining until core collapse, from core helium burning until a few days before core collapse. Represented are the evolutionary calculations with initial helium star masses of 8 (top) and 20 (bottom) M⊙, with metallicities of Z = 0.01 (left) and Z = 0.04 (right). The colour denotes the intensity of the energy generation (red) and loss (blue) rate. The hatched regions denote convective regions. The helium and carbon burning regions are indicated with arrows. The point where the products of helium burning are exposed at the surface, and the model transitions from WN to WC mass loss rates are indicated with horizontal arrows.

In the text

Fig. 4. Total lifetime of helium stars, from the start of helium core burning to core collapse, as a function of initial helium star mass. The different coloured lines represent different metallicities. The downward pointing triangles at the top indicate the lowest initial helium star mass needed for a model to transition into the WC stage before the end of its evolution. The dots on each line indicate the minimum mass at which stripped-envelope stars may be classified as WN stars, obtained via Eq. (7), and using the mass-luminosity relation of Langer (1989). Helium stars below this limit do not have optically thick WR winds, and have overestimated mass loss rates.

In the text

Fig. 5. Ratios of lifetimes of different evolutionary phases to total lifetime as a function of initial mass, for helium stars of different metallicities. Top panel: ratio of the time spent with a nitrogen-rich envelope to total lifetime. Middle panel: ratio of time spent with a carbon-rich envelope to total lifetime. Bottom panel: ratio of lifetime spent with an envelope with surface oxygen abundance higher than 0.05 after core helium burning to the total lifetime. The dashed lines in the bottom panel indicate the ratios of the lifetime spent with an envelope with surface oxygen abundance higher than 0.05 to the total lifetime, and the dot-dashed lines indicate the ratios of the lifetime after core helium depletion to the total lifetime.

In the text

Fig. 6. Evolution of surface luminosity as a function of age, between the He-ZAMS and core carbon depletion, for helium star models with different initial masses, coloured by surface helium abundance. Each panel represents a set of models with different metallicity, from 0.01 (top left) to 0.04 (bottom). Arrows on the right side of the panel with Z = 0.01 represent the luminosities of observed WC stars in the LMC, inferred from the sample of Bartzakos et al. (2001), assuming a bolometric correction of Mbol = 4.5 (e.g. Moffat 1989; Smith et al. 1994), and have typical uncertainties of 0.3 dex. Arrows on the right side of the panel with Z = 0.02 represent the luminosities of observed WC stars in the Galaxy, taken from Sander et al. (2019). The thin horizontal dashed lines represent the minimum luminosity at which carbon-rich models are formed at each metallicity, $L_{\min ,\mathrm{WC}}^{\mathrm{evo}}$. The solid black lines indicate the moment at which stellar models end core helium burning, and have surface oxygen abundance higher than 0.05. The horizontal thick dotted lines indicate the minimum luminosity at which the optical depth of nitrogen-rich stripped-envelope stars is expected to be higher than 1.45, according to Eq. (7). The horizontal thick dash-dotted lines indicate the minimum luminosity at which the optical depth of carbon-rich stripped-envelope stars is expected to be higher than 1.45, according to Eq. (8).

In the text

	Fig. 7. Minimum initial helium star mass $M_{\min ,\mathrm{WC}}^{\mathrm{evo}}$ at which a stellar model exposes its carbon- and oxygen-enriched layers (potential WC-type stars) as a function of metallicity (dots). The minimum initial helium star mass can be described as a power law, shown as a solid line (see Eq. (10)).
In the text

Fig. 8. Analytical and numerical predictions of the minimum luminosity of WC- and WO-type stars as a function of metallicity, compared with observed populations. The solid line in the left panel shows the minimum luminosity at which WC stars are predicted to be observable as a function of metallicity, according to Eq. (8). The black pentagons indicate the minimum luminosity at which our stellar models expose their carbon- and oxygen-enriched layers. The panels to the right show the observed luminosity distribution of WC and WO stars in the LMC (Hainich et al. 2014) and in the Galaxy (Sander et al. 2019). The horizontal lines compare the minimum luminosity carbon- and oxygen-rich stars produced in our models at Z = 0.01 and Z = 0.02 with the observed distributions in the LMC and the Galaxy, respectively.

In the text

	Fig. 9. Numbers of transparent-wind stripped-envelope stars, WN stars, and WC stars obtained from our population models as a function of metallicity. The numbers are normalised so that the population with metallicity of 0.01 has 25 WC stars, matching the population of the LMC (Neugent et al. 2018).
In the text

Fig. 10. Average number ratio of WC to WN stars (top panel), WO to WN stars (middle panel), and WO to WC stars (bottom panel) as a function of metallicity, as predicted from our models (solid lines). This quantity strongly depends on the minimum mass for a star to be observable as a WR. The different lines in the top and middle panels correspond to employing different values of the minimum mass of WN stars, found adopting different values of wind optical depth for the transition between WN stars and transparent wind stripped-envelope stars. The solid line corresponds to the fit to observations in Eq. (7), adopting τ = 1.45 (blue). The dashed lines are calculated adopting τ = 1 (purple) and τ = 2 (green). Data shown in triangles correspond to the observed ratios of Neugent & Massey (2019), and are lower limits due to the inclusion of hydrogen-rich WN stars. Data shown in circles for the LMC and SMC correspond to the same ratio, with WN stars with XH > 0.3 removed (according to the analyses of Hainich et al. 2014 for the LMC, and Hainich et al. 2015; Shenar et al. 2016 for the SMC).

In the text