We study in detail the evolution of two massive stars at solar metallicity Z0.014 taken from Romagnolo et al. 2024 Paper I by running evolution models for initial masses 60 and 200 Msun using MESA and GENEC. For the mass loss we adopt the selfconsistent mCAK prescription for the optically thin winds of OBstars a semiempirical formula for Hrich thick wind of WNh stars and a hydrodynamically consistent formula for the Hpoor thick wind of classical WolfRayet stars. Both codes predict different tracks across the HRD. For the 60 Msun case Genec models predict a more efficient rotational mixing and more chemically homogeneous evolution whereas MESA model predicts a large radial expansion postMS reaching the LBV phase. For the 200 Msun case differences are less relevant because their evolution is dominated by wind mass loss with a weaker dependence on internal mixing and only the treatment for superadiabacity creates an impact during the Heburning stage. The switch of the mass loss based on the proximity to the Eddington factor instead of the removal of outer layers implies the existence of WNh stars with a large mass fraction of hydrogen at the surface formed from initial masses gtrsim60 Msun. These stars are constrained in a Teff range of the HR diagram which corresponds to the MS band in agreement with the observations of Galactic WNh stars. While our models employ a fixed Gammae threshold for the switch to thick winds rather than a continuous thintothick wind model the good reproduction of observations supports the robustness of the wind model upgrades introduced in Paper I allowing its application to studies of latestage stellar evolution before core collapse.
['Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
Evolution of stars with 60 and 200 Msun predictions for WNh stars in the Milky Way
{'Evolution of stars with 60 and 200 M ⊙ : predictions for WNh stars in the Milky Way.': 'A. C. Gormaz-Matamala 1 , 2 , A. Romagnolo 1 , and K. Belczynski 1 ⋆  \n- 1 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warsaw, Poland 2 Astronomický ústav, Akademie vˇed Ceské republiky, Friˇcova 298, 251 65 Ondˇrejov, Czech Republic e-mail: alex.gormaz@asu.cas.cz', 'ABSTRACT': 'Context. Massive stars are characterised by powerful stellar winds driven by radiation, and therefore mass-loss rate plays a crucial role in their evolution.  \nAims. We study in detail the evolution of two massive stars (a classical massive star and a very massive star) at solar metallicity ( Z = 0 . 014) taken from Romagnolo et al. (2024a, Paper I); by calculating their final masses, radial expansion, and chemical enrichment, at their H-core, He-core and C-core burning stages prior to the final collapse.  \nMethods. We run evolutionary models for initial masses 60 and 200 M ⊙ , using MESA and the Geneva-evolution-code (G enec ). For the mass loss, we adopt the self-consistent m-CAK prescription for the optically thin winds of OB-type stars, a semi-empirical formula for H-rich optically thick wind of WNh stars, and a hydrodynamically consistent formula for the H-poor thick wind of classical Wolf-Rayet stars.  \nResults. As seen from Paper I, the unification of the initial setup for stellar structure and wind prescription leads into very similar M bh for both G enec and MESA codes. However, both codes predict di ff erent tracks across the HRD. For the 60 M ⊙ case, G enec models predict a more e ffi cient rotational mixing and more chemically homogeneous evolution, whereas MESA model predicts a large radial expansion post-MS reaching the LBV phase. For the 200 M ⊙ case, di ff erences between both evolution codes are less relevant because their evolution is dominated by wind mass loss with a weaker dependence on internal mixing.  \nConclusions. The switch of the mass-loss prescription based on the Eddington factor instead of the removal of outer layers, implies the existence of WNh stars with a large mass fraction of hydrogen at the surface ( X surf ≥ 0 . 3) formed from initial masses of ≳ 60 M ⊙ . These stars are constrained in a T e ff range of the HRD which corresponds to the MS band, in agreement with the observations of Galactic WNh stars at Z = 0 . 014. While our models employ a fixed Γ e,trans threshold for the switch to thick winds, rather than a continuous thin-to-thick wind model, the good reproduction of observations supports the robustness of the wind model upgrades introduced in Paper I, allowing its application to studies of late-stage stellar evolution before core collapse.  \nKey words. Stars: early-type - Stars: evolution - Stars: mass-loss - Stars: winds, outflows - Stars: Wolf-Rayet', '1. Introduction': 'The detection of the first source of gravitational waves (GW, Abbott et al. 2016, 2021) opened the door for a large number of studies dedicated to find more GW sources in the Universe, namely double compact objects (hereafter DCO) mergers. Rates of these DCO merger events depend on considered conditions for the binary system concerning the orbital separation, the formation (or not) of a common envelope (CE), as well as the initial assumptions for RLOF and their timescales (Ivanova et al. 2013; Olejak et al. 2021). Moreover, these modes are intimately linked with the single evolution of the stars in a binary system, which is crucial to determine the final remnant masses (Yungelson et al. 2008; Belczynski et al. 2010; Yusof et al. 2013; Spera & Mapelli 2017; Bavera et al. 2023; Martinet et al. 2023; Romagnolo et al. 2024a) and the maximum radial expansion that the stellar components can reach during their evolution (Agrawal et al. 2020; Romagnolo et al. 2023), reason why accurate evolution models are necessary to constrain the theoretical rate of the DCO merging events.  \nStellar black holes are the final fate of massive stars with M zams ≳ 25 M ⊙ (Heger et al. 2003). Even though the conditions  \nof the death as supernovae also constrain the mass of the final remnant (Fryer et al. 2022; Olejak et al. 2022), the stellar mass prior to the supernova explosion is explicitly dependent of the mass-loss rate (denoted as d M / d t or simply ˙ M ) carried out by the strong winds during the evolution of a massive stars (Vink 2021). Because stellar winds of hot massive stars are driven by line-force (Lucy & Solomon 1970; Castor et al. 1975; Abbott 1982), these ˙ M are directly dependent on metallicity. This explains the predominance of more massive black holes in less metallic environments (Farrell et al. 2021; Vink et al. 2021), although exceptions due to unusually weaker winds can also be theoretically described (Belczynski et al. 2020a). Moreover, the strength of the stellar winds varies along the evolution of the star, which is dominated during the main sequence (H-core burning stage) by line-driven optically thin winds for OB-type stars (Puls et al. 2008), then eventual strong outflows might stem during the LBV stage (Humphreys et al. 2016), and finally we find optically thick winds of Wolf-Rayer (WR) stars for the more evolved stages (Crowther 2007; Gräfener & Hamann 2008). Optically thick winds are naturally denser, and therefore the mass loss during the WR phase is higher in comparison with OB-type phase.  \nRecent theoretical (Krtiˇcka & Kubát 2017, 2018; GormazMatamala et al. 2019, 2022a; Björklund et al. 2021) and empir-  \nl studies (Hawcroft et al. 2021, 2024) have found theoretical values for the ˙ M of massive stars to be in the order of ∼ 3 times lower than the standard mass-loss recipe from Vink et al. (2001, hereafter V01). Such di ff erence in the value of mass loss is enough to alter the evolution of stars with M zams ≳ 20 M ⊙ (Meynet et al. 1994). Adopting new winds, massive stars will retain more mass during their main sequence phase, thus being bigger and more luminous (Gormaz-Matamala et al. 2022b; Björklund et al. 2023). Besides, massive stars keep larger rotational velocities for a more extended time, because they lose less angular momentum (Gormaz-Matamala et al. 2023). Even though these changes in the stellar evolution have been applied mostly for the early H-burning stages, they also impact the further evolution of the subsequent stages (Gormaz-Matamala et al. 2024; Josiek et al. 2024). Moreover, considering that also the mass-loss prescription for WR winds has been revisited (Bestenlehner 2020; Sander & Vink 2020) as well as its impact on the stellar evolution (Higgins et al. 2021), the entire panorama of the stellar evolution is being upgraded.  \nIn Romagnolo et al. (2024a, hereafter Paper I), we calculate the remnant mass for a set of initial masses from 10 to 300 M ⊙ for Z = 0 . 014, based on stellar evolution models adopting the state-of-the-art physical ingredients such as Tayler-Spruit dynamo, larger core-convective overshooting, and new mass-loss recipes. As a result, we found a maximum black hole mass of M bh,max = 28 . 3 M ⊙ formed from stars with M zams ≳ 250 M ⊙ , larger than the previous M bh,max ≃ 20 M ⊙ from Belczynski et al. (2010) but quite below the maximum BH masses predicted by Bavera et al. (2023, ≃ 35 M ⊙ at M zams = 75 M ⊙ ) and Martinet et al. (2023, ≃ 42 M ⊙ at M zams = 180 M ⊙ ). Even though the recently discovered largest BH in the Milky Way has a mass of M bh,max = 32 . 7 ± 0 . 70 M ⊙ (Gaia Collaboration et al. 2024), the low metallicity of its star companion suggest that this black hole was formed from a metal poor massive star and thus our predicted M bh,max remains valid.  \nMoreover, the new wind prescription introduced in Paper I entail more implications concerning the evolution of massive stars. Standard evolution models (Ekström et al. 2012; Choi et al. 2016) adopted the switch from OB-type to WR-type winds once the stellar surface removes the majority of its hydrogen. Nevertheless, theoretical studies have found that the transition from optically thin to thick winds during the evolution of massive stars is more due to the proximity to the Eddington limit rather than due to the removal of the outer envelopes (Vink et al. 2011; Gräfener et al. 2011), thus altering the expected lifetimes and total amount of mass removed at each evolutionary stage. Because of this, in Paper I we set the switch between optically thin to optically thick winds according to the proximity to the Eddington limit instead of based on the chemical composition of the stellar outer envelope. As a consequence, the most massive stars develop thick winds already during their earliest evolution stages, altering their tracks over the HR diagram and predicting the existence of WR stars with hydrogen abundance X surf ≥ 0 . 3 in their surfaces, the so-called WNh stars (Martins et al. 2008; Crowther et al. 2010).  \nIn this work, we study in detail the impact of the already mentioned upgrades on the evolution of two stars from Paper I, with special emphasis on the existence of H-rich WR stars. For that purpose, we run evolution models for the stars with initial masses of 60 and 200 M ⊙ at solar metallicity ( Z = 0 . 014), following Asplund et al. (2009). We introduce our mass-loss prescription in Sec. 2.2, together with the description of the switch from optically thin to thick winds at Γ e = 0 . 5. The resulting evolution tracks over the HR diagram for our models are detailed  \nin Sec. 3, whereas the observational diagnostics comparing with spectroscopic analyses are shown in Sec. 4. Finally, conclusions are summarised in Sec. 5.', '2.1. Stellar evolution codes': "Same as in Paper I, we perform our stellar evolution models using the codes MESA (Modules for Experiments in Stellar Astrophysics, Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) and the Geneva-evolution-code (Eggenberger et al. 2008, hereafter G enec ). For step-overshooting, we adopt a value for α ov of 0 . 5, in agreement with the most recent studies on convective cores (Martinet et al. 2021; Scott et al. 2021; Bara ff e et al. 2023). Concerning angular momentum transport, we include for both MESA and G enec codes the Tayler-Spruit dynamo (Tayler 1973; Spruit 2002; Heger et al. 2005), whereas we replace the treatment for the convective layers from Schwarzchild to Ledoux's criterion (Ledoux 1947). For the convective outer layers, MESA adopts a mixing-length parameter of α MLT = 1 . 82, whereas G enec keeps α MLT = 1 . 6. In order to reduce superadibaticity in regions near the Eddington limit we use the use\\_superad\\_reduction method, which is stated to be a more constrained and calibrated prescription than MLT ++ in MESA (Jermyn et al. 2023). In particular, this di ff erence in the superadiabacity treatment will be relevant during the VMS evolution across the He-core burning stage (Romagnolo et al. 2024b), as we observe in Sec. 3.2.", '2.2. Mass loss for massive and very massive stars': "Evolutionary tracks of massive stars adopting lower ˙ M are in rule with the observed mass loss of clumped winds (Bouret et al. 2005; Muijres et al. 2011; Šurlan et al. 2013), as well as with more detailed calculations for the wind structure of Otype stars by self-consistently solving radiative transfer equations at full NLTE conditions (Sander et al. 2017; Sundqvist et al. 2019; Gormaz-Matamala et al. 2021). However, this scenario can hardly be extended to higher masses such as VMS, whose mass-loss rates values even larger than the given ones by V01 (Vink et al. 2011). Their stellar winds are so strong that they develop WR-type spectral profile with broad emission lines before the end of the H-core burning stage (WNh stars). The limit between VMS and 'regular' massive stars is normally set around M zams ∼ 100 M ⊙ , although this value may vary by means of rotation, metallicity, and specially mass-loss input (Yusof et al. 2013; Köhler et al. 2015).  \nEvolution of VMS, revisiting the transition between formal O-type to WNh stars, has been just recenty analysed by Gräfener (2021) and Sabhahit et al. (2023). Despite this, standard evolution codes still assume the transition of the wind regime based on the surface composition (when the hydrogen mass fraction becomes X surf ≤ 0 . 3 for G enec and X surf ≤ 0 . 4 for MESA ) the massloss recipe switches from the prescription for optically thin V01 to optically thick winds according to Gräfener & Hamann (2008, hereafter GH08 1 ) in G enec , whereas MESA uses the prescription of Nugis & Lamers (2002) for WR winds. However, spectroscopic diagnostics have found WNh spectral features for VMS with surface hydrogen abundances as high as X surf ≃ 0 . 7 in the  \nMilky Way (Martins 2023) and in the Large Magellanic Cloud (Tehrani et al. 2019; Martins & Palacios 2022), thus strengthening the idea that the transition to thick winds should be based on another criteria apart from surface chemical composition.", '2.2.1. Winds in the proximity to the Eddington limit': 'Mass loss of VMS is directly dependent on the Eddington factor, and the relationship quickly escalates when the Eddington limit is approached (Vink et al. 2011; Vink & Gräfener 2012). The Eddington limit defines the limit for the hydrostatic stability of a star, and it corresponds when the outwards radiative force equals to the inwards gravitational force 2 (Langer 1997). This limit is better illustrated by the classical Eddington factor, which considers only the electron scattering for the radiative processes  \nΓ e = L ∗ κ e 4 π cGM ∗ , (1)  \nwith κ e being the electron scattering opacity. Assuming a fully ionised plasma, the Eddington factor is usually expressed as a function of stellar mass, luminosity and hydrogen abundance  \nlog Γ e = -4 . 813 + log(1 + X surf) + log( L ∗ / L ⊙ ) -log( M ∗ / M ⊙ ) , (2)  \ni.e., the classical Eddington factor depends only on the fundamental stellar parameters.  \nNumerical models found a kink in the slope of the ˙ M ∝ Γ e dependence at Γ e ≃ 0 . 7, after which the mass-loss scales with the Eddington factor in agreement with the relationship found for Gräfener & Hamann (2008) for WNL stars. Indeed, Γ e = 0 . 7 was already used as a transition criterion to switch from O-type to WR stars in evolution models (Chen et al. 2015). Posteriori observational diagnostics on very massive stars found that wind clumping and uncertainties a ff ects the accuracy to determine the kink, which can be set as low as Γ e ≃ 0 . 3 (Bestenlehner et al. 2014), even though this value for Γ e strongly depends on the mass estimations derived from the L / M relations from Gräfener et al. (2011). A more theoretically detailed study is provided by Sabhahit et al. (2023), where the switch on the wind regime is directly correlated with the wind e ffi ciency η = ˙ M 3 ∞ / ( L / c ) at di ff erent metallicities. In this case however, the mass loss of the initial O-type phase is still described by the V01 instead of more updated recipes (e.g. Krtiˇcka & Kubát 2017; Björklund et al. 2021; Gormaz-Matamala et al. 2022b), whereas the η switch is calculated on the basis of β -law velocity profiles for the 3 ∞ term.  \nInstead, we aim for a transition point between thin and thick winds, independent not only of the initial mass-loss rate adopted but also independent on the adopted velocity profile. For this, we based our transition in the upper limit of validity of our mCAK prescription (Gormaz-Matamala et al. 2022a), where both ˙ M and 3 ( r ) are self-consistently calculated together with the lineforce parameters (Gormaz-Matamala et al. 2019). Bestenlehner (2020) extended the theoretical calculation of ˙ M according to the CAK theory (Castor et al. 1975) for regimes closer to the Eddington limit, finding also a transition Γ e,trans beyond the which the slope of the relation ˙ M ∝ Γ e become steeper, in line with the observations of WNh stars exhibiting thick winds from Bestenlehner et al. (2014). This Γ e,trans ranges from ∼ 0 . 46 to ∼ 0 . 48, depending on the value of the CAK line-force parameter 3 α . For  \nα ∼ 0 . 4 -0 . 6, which are the values found from the self-consistent line-force parameters found by Gormaz-Matamala et al. (2022b) for stars with masses in the range of 25 to 120 M ⊙ and metallicities from Z = 0 . 014 (solar) to = 0 . 002 (SMC), the value of the transition Eddington factor keeps relatively constant around Γ e,trans ≃ 0 . 48. Therefore, given that our new prescription for optically thin winds is based on the m-CAK theory, we establish Γ e = 0 . 5 as the transition point between thin and thick winds for our evolution models. The thickness in this case means the dominance of the electron scattering on the radiative acceleration, above the line-driven. We will discuss the physical implications of this transition in Section 4.1.', '2.3.1. Mass loss for the H-core burning stage': 'For thin winds ( Γ e ≤ 0 . 5), we use the mass-loss prescription from Gormaz-Matamala et al. (2022b, 2023)  \nlog  \n˙  \nM  \nthin  \n=  \n-  \n40  \n+  \n1  \n.  \n.  \n314  \n+  \n15  \n0564  \ny  \n-  \nwy  \n438  \nw  \n+  \n45  \n/  \n2  \n.  \n36  \n-  \n1  \n.  \n.  \n838  \nx  \n-  \n8  \n1967  \nx y  \n+ z × " 0 . 4 + 15 . 75 M ∗ / M ⊙ # , (3)  \nwhere ˙ M thin is in M ⊙ yr -1 ; and where w , x , y , and z are defined as  \nw = log GLYPH<18> T e ff kK GLYPH<19> , x = 1 log g , y = R ∗ R ⊙ , z = log Z ∗ Z ⊙ ! .  \nAs mentioned in Gormaz-Matamala et al. (2023), the dependence between mass loss and metallicity is not constant but it intrinsically depends on the mass (and the luminosity), because the radiative acceleration for more massive and luminous stars is mostly dominated by the continuum over the line-acceleration (see also Krtiˇcka & Kubát 2018).  \nThe range of validity of Eq. 3 is for T e ff ≥ 30 kK and log g ≥ 3 . 2. Beyond these values, we return to V01. For the cases when T e ff ≤ 10 kK, mass loss is described by de Jager et al. (1988).  \nFor thick winds ( Γ e > 0 . 5), since we establish the transition from the theoretical analysis of Bestenlehner (2020), we use their ˙ M ( Γ e) recipe with the calibration performed by Brands et al. (2022)  \nlog ˙ M B22 = -5 . 19 + 2 . 69 log( Γ e) -3 . 19 log(1 -Γ e) . (4)  \nIt is easy to see that, for Γ e = 0 . 5 Eq. 4 becomes log ˙ M ≃ -5 . 04. This value is very close of the \'transition\' mass-loss rates empirically determined for the Arches cluster by Vink & Gräfener (2012) and for 30 Dor by Sabhahit et al. (2023), i.e., the value for ˙ M over the which the stars of the mentioned clusters exhibit a WNh spectrum. This implies that both approaches (Sabhahit et al. 2023, and this current study) find the transition massloss rate between thin and thick winds at roughly the same value. Moreover, because this formula was calibrated for the VMS in the R136 star cluster at the LMC, with metallicity Z = 0 . 006 (Eggenberger et al. 2021), we add an extra term with the same metallicity dependence as for thin winds. Thus  \nlog ˙ M thick = log ˙ M B22 + log Z ∗ Z LMC ! × " 0 . 4 + 15 . 75 M ∗ / M ⊙ # , (5)  \nHence, despite we are using the same Γ e,trans for all metallicities, mass loss for optically thick winds incorporates a metallicity dependence which is consistent with the intrinsic dependence  \n284  \nwx  \n.  \n.  \non mass / luminosity. Even though calibration from Brands et al. (2022) matches with the formula of Krtiˇcka & Kubát (2018) (and subsequently with their intrinsic metallicity dependence on luminosity), we adopt the intrinsic metallicity dependence on mass introduced by Gormaz-Matamala et al. (2023) to keep coherence between our adopted thin and thick winds.  \nAnalogous to Eq. 3 for thin winds, Eq. 5 is valid only for T e ff ≥ 30 kK. For thick winds at cooler temperatures, we return to V01, thus Vink et al. (2001) and their predicted bi-stability jump are still being implemented for all the temperature range between 10 and 30 kK.', '2.3.2. Mass loss for advanced evolutionary stages': 'The Eq. 5 from Bestenlehner (2020) and Brands et al. (2022) is valid for H-rich WR stars, either for the H-core or the He-core burning stages. Once the star has depleted the most part of its hydrogen ( X surf ≤ 10 -7 ), the mass-loss rate follows the formula derived from the hydrodynamically consistent models for WR stars from Sander et al. (2020), whose final formula ˙ M ∼ Γ e is given by Sander & Vink (2020)  \nlog ˙ M = a log[ -log(1 -Γ e)] -log 2 Γ e,d Γ e ! cd , b + log ˙ M o ff , (6)  \nwhere  \nΓ e,b = -0 . 324 ∗ log Z ∗ Z ⊙ ! + 0 . 244  \ncd , b = -0 . 44 ∗ log Z ∗ Z ⊙ ! + 9 . 15  \nlog ˙ M o ff = 0 . 23 ∗ log Z ∗ Z ! -2 . 61  \na = 2 . 932 , , , ⊙ .  \nWe do not incorporate the latest improvements to this formula from Sander et al. (2023), because they add extra terms dependent on a quantity T ∗ which is not equivalent to the e ff ective temperature T e ff , and thus such implementation is currently beyond the capacities of both codes MESA and G enec . The only condition of validity for Eq. 6 is the criterion X surf ≤ 10 -7 (Hpoor WR star), regardless of the Eddington factor of luminosity. Finally, we also mention that all the wind recipes considered here are thought for single stellar evolution, and for that reason we do not include wind recipes for peculiar scenarios coming from binary interaction, such as stripped low-mass helium stars (Vink 2017) or O-type stars with He-core burning (Pauli et al. 2022).', '2.4. Treatment of rotation': 'The most important di ff erence concerning rotation between MESA and G enec is the standardisation. MESA sets stellar rotation as a fraction of the critical angular rotational velocity, ω , needed to disrupt the star due to centrifugal forces (Paxton et al. 2013), i.e. 4  \nΩ MESA = ω ω crit , with ω 2 crit = GM R 3 (1 -Γ Edd) , (7)  \nwhere Γ Edd is the total Eddington factor where κ is actually the Rosseland mean opacity, calculated as a mass-weighted average in a user-specified optical depth range.  \nOn the other hand, G enec sets stellar rotation as a fraction of the linear rotational velocity ( Ω G enec = 3 / 3 crit). However, these two fractions are not equivalent because a rotating star loses its spherical shape and becomes oblate, thus meaning that the polar radius and the equatorial radius are not longer the same. Even though both codes incorporate oblate e ff ects, the formal di ff erence between linear and angular velocity is not considered in MESA (Choi et al. 2016) but only in G enec , by means of the Roche-model (Ekström et al. 2008; Georgy et al. 2011). According to this, the ratios for the critical angular and linear velocities are related in G enec as  \n3 3 crit = " ω ω crit 2 R eq R pol -1 !# 1 / 3 , with 3 2 crit = GM R e,crit , (8)  \nwith R eq and R pol being the equatorial and the polar radii respectively, and R e,crit the [maximum] equatorial radius when 3 crit is reached. The standard setup has been Ω = 0 . 4 (either ω/ω crit or 3 / 3 crit), based on the survey on Galactic BSG performed by Huang et al. (2010). For the purpose of this paper, we set for both MESA and G enec codes the initial rotation as Ω = ω/ω crit = 0 . 4 (equivalent to 3 / 3 crit ≃ 0 . 276 in G enec ).  \nThe treatment of the stellar rotation also varies between both codes, not only for the mixing for also for the enhancement of the mass loss. By default MESA follows the standard correction factor from Friend & Abbott (1986) and Langer (1998)  \n˙ M ( ω ) = ˙ M ( ω = 0) 1 1 -ω/ω crit ! 0 . 43 , (9)  \nwith the ratio ω/ω crit given by Eq. 7. MESA also incorporates a limit to the maximum mass loss following Yoon et al. (2010), to avoid large divergences of Eq. 9 due to the dependence of ω crit on the Eddington factor Γ Edd.  \nOn the contrary, G enec uses as correction factor the more complex formula provided by Maeder & Meynet (2000)  \n˙ M ( ω ) = ˙ M (0) (1 -Γ e) 1 α -1 h 1 -ω 2 2 π G ρ m -Γ e i 1 α -1 , (10)  \nwith ρ m being the internal average density, α the CAK line-force parameter, and Γ e the electron scattering Eddington factor as given by Eq. 1. Eq. 10 considers the break-up velocity not only with respect to the rotation but also with respect to the proximity to the Eddington limit ( ΩΓ -limit), together with being in rule with the so-called Ω -slow solutions (Curé 2004; Araya et al. 2018; Gormaz-Matamala et al. 2023).', '3.1. H-core burning stage': "Fig. 1 shows the evolutionary tracks across the HR diagram according to our new G enec and MESA models, for both 60 M ⊙ and 200 M ⊙ . Fig. 2 shows the evolution of the mass loss, the stellar mass, the stellar radius, and the rotational velocity, from the ZAMS until the end of the H-core burning stage.  \nWe observe that, even though the ZAMS point in the HRD is almost the same for the MESA and G enec tracks, MESA models start with a higher mass loss due to the more rigid e ff ect from Eq. 9, making a di ff erence of ∼ 0 . 2 dex in logarithmic  \n<!-- image -->  \nFig. 1: HR diagram for G enec (red) and MESA (blue) evolution models with 60 M ⊙ (left panel) and 200 M ⊙ (right panel). Black dots represent the point where Γ e = 0 . 5 and thus the wind regime transits from thin to optically thick, black crosses represent the end of the H-core burning stage, black stars represent the beginning of the He-core burning stage, cyan stars represent the switch between H-rich to H-poor helium stars, and magenta stars represent the end of the He-core burning stage.  \n<!-- image -->  \nTable 1: Properties of our evolution models during the H-core burning stage.  \nFig. 2: Left panel: Evolution of mass-loss rate, stellar mass, stellar radius, and Eddington factor during the H-core burning stage for G enec and MESA evolution models with M zams = 60 M ⊙ . Right panel: same, but for evolution models with M zams = 200 M ⊙ .  \n<!-- image -->  \nscale. However, for our M zams = 60 M ⊙ models, this di ff erence is not large enough to produce significant deviations at the beginning, and thus both tracks move upwards increasing luminosity at temperature nearly constant. Both models reach Γ e = 0 . 5 prior of reaching their maximum radial expansion in the MS, but our G enec model reaches a higher temperature because of its slightly slower mass loss, its faster rotation, and a more e ffi cient mixing. This is evidenced in Table 1, where the MESA model reaches to the switch point from thin to thick winds (black dots in Figures 1 and 2) with lower mass and lower mass fractions of helium and nitrogen, in contrast with the G enec model. Once the winds for both models become optically thick, the mass loss in both cases is described either by Eq. 5 or V01 / GH08 formulae depending on the condition of T e ff being higher or lower than 30 kK. Because our G enec model has a larger luminosity and larger Γ e, its optically thick ˙ M will be larger enough to remove more layers and therefore to end the H-core burning stage with less mass than our MESA model (42 . 1 M ⊙ vs 45 M ⊙ ). However, the most remarkable di ff erence is the final chemical composition at the H-depletion: our G enec model ends with a helium mass fraction of Y surf = 0 . 929, i.e., stellar structure is almost homogeneous, and much more compact ( ≃ 5 R ⊙ ) and hotter ( T e ff ≃ 90 kK).  \nNotice that, for both 60 M ⊙ models, the transition occurs outside the range of validity of both ˙ M thin and ˙ M thick, where V01 dominates both thick and thick wind regimes (Romagnolo et al. 2024a, Fig. 1) and the only jump in mass loss is due to the bistability jump (Vink et al. 1999), at T e ff ≃ 28 kK. Despite this, the fact of reaching Γ e = 0 . 5 during the main sequence implies that this point also represents the transition from thin to thick winds, as we will discuss later in Sec. 4.1.  \nFor our M zams = 200 M ⊙ models, the initial mass loss start around 10 -5 M ⊙ yr -1 , and again the MESA model starts with a higher ˙ M due to Eq. 9. The initial rotational velocity substantially di ff ers between both models, because both models calculate a slightly di ff erent starting point for the ZAMS in the HRD as a consequence of both the numerical uncertainties (Agrawal et al. 2022b, which tend to increase for larger masses), and the di ff erent calculation of the critical rotational velocity outlined in Sec. 2.4. Despite these corrections, both models show a very similar track across the HRD, reaching Γ e = 0 . 5 at the very beginning of their lives and subsequently the stars spend more than their 90% of their H-core burning lifetime in the optically thick wind regime. Nonetheless, in this mass regime the mass-loss escalates with metallicity only as ˙ M ∝ Z 0 . 48 (according to Eq. 5) and thus our initial mass loss for the thick wind regime is around log ˙ M ≃ -4 . 8, not large enough to produce an abrupt decrease in the luminosity, reason why both tracks keep moving upwards in the HRD diagram until log L / L ⊙ ≃ 6 . 7. Once the mass loss increases enough to stop the rise in luminosity, both models move bluewards keeping L ∗ almost constant. This time, the MESA model drops its rotational velocity slower than the G enec model, keeping a more moderate radius and ending their H-core burning stage with a larger e ff ective temperature ( ∼ 100 kK, against the ∼ 80 kK from the G enec model). Despite that our G enec expanded more redwards, it ends its H-burning with a more advanced chemical enrichment at the surface. Even though this could be partially attributed to the larger initial rotation, it is also because that G enec implements a more e ffi cient rotational mixing (Meynet & Maeder 2000; Choi et al. 2016).  \nThe correlation between the both tracks shows the most important di ff erence with respect to previous evolution studies adopting V01, where it is found that for M zams = 200 M ⊙ the mass loss starts being strong enough to make the star quickly drop in luminosity (Yusof et al. 2013; Sabhahit et al. 2022,  \n2023; Martinet et al. 2023). Di ff erences in mass loss also carries another important discrepancy with respect previous studies: our models end the H-core burning stage with final masses of ∼ 100 . 5 M ⊙ . This is larger than the M final ≃ 85 M ⊙ found by Martinet et al. (2023) for M zams = 180 M ⊙ , which is the study with the structure setup closest to ours ( α ov = 0 . 2 and Ledoux convection criterion) but still adopting old winds; but our result is even more remarkably larger than the M final ≃ 40 M ⊙ found by Sabhahit et al. (2022, 2023), who adopted a ˙ M thick from Vink et al. (2011). In contrast with these huge mass loss, our models adopting weaker (thin and thick) winds makes a M zams = 200 M ⊙ star lose 'only' one half of their mass at the end of H-core burning stage. This carries an important impact related to the spectroscopic diagnostics to be done over Galactic VMS, as we will see in Sec.4.3.", '3.2. He-core and C-core burning stages': 'Given that the star 60 M ⊙ from our MESA model ends the MS with a less homogeneous chemical structure (only a ∼ 57% of helium mass fraction, in contrast with the ∼ 93% predicted by G enec ), it quickly expands until reaching T e ff ≃ 8 000 K with a maximum radius of ≃ 740 R ⊙ . This is because stars keeping a more constant chemical composition from the core to the surface remain in the bluer part of the HR diagram (Maeder 1987). The beginning of the He-core burning stage (black star symbol in the Figures) makes the star return bluewards, thus increasing its temperature and its mass loss. On the contrary, G enec model predicts a compact star which enters into the He-core burning stage without any radial expansion. This divergence on the models indicate that M zams ≃ 60 M ⊙ represents a threshold for the Hertzsprung-gap: less massive stars will all expand post MS despite minor di ff erences in the rotational mixing, whereas more massive will not (see Fig. 2 from Romagnolo et al. 2024a). These initial di ff erences during He-core ignition later converge into the position at the HRD where both stellar models become H-poor WRstars. Hence, both models show a similar track for the rest of the He-core burning, with a large diminishing in the luminosity from log L / L ⊙ ≃ 6 . 2 to ≃ 5 . 75 at temperature almost constant T ∗ ≃ 140 kK. As a consequence, the final masses at the end of the He-core burning stage are relatively the same for MESA and G enec tracks ( ≃ 16 . 7 M ⊙ ), and thus the remaining evolution through the C-core burning is practically identical for both models.  \nFor the 200 M ⊙ case, di ff erences in both the stellar mass and the chemical composition after the H-depletion are almost negligible, and just the HRD where the He-ignition occurs is significantly relevant. Because of its proximity to the Eddington limit (with Γ e ≃ 0 . 7) and lower density, the convection at the outer envelope of a star born with M zams = 200 M ⊙ is highly superadiabatic (i.e, the temperature gradient is larger than the adiabatic temperature gradient) and therefore the energy released in the interior remains locked at the subsurface layers, thus creating numerical instabilities in the models where each evolution code adopts di ff erent strategies to handle it (Agrawal et al. 2022a). MESA reduces the superadiabacity in radiativelydominated regions by means of the use\\_superad\\_reduction method, which can be interpreted as a lowering in the opacity κ ( r ) in the sub-surface regions near the Eddington limit (Jermyn et al. 2023). On the contrary, our 200 M ⊙ G enec does not incorporate this superadiabacity rescaling and thus excess in the temperature gradient leads into a considerable inflation when the last traces of hydrogen are removed from the surface and the opacity drastically changes (Yusof et al. 2013; Martinet et al. 2023),  \nuntil this excess of energy is released due to Supra-Eddington mass loss (Ekström et al. 2012) and the star returns to its initial radius. During this expansion, only less than the 0 . 1% of the total mass is inflated, keeping the remaining 99 . 9% of the total mass constrained into the same stellar radius, reason why this expansion should not be considered for any binary interaction Belczynski et al. (2022). Indeed, we show this expansion due to superadiabatic instability as a sudden expansion redwards in the HRD of Fig. 1, but we do not consider it for our plot of radial evolution during the He-core burning (Fig. 3) nor for the maximum radius in Table 2. Moreover, the strategies to deal with these superadiabatic instabilities are approached from a numerical point of view, and possibilities of a physical phenomenon have been barely considered. After this phenomenon, the star will decrease its luminosity at almost constant temperature until log L / L ∗ ≃ 6 . 0 according to both MESA and G enec models, then the He-core burning stage ends with a stellar mass of ≃ 25 . 5 M ⊙ being the star a ball of ≃ 94% carbon and oxygen. We ultimately reach the end of the C-core burning stage with final masses of 24 . 4 M ⊙ according to G enec and 24 . 8 M ⊙ according to MESA .', '4. Theoretical and observational diagnostics': '4.1. Self-consistent mass loss close to Γ e,trans  \nAs shown in Sec. 2.2.1, our evolution models use as a criterion to switch from optically thin winds of OB-type stars to thick winds of WR-type stars during the main sequence, the Eddington factor instead of the surface hydrogen abundance. In particular we use the transition defined by Bestenlehner (2020), which corresponds to the Γ e where the line-driven contribution to the CAK mass-loss equalises with the continuum-driven contribution. Below Γ e,trans, the wind is predominantly line-driven and thus ˙ M it can be described by our self-consistent m-CAK prescription, whereas over Γ e,trans the continuum radiation by electron scattering dominates.  \nTherefore, a question that arises is how the transition towards thick winds approaches from self-consistent mass-loss rates at optically thin regime towards Γ e,trans. To do so, we plot in Fig.4 the ˙ M sc as a function of the Eddington factor coming from the models introduced in Gormaz-Matamala et al. (2023), plus the G enec models introduced in this work. We can see that the ˙ M -Γ e easily correlates with a linear fit in the log-log plane, as  \nlog ˙ M sc ≃ -4 . 547 + 2 . 709 log Γ e (11)  \nthis rough fit of Eq. 3 agrees reasonably well with the ˙ M ∝ Γ 2 . 73 ± 0 . 43 e from Bestenlehner et al. (2014) and with the ˙ M ∝ Γ 2 . 69 e from (Brands et al. 2022), for the low Eddington factor regime. In other words, ˙ M sc behaves like Eq. 5 when Γ e ≲ 0 . 5, which is appreciated in Fig. 4, where the both fits run in parallel for Γ e ≪ 1. It is important to mention that even though the ˙ M -Γ e correlation is a very useful tool to describe the escalation of the mass loss in the proximity of the Eddington limit, it is not a formal function. The value of ˙ M may vary for a same value of Γ e, if the stellar mass or the luminosity or the hydrogen abundance are di ff erent. In particular, the linear behaviour of ˙ M sc dissipates as we approach Γ e ≃ 0 . 5, thus suggesting that the m-CAK prescription stops being valid for regimes with larger Eddington factor (where the contribution of the electron scattering to the radiative acceleration becomes more dominant) and a prescription considering optically thick winds needs to be adopted.  \nWe also observe in Fig. 4 that Eq. 11 and Eq. 5 intersect at log Γ e ≃ -0 . 62 ( Γ e ≃ 0 . 24), very close to the empirical kink found by Bestenlehner et al. (2014) at log Γ e ≃ -0 . 58. After that  \nTable 2: Properties of our evolution models during the He-core and C-core burning stages.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 3: Left panel: Evolution of mass-loss rate, stellar mass, stellar radius, and Eddington factor during the He-core and C-core burning stages for G enec and MESA evolution models with M zams = 60 M ⊙ . Right panel: same, but for evolution models with M zams = 200 M ⊙ .  \n<!-- image -->  \nTable 3: Lifetime (in Myr) predicted by our G enec and MESA models in the spectroscopic phases according to stellar evolution. The percentages with respect to the total lifetime are shown in parenthesis.  \npoint, formula for ˙ M thick winds clearly exceeds ˙ M thin, finishing with a di ff erence of ≃ 0 . 5 dex when Γ e = 0 . 5, as shown in Fig. 2. Hence, we could have a smooth transition from ˙ M thin to ˙ M thick by simply setting the transition at the point where both formulae intersect. However, we decide for the purposes of this paper to employ Γ e,trans as originally set by Bestenlehner (2020), to preserve the physical meaning and exploring the caveats.  \nThe rough fit introduced in Eq. 11 is set for solar metallicity, whereas the fit introduced in Eq. 5 corresponds to the calibration done by Brands et al. (2022) at LMC metallicity ( Z = 0 . 006) plus an extra metallicity-dependence prescribed by us. Actually, the general expression derived by Bestenlehner (2020) is  \nlog ˙ M = log M 0 + 1 α + 0 . 5 ! log( Γ e) -1 -α α + 2 ! log(1 -Γ e) , (12)  \nwith α being the CAK line-force parameter and ˙ M 0 the mass-loss rate at Γ e,trans. In principle, we could merge both Eq. 11 and Eq. 5 into one single expression as Eq. 12, by using the tabulated α val-  \nfrom Gormaz-Matamala et al. (2023). In such case, the transition from thin to thick winds would be smooth, given that the continuum already contributes to the wind mass-loss rate even below Γ e,trans = 0 . 5. However, that merging would require the analysis of di ff erent metallicities, and the present study is focused on the Milky Way. In a forthcoming paper we aim to explore the possibility of adopting one unique Eq. 12 for both thin and thick winds, together with more self-consistent calculations of the line-force parameters (Gormaz-Matamala et al. 2022a) for stellar models around Γ e,trans.', '4.2. Eddington factor and wind efficiency': 'The transition from thin to thick winds based on the proximity of the Eddington factor has been recently studied by Gräfener (2021) and Sabhahit et al. (2023), where the switch from OB to WNh stars is set by using the relationship between the optical depth at the sonic point of the wind structure ( τ s) and the wind e ffi ciency number ( η ≡ ˙ M 3 ∞ / ( L ∗ / c )) from Gräfener et al. (2017)  \nFig. 4: ˙ M -Γ e correlations taken from the G enec self-consistent evolution models from Gormaz-Matamala et al. (2023) and this work. The dashed line represents the best linear fit for these M sc, whereas the gray dot-dashed line represents the fit of Eq. 5 adopted for Z = 0 . 014.  \n<!-- image -->  \nτ s ≃ ˙ M 3 ∞ L ∗ / c 1 + 3 2 esc 3 2 ∞ ! , (13)  \nwith 3 ∞ and 3 esc being the wind terminal and the escape velocities respectively.  \nSince η denotes the ratio between the mechanical wind momentum and the total momentum of the radiatively-driven wind, the value η = 1 corresponds to the single scattering limit. If η ≥ 1, it means that photons can be absorbed and re-emitted through the wind structure more than once (Lucy & Abbott 1993). The condition η = τ = 1 at the transition point between thin and thick winds, relies in the assumption that the radiative acceleration largely exceeds gravity at larger distances. Posteriori calibrations found η ≃ 0 . 6 as the transition point, based on the observed mass-loss rates for Galactic and LMC WNh stars Sabhahit et al. (2022), plus finding that this η switch decreases for lower metallicities (Sabhahit et al. 2023). These calibrations make the value of this transition mass-loss rate ˙ M trans being independent of factors such as clumping. However, the temporal location of this thin-to-thick transition depends on the evolution of η , which in turns depends on the mass-loss recipe adopted a priori for the initial optically thin regime (V01 for this case, even for M zams ≃ 40 -60 M ⊙ ), besides of the 3 ∞ / 3 esc ratio assumed for β -law wind velocity profiles.  \nOn the contrary, the proposed transition towards thicker winds at Γ e = κ e L ∗ / 4 π cGM ∗ = 0 . 5 is explicitly independent on the mass-loss regime adopted for the optically thin regime, together with not needing any β -law assumption for the velocity profile (Gormaz-Matamala et al. 2022a). In this case, we enter into the optically thick wind regime because the continuum contribution to the radiative-driving mechanism becomes more relevant as we approach the Eddington limit, reason why the Γ e,trans is almost constant for di ff erent values of α CAK parameter. Hence, this transition di ff ers from the transition based on the single-line scattering limit established by Sabhahit et al. (2022, 2023), where the kink in the ˙ M -Γ e is due to the dominance  \nFig. 5: η -Γ e correlations taken from theoretical models tabulated by Krtiˇcka & Kubát (2017), Björklund et al. (2021), and GormazMatamala et al. (2023).  \n<!-- image -->  \nof multi-scattering for the line-driven mechanism. This means, both transitions are connected to a di ff erent physical mechanisms at di ff erent wind conditions. The existence of the transition at Γ e,trans = 0 . 5 where continuum-driven wind becomes predominant does not contradict the existence of a limit where the line-driving needs to adopt multi-scattering. Indeed, the trend shown in Fig. 5 suggests that the condition η ≃ 0 . 6 should be satisfied at some point Γ e ≫ 0 . 5, even though it is not possible to assure it properly because of the lack of studies adopting new ˙ M recipes (Krtiˇcka & Kubát 2017; Björklund et al. 2021; Gormaz-Matamala et al. 2023) for Γ e > 0 . 4, whereas the transition Eddington factor determined by Sabhahit et al. (2022) is around ∼ 0 . 4 (see their Table 3). This implies, that either new winds predict a transition wind e ffi ciency lower than the η switch found by keeping V01 for thin winds (thus meaning that the η transition depends on the mass-loss rate adopted for low Γ e), or the transition to thick winds might be better represented with the predominance of electron scattering for the radiative acceleration when Eddington factor goes around ∼ 0 . 4 -0 . 5. As discussed in Sec. 4.1, this is exactly the range of Γ e where self-consistent mass-loss rate shows more dispersion.  \nAnother aspect that should be explored, is the metallicity dependence. Sabhahit et al. (2023) finds that the η switch value where the wind becomes optically thick varies with the metallicity, given the direct dependence on ˙ M and terminal velocity for the wind e ffi ciency number (Vink & Gräfener 2012), and the adopted scaling for the terminal velocities a di ff erent Z . From our side, the transition based on Γ e,trans = 0 . 5 suggests to be invariant for di ff erent metallicities, because of its invariance for di ff erent values of the line-force α (which it is metallicity dependent). However, as seen in Sec. 4.1, a deeper analysis of the m-CAK around Γ e > 0 . 4 will be required at various values of Z , both to evaluate the possibilities of some metallicity dependence for Γ e,trans but also to get ˙ M thick( Z ). Up to now, we just have the calibration done by Brands et al. (2022) plus the ˙ M -Z dependence from Gormaz-Matamala et al. (2023), fitted only for the mass range of 25 -120 M ⊙ and thus extrapolated up to M zams = 200 M ⊙ in this work.  \nFor a more general prospective, it is important to mention the increase of the Eddington factor as the star evolves through the main sequence, with weaker winds boosting this tendency for the cases with M zams ≳ 70 M ⊙ (Gormaz-Matamala et al. 2023, Fig. A7). The more massive the star, the faster it reaches Γ e,trans = 0 . 5 and therefore the shortest its fraction of MS lifetime adopting ˙ M sc. In the case of the evolution models introduced in this work, a 60 M ⊙ star will spend ∼ 92% of its MS lifetime exhibiting optically thick winds ( ∼ 92% for the MESA model), in contrast with the ∼ 7% of the MS lifetime for a 200 M ⊙ star ( ∼ 4%for the MESA model). These percentages of the lifetime are given by the location and track across the HRD and not by their ˙ M adopted during the optically thin regime. Therefore, the chosen stellar masses for the analysis in this work represent the approximated limits where optically thin and thick winds can coexist during the main sequence stage. For stars with M zams < 60 M ⊙ the limit Γ e = 0 . 5 may not be reached during the MS and therefore the stars will just become classical WR after the depletion of hydrogen, as in traditional evolution models. On the other limit, for stars with M zams ≥ 200 M ⊙ it is satisfied that Γ e = 0 . 5 from the ZAMS and therefore the wind is fully optically thick, as also found by Sabhahit et al. (2023). In the intermediate mass range, stars will spend one fraction of their lifetime with thin OB-type winds and reduced mass-loss as reported by Gräfener (2021), and the other fraction with the enhanced ˙ M if WNh stars. Summarising, the new weaker winds from Gormaz-Matamala et al. (2022b, 2023) can be applied for all initial masses for the earliest stages, even for VMS, for the segment of time when Γ e < 0 . 5 is satisfied.', '4.3. Spectroscopic analysis': "The new mass-loss prescription, plus the new transition from thin to thick winds, modify the evolutionary tracks across the HRdiagram and thus the spectroscopic phases. According to the output of our evolution models, we can classify the spectroscopic phases of our evolution models according to  \n- -O-type: optically thin wind and T e ff ≥ 26 . 3 kK.\n- -B-type: optically thin wind and T e ff < 26 . 3 kK.\n- -eWNL: WR star with X surf ≥ 0 . 05.\n- -eWNE: WR star and X surf < 0 . 05 and ( 12 C + 16 O) / 4 He ≤ 0 . 03.\n- -WC: WR star wind and ( 12 C + 16 O) / 4 He > 0 . 03.\n- -WO: WR star with log T e ff ≥ 5 . 25.  \nThese WR boundaries are derived from the conditions settled by Yusof et al. (2013), Aadland et al. (2022), and Martinet et al. (2023); whereas the distinction between O and B-type is taken from Groh et al. (2014). It is important to note, that this classification is not based on spectral analysis but on the output of our theoretical models, and thus discrepancies may be observed. For example, we separate the evolutionary eWNL and eWNE Foellmi et al. (where 'e' stands from evolutionary, following 2003) according to the hydrogen abundance as done by Meynet & Maeder (2003), but the separation between WNL and WNE spectral subtypes is based on the ionisation levels of the surface nitrogen (Smith et al. 1996). This could suggest that all WNE stars are H-poor, albeit there are evidence of few H-rich WNEstars in the Milky Way (Hamann et al. 2019). These H-rich WNE stars however, have luminosities below log L ∗ / L ⊙ ≃ 5 . 5 which suggest their progenitors were stars with M zams ≃ 40 M ⊙ (Georgy et al. 2012), thus being out of the range of our analysis. The results are shown in Table 3, whereas the location of these spectral phases across the HRD is shown in Fig. 6.  \nEven though is not a formal spectroscopic phase, we classify our star as a luminous blue variable (LBV) if the HD-limit ( L > 6 × 10 5 and 10 -5 L 1 / 2 R > 1) is exceeded (Humphreys & Davidson 1994; Hurley et al. 2000), although recent works have found a lower luminosity limit of L > 3 × 10 5 (Davies & Beasor 2020). This is particularly evident for our 60 M ⊙ MESA model, where the post MS expansion exceeds the HD limit and therefore the star spends 12 kyr (0 . 26% of its total lifetime) as a LBV. Despite this short period, the star had a mass of M pre,lbv ≃ 44 . 1 M ⊙ and X H ≃ 0 . 4, to finish with M post,lbv ≃ 41 . 9 M ⊙ and X H ≃ 0 . 24, i.e., it experienced an average mass loss of ˙ M lbv ≃ 1 . 8 × 10 -4 M ⊙ yr -1 (log ˙ M lbv ≃ -3 . 7). Such value is close to the adopted ˙ M lbv by Belczynski et al. (2010), but still it comes from the application of the wind recipes of Vink et al. (2001) and de Jager et al. (1988) in a region of the HRD beyond the HD limit. Models for eruptive mass loss (Quataert et al. 2016; Cheng et al. 2024), more in agreement with the nature of LBVs, would be more require to better represent the evolution during this phase. In parallel, our 60 M ⊙ G enec does not become a LBV because its more chemically homogeneous structure avoids any radial expansion after the H-depletion. Therefore, despite that the final masses tabulated for both codes are quite similar, the prediction about the existence or not of the LBV phenomenon strongly depends on the calculated mixing during the main sequence stage. In opposition, our 200 M ⊙ G enec evolution track moves redwards, crossing the HD limit, between its phases WNE and WC. This reinforces the uncertain relation of the LBV phenomenon with the evolution of massive stars because, if we consider that redwards expansion as a physical phenomenon instead of a numerical artefact, it would mean that LBV events could occur even after the star becomes a hydrogen depleted WN star. However, this has never been observationally found, and then the idea of a numerical artefact due to a drastic change in the opacity due to the removal of the last traces of hydrogen looks more compelling.  \nBesides the standard classification for WR stars, in Table 3 and Fig. 6 we also separate the WNL stars with X H > 0 . 3. Given that typical stellar evolution models enter into the WR phase when either X H < 0 . 3 (Ekström et al. 2012) or X H < 0 . 4 Brott et al. (2011), this distinction is made to incorporate the WNh stars whose optically thick winds are because of the proximity to the Eddington limit before the H-depletion, according to the criterion Γ e ≥ 0 . 5 established in Sec. 2.2.1. Thus we discern from the WNL stars which are result of the 'classical' evolution, where the star becomes a WR due to the drift bluewards after the LBV episodes during the He-core burning stage (Conti 1975; Groh et al. 2014). Indeed, WNh stars with X H ≥ 0 . 3 are constrained to a very narrow band of temperature in the HRD, log T e ff ≃ 4 . 5 -4 . 7 (Martins et al. 2008; Crowther et al. 2010; Hamann et al. 2019), whereas WNL stars with X H ≤ 0 . 3 are spread across a wider range of temperatures (Martins 2023). To illustrate better this point, we select the Galactic WNh stars (with X H ≥ 0 . 3) from the Arches cluster (Martins et al. 2008), from the NGC 3603 cluster (Crowther et al. 2010), and from the local neighbourhood (Hamann et al. 2019), and we tabulated them in Table 4. In addition, we plot our G enec and MESA models for 60 and 200 M ⊙ , highlighting the segments where Γ e ≥ 0 . 5 and X H ≥ 0 . 3. The results are shown in Fig. 7. Despite the limited di ff erences between the evolution tracks, our prescription adequately describes the narrow band where WNh exist, thus meaning that WNh are the product of stars born with M zams ≳ 60 M ⊙ which produce optically thick winds before the H-depletion due to their proximity to the Eddington limit. In the middle we append models for 100 M ⊙ , showing that the thick band also fits the WNh of our sample.  \nFig. 6: Distribution across the HRD of the evolutionary spectral types introduced in Sec. 4.3. The red and blue dots represents time intervals of 20 000  \n<!-- image -->  \nFor the 200 M ⊙ cases the observational diagnostics are more uncertain. For instance, the stars tabulated in Table 4 are also the brightest WNh stars in the Milky Way, and none of them reaches the luminosity of log L ∗ / log L ⊙ ≳ 6 . 55. The early approaching to the Eddington limit from the very beginning of the main sequence implies that the O-type phase in our models is very short, and therefore the star will spend the largest percentage of their lifetime (almost a ∼ 90% for both models) as a WNh. This is a remarkable di ff erence with respect to previous studies, where very massive stars still spent the majority of their lifetime as Otype stars until they remove more than the half of their initial surface H abundance (Martinet et al. 2023). In other words, our models predict that extremely luminous Galactic stars have more probabilities of being WNh than O-type. The star W49-#2, with log L ∗ / L ⊙ ≃ 6 . 64 ± 0 . 25 (one of the brightest stars known in the Milky Way) is classified as O2-3.5If* (Wu et al. 2016) but quite close of being a Of / WN star based on the spectral classification of Crowther & Walborn (2011). In any case, the high stellar mass ( ∼ 250 M ⊙ ) estimated by Wu et al. (2016) implies Γ e ≃ 0 . 35 -0 . 46 (depending on the X H to be considered in Eq. 2), thus not being close to the Eddington limit enough to develop optically thick winds.  \nBeyond WNh stars, our new evolution models for 200 M ⊙ also predicts the existence of H-poor WR stars with luminosities larger that log L ∗ / L ⊙ > 6 . 0, same as the rotating models of 180 M ⊙ and 250 M ⊙ from Martinet et al. (2023). If we sum up the lifetime of eWNE, WC and WO spectroscopic phases from Table 3, we find that the our 200 M ⊙ models will spent a ∼ 14% of their lifetime according to the G enec model ( ∼ 13% in the case of MESA ), very close to the ∼ 11% for the same phases predicted by them. Such small discrepancy can be easily explained by our choice of setting eWNE for X surf ≤ 0 . 5 (Yusof et al. 2013, from) instead of ≤ 10 -5 , which implies that eWNE will form around the ∼ 3 . 5% of the total lifetime (either H-rich or H-poor) in contrast with the absence of eWNE phase for 200 M ⊙ from Martinet et al. (2023). Nonetheless, only two WN (WR 18 and WR 37) stars with X H = 0 satisfies log L ∗ / L ⊙ > 6 . 0 in the Galactic catalog of Hamann et al. (2019); whereas for WC type, WR 126 is the only WN / WC star with log L ∗ / L ⊙ > 6 . 0 in the Milky Way  \n(Sander et al. 2019). However, the trace of the red and blue dots from Fig. 6 confirms that eWNE is a very short period, and these stars quickly drops their luminosities from log L ∗ / L ⊙ ≥ 6 . 6 to ≃ 6 . 0 in just ∼ 200 kyr. Even though the mentioned catalogs still cover a small region of the Milky Way, the absence of WNE and WC stars with log L ∗ / L ⊙ ≥ 6 . 1 (which should be evolved from models born with masses ≃ 60) should suggest that Eq. 5 might be underestimating the actual value of the mass-loss rate for the most massive H-rich WN stars, as a consequence of the extrapolation up to 200 M ⊙ mentioned in Sec. 4.2. Even though it is still possible that we might have been underestimating the actual value of the mass-loss rate for the most massive H-rich WN stars as a consequence of the extrapolation up to 200 M ⊙ mentioned in Sec. 4.2 (and thus overpredicting the existence of Hpoor WR stars at larger luminosities), it is important to mention that our Galactic WR catalogs are still quite constrained to our solar neighbourhood. The study and discoveries of new sources at the edge of the known stellar luminosities, together with less uncertain data, is beyond the scope of this paper.", '5. Summary and conclusions': "In Paper I (Romagnolo et al. 2024a) we explored the expected final masses prior to the core collapse, when we calculate stellar evolution updating the mass-loss prescription. In this work, we confirm the importance of these new wind prescriptions, by analysing its impact over the expected spectroscopic stages (OBtype, WNh and classical WR).  \nRegardless of the unification of the overshooting ( α ov = 0 . 5), the rotation ( ω/ω crit = 0 . 4), the angular momentum transport (Tayler-Spruit dynamo), the convective boundary criterion (Ledoux), and the line-driven mass-loss prescriptions for both codes G enec and MESA , there are important di ff erences between the calculated evolutionary tracks. For instance, G enec computes a more prominent rotational mixing during the main sequence phase, leading to a more homogeneous stellar structure at the Hcore depletion in comparison with MESA . This is particularly evident for our 60 M ⊙ models, where our less homogeneous model largely expands and crosses the HD limit before the He-ignition.  \nO  \nB  \nWN  \neW  \neW  \nWC  \nWO  \nLBV  \nTable 4: WNh stars included in Fig. 7.References. (1) Martins et al. (2008); (2) Crowther et al. (2010); (3) Hamann et al. (2019).  \nFig. 7: Evolutionary G enec and MESA tracks for 60 and 200 M ⊙ from Fig.1 and additional 100 M ⊙ models, where the thick solid lines represent the segments where the wind is optically thick ( Γ e ≥ 0 . 5) and X surf ≥ 0 . 3. The gray stars correspond to the sample of WNh stars from Martins et al. (2008), the purple stars to the sample of Crowther et al. (2010) and the green stars to the sample of Hamann et al. (2019), tabulated in Table 4. Yellow shadowed area correspond to the HRD region where the HD limit is exceeded.  \n<!-- image -->  \nSuch discrepancy between both models can be interpreted as a boundary in the mass range below which stars normally expands, and over which stars keep a moderate size because of the strong outflows (Romagnolo et al. 2023). For the 200 M ⊙ models, due to the large size of the convective core the di ff erences between  \nboth evolution codes are less relevant, and only the treatment for superadiabacity during the H-core burning phase is a matter of debate.  \nConcerning stellar masses, both G enec and MESA codes show similar results at the end of each core burning process, with the most remarkable mass losses during the post-MS stages. Models with M zams = 60 M ⊙ end the H-core burning with final masses of 42 . 1 and 44 . 5 M ⊙ respectively, subsequently experiencing strong outflows during He-core burning stages to finish with 16 . 1 M ⊙ at the C-core depletion. This is a striking discrepancy with the M bh ≃ 30 M ⊙ found at M zams = 60 M ⊙ by Bavera et al. (2023), where the rotational e ff ects were not considered, or with Vink et al. (2024), where the WR stage is suppressed for M zams = 40 M ⊙ . For our 200 M ⊙ case, where the main sequence is almost fully dominated by thick winds, both G enec and MESA models predict a final mass of ∼ 100 M ⊙ at the end of the H-core burning stage, finishing later with M fin ≃ 24 . 6 M ⊙ at the C-core depletion. This result is also largely discrepant with other studies such as Yusof et al. (2013) or Sabhahit et al. (2022, 2023), who found that a star with M zams = 200 M ⊙ would loss almost the ∼ 80% of its mass during the H-burning stage only. Despite the discussion introduced in Section 4.1, such di ff erence in the calculated mass remnants is not produced due to the approach of the transition from thin to thick winds close to the Eddington factor, but due to the lower recipes for ˙ M for both thin -(GormazMatamala et al. 2023) instead of V01- and thick -Bestenlehner (2020) instead of Vink et al. (2011)- wind regimes.  \nHowever, the most remarkable result is that our transition between the optically thin winds of OB-type stars and the thick winds of WNh stars based on the dominance of electron scattering over line-driving, results in the prediction of a band in the HRD where stars are expected to develop optically thick winds but with large hydrogen fractions at the surface ( X surf ≫ 0 . 3). This band corresponds to e ff ective temperatures between log T e ff ≃ 4 . 5 -4 . 7, in agreement with the spectroscopic diagnostics of Martins et al. (2008), Crowther et al. (2010), and Hamann et al. (2019). The basis for the proposed switch of wind regime at Γ e,trans = 0 . 5 follows the formulation of Bestenlehner (2020) based on CAK theory, and it adequately considers that at early evolutionary stages the ˙ M for regular massive stars (i.e., not VMS) is ∼ 3 times lower than the traditionally assumed for V01. This statement is compatible with the idea of having WNh stars born from M zams ≃ 60 M ⊙ , whereas we find reduced ˙ M from hydrodynamic analysis on winds for the same mass range (Sander et al. 2017; Gormaz-Matamala et al. 2021). As we escalate in the mass range, the fraction of MS time at which the stellar wind is optically thin becomes shorter, and therefore the wind of very massive stars posses enhanced values of mass-loss rates because they spent larger fractions of their MS dominated by thick winds. Stars at M zams ≳ 200 M ⊙ are already dominated by thick winds from the very beginning of their lifetimes, similar as found by Sabhahit et al. (2023).  \nThis work focuses on calculations applicable to stars with solar metallicity and in earlier stages of evolution. However, recent observations reveal a significant population of Wolf-Rayet stars (WNh) in young and metal-poor clusters like R136 (Brands et al. 2022) and 30 Dor Ramírez-Agudelo et al. (2017), both located in the LMC. This necessitates extending our procedures and diagnostics to encompass very massive stars with metallicities as low as Z ≲ 0 . 006, as suggested by Martins & Palacios (2022). Additionally, these upgrades in mass-loss prescriptions need to be applicable to stars with T e ff ≤ 30 kK. Evaluating approaches like those proposed by Krtiˇcka et al. (2024) for B-supergiants is a promising avenue for future work, together with incorporating  \nmore detailed prescriptions for the stellar interiors based on the state-of-the-art asteroseismology (Pedersen et al. 2021).  \nAcknowledgements. We thank to the anonymous referee for their valuable comments and feedback, together with to Joachim Bestenlehner, Sylvia Ekström, and Michel Curé for the fruitful discussions. Authors acknowledge support from the Polish National Science Center grant Maestro (2018 / 30 / A / ST9 / 00050). The authors also thank the G enec and MESA communities at large for their helpful feedback on the models creation. Computations for this article have been performed using the computer cluster at CAMK PAN. ACGM thanks the support from project 10108195 MERIT (MSCA-COFUND Horizon Europe). We dedicate this paper to Krzysztof Belczy'nski, who contributed to this research before his untimely passing on 13th January 2024.", 'References': '- Aadland, E., Massey, P., Hillier, D. J., et al. 2022, ApJ, 931, 157\n- Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102\n- Abbott, D. C. 1982, ApJ, 259, 282\n- Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 913, L7\n- Agrawal, P., Hurley, J., Stevenson, S., Szécsi, D., & Flynn, C. 2020, MNRAS, 497, 4549\n- Agrawal, P., Stevenson, S., Szécsi, D., & Hurley, J. 2022a, A&A, 668, A90  \nAgrawal, P., Szécsi, D., Stevenson, S., Eldridge, J. J., & Hurley, J. 2022b, MN-  \nRAS, 512, 5717  \n- Araya, I., Curé, M., ud-Doula, A., Santillán, A., & Cidale, L. 2018, MNRAS, 477, 755\n- Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 Bara ff e, I., Clarke, J., Morison, A., et al. 2023, MNRAS, 519, 5333\n- Bavera, S. S., Fragos, T., Zapartas, E., et al. 2023, Nature Astronomy, 7, 1090\n- Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217\n- Belczynski, K., Hirschi, R., Kaiser, E. A., et al. 2020a, ApJ, 890, 113\n- Belczynski, K., Klencki, J., Fields, C. E., et al. 2020b, A&A, 636, A104  \nBelczynski, K., Romagnolo, A., Olejak, A., et al. 2022, ApJ, 925, 69  \n- Bestenlehner, J. M. 2020, MNRAS, 493, 3938\n- Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38  \nBjörklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, 648, A36  \n- Björklund, R., Sundqvist, J. O., Singh, S. M., Puls, J., & Najarro, F. 2023, A&A, 676, A109\n- Bouret, J. C., Lanz, T., & Hillier, D. J. 2005, A&A, 438, 301\n- Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36\n- Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115\n- Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157\n- Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068\n- Cheng, S. J., Goldberg, J. A., Cantiello, M., et al. 2024, arXiv e-prints, arXiv:2405.12274\n- Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102\n- Conti, P. S. 1975, Memoires of the Societe Royale des Sciences de Liege, 9, 193\n- Crowther, P. A. 2007, ARA&A, 45, 177\n- Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731\n- Crowther, P. A. & Walborn, N. R. 2011, MNRAS, 416, 1311\n- Curé, M. 2004, ApJ, 614, 929\n- Curé, M. & Rial, D. F. 2004, A&A, 428, 545\n- Davies, B. & Beasor, E. R. 2020, MNRAS, 493, 468\n- de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259 Eggenberger, P., Ekström, S., Georgy, C., et al. 2021, A&A, 652, A137  \nEggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43  \n- Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146  \nEkström, S., Meynet, G., Maeder, A., & Barblan, F. 2008, A&A, 478, 467  \n- Farrell, E., Groh, J. H., Hirschi, R., et al. 2021, MNRAS, 502, L40\n- Foellmi, C., Mo ff at, A. F. J., & Guerrero, M. A. 2003, MNRAS, 338, 1025 Friend, D. B. & Abbott, D. C. 1986, ApJ, 311, 701\n- Fryer, C. L., Olejak, A., & Belczynski, K. 2022, ApJ, 931, 94\n- Gaia Collaboration, Panuzzo, P., Mazeh, T., et al. 2024, arXiv e-prints, arXiv:2404.10486\n- Georgy, C., Ekström, S., Meynet, G., et al. 2012, A&A, 542, A29\n- Georgy, C., Meynet, G., & Maeder, A. 2011, A&A, 527, A52\n- Georgy, C., Saio, H., & Meynet, G. 2014, MNRAS, 439, L6\n- Gormaz-Matamala, A. C., Cuadra, J., Ekström, S., et al. 2024, arXiv e-prints, arXiv:2404.07570\n- Gormaz-Matamala, A. C., Cuadra, J., Meynet, G., & Curé, M. 2023, A&A, 673, A109\n- Gormaz-Matamala, A. C., Curé, M., Cidale, L. S., & Venero, R. O. J. 2019, ApJ, 873, 131\n- Gormaz-Matamala, A. C., Curé, M., Hillier, D. J., et al. 2021, ApJ, 920, 64\n- Gormaz-Matamala, A. C., Curé, M., Lobel, A., et al. 2022a, A&A, 661, A51\n- Gormaz-Matamala, A. C., Curé, M., Meynet, G., et al. 2022b, A&A, 665, A133\n- Gräfener, G. 2021, A&A, 647, A13\n- Gräfener, G. & Hamann, W. R. 2008, A&A, 482, 945\n- Gräfener, G., Owocki, S. P., Grassitelli, L., & Langer, N. 2017, A&A, 608, A34\n- Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56\n- Groh, J. H., Meynet, G., Ekström, S., & Georgy, C. 2014, A&A, 564, A30\n- Hamann, W. R., Gräfener, G., Liermann, A., et al. 2019, A&A, 625, A57\n- Hawcroft, C., Mahy, L., Sana, H., et al. 2024, arXiv e-prints, arXiv:2407.06775\n- Hawcroft, C., Sana, H., Mahy, L., et al. 2021, A&A, 655, A67\n- Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288\n- Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350\n- Higgins, E. R., Sander, A. A. C., Vink, J. S., & Hirschi, R. 2021, MNRAS, 505, 4874\n- Huang, W., Gies, D. R., & McSwain, M. V. 2010, ApJ, 722, 605\n- Humphreys, R. M. & Davidson, K. 1994, PASP, 106, 1025\n- Humphreys, R. M., Weis, K., Davidson, K., & Gordon, M. S. 2016, ApJ, 825, 64\n- Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543\n- Ivanova, N., Justham, S., Chen, X., et al. 2013, A&A Rev., 21, 59  \nJermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15  \n- Josiek, J., Ekström, S., & Sander, A. A. C. 2024, arXiv e-prints, arXiv:2404.14488\n- Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71\n- Krtiˇcka, J. & Kubát, J. 2017, A&A, 606, A31\n- Krtiˇcka, J. & Kubát, J. 2018, A&A, 612, A20\n- Krtiˇcka, J., Kubát, J., & Krtiˇcková, I. 2024, A&A, 681, A29\n- Langer, N. 1997, in Astronomical Society of the Pacific Conference Series, Vol.  \n120, Luminous Blue Variables: Massive Stars in Transition, ed. A. Nota &  \n- H. Lamers, 83\n- Langer, N. 1998, A&A, 329, 551\n- Ledoux, P. 1947, ApJ, 105, 305\n- Lucy, L. B. & Abbott, D. C. 1993, ApJ, 405, 738\n- Lucy, L. B. & Solomon, P. M. 1970, ApJ, 159, 879\n- Maeder, A. 1987, A&A, 178, 159\n- Maeder, A. & Meynet, G. 2000, A&A, 361, 159\n- Martinet, S., Meynet, G., Ekström, S., Georgy, C., & Hirschi, R. 2023, arXiv e-prints, arXiv:2309.00062\n- Martinet, S., Meynet, G., Ekström, S., et al. 2021, A&A, 648, A126  \nMartins, F. 2023, A&A, 680, A22  \n- Martins, F., Hillier, D. J., Paumard, T., et al. 2008, A&A, 478, 219\n- Martins, F. & Palacios, A. 2022, A&A, 659, A163\n- Meynet, G. & Maeder, A. 2000, A&A, 361, 101\n- Meynet, G. & Maeder, A. 2003, A&A, 404, 975\n- Meynet, G., Maeder, A., Schaller, G., Schaerer, D., & Charbonnel, C. 1994, A&AS, 103, 97\n- Muijres, L. E., de Koter, A., Vink, J. S., et al. 2011, A&A, 526, A32\n- Nugis, T. & Lamers, H. J. G. L. M. 2002, A&A, 389, 162\n- Olejak, A., Belczynski, K., & Ivanova, N. 2021, A&A, 651, A100\n- Olejak, A., Fryer, C. L., Belczynski, K., & Baibhav, V. 2022, MNRAS, 516, 2252\n- Pauli, D., Oskinova, L. M., Hamann, W. R., et al. 2022, A&A, 659, A9\n- Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3\n- Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4\n- Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15\n- Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34\n- Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10\n- Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nature Astronomy, 5, 715  \nPuls, J., Vink, J. S., & Najarro, F. 2008, Astronomy and Astrophysics Review,  \n16, 209  \n- Quataert, E., Fernández, R., Kasen, D., Klion, H., & Paxton, B. 2016, MNRAS, 458, 1214\n- Ramírez-Agudelo, O. H., Sana, H., de Koter, A., et al. 2017, A&A, 600, A81\n- Romagnolo, A., Belczynski, K., Klencki, J., et al. 2023, MNRAS, 525, 706\n- Romagnolo, A., Gormaz-Matamala, A. C., & Belczynski, K. 2024a, ApJ, 964, L23\n- Romagnolo, A., Klencki, J., Vigna-Gomez, A., & Belczynski, K. 2024b, arXiv e-prints, arXiv:2410.17315\n- Sabhahit, G. N., Vink, J. S., Higgins, E. R., & Sander, A. A. C. 2022, MNRAS, 514, 3736\n- Sabhahit, G. N., Vink, J. S., Sander, A. A. C., & Higgins, E. R. 2023, MNRAS, 524, 1529\n- Sander, A. A. C., Hamann, W.-R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86\n- Sander, A. A. C., Hamann, W. R., Todt, H., et al. 2019, A&A, 621, A92\n- Sander, A. A. C., Lefever, R. R., Poniatowski, L. G., et al. 2023, A&A, 670, A83 Sander, A. A. C. & Vink, J. S. 2020, MNRAS, 499, 873\n- Sander, A. A. C., Vink, J. S., & Hamann, W. R. 2020, MNRAS, 491, 4406\n- Scott, L. J. A., Hirschi, R., Georgy, C., et al. 2021, MNRAS, 503, 4208\n- Sibony, Y., Georgy, C., Ekström, S., & Meynet, G. 2023, A&A, 680, A101  \nSmith, L. F., Shara, M. M., & Mo  \nff  \nat, A. F. J. 1996, MNRAS, 281, 163  \n- Spera, M. & Mapelli, M. 2017, MNRAS, 470, 4739\n- Spruit, H. C. 2002, A&A, 381, 923  \nSundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126 Tayler, R. J. 1973, MNRAS, 161, 365  \n- Tehrani, K. A., Crowther, P. A., Bestenlehner, J. M., et al. 2019, MNRAS, 484, 2692  \nVenero, R. O. J., Curé, M., Puls, J., et al. 2024, MNRAS, 527, 93  \n- Vink, J. S. 2017, A&A, 607, L8\n- Vink, J. S. 2021, arXiv e-prints, arXiv:2109.08164\n- Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181\n- Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574\n- Vink, J. S. & Gräfener, G. 2012, ApJ, 751, L34\n- Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146\n- Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132\n- Vink, J. S., Sabhahit, G. N., & Higgins, E. R. 2024, arXiv e-prints, arXiv:2407.07204  \nŠurlan, B., Hamann, W. R., Aret, A., et al. 2013, A&A, 559, A130  \nWu, S.-W., Bik, A., Bestenlehner, J. M., et al. 2016, A&A, 589, A16  \nYoon, S. C., Woosley, S. E., & Langer, N. 2010, ApJ, 725, 940  \n- Yungelson, L. R., van den Heuvel, E. P. J., Vink, J. S., Portegies Zwart, S. F., & de Koter, A. 2008, A&A, 477, 223  \nYusof, N., Hirschi, R., Meynet, G., et al. 2013, MNRAS, 433, 1114', 'Appendix A: Effects of TS dynamo for angular momentum transport': 'For our G enec models, we include Tayler-Spruit dynamo (Tayler 1973; Spruit 2002; Heger et al. 2005) for the treatment of the angular momentum transport, besides adopting Ledoux instead of Schwarzschild criterion for the treatment of convective layers (Martinet et al. 2023). The evolution during the MS phase is basically not a ff ected by di ff erences in convection treatment, and just for the H-core burning the Ledoux criterion improves the agreement between models and observations for the surface abundances of BSGs evolving back from RSGs (Georgy et al. 2014). For stars with M zams ≥ 60 M ⊙ (as in this study), the impact of the di ff erent convective criterion is almost negligible because mass loss become the main driver of the evolution (Sibony et al. 2023).  \nIn Fig. A.1 we show the evolution tracks of our models shown in Section 3, compared to models with the same input but adopting the original inner structure of G enec . The initial rotation this time follows the framework Ω = 3 / 3 crit = 0 . 4 (equivalent to ω/ω crit ≃ 0 . 58, see Eq. 8). The evolution of the rotational speed during the main sequence is shown in Fig. A.2, whereas the final parameters between both cases are displayed in Table A.1.  \nAs expected, the di ff erences between adopting the new or the old inner structure are minimal, and are reduced to braking in the stellar rotation since the ZAMS to the TAMS. TS dynamo produces a more e ffi cient transport of angular momentum from the core to the surface, and therefore the stellar surface will rotate faster than in the case without TS. This impact is stronger for the 60 M ⊙ and almost negligible for the 200 M ⊙ models, because of the larger convective core of VMS.  \nRegardless of the HRD tracks, TS dynamo actually a ff ects the stellar inner structure. Fig. A.3 shows how the specific angular momentum jr = ω iR 2 at the end of the He-core burning stage is more depleted when TS dynamo is incorporated, due to the more e ffi cient angular momentum transport. Again, the influence is more prominent for the case with the smaller convective core. Even though, adoption of TS dynamo have an impact on the expected spin of the resulting black hole (Belczynski et al. 2020b), for the analysis of the evolution of stars with M zams ≥ 60 M ⊙ its impact is not relevant enough.  \nTable A.1: Table of final parameters for our G enec evolution models, comparing the new inner structure (Tayler-Spruit dynamo and Ledoux convective criterion) with the old inner structure (no-TS dynamo and Schwarzchild convective criterion).  \n<!-- image -->  \nFig. A.1: HR diagram for our G enec evolution models, comparing new inner structure (Tayler-Spruit dynamo and Ledoux convective criterion, red solid lines) with old inner structure (no-TS dynamo and Schwarzchild convective criterion, dashed brown lines).  \n<!-- image -->  \n<!-- image -->  \nFig. A.2: Rotational velocity during the MS for our G enec evolution models, comparing new inner structure (Tayler-Spruit dynamo and Ledoux convective criterion, red solid lines) with old inner structure (no-TS dynamo and Schwarzchild convective criterion, dashed brown lines).  \n<!-- image -->  \nGormaz-Matamala et al.: Evolution of stars with 60 and 200 M ⊙ : predictions for WNh stars in the Milky Way.  \n<!-- image -->  \nFig. A.3: Specific angular momentum at the end of the He-core burning for our G enec evolution models, comparing new inner structure (TaylerSpruit dynamo and Ledoux convective criterion, blue) with old inner structure (no-TS dynamo and Schwarzchild convective criterion, orange).  \n<!-- image -->'}