© The Authors 2025

1. Introduction

With a metallicity of Z ∼ 0.0025, the Small Magellanic Cloud (SMC) is a perfect laboratory to study the evolution of young, massive, and relatively metal-poor stars in the backyard of the Milky Way (Vink et al. 2023). The SMC hosts 12 known Wolf-Rayet (WR) stars, i.e., hot (T_eff > 40 kK) massive stars whose spectrum is dominated by broad emission lines. At low metallicity, the formation of WRs is usually attributed to envelope stripping by a companion star (Vanbeveren et al. 1997) because stellar winds are deemed to be too inefficient to (almost) completely remove the H-rich envelope (Shenar et al. 2020). However, observational searches have failed to identify a companion for seven out of the 12 known WRs in the SMC (Westerlund & Smith 1964; Smith 1968; Sanduleak 1968, 1969; Breysacher & Westerlund 1978; Azzopardi & Breysacher 1979; Vanbeveren & Conti 1980; Moffat 1982; Moffat et al. 1985; Morgan et al. 1991; Bartzakos et al. 2001; Massey & Duffy 2001; Massey et al. 2003; Foellmi et al. 2003; Foellmi 2004; Hainich et al. 2015; Shenar et al. 2016, 2018; Neugent et al. 2018; Schootemeijer et al. 2024). In particular, with a 95% confidence level, Schootemeijer et al. (2024) rule out the presence of companions more massive than 5 M_⊙ and orbital periods shorter than one year for all seven WRs. This result implies that massive single stars at SMC metallicity lose the majority of their hydrogen-rich envelope without requiring a binary companion.

Stellar-structure codes are currently not able to reproduce the observed properties of single WRs in the SMC (Hainich et al. 2015; Schootemeijer et al. 2024). Without invoking binary interactions, reproducing the population of WRs is a difficult task, even at solar metallicity, as the predicted Hertzsprung-Russell (HR) diagram positions are commonly off with respect to their effective temperatures and luminosities. A temperature discrepancy arises as evolutionary models tend to produce WRs that are too hot (T_eff > 100 kK) compared to the observed ones (40 < T_eff/kK < 100). This is also referred to as the WR radius problem (e.g., Hillier 1987; Langer et al. 1988; Hamann et al. 2006; Grassitelli et al. 2018). Detailed atmosphere models taking into account the effect of the (hot) iron opacity bump on the wind dynamics (Gräfener & Hamann 2005; Sander et al. 2020, 2023) could demonstrate that for stars close to the Eddington limit, strong winds with high mass-loss rates can be launched, which effectively cloak the hydrostatic layers of the star. In these cases, the effective temperature referring to the hydrostatic layers is significantly higher than the observed effective temperature, as the wind can remain optically thick out to many stellar radii. In regimes where the iron bump is not able to launch a wind (Moens et al. 2022; Sander et al. 2023), the proximity of the star to the Eddington limit could still lead to a considerable inflation of the envelope (e.g., Gräfener et al. 2011), which also effectively makes the star look cooler. Hence, the temperature or radius discrepancy can – at least partially – be attributed to a problem of the structure models in predicting accurate effective temperatures for stars with optically thick winds.

The luminosity discrepancy, instead, comes from the fact that theoretical models, even at solar metallicity, are able to reproduce only the most luminous WRs (log L/L_⊙ > 5.2), whereas they struggle with the low-luminosity ones. This happens because they are able to peel off the envelope only for the most massive stars M > 40 M_⊙, which have stronger winds, and not for initial masses in the range 20 M_⊙ < M < 40 M_⊙. Some codes, such as FRANEC (Chieffi & Limongi 2013), succeed in peeling off the H-rich envelope for all the stars with initial masses of M > 20 M_⊙ thanks to rotation, but only at Galactic metallicity (Z ∼ 0.009 − 0.015). Other codes struggle even at Galactic metallicity (see Sander et al. 2019; Hamann et al. 2019). With significantly weaker winds at SMC metallicity, reproducing observed WRs with single stellar evolution seems almost impossible.

Another aspect of the same problem is that current evolutionary models overpredict the number of over-luminous supergiants falling inside the Humphreys-Davidson (HD) limit (Humphreys & Davidson 1979). Gilkis et al. (2021) demonstrated that current stellar evolution models overproduce very luminous supergiant stars, even when varying parameters such as rotation, semi-convection, and overshooting. Such an overproduction of luminous supergiant stars in theoretical models might be due to the mass-loss prescription adopted for massive stars. In particular, if stellar winds are too weak, the star is not able to lose its envelope, and it inflates even during the core hydrogen burning phase and remains inside the HD limit until the end of its life (see e.g. Sabhahit et al. 2022, and references therein).

Mass-loss plays a key role during the evolution of massive stars and can determine their fate. Hot massive stars lose mass through radiation-driven stellar winds. The CAK model provides the traditional description of radiation-driven winds (Castor et al. 1975). This model is able to match the fundamental properties of OB star winds (e.g. Friend & Abbott 1986; Pauldrach et al. 1986), but, among other issues, it struggles at explaining the properties of the most massive stars, especially if they are close to their Eddington limit (Gräfener & Hamann 2008; Gräfener et al. 2011; Vink et al. 2011). While an increase of the mass-loss rate very close to the Eddington limit, i.e., Γ_e ≈ 1 (where Γ_e ≡ L/L_e is the ratio between the stellar luminosity L and the Eddington luminosity L_e), in principle is already expected from CAK, Vink et al. (2011) found a significant change of slope, a kink, for high Γ_e with a very steep trend of $\dot{M}\propto {\mathrm{\Gamma }}_{\text{e}}^{4.77}$. Observationally, there is also a change of spectral appearance from Of to WNh-type, even on the main sequence (e.g., de Koter et al. 1997; Crowther & Walborn 2011), meaning that the stars are eventually classified as Wolf-Rayet stars due to their emission-line spectra independent of their evolutionary status. Bestenlehner et al. (2014) observationally studied this transition regime, confirming the findings from Vink et al. (2011) based on Monte Carlo wind models.

As we discuss further in Sect. 2.2, the change in the Ṁ(Γ_e) behavior can be related to a switch from an optically thin to an optically thick wind regime. Sabhahit et al. (2022) developed a formalism to account for this switch and applied it to reconcile the temperature evolution of very massive stars (M > 100 M_⊙) to observations in the Milky Way and the Large Magellanic Cloud (LMC). In particular, they showed that by switching to the optically thick regime, very massive stars avoid envelope inflation and spend their main sequence evolution in a narrow range of temperatures that are in agreement with observations. Sabhahit et al. 2023 (2023, hereafter S23) extended this formalism to lower metallicities (down to Z_⊙/100) to study the threshold of pair-instability supernovae.

In this work, we show that the optically thick wind model of S23 is able to explain the self-stripping of single WRs at SMC metallicities when combined with fast stellar rotation. To avoid the luminosity discrepancy, we applied the S23 formalism to stars with a lower mass (M ≥ 20 M_⊙). The paper is organized as follows: Section 2 presents our methods, describing the code used for our simulations and quickly summarizing the S23 formalism for the optically thick wind regime. In Section 3.1, we show our main results and their dependence on metallicity and rotation. Section 3.2 discusses additional observable features of our simulated stars (e.g., surface abundances). In Section 4, we discuss some possible modifications and implications of our model. Finally, Section 5 presents a summary of our work.

2. Methods

2.1. General setup

We used the Modules for Experiments in Stellar Astrophysics code (MESA; version r12115, Paxton et al. 2011, 2013, 2015, 2018, 2019) to evolve stars with different zero-age main sequence (ZAMS) masses M_ZAMS/M_⊙ = 20, 25, 30, 40, 50, 60, 70, and 80. We simulated six different metallicities that might be representative of different regions of the SMC Z = 0.002, 0.0025, 0.003, 0.0035, 0.004, 0.0045, and eight initial rotation speeds Ω/Ω_c = 0, 0.3, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, with critical velocity Ω_c² = (1 − Γ_e) G M/R³, where Γ_e ≡ L/L_e = χ_e L/(4π G c M) is the electron scattering Eddington parameter (L_e and χ_e being respectively the Eddington luminosity and the electron scattering opacity), G is the gravity constant, M the stellar mass, and R the stellar radius. The rotational mixing efficiency is set by the ratio of the turbulent viscosity to the diffusion coefficient f_c (the amDmixfactor parameter of MESA) and the ratio between the actual molecular weight gradient and the value used for computing the mixing coefficients f_μ (amgradmufactor), calibrated as in Heger et al. (2000): f_c = 1/30, f_μ = 0.05.

We used the Grevesse & Sauval (1998) initial composition and opacity tables. For the equation of state (EOS), we used the default tables available in MESA, that are mainly based on the OPAL EOS tables (Rogers et al. 1996). Convective regions are defined by the Ledoux criterion, with mixing length parameter α_MLT = 1.5 and semi-convective diffusion coefficient α_sc = 1.

We divided our simulations into two steps. The first roughly corresponds to the main sequence and lasts until the central hydrogen fraction drops below 10⁻², and the second represents a more advanced phase and lasts until the central temperature reaches the value log(T_c/K) = 9.55. With this choice, we ran all of our stellar models at least until the onset of O-burning. During the first phase, we adopt the default nuclear reaction network in MESA (basic net), whereas in the second phase we use the coburn net (parameter autoextendnet = true), which includes ²⁸Si, in order to have a more complete and extended network for C- and O-burning and α-chains with respect to S23.

As in S23, in the first evolutionary step, the MLT++ prescription is switched off, while in the second step it is switched on. MLT++ is a numerical patch to standard MLT used in MESA to artificially reduce the superadiabatic gradient in radiation-dominated, convective envelopes. Following S23, the diffusion coefficient in the overshooting region is described by the exponential criterion with overshooting coefficient f_ov = 0.03 in the first step. The choice of f_ov = 0.03 during the core hydrogen burning phase is motivated by several recent studies, which indicate that massive stars require increased convective boundary mixing than low-mass stars in order to match observations¹. In Section 4 we discuss how variations of f_ov affect the stellar evolution models.

Mixing in later evolutionary stages, instead, is more difficult to constrain and remains highly uncertain. The He-burning core develops a much steeper composition gradient compared to the H-burning core. This strong chemical discontinuity above He core can restrict the growth of the convective region (Langer 1991). For this reason we choose a smaller value of f_ov = 0.01 during the advanced phase, as done by S23. We have verified that this choice produces only minor differences in the evolution of our stellar tracks.

2.2. Wind mass-loss

For the wind mass-loss, we closely followed the model implemented by S23, of which we give a brief recap here. S23 implement a consistent procedure to find the switch point between optically thin and thick winds Γ_e,switch as a function of stellar parameters, and adopt the steep Γ_e scaling for optically thick winds found by Vink et al. (2011). Specifically, S23 follow the formalism by Vink & Gräfener (2012), who found a model-independent way to characterize the transition from optically thin winds of O-type stars to optically thick winds of WNh stars, irrespective of the assumptions on the wind clumping factor. They derived a relation between the wind efficiency parameter, η, and the optical depth parameter, τ:

$$\begin{array}{c}\hfill \eta \equiv \frac{M{v}_{\infty }}{L/c}=f\tau .\end{array}$$(1)

The switch between optically thin and thick winds occurs when τ = 1 and corresponds to a wind efficiency parameter η_switch = f ≃ 0.6.

Vink & Gräfener (2012) have prescribed a way to derive the transition mass-loss rate, Ṁ_switch, by combining this equation with empirical observations of luminosities and terminal velocities of slash (i.e., Of/WNh) stars. If applied to slash stars in the Arches cluster, at Galactic metallicity, the derived Ṁ_switch agrees well with the empirically determined one assuming a clumping factor of D = 10 (Martins et al. 2008) and with the optically thin mass-loss recipe from Vink et al. (2001), confirming that such a model for optically thin winds correctly predicts the mass-loss rates at the transition.

S23 employed a set of stellar atmosphere models computed with the dynamically consistent branch of the Potsdam Wolf–Rayet (PoWR) code (Sander et al. 2017, 2020) to study the dependence of f on various stellar parameters. They find that f mainly depends on the ratio v_esc/v_∞, where v_esc ≡ (2 G M/R)^1/2, and $v_{\infty }\equiv a{[2GM(1-{\mathrm{\Gamma }}_{\text{e}})/R]}^{1/2}{(Z/Z_{\odot })}^{0.2}$ (Lamers et al. 1995), with a = 2.6 above the bi-stability jump (T_eff > 25 kK) and a = 1.3 below it (T_eff < 25 kK). Specifically, S23 found an equation for the wind efficiency parameter at the switch:

$$\begin{array}{c}\hfill {\eta }_{\mathrm{switch}}=f=0.75{(1+\frac{v_{\text{esc}}^2}{v_{\infty }^2})}^{-1}.\end{array}$$(2)

This equation was initially derived for very massive stars, M > 50 M_⊙, and is in agreement with the f values reported by Vink & Gräfener (2012) for the stratification-wind parameter β = 1.5 (see their Table 2), which is an appropriate value for stars of this mass. Since we are interested in applying the S23 formalism to stars with M > 20 M_⊙, we checked what happens if we rescale the value 0.75 in Eq. (2) to be in agreement with the values of f reported by Vink & Gräfener (2012) for β = 1 at v_∞/v_esc = 1.5, 2.5. We found negligible differences, and thus we use η_switch in our work.

By means of Eq. (2), we can calculate the value of the Eddington factor Γ_e,switch at which the transition from optically thin to optically thick winds takes place η(Γ_e,switch) = η_switch(Γ_e,switch), where the mass-loss at the transition is computed via the optically thin prescription by Vink et al. (2001) $\dot{M}={\dot{M}}_{\text{V01}}(L,M,T,v_{\infty },Z)$, M is evaluated from the mass-luminosity relations under homogeneity assumption, and the radius is computed from the Stefan-Boltzmann law. The value of Γ_e,switch is therefore determined by the functions η and η_switch. An increase in η or a decrease in η_switch results in a lower Γ_e,switch and vice versa.

As in S23, in the case of rotation, we multiplied the wind mass-loss rate by a factor (Maeder & Meynet 2000):

$$\begin{array}{c}\hfill \frac{\dot{M}(\mathrm{\Omega })}{\dot{M}(\mathrm{\Omega }=0)}=\frac{{(1-\mathrm{\Gamma })}^{\frac{1}{\alpha }-1}}{{(1-\mathrm{\Gamma }-\frac{4}{9}{\left(\frac{v}{v_{\text{crit}}}\right)}^2)}^{\frac{1}{\alpha }-1}},\end{array}$$(3)

where Γ ≡ χ L/(4π G c M) is the total Eddington parameter (χ being the total opacity), α = 0.52, and $v_{\text{crit}}={[2GM/(3R_{\text{pb}})]}^{1/2}$, where R_pb is the polar radius at breakup, which we assumed does not change with rotation.

As in S23, we used the following recipes for wind mass-loss:

• For temperatures in the range 4000 K < T_eff < 10⁵ K, in the optically thin regime (η < η_switch) the standard Vink et al. (2001) prescription is used, while in the optically thick regime (η > η_switch), the Vink et al. (2011) prescription is used.

• For temperatures T_eff < 4000 K, we used the red supergiants prescription by de Jager et al. (1988).

• For temperatures T_eff > 10⁵ K, we used the maximum between the WR wind prescription by Sander & Vink (2020) and Vink (2017). We call these winds WR-type winds in the rest of the paper.

Throughout the paper, we compare the S23 wind model with a model where only optically thin winds with the recipe of Vink et al. (2001) are considered. We call this case V01.

3. Results

3.1. Self-stripping of WRs at SMC metallicities

Figure 1 shows our stellar tracks at metallicity Z = 0.0025, which corresponds to the SMC baseline metallicity computed from the element abundances reported in Vink et al. (2023). This figure compares stellar tracks in the V01 case, where only optically thin winds are used (left), with the S23 case, where optically thick winds can be activated (right) if the condition Γ_e > Γ_e,switch is met. We also compare models without and with rotation. For the cases with rotation, we assume initial rotation Ω/Ω_c = 0.6, corresponding to an initial rotational velocity of v_rot,i ≃ 450 km/s.

Fig. 1. Stellar tracks on the HR diagram for initial masses M = 25,  30,  40,  60, and 80 M⊙ at metallicity Z = 0.0025. In the left-hand panels, the V01 model with only optically thin winds is implemented. In the right-hand panels, the S23 model is enforced, with the possibility to activate optically thick winds. The upper panels show the case with no rotation, while the lower panels the Ω/Ωc = 0.6 case. Different line styles represent different wind regimes: dotted for optically thin winds, solid for optically thick winds, dashed for WR-type winds. Colors represent different burning stages, defined as the element whose burning generates most of the energy of the star. Blue is for H-burning through CNO, orange for He-burning through triple α, red for carbon burning, and green for oxygen burning. Green circles are observations of single WRs (Hainich et al. 2015), filled for the hottest, open for the coldest. The lower-right plot shows that the combination of optically thick winds and high initial rotation is able to self-strip stars even at SMC metallicity. Stars in the lower-right plot (S23, Ω/Ωc = 0.6) avoid the HD limit and transition from main-sequence stars to WRs.

Figure 1 shows that optically thin winds at the metallicity of the SMC (top left) are not able to peel off the stellar envelope and the stars evolve to the cool side of the HR diagram after the main sequence, entering and remaining inside the HD limit if they are massive enough (M ≳ 40 M_⊙).

A rotation as high as v_rot,i ≃ 450 km/s (bottom left) changes the main sequence evolution of the star, which enters the chemically homogeneous evolution regime, becoming more luminous and more compact, but it is not enough to peel off the envelope. Even at this high initial rotation speed, if only thin winds are allowed, all stars end up being cool supergiants, entering the HD limit at the end of their life.

Contrariwise, in the S23 model, some stars are able to switch to optically thick winds, peel off their hydrogen-rich envelope, move to the hot side of the HR diagram, and avoid the HD limit (Romagnolo et al. 2024). In the case with no rotation, this happens only for the most massive stars (M ≃ 800 M_⊙), as already shown by S23. However, for rotating stars with Ω/Ω_c = 0.6, this behavior extends to much lower initial masses (M ≃ 25 M_⊙).

This happens because of several concurring factors. The most straightforward one is the boost factor multiplying the wind for rotating stars (Eq. 3). This increases the strength of the wind and the wind efficiency parameter η and, consequently, reduces Γ_e,switch, as also shown in Fig. 7 by S23.

However, the most important factors are related to the evolution of the star during the main sequence. Fast rotating stars, as the ones shown in Figure 1, bottom panels, evolve by chemically homogeneous evolution and become more luminous and more compact than non-rotating stars. The increase in luminosity has a twofold effect: (i) Γ_e increases, (ii) the wind strength Ṁ and wind efficiency parameter η increase, reducing Γ_e,switch. Furthermore, the smaller radius of rotating stars contributes to the increase of v_∞ and consequently to a further increase of η, which reduces Γ_e,switch. All these effects concur in facilitating the activation of optically thick winds, even for less massive stars. If optically thick winds are activated before the onset of core He-burning, the star is able to peel off its hydrogen-rich envelope and becomes hot and compact, reaching temperatures > 10⁵ K and avoiding the HD limit.

Figure 1 demonstrates that the combination of optically thick winds and rotation might be the key to reproduce single WR observations in the SMC, especially the five hottest single WR stars highlighted by the solid green circles in the Figure². In our models, we overcome the claimed luminosity discrepancy as stars with initial mass as low as 25 M_⊙ can self-strip their envelope. We can reproduce the empirical HR diagram location of the single WRs with 5.6 < log(L/L_⊙) < 6 with rapidly rotating stars with initial masses in the range 25 < M/M_⊙ < 40formula>.

The ability of stars in this mass range to lose their envelope also helps to address another issue: the observational lack of massive WR progenitors in the SMC. Current evolutionary models predict that luminous single WRs at low Z originate from very massive O-type stars (> 40 M_⊙). However, observations (Ramachandran et al. 2019; Schootemeijer et al. 2021) reveal a scarcity of such massive O-type stars in the SMC. Furthermore, even some of the seemingly luminous O-type stars in the SMC are actually stripped and have lower masses (Pauli et al. 2022). Given that the helium burning WR phase constitutes only about ∼10% of a star’s lifetime, we would expect to observe hundreds of massive O-type stars to account for the current WR population. This discrepancy suggests that WRs in the SMC likely evolved from less massive O-type stars, consistent with our findings that enhanced envelope stripping via thick winds can enable such evolutionary pathways.

Rapidly rotating stars with initial velocity Ω/Ω_c ≳ 0.6, corresponding to v_rot ≳ 450 km/s, are not representative of the bulk of the O-type stars’ rotational velocity distribution in the SMC. We compiled v sin i values estimated for SMC O-type stars based on the recent literature from the BLOeM survey (Bestenlehner et al. 2025), which covers the main bar of the SMC, combined with data from (Ramachandran et al. 2019; Rickard et al. 2024), that primarily cover the SMC wing and the core of the NGC346 star cluster. This combined dataset resulted in a total of 171 O-type stars, and the distribution of their observed v sin i and true rotational velocities (v_rot) are shown in Figure 2. To determine the intrinsic rotational velocity distribution, we employed the Lucy-Richardson deconvolution method, an iterative technique that converges to a maximum likelihood solution, as demonstrated by Lucy (1974). Our analysis indicated that O-type stars with v_rot ≳ 450 km/s approximately represent ∼10% of the population in the SMC. Most of the observed O-type stars are no longer on the ZAMS, implying they may have already spun down. Therefore, even a small number of rapidly rotating O-type stars might be enough to explain the observed seven single WRs in the SMC.

Fig. 2. Velocity distribution of O-type stars in the SMC. The plot shows the normalized observed v sin i distribution (blue) and the reconstructed initial rotational velocity (vrot, red) for 171 O-type stars. The dataset is a compilation from Ramachandran et al. (2019), Rickard et al. (2024), and Bestenlehner et al. (2025). About 10% of the O-type stars feature vrot > 450 km/s.

3.1.1. The dependence on metallicity and rotation

In this section, we aim to describe the complex dependence of the activation of thick winds on metallicity and initial rotation. We focus on stars with mass M < 40 M_⊙, which are the most important to reproduce the observed WRs in the luminosity range 5.6 < log(L/L_⊙) < 6. In principle, one would expect weaker winds at a lower metallicity. While this remains generally true, the activation of the optically thick wind regime follows a more complex behavior.

Stars at a lower metallicity are generally more compact and luminous than stars at a higher metallicity (Kippenhahn et al. 2013; Klencki et al. 2020; Gilkis et al. 2021). This has two effects: On the one hand, a higher luminosity means a higher Γ_e. On the other hand, a higher luminosity and compactness increase the wind efficiency parameter, η, reducing Γ_{e, switch} and facilitating the activation of thick winds. This scenario is further complicated by the bi-stability jump (Vink et al. 1999, 2000, 2001), which enhances wind strength for stars cooler than 25 kK.

To qualitatively describe the impact of all these effects, we refer to Figure 3. While most of the very massive stars simulated in S23 enter the optically thick regime during the main sequence thanks to their high luminosities, lower mass stars (as the ones simulated in this work) can enter the optically thick wind regime only at the end of the main sequence. When core hydrogen burning stops, the star undergoes a rapid phase of contraction, usually referred to as ’hook’ in the HR diagram. This rapid contraction phase leads to an increase in the terminal velocity v_∞, with a consequent increase in the wind efficiency parameter η and a reduction of Γ_{e, switch}. Moreover, if the star is rapidly rotating, this phase is followed by a quick increase in luminosity and Γ_e. For these reasons, this is the first good spot for a star with mass M < 100 M_⊙ to enter the optically thick regime, provided that its luminosity is high enough and the star is sufficiently compact. We refer to this activation scenario as channel 1.

Fig. 3. Illustration depicting the possible evolutionary paths of a chemically homogeneous star with initial mass ≳ 20 M⊙ at SMC metallicity. After core H-burning (yellow ball with blue center), the star contracts, Γe,switch decreases and Γe increases. If Γe > Γe,switch the star enters the optically thick winds regime (left arrow) and ends up being a WR (channel 1). If Γe < Γe,switch (right arrow), winds are still optically thin during shell H-burning (yellow ball with blue shell) and the star expands and cools, increasing Γe,switch. At this stage, there are two possible outcomes: (i) If the star cools enough to cross the bi-stability jump Teff < 25 kK (right arrow), it may enter the optically thick wind regime and become a WR (channel 2). (ii) If the star starts core He-burning at Teff > 25 kK (left arrow), Teff rises again, and the star does not cross the bi-stability jump. Optically thick winds are not activated (or are activated too late) and the star ends up being a cool supergiant.

If the star does not enter the optically thick regime through channel 1, the contraction phase is followed by the stellar expansion during shell hydrogen burning. During the expansion, both the terminal velocity and the rotation boost factor decrease. Consequently, η decreases and Γ_e,switch increases, and the star remains in the optically thin regime. At this point two things may occur: (1) Core He-burning begins when the star is still above the bi-stability jump. (2) The star becomes cold enough (T_eff < 25 kK) to cross the bi-stability jump before the onset of core He-burning. In case (1), core He-burning stops the expansion of the star, thus preventing it from crossing the bi-stability jump. In this case, winds remain optically thin for all the core He-burning phase. Afterward, some stars are able to activate optically thick winds, but that happens too late in their evolution to remove a substantial fraction of the envelope. These stars will end their lives as cool supergiants.

In case (2), instead, the bi-stability jump is crossed before core He-burning. The bi-stability jump reduces the terminal velocity by a factor of ∼2 (Vink et al. 2001). This has a twofold effect. First, since $\eta \propto {\dot{M}}_{\text{V01}}{v}_{\infty }\propto v_{\infty }^{-0.226}$, a reduction in v_∞ leads to a slight increase in η. Second and most important, the ratio $v_{\text{esc}}^2/v_{\infty }^2$ increases by a factor of ∼4, leading to a reduction of η_switch and, consequently, of Γ_{e, switch}. This drastic reduction of Γ_{e, switch} can lead to the activation of optically thick winds if the star is luminous enough (see also S23). This is the second good spot for the star to enter the optically thick wind regime. We call this channel 2. This general behavior is depicted in Figure 3.

To see these effects in more detail, the upper panel of Figure 4 shows three different tracks for a star with initial mass M = 25 M_⊙ and initial rotation Ω/Ω_c = 0.6 at three different metallicities: Z = 0.002, 0.003, and 0.004. As said above, lower metallicity tracks have higher luminosities and temperatures. The Z = 0.002 track is luminous enough to activate optically thick winds right after the contraction phase at the end of the main sequence (channel 1).

Fig. 4. Stellar tracks for a 25 M⊙ star. The upper panel shows a rotating star with Ω/Ωc = 0.6 for three different metallicities Z = 0.002,  0.003,  0.004. The lower panel shows a star with metallicity Z = 0.003 for three different values of the initial rotation speed Ω/Ωc = 0.5,  0.6,  0.7. The linestyles and color code are the same as in Figure 1. The vertical dotted black line represents the temperature below which we expect the bi-stability jump. This Figure shows the non-monotonic dependence of mass-loss on metallicity and rotation. Low metallicity (Z = 0.002) favors the activation of optically thick winds through channel 1, because the star is more luminous and compact. Higher metallicity (Z = 0.004) can also activate optically thick winds through channel 2. Stars in the intermediate metallicity case (Z = 0.003), instead, fail to activate optically thick winds soon enough and end up as cool supergiant stars. The same trend happens for rotation.

The Z = 0.003 track, instead, is not luminous enough to activate optically thick winds at the end of core H-burning and evolves toward cooler temperatures. At log(T_eff/K)≃4.5, core He-burning starts and the star contracts again. It is only later, during carbon burning, that the star crosses the bi-stability jump and activates optically thick winds. This is too late to peel off the hydrogen-rich envelope, and the star ends up as a cool supergiant, eventually crossing the HD limit.

The Z = 0.004 track is even cooler and not luminous enough to activate optically thick winds through channel 1. It evolves toward cooler temperatures and crosses the bi-stability jump before core He-burning. In our model, this leads to the activation of optically thick winds through channel 2, and the star moves to the hot region of the HR diagram.

Therefore, the metallicity dependence is not monotonic. At lower metallicities, optically thick winds activate because of the high luminosity and hot temperatures (channel 1). At higher metallicities, optically thick winds activate because cooler stars can reach the bi-stability jump before core He-burning (channel 2). At intermediate metallicities, instead, the stars struggle to enter the optically thick wind regime, as they are not luminous enough to enter the optically thick wind regime at the end of the main sequence, but they are too hot to cross the bi-stability jump before the onset of core He-burning.

A similar behavior applies to rotation, with higher rotations corresponding to more luminous and hotter stars and lower rotations to less luminous and cooler ones. This is shown in Figure 4 (lower panel) where we display three different tracks for a star with initial mass M = 25 M_⊙ at metallicity Z = 0.003, with initial rotation Ω/Ω_c = 0.5,  0.6, and 0.7. The fastest rotating star in our example (Ω = 0.7 Ω_c) activates optically thick winds through channel 1, while the one rotating more slowly (Ω = 0.5 Ω_c) becomes a WR through channel 2. The star with intermediate rotation velocity (Ω = 0.6 Ω_c), instead, is not able to activate any of the two channels. However, we note that an even slower rotation (Ω < 0.5 Ω_c in our example) will not activate optically thick winds even if the star crosses the bi-stability jump before core He-burning. This happens because the star needs enough rotation to enter the chemically homogeneous evolution regime; otherwise its luminosity would be too low to activate optically thick winds.

A similar dependence is found for the overshooting parameter, which has a very similar effect to rotation (see Section 4.1). A more detailed analysis of the behavior of Γ_e, ${\mathrm{\Gamma }}_{\text{e, switch}}$, η and η_switch can be found in Appendix A. The outcome of all the simulations with different masses, metallicities, and initial rotations can be found in Table B.1 in Appendix B.

As a final remark, we note that while channel 1 implies that the star enters the optically thick regime due to its compactness and luminosity, channel 2 heavily relies on the bi-stability jump and thus has to be considered more uncertain. The existence and strength of the bi-stability jump have intensely been debated over the last few years. We refer to Appendix E for additional discussion about the bi-stability jump.

3.2. WR evolution and their properties

We have shown that even stars with mass down to M ∼ 20 M_⊙ can develop optically thick winds if their initial rotation is high enough. Winds peel off the stellar envelope and make the star move to the hot part of the HR diagram. We now describe additional properties of these stars and assess whether they are compatible with single WRs in the SMC. We consider different properties of our simulated stars: the time they spend in a given region of the HR diagram, their surface abundances, their average rotation velocity, and the transformed mass-loss rate Ṁ_t. Given the uncertainties on the bi-stability jump mentioned above, we take as reference the case at metallicity Z = 0.002 and Ω/Ω_c = 0.6, where all stars become WRs through channel 1.

3.2.1. Time evolution on the HR diagram

To calculate the time spent by our simulated stars in each region of the HR diagram, we first defined the quantity

$$\begin{array}{c}\hfill v_{\mathrm{HRD}}=\sqrt{{\left|\frac{dlogL}{\mathit{dt}}\right|}^2+{\left|\frac{dlogT}{\mathit{dt}}\right|}^2},\end{array}$$(4)

which represents the “transition” velocity of the star on the HR diagram and has units of dex yr⁻¹. The inverse of v_HRD, ${\tau }_{\mathrm{HRD}}\equiv 1/v_{\mathrm{HRD}}$, with units of yr dex⁻¹, measures the time spent by the star to move by 1 dex in the HR diagram.

The upper panel of Figure 5 shows stellar tracks on the HR diagram for Z = 0.002 and Ω/Ω_c = 0.6 colored by τ_HRD. The lower panel displays the same tracks with a cross marker every ∼ 5 × 10⁴ yr, colored by burning stage. It is clear that all stars spend most of their post main sequence lifetime in two regions: (1) near the turn-off point, when they start losing their envelope and moving blueward, and (2) during the hot WR stage. The exact effective temperatures of these stages vary with stellar mass.

Fig. 5. Time spent on different regions of the HR diagram for stars with mass M = 20,  30,  40,  60, and 80 M⊙ at Z = 0.002 and Ω/Ωc = 0.6. In the upper panel, the color code represents the quantity τHRD, in yr dex−1, which is the time spent to move by 1 dex in the HR diagram. Green circles are WRs in the SMC ,and blue squares are O-type stars in the SMC with v sin i > 200 km/s from (Mokiem et al. 2006; Bouret et al. 2013, 2021; Ramachandran et al. 2019; Dufton et al. 2019; Rickard et al. 2024; Backs et al. 2024). In the lower panel, the color code represents different burning stages: core H-burning (blue), He-burning (orange), carbon burning (red), and oxygen burning (green), while the line styles represent different wind regimes: dotted for optically thin winds, solid for optically thick winds, dashed for WR-type winds. Cross markers are spaced by ∼ 5 × 104 yr each. After the main sequence, stars spend most of their time at the turn-off point, where they start to peel off their H-rich envelope and move blueward, and in the hot WR phase, at 4.9 < log(Teff/K) < 5.2.

For lower mass tracks (20−40 M_⊙), region (1), corresponding to early core He burning, is at log(T_eff/K) ≳ 4.6 , while region (2) is at log(T_eff/K) ∼ 5−5.1 during mid/late core He-burning. We accurately reproduce the position of the lowest luminosity observed WR at log(T_eff/K) ∼ 5, where the M = 20 M_⊙ track spends most of its post main sequence time. Tacks with 30−40 M_⊙, instead, tend to have slightly higher effective temperatures than observed stars, by ∼ 0.1 dex.

The 20−40 M_⊙ tracks intersect the HR positions of the five observed hot WRs during mid-to-late core He burning. Although the “WNh” label is often associated with core hydrogen-burning stars having a WN-type appearance, we can exclude this possibility here. Extending the main sequence to effective temperatures as high as log(T_eff/K) ∼ 5 would require fully chemically homogeneous evolution, which in turn would imply very low hydrogen abundances throughout the star, including the surface. This is at odds with their observed WNh classification. Alternatively, assuming high hydrogen surface abundances would reflect a large, hydrogen-rich shell, and would imply a stellar mass that is so high that the resulting L/M-ratio would not yield a WR-type spectrum (cf. Shenar et al. 2020; Sander 2024a,b).

Overall, we interpret the five hot WNh in the SMC as being in a similar evolutionary stage as WN-early stars observed in the Milky Way or the Large Magellanic Cloud except that they retain some surface hydrogen, likely due to weaker winds, which preserves their WNh spectral features. Instead, the two cooler WNh stars can be interpreted as late main sequence or early core He burning objects. From Figure 5, the 2:5 ratio of cold to hot WNh roughly corresponds to the relative time stars spend in these two HR regions. However, we warn that this ratio may substantially depend on stochastic fluctuations and on the details of the star formation history of the SMC, which we do not consider in this work. Thus, we refrain from drawing further conclusions based on such number counts.

For higher mass tracks (60−80 M_⊙), region (1) is at log(T_eff/K) ≳ 4.4 and region (2) at log(T_eff/K) ∼ 4.9. These high-mass stars also spend some time at higher temperatures, log(T_eff/K) ∼ 5.2, corresponding to carbon burning stage. The absence of observed WRs in the SMC in this mass regime is not surprising given the general deficiency of very luminous massive stars in the SMC and it is likely the result of its star formation history (Ramachandran et al. 2019; Schootemeijer et al. 2021).

In the upper panel of Figure 5, we also report data for O-type stars in the SMC with high projected rotational velocity (v sin i > 200 km/s). Most of these stars cluster around the main-sequence evolutionary tracks for 20, 30, 40 M_⊙ stars, where they spend most of their lifetime. However, one rapidly rotating O-type star lies right after the hook, at log(L/L_⊙)≲6, where the star temporarily stalls before evolving toward the hotter region of the HR diagram.

The results shown here are for Z = 0.002, where we find WRs to be mainly produced by channel 1. In Appendix C, we present the same analysis for Z = 0.004, where channel 2 is the primary WR formation channel. From Figure C.1, it is evident that hot WRs formed through channel 2 are more evolved than those formed through channel 1. Specifically, the channel 2 tracks are already in the core C burning phase when they reach the observed HRD location.

3.2.2. Surface abundances

Figure 6 shows the surface abundances at Z = 0.002 and rotation Ω/Ω_c = 0.6. The upper panel shows the evolution on the HR diagram and the color represents the WR type. We define as WNh (or slash) all stars with surface H abundance X_s in the range 5 − 40% and in the optically thick winds regime. The star is a WN if the surface H abundance drops below 5% and the ¹²C or ¹⁶O surface abundances are < 20%. If the ¹²C (¹⁶O) surface abundance is > 20%, we label the star as WC (WO). WN, WC, and WO stars are also called classic WRs.

Fig. 6. Upper panel: Classification of the star along the HR diagram: O-type (black), WNh (blue), WN (orange), and WC and WO (red). Lower panel: Surface abundances as a function of the effective temperature, Teff. Line styles are for different elements: hydrogen (dotted line), 14N (solid line), 12C (dashed line). Different colors are for different initial masses: 20 (blue), 30 (orange), 40 (red), and 50 M⊙ (green). Green triangles (circles) represent hydrogen (nitrogen) abundances for the five hottest WRs observed (Hainich et al. 2015). The gray horizontal band covers the ∼5 − 40% range; if the hydrogen abundance is inside this region, the star is considered a WNh.

The WNh phase begins at the end of the main sequence, when thick winds are activated, and lasts until log(T_eff/K)≳5.1. Afterward, surface hydrogen is completely removed by the wind, and the star has a WN phase before becoming carbon or oxygen dominated. The classic WR phase occurs at very high temperatures. Lower mass stars M ≤ 30 M_⊙, do not reach the carbon or oxygen dominated phase by the end of our simulation.

The bottom panel shows the evolution of surface abundances of hydrogen, ¹⁴N and ¹²C for different stellar masses, as a function of the effective temperature. Since the evolution of T_eff is not monotonic, we plot only the log(T_eff/K)≳4.9 regime, for the sake of clarity. In this regime, T_eff always increases with time. All stars feature ∼20% surface hydrogen for a huge part of their evolution, until log(T_eff/K)∼5.1, corresponding to ∼97% of the star lifetime. During this phase, they are classified as WNh, in agreement with observational data, although the observed WNh in the SMC show a wider range of surface hydrogen abundances, ∼20 − 50%.

This discrepancy was already discussed by Schootemeijer & Langer (2018), who showed that reproducing the effective temperatures of the hot SMC WRs through chemical homogeneous evolution leads to underestimation of their surface hydrogen content. In our scenario, chemical homogeneous evolution is instrumental for the activation of thick winds, but it is not the main driver of the transition to the hot side of the HR diagram. Because lower initial rotations are sufficient to trigger optically thick winds and envelope removal, our models retain somewhat higher surface hydrogen fractions than in a purely homogeneous scenario. While they still fall short of reproducing surface H abundances as high as ∼50%, this is an improvement compared to a purely chemically homogeneous evolution scenario. Moreover, Appendix C shows that models at Z = 0.004, which evolve through channel 2, are in better agreement with the observed surface hydrogen abundance. At log(T_eff/K)≳5.1, hydrogen is rapidly eroded by winds, which progressively expose the internal structure and, when surface hydrogen drops below 0.05%, we classify the star as a WN.

The nitrogen abundance increases quite rapidly during the main sequence as a result of the CNO cycle, rising from an initial value of ∼10⁻⁴ to an equilibrium value of ∼10⁻³ (see Appendix D). This enhancement is distributed nearly homogeneously throughout the star owing to rotational mixing. Once core He-burning begins, nitrogen in the core is quickly depleted, while the surface abundance remains unchanged. At log(T_eff/K)≳5.1, there is a further enhancement of surface nitrogen. This is due to the combined effect of envelope hydrogen burning and winds exposing the interior of the star. Afterward, the nitrogen surface abundance declines as it is eroded by winds, which expose the nitrogen depleted core of the star. The rise and fall of surface nitrogen occur more rapidly in more massive stars, due to their stronger winds. The predictions of our models are in agreement with the two data points at log(T_eff/K) < 5, featuring a surface ¹⁴N ∼ 2 × 10⁻³.

We note that here we show only the lowest metallicity case Z = 0.002. For an initial metallicity of Z = 0.004, the models predict an equilibrium nitrogen abundance of ∼2 × 10⁻³, in good agreement with the data (see Appendix C). Observed WRs in the SMC, however, display surface nitrogen enrichment already at log(T_eff/K)≳5, whereas in our models the enrichment becomes visible only beyond log(T_eff/K)≳5.1. Nevertheless, our models do produce a nitrogen “bump” around log(T_eff/K)∼5, but this enhancement is confined to layers just beneath the surface and only becomes observable once stellar winds peel off the outermost material (Appendix D).

Overall, we attribute this discrepancy in the timing of surface nitrogen enrichment to the poor understanding of the rotational mixing processes. As discussed in Section 2, we adopt the Heger et al. (2000) calibration of the mixing efficiency parameters f_c and f_μ, based on nitrogen abundances in solar metallicity terminal-age main-sequence stars. Alternative calibrations (e.g., Brott et al. 2011; Costa et al. 2019) may be more appropriate for SMC metallicities. However, a detailed tuning of these parameters is beyond the scope of the present work.

Another possible explanation for the mismatch in the effective temperature of the nitrogen enrichment could be traced back to the radius problem. Since simulated WRs tend to be at a higher T_eff than the observed ones by ∼0.1 dex (see also Figure 5), a shift of the data to higher temperatures by this amount would allow the effective temperature of the nitrogen enrichment to be matched.

Carbon surface abundance follows a trend similar to that of nitrogen. At log(T_eff/K) < 5.1, it remains nearly constant at ∼2 × 10⁻⁵, independent of stellar mass. Beyond this point, the surface carbon abundance increases as stellar winds expose deeper layers. Oxygen features a very similar trend as carbon, but it is not shown in the figure for clarity. By the end of the evolution, more massive stars (40−50 M_⊙) reach ≳80% in ¹²C + ¹⁶O abundance, being classified as WC or WO. In contrast, the 20 and 30 M_⊙ tracks have not reached the WC or WO phase even at temperatures as high as log(T_eff/K)∼5.3, with their surface carbon abundance still increasing with temperature.

3.2.3. Rotation

Figure 7 shows the evolution of the rotation velocity of the star in the HR diagram (top) and as a function of T_eff (bottom). The star needs a high initial rotation to develop chemically homogeneous evolution (Brott et al. 2011). However, after approximately half of the stellar lifetime, the velocity drops below 400 km/s, reaching values more similar to those observed for the bulk of the O-type star population. The drop is faster for more massive stars, with the 50 M_⊙ star that can reach a rotation speed as low as ∼100 km/s during the main sequence. After the end of core H-burning the star contracts and the rotational velocity increases, before dropping to almost zero due to the activation of optically thick winds which remove angular momentum. Finally, toward the end of the evolution, when log(T_eff/K)≳5.05 (≳97% of the star lifetime), the rotational velocity increases again to values ∼100 − 250 km/s due to winds exposing the rotating core of the star.

Fig. 7. Upper panel: Average rotational velocity along the HR diagram (color code). Lower panel: Evolution of the average rotational velocity as a function of the Teff. The color code is the same as in Figure 6. Dashed lines represent the main sequence evolution, solid lines the advanced phases. The velocity drops during the main-sequence evolution. This is followed by the hook at the end of the main sequence, when the velocity quickly rises due to post main sequence contraction, and by a subsequent drop. In the late phases of the evolution log(Teff/K)≳5.05, the velocity rises again due to winds peeling-off the outer layers of the star and exposing the rotating core.

While we do not have measures of the rotational velocity of WRs, there are some broad upper limits on v sin i reported in the literature. From the shape of narrow absorption/emission lines Martins et al. (2009) find v sin i < 50 km/s for WRs in the SMC. Hainich et al. (2015) revised this estimate by considering that spectral lines of WRs are formed in the stellar wind, obtaining less stringent velocity constraints of v sin i < 200 km/s. Using spectropolarimetry, Vink & Harries (2017) found a quite low ∼10% incidence of line effects on WRs in the LMC and detected only one line effect in the SMC WRs. This implies that observed WRs should feature moderate rotational velocities. Our tracks are compatible with these upper limits. They start to rise to v_rot > 100 km/s when the temperature log(T_eff/K)≳5.05 due to the exposition of the core. Since all observed WRs in the SMC have log(T_eff/K) < 5.05, they are still in the phase with v_rot < 50 km/s, consistent even with the most stringent upper limits.

3.2.4. Mass-loss

Figure 8 shows the evolution of mass-loss and effective temperature as a function of the normalized stellar age, defined as the fraction of the total stellar lifetime. During the main sequence, mass-loss is driven by optically thin winds, whose strength gradually increases as luminosity rises. At the end of hydrogen burning, the star contracts, both effective temperature and luminosity increase, and optically thick winds are activated through channel 1. This correspond to an enhancement of mass-loss of about one order of magnitude. In the optically thick-wind regime, mass-loss remains approximately constant or shows a slight decline due to the modest luminosity decrease. At log(T_eff/K) > 5 WR-type winds are activated, and the mass-loss rate can exhibit sharp variations owing to its steep dependence on Γ_e.

Fig. 8. mass-loss, Ṁ, (upper panel) and log(Teff/K) (lower panel) as a function of normalized stellar age. The color code is the same as in Figure 6. Line styles represent different wind regimes: dotted for optically thin winds, solid for optically thick winds, dashed for WR-type winds.

We computed the “transformed mass-loss rate” (Gräfener & Vink 2013), which is defined as

$$\begin{array}{c}\hfill {\dot{M}}_{\text{t}}=\dot{M}\sqrt{D}\left(\frac{1000\mathrm{km}{\mathrm{s}}^{-1}}{v_{\infty }}\right){\left(\frac{{10}^6{\mathrm{L}}_{\odot }}{L}\right)}^{3/4},\end{array}$$(5)

where D is the clumping factor assumed to be D ≃ 10 when calculating it from the evolution models. The quantity Ṁ_t provides the (unclumped) mass-loss rate a star would have if it had a luminosity of 10⁶ L_⊙ and a wind terminal velocity of 1000 km s⁻¹. Stars with the same temperature and Ṁ_t have approximately the same normalized line spectrum (see also Schmutz et al. 1989). As the “raw” mass-loss rate increases considerably with the luminosity, Ṁ_t is a better quantity to compare the wind strength of stars with different luminosities and to judge whether our tracks not only lie in the same HR diagram position as the single WRs of the SMC, but also show a mass-loss rate approximately compatible with the observed strength.

Figure 9 shows the evolution of Ṁ_t in the HR diagram (top) and as a function of the effective temperature (bottom). With the activation of the optically thick winds, Ṁ_t rises to values ∼10^−4.5 M_⊙/yr. This is followed by a decrease, more evident for less massive stars, due to a drop in the applied mass-loss recipe. Afterward, when T_eff > 10⁵ K, Ṁ_t rises again due to the activation of WR-type winds (Vink 2017; Sander & Vink 2020). This is then followed by a shallow decrease, due to a decline in Γ_e, and by another increase at log(T_eff/K)≳5.1 due to the activation of core C burning, which increases the luminosity and thus Γ_e. The decline and rise at log(T_eff/K) = 5 is very strong for the 20 M_⊙ star because of two reasons: (i) the WR-type winds are extremely sensitive to Γ_e (see Sander et al. 2020), and thus a slightly lower Γ_e implies much weaker winds. (ii) The activation of core C burning for the 20 M_⊙ star occurs exactly at log(T_eff/K) = 5, causing the steep rise. We note that we do not include the temperature dependence on Ṁ from Sander et al. (2023), which could soften this mass-loss increase as Ṁ is expected to decline for more compact stars.

Fig. 9. Upper panel: Transformed mass-loss, Ṁt, along the HR diagram track (color code). Lower panel: Evolution of Ṁt as a function of Teff. The color code is the same as in Figure 6. The green points are computed from the observed Ṁ, L, and v∞ from Hainich et al. (2015).

The observed Ṁ_t-values are computed from Eq. (5) using the values of Ṁ, L, and v_∞ reported by Hainich et al. (2015). Our tracks, especially the ones for M ≥ 30 M_⊙, overestimate the value of Ṁ_t by a factor ∼10, predicting a mass-loss that is too high, at least during the WR phase. However, as already mentioned above, WR-type winds feature a steep dependence on Γ_e. At the observed T_eff, our tracks have Γ_e ∼ 0.6, while a value of Γ_e ∼ 0.5 would be more than sufficient to be in agreement with observations. Still, such an overprediction of the mass-loss should be carefully regarded for future refinements of our model, and points toward the direction of having less aggressive winds at least during the WR phase. This would also help in reproducing the higher surface hydrogen abundances of the WNh stars in the SMC (see Figure 6).

4. Discussion

4.1. Overshooting

The results shown in Sections 3.1 and 3.2, are obtained with an overshooting parameter f_ov = 0.03 during the first step of the simulation, corresponding to the main sequence, and f_ov = 0.01, during the advanced phase. We choose these values according to S23. However, while these values are expected for stars with mass ≲ 5 M_⊙ (Herwig 2000), there are some hints that the overshooting parameter might be higher ≳0.05 for stars with larger masses (Vink et al. 2010; Higgins & Vink 2019). Since overshooting enhances mixing, it has an effect similar to rotation, making the star more luminous and compact, and facilitating the activation of optically thick winds. In particular, for higher f_ov values a lower initial rotation is required to enter chemically homogeneous evolution and to activate thick winds. This is shown in Figure 10. The three panels represent the stellar tracks of a 20 M_⊙ star at Z = 0.0025, with three initial rotation values. The different tracks are for f_ov = 0.03,  0.045,  0.07 during the main sequence phase, while we leave f_ov = 0.01 for the advanced stage. A 20 M_⊙ star does not reach the optically thick wind regime for Ω/Ω_c = 0.55 and f_ov = 0.03. In contrast, by slightly increasing the overshooting parameter to f_ov = 0.045, the same star enters the thick wind regime and self-strips its envelope. The same occurs at Ω/Ω_c = 0.5 and f_ov = 0.07.

Fig. 10. Stellar tracks for a star of initial mass M = 20 M⊙ at Z = 0.0025 for initial rotation speed Ω/Ωc = 0.55 (left), 0.5 (middle), and 0.45 (right) and overshooting parameter fov = 0.03, 0.045, 0.07, and 0.1 (right panel) during the main sequence. The effect of overshooting is almost degenerate with that of rotation. A slight increase of fov corresponds to a reduction of the rotation threshold needed for self-stripping.

The right-hand panel shows that for very large overshooting values f_ov = 0.109, even stars with moderate initial rotation Ω/Ω_c = 0.45 can activate optically thick winds. Variations of the overshooting parameter help in reproducing observations of single WRs with lower initial rotations. An accurate investigation of the effect of overshooting and its relation with rotation and metallicity should be carried out with a population study and is deferred to a future work.

4.2. Thick winds and the HD limit

We have shown that a few rapidly rotating stars can evolve into single WRs at the metallicity of the SMC as an effect of optically thick winds. We now discuss whether the interplay of enhanced stellar winds and rotation may possibly address the problem of over-luminous cool supergiant stars crossing the HD limit. Answering this question goes beyond the goal of this work and requires a population study, which we will address in a follow-up paper. However, here we briefly speculate how optically thick winds might address the HD limit problem as well.

Figure 11 shows stellar tracks with different initial rotation velocities and with metallicity Z = 0.002 and Z = 0.004 to encompass the metallicity spread in the SMC. Since stars with initial mass ≤30 M_⊙ barely enter the HD limit, we focus on higher initial masses (M = 40, 50, 60 M_⊙).

Fig. 11. Stellar tracks for stars with mass M = 40,  50, and 60 M⊙ at Z = 0.002 (upper panels) and Z = 0.004 (lower panels), encompassing the metallicity range of the SMC. Columns represents different initial rotations: Ω = 0 (left), Ω/Ωc = 0.3 (middle), Ω/Ωc = 0.6 (right), corresponding to the minimum, average, and maximum observed rotation velocities of O-type stars in the SMC. The color code shows the time spent in different regions of the HR diagram τHRD.

For high initial rotation velocity Ω/Ω_c = 0.6, stars at the lowest considered metallicity (Z = 0.002) barely enter the HD limit, because of the activation of optically thick winds through channel 1. At higher metallicity, the situation is slightly worse. Indeed, at Z = 0.004, stars with mass M ≥ 40 M_⊙ are able to self-strip their envelope and become WRs. However, since they activate optically thick winds through channel 2, all of them enter the HD limit and spend there a fraction ∼2 − 2.5% of their total lifetime, before moving to the hot part of the HR diagram, where they spend the remaining ∼4% of their lifetime.

Reducing the initial rotation worsens the situation significantly for Z = 0.002 and only slightly for Z = 0.004. At Ω/Ω_c = 0.3, representing the rotation speed of most O-type stars in the SMC, only stars with M ≥ 60 M_⊙ exit the HD limit at Z = 0.002. This means that all stars with initial mass 40 < M/M_⊙ < 60 enter and die inside the HD limit. This threshold is lower (M ≥ 50 M_⊙) for Z = 0.004, but the amount of time spent inside the HD region increases to ∼3 − 6% of the total lifetime, which becomes comparable to the time spent as WR. In the case with no rotation, all the stars with M ≤ 80 M_⊙ enter and die inside the HD limit at Z = 0.002, while the threshold stays at M ≥ 50 M_⊙ for Z = 0.004.

Overall, the situation is complex: For high rotation velocities, the HD limit is almost empty at low metallicity (Z = 0.002), whereas for low or no rotation the HD limit is less populated at high than low metallicity. In any case, stars with initial mass 40 < M/M_⊙ < 60 tend to enter and remain in the HD limit if they are not rotating fast enough. As we have shown in Section 4.1, the situation could be improved by increasing the overshooting parameter, which facilitates activation of optically thick winds through channel 1. A population study such as the one performed by Gilkis et al. (2021), but with the inclusion of optically thick winds, is required to analyze the situation in more detail, changing metallicity and rotation on a finer grid and testing the effect of different overshooting parameters.

4.3. Consequences for black hole masses

The onset of optically thick winds has severe consequences also for black hole masses, since optically thick winds drastically change the evolution of the stellar mass, radius, and, consequently, of the envelope compactness parameter defined as ξ ≡ M/R (total mass over radius, Fernández et al. 2018). Figure 12 displays the evolution of these three quantities as a function of time at Z = 0.0025 and Ω/Ω_c = 0.6. We can see that after the activation of optically thick winds both the stellar mass and radius drop drastically compared to the optically thin wind case. At the end of the evolution, the stellar mass (radius) can be a factor of ∼2 (∼10³) smaller in the S23 case compared to the V01 case. This huge variation of the radius is reflected in the compactness parameter, which is extremely high (∼100) in the S23 case. In other words, the WR stars we obtain in the S23 scenario are much more compact (i.e., have a larger value of ξ = M/R, Fernández et al. 2018) compared to the giant stars we obtain as a result of the optically thin wind model.

Fig. 12. Evolution of stellar mass (top), radius (middle), and envelope compactness (bottom) as a function of the stellar age for Z = 0.0025 and Ω/Ωc = 0.6. Black dotted lines: V01 case. Red lines: S23 case. Dotted lines are for optically thin winds, solid lines for optically thick winds, and dashed lines for WR-type winds.

If the star directly collapses to a black hole at the end of its life, the large final compactness of our WRs implies that the mass of the black hole will be very similar to the total mass of the star, because neutrino emission during a failed supernova cannot unbind such a compact envelope (Fernández et al. 2018; Costa et al. 2022). In contrast, the low final compactness of the models with optically thin winds implies that a (large) fraction of the envelope might get unbound even in the case of a failed supernova explosion, because of the shock triggered by neutrino emission and propagating through the loosely bound envelope.

The main uncertainty is whether (and which of) our WRs collapse directly into black holes or undergo a successful core-collapse supernova. Here, we estimate the explodability of our stellar models with the delayed explosion formalism by Fryer et al. (2012). According to this simplified model, the final mass of the compact object only depends on the final CO-core mass and total mass of the progenitor star (Figure 13). We will perform a more detailed explodability study in a follow-up paper.

Fig. 13. Relation between MZAMS and final mass of the star (Mfin, thin solid line), He core mass (MHe, thick dotted line), CO core mass (MCO, thick dot-dashed line), and black hole mass (MBH, thick solid line), for optically thin winds (V01, magenta) and for the S23 model (S23, black). Upper panel: Ω/Ωc = 0; lower panel: Ω/Ωc = 0.6. Gray shaded area: range of MHe in which we expect a pulsational pair instability supernova.

Figure 13 clearly shows a difference between fast rotating and non-rotating models. For non-rotating models (upper panels), the final stellar mass (M_fin), final He core mass (M_He), and final CO core mass (M_CO) of model S23 are identical to those of model V01 until M_ZAMS < 60 M_⊙. In fact, for the metallicity shown in the figure (Z = 0.0025), non-rotating stars enter the optically thick wind regime only if M_ZAMS ≥ 60 M_⊙. As a consequence, the black hole mass M_BH is identical in model S23 compared to V01 up to M_ZAMS = 60 M_⊙. The differences in the black hole mass are minimal even for larger ZAMS masses, because pulsational pair instability kicks in at M_He ≥ 32 M_⊙ for the V01 model, removing the residual envelope mass for both model V01 and S23. The S23 model never enters pulsational pair instability (Simonato et al. 2025; Shepherd et al. 2025; Torniamenti et al. 2025), but it has already lost all of the H-rich envelope mass.

In the fast rotating case (Ω = 0.6 Ω_c, lower panel of Figure 13), the differences in the final masses (M_fin, M_He, M_CO) between models V01 and S23 are large already for M_ZAMS = 25 M_⊙. In fact, if the star rotates fast, the optically thick winds activate already at M_ZAMS = 25 M_⊙, peeling-off the H-rich envelope. Indeed, the lower panel of Figure 13 shows that the final total mass (M_fin) of a fast rotating star in the S23 model is almost always the same as its final He core mass (M_He). In the rotating case, both models V01 and S23 tend to develop more massive He and CO cores because of rotational mixing. In model S23, the growth of the He core is limited by core erosion from stellar winds. This explains why the black hole mass of model V01 is larger than that of model S23 for the entire considered mass range. Pulsational pair instability kicks in at M_ZAMS ≥ 40 M_⊙ and M_ZAMS ≥ 70 M_⊙ for models V01 and S23, respectively³.

The mass and radius evolution of our new models potentially has a large impact on the formation of binary compact objects. On the one hand, stronger winds lead to less bound binaries. This is due to (1) the rapid drop in stellar mass, which increases the binary separation, and (2) the radius evolution, which prevents the star from filling the Roche lobe and reduces the probability of a common envelope. On the other hand, a milder stellar expansion may prevent the merger of close stellar binaries at the end of the main sequence phase, increasing the number of close objects that survive to form a tight compact object binary. The balance of these two effects will be studied in more detail in a future work.

4.4. Optically thick winds at lower metallicity

To explore whether single WRs could be produced by self-stripping and the activation of the thick wind regime at metallicities even lower than that of the SMC, we simulate the evolution of a 20 M_⊙ star with initial rotation velocity Ω/Ω_c = 0.6 at Z = 0.001,  0.0006,  0.0003, corresponding roughly to Z_⊙/15, Z_⊙/25, and Z_⊙/50. The latter corresponds to the metallicity of IZw18, for which WR signatures have been observed in the past (Izotov et al. 1997) and chemically homogeneous evolution caused by strong initial rotation has been suggested as a possible channel to create these (e.g., Szécsi et al. 2015, 2025; Kubátová et al. 2019). Figure 14 shows that even at these low metallicities, rotating stars are indeed able to activate the thick wind regime according to the S23 prescription and thus would be able strip their envelope. This happens because stars at lower metallicities evolve on more luminous and compact tracks, facilitating the activation of thick winds through channel 1. Notably, the different wind prescription leads to a distinctively different shape of our tracks compared to early studies simulating chemically homogeneous evolution caused by high initial rotation (e.g. Maeder 1987; Yoon & Langer 2005; Szécsi et al. 2015). The characteristic hook-like feature when leaving the main sequence is not observed in older studies as they tend to assume a generally stronger wind mass-loss, while with the new treatment from S23, only the activation of optically thick winds keeps the stars eventually compact. Nonetheless, assuming the high amount of initial rotation assumed in our models is realized in nature, our models confirm a potential route for generating single WRs in extremely metal-poor galaxies such as IZw18 (Legrand et al. 1997; Aloisi et al. 1999; Schaerer et al. 1999; Shirazi & Brinchmann 2012; Kehrig et al. 2013; Szécsi et al. 2015, 2025; Kubátová et al. 2019).

Fig. 14. Tracks of a 20 M⊙ star with initial rotation Ω/Ωc = 0.6 and metallicity Z = 0.001, 0.0006, and 0.0003. The line styles and color code are the same as in Figure 1. We plot the tracks only until the end of C burning; afterward, in the very last stages of their evolution, the Teff and luminosity oscillate widely in our MESA models, and we did not plot them in order to have a cleaner image (see, e.g., Costa et al. 2022). Even at these low metallicities, stars are able to activate thick winds and peel off their envelopes.

5. Summary

Motivated by the recent confirmation that there is a group of companionless WR stars in the SMC (Schootemeijer et al. 2024), we have created a novel set of stellar evolution models to provide a scenario that can explain the self-stripping of these and other stars at low metallicity. This scenario requires the stars to reach a regime of optically thick winds, which we describe with the predictions of S23, through high initial rotation. Our main results are the following:

• We have shown that optically thick winds (S23) combined with high initial rotation, Ω/Ω_c ∼ 0.6, are able to strip single metal-poor O-type stars and transform them into WNh stars even at the metallicity of the SMC.

• The dependence on metallicity and rotation is not trivial. A metallicity Z ≤ 0.002 and high rotation Ω ≳ 0.6 Ω_c enable chemically homogeneous evolution, as the star is more luminous and compact than a slowly rotating star. This facilitates the activation of thick winds by the end of the core H-burning phase (channel 1). In contrast, a higher metallicity (Z ≥ 0.004) and lower rotation velocity lead to lower luminosities and effective temperatures and may lead to the activation of thick winds through the bi-stability jump (channel 2). Intermediate metallicity (Z ∼ 0.003) and rotation, instead, can lead to the formation of cool supergiant stars.

• We find a reasonable agreement between our models and the observed parameters of the hot SMC WNh stars with respect to the time spent in the different regions of the HR diagram, their surface abundance, and rotation velocity. We therefore consider our scenario a plausible path for the formation of single WRs at a lower metallicity.

• Within our model setup for the stripping of the stars, overshooting has a similar effect as rotation. For higher values of the overshooting parameter, f_ov, lower initial rotations are required to self-strip the envelope.

• The activation of optically thick winds combined with rotation also can help in alleviating the problem of overluminous giants inside the HD limit.

• Within our framework, the final stellar masses, core masses, and radii are smaller if the optically thick winds are enabled. This can impact the supernova explosion scenario, the final black hole masses, and the formation of binary compact objects.

• Fast rotating massive O-type stars develop thick winds even at a lower metallicity than the SMC (down to Z = 0.0003 ∼ Z_⊙/50) and can self-strip their envelope at least until carbon burning. This provides a possible explanation for the presence of WRs in extremely metal poor galaxies as IZw18.

Acknowledgments

We thank the anonymous referee for helpful and constructive comments. We thank Filippo Simonato, Marco Dall’Amico, and Sandro Bressan for useful discussions. MM acknowledges financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. LB, MM and ST acknowledge financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) STRUCTURES. AACS and VR acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander). This project was co-funded by the European Union (Project 101183150 – OCEANS). ST acknowledges financial support from the Alexander von Humboldt Foundation for the Humboldt Research Fellowship. EK and MM acknowledge support from the PRIN grant METE under contract No. 2020KB33TP. EK acknowledges financial support from the Fondazioni Gini for the Gini grant. GNS ad JSV are supported by STFC funding under grant number ST/Y001338/1. We acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grants INST 35/1597-1 FUGG and INST 35/1503-1 FUGG. We made use of the MESA software (https://docs.mesastar.org/en/latest/) version r12115; (Paxton et al. 2011, 2013, 2015, 2018, 2019). We made use of the MESA inlists https://github.com/Apophis-1/VMS_Paper1, and https://github.com/Apophis-1/VMS_Paper2 from Sabhahit et al. (2023). This research made use of NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), and MATPLOTLIB (Hunter 2007).

References

Appendix A: Evolution of Γ_e and Γ_{e, switch}

Here, we show the evolution of Γ_e and Γ_e,switch, as well as of the wind efficiency parameter η and η_switch to better understand the activation of thick winds through different channels. Figures A.1, A.2 and A.3 are for a star of 30 M_⊙ at Z = 0.002,  0.003,  0.004 respectively.

At Z = 0.002 (Figure A.1), at the end of the hydrogen burning phase, the star contracts, T_eff increases, luminosity and hence Γ_e increase, while the threshold for thick winds activation Γ_{e, switch} decreases. The evolution of Γ_{e, switch} can be understood from the bottom panel which compares η and η_switch as a function of Γ_e at different times. Γ_{e, switch} corresponds to the Γ_e where η = η_switch. We see that, during the contraction phase at the end of the main sequence, the increase of the wind efficiency parameter η, due to the growth of v_∞, implies a decrease of Γ_{e, switch}. Optically thick winds are activated during this contraction phase, when Γ_e > Γ_{e, switch} (middle panel), through channel 1.

At Z = 0.003 (Figure A.2) the behavior of Γ_{e, switch} is similar, but Γ_e is lower, due to lower luminosity, and never oversteps Γ_{e, switch} (middle panel). The contraction is then followed by an expansion, where Γ_{e, switch} rapidly increases due to the drop of v_∞ and of the rotation boost factor. The activation of helium burning stops the expansion and prevents the star from crossing the bi-stability jump at log(T_eff/K)≃4.4. In this case optically thick winds are not activated and the star becomes a cool supergiant.

At Z = 0.004 (Figure A.3) luminosity and Γ_e are even lower and Γ_e < Γ_{e, switch} during the contraction phase. However, after only ∼ 10⁴ yr, the star crosses the bi-stability jump, Γ_{e, switch} drops and optically thick winds are activated. As can be seen from the bottom panel, the effects of the bi-stability jump are: (i) increase of η due to an increase in Ṁ, (ii) decrease of η_switch due to the increase of the ratio $v_{\text{esc}}^2/v_{\infty }^2$ (see Eq. (2)). The combination of these two effects leads to a drop in Γ_{e, switch} from ∼0.7 to ∼0.5. Optically thick winds are activated through channel 2 (middle panel).

Fig. A.1. Evolution of log(Teff/K) (top panel), Γe and Γe, switch (middle panel, solid and dashed lines) as a function of stellar age for a star with M = 30 M⊙, Ω/Ωc = 0.6, at Z = 0.002. The color of the solid line indicates the burning stage: blue for hydrogen burning, orange for helium burning. The bottom panel shows η and ηswitch as a function of Γe at the times indicated by the vertical dotted lines in the top and middle panels. The lines’ color goes from blue to red increasing the time. The point where η = ηswitch is marked and represents the threshold for thick winds activation Γe, switch.

Fig. A.2. Same as Figure A.1 but for Z = 0.003.

Fig. A.3. Same as Figure A.1 but for Z = 0.004.

Appendix B: Outcomes for different metallicities and rotation

Fig. B.1. Tables showing the outcome of the simulation for a star with 20 M⊙ (top), 25 M⊙ (middle), 30 M⊙ (bottom), at different metallicities (columns) and initial rotations (rows). The color code represents the outcome: blue for WR formation through channel 1, yellow for WR formation through channel 2, red for cool supergiant formation.

Figure B.1 displays the outcome of the simulations for different metallicities and rotations for stars with 20, 25, 30 M_⊙. We see that increasing the initial stellar mass the WR outcome is more and more probable through both channel 1 and 2. For M = 20 M_⊙ the minimum rotation required to form a WR is Ω/Ω_c ∼ 0.6, for M = 25 M_⊙ this reduces to ∼0.5 and for M = 30 M_⊙ to ∼0.45. The dependence on metallicty and rotation is as we described in Section 3.1.1. At fixed rotation lower metallicity favors channel 1 activation, high metallicity favors channel 2, and intermediate metallicity leads to a cool supergiant. The same trend can be seen at fixed metallicity and decreasing rotation.

Appendix C: Results for Z = 0.004, channel 2

In the main text, we have shown most of the results for Z = 0.002, where channel 1 activation is favored. In this Appendix, we show the same results for Z = 0.004, where optically thick winds are activated through channel 2, thanks to the bi-stability jump. Figure C.1 shows the time spent on the HR diagram for stars with 20, 30, 40, 60 M_⊙. As described in Section 3.1.1, after hydrogen burning, these stars expand until reaching log(T_eff/K)∼4.4, where they enter the optically thick wind regime thanks to the bi-stability jump. Most of their He-burning phase is spent at lower temperatures log(T_eff/K) < 4.6, and they reach the hot WR stage during core carbon burning.

Fig. C.1. Same as Figure 5 but for Z = 0.004. Stars in the hot WR phase are generally more evolved with respect to Z = 0.002 and already in the carbon burning phase.

Figure C.2 shows element surface abundances. All the stellar tracks are still compatible with our definition of WNh stars when they pass through the data points. Their are less hydrogen depleted then the Z = 0.002 case, with surface hydrogen in the range 20 − 30%. Nitrogen is also slightly more abundant, with equilibrium value ∼2 × 10⁻³, in great agreement with data at log(T_eff/K) < 5. However, even in this case, further nitrogen surface enrichment occurs at higher temperatures with respect to observations.

Fig. C.2. Same as Figure 6 but for Z = 0.004.

Stellar rotation is shown in Figure C.3. The general behavior is similar to the Z = 0.002 case, but the drop during the main sequence is more pronounced. While stars expand to lower temperatures, their rotation drops further, and when they transition to the hot WRs phase they are practically non-spinning. Rotation starts to rise again at log(T_eff/K) > 4.9 but in a shallower way with respect to the Z = 0.002 case. Our models predicts v_rot ∼ 50 km/s at log(T_eff/K)∼5, still consistent with the upper limits mentioned in the text.

Fig. C.3. Same as Figure C.3 but for Z = 0.004.

The mass-loss and transformed mass-loss are displayed in Figures C.4 and C.5. Also in this case, the predicted transformed mass-loss is larger than observational data.

Fig. C.4. Same as Figure 8 but for Z = 0.004.

Fig. C.5. Same as Figure 9 but for Z = 0.004.

Appendix D: Stellar interiors: Abundance profiles

Figure D.1 shows the evolution of the abundance profiles of our simulated stars. We take as an example a star with initial mass M = 30 M_⊙, representative of those passing through observational data, for Z = 0.002 (left) and Z = 0.004 (right), during main sequence (top) and during more advanced phases (bottom). For the main sequence evolution we plot the initial configuration at ZAMS, the final configuration, when central hydrogen drops below 10⁻² (8.8 Myr for Z = 0.002, 8.49 Myr for Z = 0.004), and the configuration at half main sequence time.

The evolution during the first half of the main sequence is almost purely homogeneous, with helium and nitrogen increasing and hydrogen and carbon decreasing throughout the whole the star. Nitrogen abundance has already reached equilibrium at 10⁻³ and 2 × 10⁻³ for the two cases. Afterward, pure homogeneity is broken, due to rotational velocity reduction, and some hydrogen remains on the surface ∼20% for Z = 0.002, and ∼30% for Z = 0.004. The helium core, defined as the region where H abundance < 10⁻² and He abundance > 10⁻², is slightly larger in the low metallicity case, ∼ 18 M_⊙ against ∼ 15 M_⊙.

As for the post main sequence phase, we show 3 times, corresponding to:

1. the time when the effective temperature is at its minimum, i.e., the star is in the coldest region of the HR diagram, soon after the end of hydrogen burning, corresponding to log(T_eff/K)∼4.65 (age 8.93 Myr) and log(T_eff/K)∼4.38 (age 8.58 Myr);

2. the time when log(T_eff/K)∼5, i.e., most of the observed SMC WRs are found, corresponding to 9.13 Myr and 8.91 Myr;

3. the time when log(T_eff/K)∼5.15, corresponding to very late evolutionary stages (9.32 Myr and 8.96 Myr).

First, we notice that stellar mass is fastly reduced by winds, with both stars loosing ∼ 7 M_⊙ in ∼ 0.4 Myr. Second, we see that hydrogen is still present on the surface at log(T_eff/K)∼5, making the star appear as a WNh, but it is rapidly eroded by winds. Third, we notice that core helium burning quickly depletes nitrogen in the core. However, in the outer layers, the CNO equilibrium value of nitrogen remains constant even at log(T_eff/K)∼5, in tension with observational data that show nitrogen enhancement at these temperatures.

Both cases present a small nitrogen bump just beneath the surface, but probably rotation is too low to bring it on the surface. In our models, nitrogen enhancement occurs later, at log(T_eff/K)∼5.15, when winds reveal the inner nitrogen bump.

Fig. D.1. Stellar profiles for a 30 M⊙ star. Abundances of hydrogen (blue), helium (orange), 12C (red), and 14N (black) as a function of the mass enclosed in a given shell. Left panels are for Z = 0.002, right for Z = 0.004. Top panels: main sequence evolution. Line styles are for different times: dotted for ZAMS, dashed for half main sequence, solid for end of the main sequence. Bottom panels: post main sequence phases. Line styles are for different times: dotted for the time corresponding to minimum temperature, dashed for log(Teff/K)∼5, solid for log(Teff/K)∼5.15.

Appendix E: Discussion on the bi-stability jump

The existence and strength of the mass-loss increase found by Vink et al. (1999) has been an object of investigation over the last 25 years, with very different results based on the adopted stellar atmosphere model. While studies using a version of the FASTWIND code aiming for a locally consistent solution of the wind dynamics find a steady decline of mass-loss rate with T_eff with no bi-stability jump at T_eff ∼ 25 kK (Björklund et al. 2021, 2023), Petrov et al. (2016), using a globally consistent approach with the CMFGEN code, confirm the existence of a jump, albeit at a somewhat lower T_eff ∼ 20 kK. Using a different code (METUJE), Krtička et al. (2021, 2024) do find an increase of the mass-loss at cooler T_eff ∼ 19 kK, but less pronounced than Vink et al. (1999).

Observationally, changing wind parameters and specific line profiles in individual luminous blue variable stars provides evidence for a mass-loss rate change (Groh et al. 2011; Groh & Vink 2011). This could also be found in recent hydrodynamical evolution models Grassitelli et al. (2021) implementing a specific formalism of the bi-stability jump mass-loss description Smith et al. (2004). In contrast, empirical evidence for a mass-loss rate increase in normal OB supergiants is much less clear. Already Vink et al. (2000) noted a discrepancy between their theoretical predictions and some of the measurements by Kudritzki et al. (1999) in the 12500 < T_eff/K < 22500 temperature regime. Also Trundle et al. (2004) confirmed the discrepancy between theoretical mass-loss rates and empirical modeling. Recent investigations of blue supergiants in the Galaxy (Bernini-Peron et al. 2023; de Burgos et al. 2024), LMC (Verhamme et al. 2024; Alkousa et al. 2025), and SMC (Bernini-Peron et al. 2024) could not find a clear mass-loss enhancement when transitioning to cooler T_eff, but in particular the detailed analysis efforts of individual targets (Bernini-Peron et al. 2023, 2024; Alkousa et al. 2025) show that there is no clear monotonic downward trend in the mass-loss rates predicted by Björklund et al. (2023). Among the available descriptions the results align reasonably well with the predictions from Krtička et al. (2021, 2024) including a shallow mass-loss increase, albeit with some scatter depending on the sample.

However, in the context of the S23-formalism, it is not the direct Ṁ-jump in the Vink et al. (2001) description, but actually the reduction in the η_switch-value that is the major driving factor for the activation of the thick wind regime. This reduction is caused by a drop of the v_∞/v_esc ratio. Despite some uncertainty in the stellar masses, which can affect the determination of v_esc, a decline of v_∞/v_esc is well observed, also across different metallicities (e.g., Lamers et al. 1995; Crowther et al. 2006; Markova & Puls 2008; Bernini-Peron et al. 2024; Alkousa et al. 2025). Depending on the sample it is currently unclear whether there is a sharp drop in this ratio or a more gradual decline, but in each case the S23-condition would be met. While this does not imply that channel 2 must occur in nature, its activation criteria is not in conflict with observations and thus justified within our framework.

All Figures

Fig. 1. Stellar tracks on the HR diagram for initial masses M = 25,  30,  40,  60, and 80 M⊙ at metallicity Z = 0.0025. In the left-hand panels, the V01 model with only optically thin winds is implemented. In the right-hand panels, the S23 model is enforced, with the possibility to activate optically thick winds. The upper panels show the case with no rotation, while the lower panels the Ω/Ωc = 0.6 case. Different line styles represent different wind regimes: dotted for optically thin winds, solid for optically thick winds, dashed for WR-type winds. Colors represent different burning stages, defined as the element whose burning generates most of the energy of the star. Blue is for H-burning through CNO, orange for He-burning through triple α, red for carbon burning, and green for oxygen burning. Green circles are observations of single WRs (Hainich et al. 2015), filled for the hottest, open for the coldest. The lower-right plot shows that the combination of optically thick winds and high initial rotation is able to self-strip stars even at SMC metallicity. Stars in the lower-right plot (S23, Ω/Ωc = 0.6) avoid the HD limit and transition from main-sequence stars to WRs.

In the text

	Fig. 2. Velocity distribution of O-type stars in the SMC. The plot shows the normalized observed v sin i distribution (blue) and the reconstructed initial rotational velocity (vrot, red) for 171 O-type stars. The dataset is a compilation from Ramachandran et al. (2019), Rickard et al. (2024), and Bestenlehner et al. (2025). About 10% of the O-type stars feature vrot > 450 km/s.
In the text

Fig. 3. Illustration depicting the possible evolutionary paths of a chemically homogeneous star with initial mass ≳ 20 M⊙ at SMC metallicity. After core H-burning (yellow ball with blue center), the star contracts, Γe,switch decreases and Γe increases. If Γe > Γe,switch the star enters the optically thick winds regime (left arrow) and ends up being a WR (channel 1). If Γe < Γe,switch (right arrow), winds are still optically thin during shell H-burning (yellow ball with blue shell) and the star expands and cools, increasing Γe,switch. At this stage, there are two possible outcomes: (i) If the star cools enough to cross the bi-stability jump Teff < 25 kK (right arrow), it may enter the optically thick wind regime and become a WR (channel 2). (ii) If the star starts core He-burning at Teff > 25 kK (left arrow), Teff rises again, and the star does not cross the bi-stability jump. Optically thick winds are not activated (or are activated too late) and the star ends up being a cool supergiant.

In the text

Fig. 4. Stellar tracks for a 25 M⊙ star. The upper panel shows a rotating star with Ω/Ωc = 0.6 for three different metallicities Z = 0.002,  0.003,  0.004. The lower panel shows a star with metallicity Z = 0.003 for three different values of the initial rotation speed Ω/Ωc = 0.5,  0.6,  0.7. The linestyles and color code are the same as in Figure 1. The vertical dotted black line represents the temperature below which we expect the bi-stability jump. This Figure shows the non-monotonic dependence of mass-loss on metallicity and rotation. Low metallicity (Z = 0.002) favors the activation of optically thick winds through channel 1, because the star is more luminous and compact. Higher metallicity (Z = 0.004) can also activate optically thick winds through channel 2. Stars in the intermediate metallicity case (Z = 0.003), instead, fail to activate optically thick winds soon enough and end up as cool supergiant stars. The same trend happens for rotation.

In the text

Fig. 5. Time spent on different regions of the HR diagram for stars with mass M = 20,  30,  40,  60, and 80 M⊙ at Z = 0.002 and Ω/Ωc = 0.6. In the upper panel, the color code represents the quantity τHRD, in yr dex−1, which is the time spent to move by 1 dex in the HR diagram. Green circles are WRs in the SMC ,and blue squares are O-type stars in the SMC with v sin i > 200 km/s from (Mokiem et al. 2006; Bouret et al. 2013, 2021; Ramachandran et al. 2019; Dufton et al. 2019; Rickard et al. 2024; Backs et al. 2024). In the lower panel, the color code represents different burning stages: core H-burning (blue), He-burning (orange), carbon burning (red), and oxygen burning (green), while the line styles represent different wind regimes: dotted for optically thin winds, solid for optically thick winds, dashed for WR-type winds. Cross markers are spaced by ∼ 5 × 104 yr each. After the main sequence, stars spend most of their time at the turn-off point, where they start to peel off their H-rich envelope and move blueward, and in the hot WR phase, at 4.9 < log(Teff/K) < 5.2.

In the text

Fig. 6. Upper panel: Classification of the star along the HR diagram: O-type (black), WNh (blue), WN (orange), and WC and WO (red). Lower panel: Surface abundances as a function of the effective temperature, Teff. Line styles are for different elements: hydrogen (dotted line), 14N (solid line), 12C (dashed line). Different colors are for different initial masses: 20 (blue), 30 (orange), 40 (red), and 50 M⊙ (green). Green triangles (circles) represent hydrogen (nitrogen) abundances for the five hottest WRs observed (Hainich et al. 2015). The gray horizontal band covers the ∼5 − 40% range; if the hydrogen abundance is inside this region, the star is considered a WNh.

In the text

Fig. 7. Upper panel: Average rotational velocity along the HR diagram (color code). Lower panel: Evolution of the average rotational velocity as a function of the Teff. The color code is the same as in Figure 6. Dashed lines represent the main sequence evolution, solid lines the advanced phases. The velocity drops during the main-sequence evolution. This is followed by the hook at the end of the main sequence, when the velocity quickly rises due to post main sequence contraction, and by a subsequent drop. In the late phases of the evolution log(Teff/K)≳5.05, the velocity rises again due to winds peeling-off the outer layers of the star and exposing the rotating core.

In the text

	Fig. 8. mass-loss, Ṁ, (upper panel) and log(Teff/K) (lower panel) as a function of normalized stellar age. The color code is the same as in Figure 6. Line styles represent different wind regimes: dotted for optically thin winds, solid for optically thick winds, dashed for WR-type winds.
In the text

	Fig. 9. Upper panel: Transformed mass-loss, Ṁt, along the HR diagram track (color code). Lower panel: Evolution of Ṁt as a function of Teff. The color code is the same as in Figure 6. The green points are computed from the observed Ṁ, L, and v∞ from Hainich et al. (2015).
In the text

Fig. 10. Stellar tracks for a star of initial mass M = 20 M⊙ at Z = 0.0025 for initial rotation speed Ω/Ωc = 0.55 (left), 0.5 (middle), and 0.45 (right) and overshooting parameter fov = 0.03, 0.045, 0.07, and 0.1 (right panel) during the main sequence. The effect of overshooting is almost degenerate with that of rotation. A slight increase of fov corresponds to a reduction of the rotation threshold needed for self-stripping.

In the text

Fig. 11. Stellar tracks for stars with mass M = 40,  50, and 60 M⊙ at Z = 0.002 (upper panels) and Z = 0.004 (lower panels), encompassing the metallicity range of the SMC. Columns represents different initial rotations: Ω = 0 (left), Ω/Ωc = 0.3 (middle), Ω/Ωc = 0.6 (right), corresponding to the minimum, average, and maximum observed rotation velocities of O-type stars in the SMC. The color code shows the time spent in different regions of the HR diagram τHRD.

In the text

	Fig. 12. Evolution of stellar mass (top), radius (middle), and envelope compactness (bottom) as a function of the stellar age for Z = 0.0025 and Ω/Ωc = 0.6. Black dotted lines: V01 case. Red lines: S23 case. Dotted lines are for optically thin winds, solid lines for optically thick winds, and dashed lines for WR-type winds.
In the text

Fig. 13. Relation between MZAMS and final mass of the star (Mfin, thin solid line), He core mass (MHe, thick dotted line), CO core mass (MCO, thick dot-dashed line), and black hole mass (MBH, thick solid line), for optically thin winds (V01, magenta) and for the S23 model (S23, black). Upper panel: Ω/Ωc = 0; lower panel: Ω/Ωc = 0.6. Gray shaded area: range of MHe in which we expect a pulsational pair instability supernova.

In the text

Fig. 14. Tracks of a 20 M⊙ star with initial rotation Ω/Ωc = 0.6 and metallicity Z = 0.001, 0.0006, and 0.0003. The line styles and color code are the same as in Figure 1. We plot the tracks only until the end of C burning; afterward, in the very last stages of their evolution, the Teff and luminosity oscillate widely in our MESA models, and we did not plot them in order to have a cleaner image (see, e.g., Costa et al. 2022). Even at these low metallicities, stars are able to activate thick winds and peel off their envelopes.

In the text

Fig. A.1. Evolution of log(Teff/K) (top panel), Γe and Γe, switch (middle panel, solid and dashed lines) as a function of stellar age for a star with M = 30 M⊙, Ω/Ωc = 0.6, at Z = 0.002. The color of the solid line indicates the burning stage: blue for hydrogen burning, orange for helium burning. The bottom panel shows η and ηswitch as a function of Γe at the times indicated by the vertical dotted lines in the top and middle panels. The lines’ color goes from blue to red increasing the time. The point where η = ηswitch is marked and represents the threshold for thick winds activation Γe, switch.

In the text

	Fig. A.2. Same as Figure A.1 but for Z = 0.003.
In the text

	Fig. A.3. Same as Figure A.1 but for Z = 0.004.
In the text

	Fig. B.1. Tables showing the outcome of the simulation for a star with 20 M⊙ (top), 25 M⊙ (middle), 30 M⊙ (bottom), at different metallicities (columns) and initial rotations (rows). The color code represents the outcome: blue for WR formation through channel 1, yellow for WR formation through channel 2, red for cool supergiant formation.
In the text

	Fig. C.1. Same as Figure 5 but for Z = 0.004. Stars in the hot WR phase are generally more evolved with respect to Z = 0.002 and already in the carbon burning phase.
In the text

	Fig. C.2. Same as Figure 6 but for Z = 0.004.
In the text

	Fig. C.3. Same as Figure C.3 but for Z = 0.004.
In the text

	Fig. C.4. Same as Figure 8 but for Z = 0.004.
In the text

	Fig. C.5. Same as Figure 9 but for Z = 0.004.
In the text

Fig. D.1. Stellar profiles for a 30 M⊙ star. Abundances of hydrogen (blue), helium (orange), 12C (red), and 14N (black) as a function of the mass enclosed in a given shell. Left panels are for Z = 0.002, right for Z = 0.004. Top panels: main sequence evolution. Line styles are for different times: dotted for ZAMS, dashed for half main sequence, solid for end of the main sequence. Bottom panels: post main sequence phases. Line styles are for different times: dotted for the time corresponding to minimum temperature, dashed for log(Teff/K)∼5, solid for log(Teff/K)∼5.15.

In the text