Evolution of rotating massive stars with new hydrodynamic wind models

Context. Mass loss due to radiatively line-driven winds is central to our understanding of the evolution of massive stars in both single and multiple systems. This mass loss plays a key role in modulating the stellar evolution at di ﬀ erent metallicities, particularly in the case of massive stars with M ∗ ≥ 25 M (cid:12) . Aims. We extend the evolution models introduced in Paper I, where the mass-loss recipe is based on the simultaneous calculation of the wind hydrodynamics and the line acceleration, by incorporating the e ﬀ ects of stellar rotation. Methods. As in Paper I, we introduce a grid of self-consistent line-force parameters ( k ,α,δ ) for a set of standard evolutionary tracks using G enec . Based on this grid, we analysed the e ﬀ ects of stellar rotation, CNO abundances, and He / H ratio on the wind solutions to derive additional terms for the recipe with which we predict the self-consistent mass-loss rate, ˙ M sc . With this, we generated a new set of evolutionary tracks with rotation for M ZAMS = 25, 40, 70, and 120 M (cid:12) , and for metallicities Z = 0 . 014 (Galactic) and 0.006 (Large Magellanic Cloud). Results. In addition to the expected correction factor due to rotation, the mass-loss rate decreases when the surface becomes more helium rich, especially in the later moments of the main-sequence phase. The self-consistent approach gives lower mass-loss rates than the standard values adopted in previous G enec evolution models. This decrease strongly a ﬀ ects the tracks of the most massive models. Weaker winds allow the star to retain more mass, but also more angular momentum. As a consequence, weaker wind models rotate faster and show a less e ﬃ cient mixing in their inner stellar structure at a given age. Conclusions. The self-consistent tracks predict an evolution of the rotational velocities through the main sequence that closely agrees with the range of v sin i values found by recent surveys of Galactic O-type stars. As subsequent implications, the weaker winds from self-consistent models also suggest a reduction of the contribution of the isotope 26 Al to the interstellar medium due to stellar winds of massive stars during the MS phase. Moreover, the higher luminosities found for the self-consistent evolutionary models suggest that some populations of massive stars might be less massive than previously thought, as in the case of Ofpe stars at the Galactic centre. Therefore, this study opens a wide range of consequences for further research based on the evolution of massive stars.


Introduction
Evolution of massive stars is an important topic in stellar astrophysics.Indeed some of them are the progenitors of corecollapse supernova events, and give birth to neutron stars and black holes (Heger et al. 2003).Massive stars are also important for the study of nucleosynthesis, production of ionising flux, feedback due to wind momentum, studies of star formation history, and galaxy evolution.
Massive stars with M * 25 M are characterised by strong stellar winds and large mass loss rates ( Ṁ), which determine the evolution of the star such as the phases involving Luminous Blue Variables (LBV) and Wolf-Rayet (WR) stars (Groh et al. 2014).The wind properties through all these stages constrain which kind of core-collapse process the star will experience at the end of its life, and what will be the final mass of the remnant (Belczynski et al. 2010).This last item is relevant, because theoretical masses of stellar black holes need to be in agreement with those measured by gravitational waves, in merger events detected by the Advanced-LIGO interferometers (Abbott et al. 2016).The study of the stellar winds has further implications that reach far beyond stellar astrophysics.For instance, the stellar winds from numerous evolved massive stars collide to produce the plasma that fills up the interstellar medium in the Galactic Centre (e.g., Cuadra et al. 2008Cuadra et al. , 2015;;Ressler et al. 2018Ressler et al. , 2020;;Calderón et al. 2020b).The properties of the stellar winds determine whether the plasma is able to cool and form clumps, which affects the accretion process on to the Galactic supermassive black hole, Sgr A* (Cuadra et al. 2005;Calderón et al. 2016Calderón et al. , 2020a)).Therefore, a proper determination of the stellar winds is needed to correctly interpret the images of the black hole silhouette obtained by the Event Horizon Telescope Collaboration et al. (2022).These examples highlight that theoretical calculations for the stellar winds can even influence the field of relativistic astrophysics.
In the last decade, computational codes such as Mesa (Paxton et al. 2019) and the Geneva evolution code (Eggenberger et al. 2008, hereafter Genec) have been used to study stellar evolution.The recipe for theoretical mass loss rate varies, depending on the stellar evolution stage and the region of the HR dia-gram (HRD).For O-type main sequence (MS) stars, the recipe most commonly adopted has been the so-called Vink's formula (Vink et al. 2001).However, diagnostics of mass loss rates performed in recent years consider that values from Vink's formula are overestimated by a factor of ∼ 3 (Bouret et al. 2012;Šurlan et al. 2013;Vink 2021), and therefore the quest has been the development of evolution models adopting more updated and accurate recipes for mass loss rate.In that direction, we mention the studies from Björklund et al. (2022) using Mesa, and Gormaz-Matamala et al. (2022b, hereafter Paper I) using Genec.
In Paper I, we developed new evolution models for stars born with masses M zams ≥ 25 M , introducing a new recipe for the theoretical mass loss rate based on the self-consistent m-CAK wind solutions from Gormaz-Matamala et al. (2019, 2022a).These new self-consistent tracks show that stars retain more mass through their evolution, therefore remaining larger and more luminous compared with models for massive stars using Vink's formula for Ṁ (Ekström et al. 2012;Georgy et al. 2013;Eggenberger et al. 2021).A similar result is obtained by Björklund et al. (2022), where evolution models adopting new mass loss rates are shifted to higher luminosities in HRD.It is important to remark that, to our knowledge, both Björklund et al. (2022) and Gormaz-Matamala et al. (2022b) are the first studies on developing tracks for stellar evolution of massive stars adopting a self-consistent treatment for the stellar wind.
In this work, we extend the self-consistent evolutionary tracks introduced in Paper I, including rotation to our models.Rotation is important, because it modifies the mass loss rate and the surface abundances due to the rotational mixing (Maeder & Meynet 2010), and because it makes the star not only mass but also losing angular momentum through the stellar wind (Georgy et al. 2011;Keszthelyi et al. 2017).Because of the extra details to be considered into the analysis of our results, this work includes only Galactic (Z = 0.014) and LMC (Z = 0.006) metallicities, leaving the supra-solar (Z > 0.014) and SMC and lowermetallicity (Z ≤ 0.002) cases to forthcoming studies.
This paper is organised as follows.Section 2 summarises the most relevant results and conclusions obtained from our self-consistent evolutionary tracks in Paper I. Section 3 introduces the updates on Ṁsc to account for the effects of rotation.Section 4 presents the new evolution models adopting selfconsistent winds, whereas their comparison with observational diagnostics are analysed in Section 5.In Section 6, we deliberate about the implications of these new evolution models in a Galactic scale.Finally, conclusions of this work are summarised in Section 7.

Self-consistent evolution models
We call self-consistent evolutionary tracks the evolution models that adopt a self-consistent recipe for the mass loss rate.Selfconsistency implies that the wind hydrodynamics (i.e., velocity and density profiles for the wind) are calculated in agreement with the radiative acceleration (i.e., the acceleration over the wind particles due to the stellar radiation).As a consequence, mass loss rate is a direct solution for the wind under a specific set of stellar parameters (temperature, radius, mass, abundances) without assuming any extra condition a priori for the wind (such as a velocity law or a ∞ / esc ratio).For that purpose, we simultaneously calculate the line-acceleration g line by means of the socalled line-force parameters (k, α, δ) from m-CAK theory (Castor et al. 1975;Abbott 1982;Pauldrach et al. 1986), whereas we calculate the hydrodynamics by solving the stationary equation of motion for the wind using the code HydWind (Curé 2004), using ṀVink and Z/Z = 1.0 covering the main sequence stage, hereafter original tracks.Circles represent the location of the selected stellar models to be tabulated in Table 1 of Section 3.
for different sets of stellar parameters for massive stars.This iterative procedure is hereafter named m-CAK prescription, it was introduced in Gormaz- Matamala et al. (2019) and it has demonstrated to provide values for Ṁ of the same order of magnitude than those obtained with self-consistent NLTE studies in the comoving frame such as Krtička & Kubát (2017, 2018) or Gormaz-Matamala et al. (2021).Moreover, the wind hydrodynamics used as input in the NLTE radiative transport code is able to reproduce reliable synthetic spectra, in agreement with the observations (Gormaz-Matamala et al. 2022a).
However, although the calculation of the line-force parameters has been extended for even cooler temperatures (Lattimer & Cranmer 2021;Poniatowski et al. 2022), in Gormaz-Matamala et al. (2019) our m-CAK prescription is restricted into the range of stellar parameters where (k, α, δ) can be assumed as constants.This is an important consideration because, even though lineforce parameters have been commonly treated as constants (Shimada et al. 1994;Noebauer & Sim 2015), there are studies which suggest that the values for k, α and δ may radially change with the wind (Schaerer & Schmutz 1994;Puls et al. 2000).About this, Gormaz-Matamala et al. (2022a, Sec. 2) performed an quantitative analysis on the standard deviation and numerical errors in the calculation of our line-force parameters, and by extension on the error margins for the self-consistent wind parameters Ṁ and ∞ , finding that uncertainties for (k, α, δ) can be neglected for stars with T eff ≥ 30 kK and log g ≥ 3.2.Therefore, we set these thresholds as the range of validity for our m-CAK prescription.
Thus, self-consistent mass loss rate ( Ṁsc ) derived from the m-CAK prescription is parametrised as a function of stellar parameters as shown in Eq. 7 from Paper I: log Ṁsc,predicted = − 40.
Table 1: Self-consistent line-force parameters (k, α, δ) calculated for the stellar parameters corresponding to the positions indicated by circles in Fig. 1, together with their resulting terminal velocities and mass loss rates ( Ṁsc ).Ratios showing the enhancement of mass loss rate due to rotation, adopting either our self-consistent procedure or MM00, are also shown in the last columns.where w, x, y, and z are defined as:

Input stellar parameters
This intervariable fitting, besides minimising the separation with respect to the real calculated wind solution by line-force parameters, allows a deeper examination of the dependance of mass loss rate on stellar parameters.As an example, we found from Eq. 1 that dependency in metallicity does not follow the classical exponential law Ṁ ∝ Z a with a constant exponent, but with a ranging from ∼ 0.53 for M * ∼ 120 M to ∼ 1.02 for M * ∼ 25 M .In other words, we can approximate this Z-dependence for mass loss rate as an intrinsic relationship of stellar mass, as Ṁsc ∝ Z 0.4+ 15.75  (M * /M ) . (2) An analogous result, but for an intrinsic luminositydependence, was found by Krtička & Kubát (2018) for their global wind models; where the exponent for the metallicity was less steep for higher stellar luminosity (and subsequently, stellar mass).Explanations of this intrinsic mass-dependence may rely on the fact that the more massive the star, the nearer it is from the Eddington limit.In that case, the continuum should become increasingly important in contributing to the radiative acceleration g rad , in decline of the line-driving (see Eq. 6 from Gormaz-Matamala et al. 2021).The effect of the continuum is less dependent on metallicity because it does not involve interaction between the radiative field and absorption lines.This condition for the mass loss rate was previously pointed out by Gräfener & Hamann (2008) for WNL stars, where the Z-dependence was weaker closer to the Eddington limit (i.e., more mass).All in all we can state that Eq. 2 is an important step to understand the dependence on metallicity for the mass loss rates at the earlier evolutionary stages of massive stars.
Our new recipe for mass loss rate is implemented inside the Genec code, which runs stellar evolutionary models to explain and predict the physical properties of massive stars (Eggenberger et al. 2008).We used the same prescriptions as described in Ekström et al. (2012); that is, solar abundances from Asplund et al. (2005Asplund et al. ( , 2009)), opacities from Iglesias & Rogers (1996), and overshoot parameter α ov 0.10, except for the mass loss rates.The self-consistent m-CAK prescription has been already incorporated inside Genec for stars satisfying T eff ≥ 30 kK and log g ≥ 3.2.Below such thresholds, the recipe for mass loss rate reverts to Vink's formula.
Because self-consistent mass loss rates are a factor of ∼ 3 lower with respect to Vink's formula at the beginning of MS (i.e., for T eff 40 kK and log g ∼ 4.0, see Gormaz-Matamala et al. 2022a), self-consistent evolution models predict that stars will retain more mass during their H-core burning stage and therefore they are larger and more luminous.These changes are proportional to the initial mass of the star; from important differences in the final stellar properties at the end of MS for M zams = 120 M , to almost negligible differences for M zams = 25 M .Also, evolution models adopting Ṁsc predict a smaller production of the isotope Aluminium-26, which means that massive stars contribute in a lesser proportion to feed the interstellar medium with this isotope, compared to other sources (Palacios et al. 2005;Wang et al. 2009).
Such results are important, but still represent only a first glance to the implications that the self-consistent wind models have for stellar evolution.The next step is the incorporation of rotation.Stellar rotation is important for stellar evolution because of two reasons.First, because the angular momentum in the core at the time of core collapse may strongly impact the final stages of a massive star (Yoon et al. 2006).Powerful stellar winds not only make stars lose mass, but also angular momentum, subsequently affecting the rotation and the evolution of the star (Georgy et al. 2011;Keszthelyi et al. 2017).And second, because rotational mixing (Maeder & Meynet 2010) changes the distribution of the abundances of the elements inside the star and in particular at the surface (see, e.g., Brott et al. 2011;Ekström et al. 2012).The changes of the abundances at the surface is particularly interesting, because modification in the chemical composition of the atmosphere induces variations in the lineacceleration and therefore in the self-consistent mass loss rate, something that it has been either not considered or neglected by previous authors (Björklund et al. 2022).For that reason, before running new evolution models with rotation adopting Ṁsc , it is necessary to study the influence of rotation over the m-CAK prescription.

Enhancing mass loss due to rotational effects
Models by Genec are computed from the formula by Vink (Vink et al. 2001) with the correcting factor induced by rotation as given by Maeder & Meynet (2000, hereafter MM00), that is where Ṁ(0) is the mass loss rate without rotation, α is the lineforce parameter from CAK theory, ω is the angular velocity at the stellar surface, and Γ Edd is the Eddington factor with κ being the total opacity.Notice that Eq. 3 introduces a correction factor independent of the adopted recipe for mass loss rate without rotation; in other words, Ṁ(0) can be computed either from Vink's formula or from Eq. 1. Actually, Ṁ(ω) results from an integration over the surface and the integration accounts for the change of the local mass loss rate with the latitude.Equation 3 incorporates the impacts on mass loss rate produced not only by rotation but also the radiative acceleration due to the continuum (hence the presence of the Eddington factor), through the so-called ΩΓ-limit.This formulation implies that the stellar break-up or Ω-limit (i.e., when the centrifugal forces compensate gravity in a rotating star) is reached for reduced rotation velocities if Γ is large enough.The velocity associated to this break-up is the so called critical rotational speed, and hence, the rotation of the star is commonly expressed as a dimensionless ratio between the rotational velocity and this critical value as On the other hand, the m-CAK prescription incorporates rotational effects in the hydrodynamic solution calculated by Hyd-Wind (Curé 2004), also using Ω as an input.Hydwind solves the equation of motion for the stationary standard line-driven theory, by finding the eigenfunction that satisfies the velocity profile, = (r), and the eigenvalue, which is proportional to the mass loss rate, Ṁ (Curé & Rial 2007).The rotation factor Ω is implicitly included in Eq. 7 by means of φ = rot R * /r, whereas the additional terms are the same as described by Araya et al. (2017Araya et al. ( , 2021) ) or Gormaz-Matamala et al. (2022a).In this work, we use moderate to intermediate values of Ω (= 0.4), therefore all our wind solutions are in the fast wind regime and not in the Ω-slow wind regime which starts approximately from Ω 0.75 (Araya et al. 2018).
In order to evaluate the effect of the rotation over the selfconsistent wind solutions, we repeat the methodology of Paper I and we run evolutionary rotating models in Genec for initial masses of 120, 70, 40, and 25 M adopting Vink's formula for these tracks, hereafter original tracks, as shown in Fig. 1.Then we select a set of representative points over these original tracks, and based on their stellar parameters (T eff , log g, R * /R , Ω) we calculate the line-force parameters (k, α, δ) and their respective self-consistent wind parameters Ṁsc and ∞,sc .Initial rotational velocity for all evolution models is Ω = 0.4, as in Ekström et al. (2012).The results are tabulated in Table 1 where, besides the former columns already included in the Tables 2, 3 and 4 from Paper I, we add the rotation factor Ω as an extra stellar parameter for the input.For each one of the points, we calculate an independent solution for (k, α, δ) for the tabulated Ω but also for Ω = 0.0 (i.e., without rotation), with their respective self-consistent solution for terminal velocity and mass loss rate.The difference in the self-consistent mass loss due to the incorporation of rotation is expressed as the ratio Ṁ(Ω)/ Ṁ(0), which is compared with the correction factor predicted by MM00 (i.e., Eq. 3).Here, we found that both approaches for the rotational effects over Ṁ produce similar results, which is remarkable considering the error bars associated with the self-consistent mass loss rate.Based on this outcome, we decide to perform our new evolutionary tracks keeping Eq. 3 as correction factor for Ṁsc .
The other physical ingredients related to rotation are the same as in Ekström et al. (2012), Georgy et al. (2013) and Eggenberger et al. (2021).Angular momentum is transported by an advective equation as described in Zahn (1992).The horizontal and vertical shear diffusion coefficients are from respectively Zahn (1992) and Maeder (1997).

Variation in abundances
Besides the stellar parameters such as temperature, surface gravity, radius and rotation, we included in Table 1 the values for the abundances of helium (expressed as a fraction of the hydrogen abundance by number) and the individual abundances of carbon, nitrogen and oxygen (expressed as a fraction of the solar abundances).Due to the nature of the CNO-cycle, the abundances Self-consistent tracks are marked with two black dots which represent the zero-age main sequence and the switch of the mass loss recipe from self-consistent to Vink's formula (when either T eff = 30 kK or log g = 3.2), whereas the end of the main sequence (terminal-age main-sequence, TAMS) is represented with a black cross.
of helium and nitrogen increase over time, whereas carbon and oxygen decrease (Przybilla et al. 2010).The chemical composition of the surface is enriched in material processed in the core by rotational mixing in the outer radiative envelope.Since we assume that chemical composition in the stellar wind is exactly the same as in the photosphere, the chemical composition of the stellar wind is altered.The question that arises here is, do these alterations in the wind chemical abundances affect the lineacceleration and afterwards the self-consistent mass loss rate?
To solve that question, we artificially alter the individual abundances for the He and the CNO elements from our stellar models of Table 1, and then calculate the respective line-force with Ṁsc the mass losses tabulated in Table 1, and Ṁsc, the selfconsistent mass loss rate calculated if the abundance of the evaluated element were solar.
Results of these modifications over the Ṁsc are displayed in Fig. 2a for changes in the He/H ratio, and in Fig. 2b for modifications of the CNO elements.We see that the variation of each one of the CNO elements even by a factor of ∼ 10 barely affects the Ṁsc , keeping the difference ∆ log Ṁsc close to zero.Indeed, differences in mass loss rate are just of the order of ±0.05 in log scale (∼ 12% of Ṁsc ), smaller than the error bars shown in Table 1.On the contrary, the increase in the He/H ratio makes Ṁsc decrease in larger magnitudes.
Results from Fig. 2 are quickly explained due to their absolute numbers of the evaluated abundances.Even though the changes in CNO abundances are by a factor of ∼ 10 (compared with the initial values at ZAMS), they still represent a small percentage of the total composition of the star.On the contrary, the abundance of helium plays a major role because it grows (in mass fraction) from ∼ 26% up to ∼ 70%, thus significantly changing the mean molecular weight of the wind.Also, helium abundance is involved in the computation of the ratio N e /ρ and in the radiative acceleration due to Thomson scattering 1 where X He is the helium abundance by number, N e is the density of free electrons, and m p is the proton mass (the meaning of the other symbols are detailed in Curé 2004, Sec. 2).Hence, when the helium abundance increases, the number of free electrons is reduced and then Γ decreases.As a consequence, the effective mass M eff = M * (1 − Γ) becomes larger in Eq. 7, and accordingly effective gravity, meaning that the line-acceleration now is able to remove less material as stellar wind (i.e., mass loss rate decreases).However, it is important to mention that this theoretical analysis is made upon the basis of the isolated alteration of the He/H ratio, and therefore this decreasing in the resulting Ṁ does not imply any modifications in temperature or in any other stellar parameter.For this reason, it differs from the Ṁ ∝ X He relationship found for Hamann et al. (1995), which was the result of a global analysis over the wind regime of WR stars.Thus we decide to introduce an extra variable to the recipe of Eq. 1, based on a simple fit of Fig. 2a log Ṁsc = log Ṁsc,0 − 0.559 × (He/H) * − (He/H) 0 , where (He/H) 0 is the He/H ratio at ZAMS (0.085 from Asplund et al. 2009).We see from Table 1 that (He/H) * may grow up to ∼ 0.3 − 0.45 at the end of the self-consistent regime, thus leading to a reduction of the mass loss rate of around ∼ 30−40%.However, same as in Paper I, this modification applies only when the self-consistent mass loss rate is adopted in the evolutionary tracks, for T eff ≥ 30 kK and log g ≥ 3.2, and hence it is still unknown the magnitude of the decreasing factor in Ṁ prior to the end of the MS, when hydrogen abundance is X = 0.3 at the surface.
Finally, details of the formulae used for both old and new evolution models are summarised in Table 2.

Results
Evolutionary models with rotation, both original and selfconsistent, are shown in Fig. 3 in the HR diagram, whereas the rotation and equatorial velocities for tracks adopting both mass loss recipes are shown in Fig. 4. Additional plots such as the comparison with self-consistent non-rotating models from Paper I, the evolutionary tracks in the log g vs. T eff plane, evolution of the mass loss rate, stellar mass and radii, He/H ratio, Eddington factor, and convective core mass are shown in Appendix A. Also, hereafter all plots comparing new and old recipes of mass loss rate will use solid lines for models with Ṁsc and dashed lines for models with ṀVink .Major features of the stellar models at the end of the MS, comparing self-consistent tracks with the tracks adopting Vink's formula (Ekström et al. 2012;Eggenberger et al. 2021), are tabulated in Table 3.
Overall, compared with Paper I (see Table 8 from Gormaz-Matamala et al. 2022b) we find that the lifetime during the main sequence is ∼ 18 − 22% longer than for the non-rotating models, regardless of the adopted mass loss recipe.The difference with 1 Here, Γ is the Eddington factor when only free electron scattering opacity is accounted for (the so called classical Eddington factor).
non-rotating models is also appreciated in when Ṁsc is adopted.Moreover, besides the straightforward differences in the final masses and the final radii, we also observe that the less strong winds from the self-consistent prescription predict in general weaker surface enrichments of helium and nitrogen at the surface, as seen in the last columns of Table 3 and in the evolution of the He/H ratio (Fig. A.6).The only noticeable exception is the model with M zams = 25 M and Z = 0.006, described in detail in the subsection below.
Concerning rotation, self-consistent evolutionary tracks have higher velocities for all models as appreciated in Fig. 4, where we compute equatorial velocities and Ω = eq / crit as a function of the core hydrogen abundance from ZAMS (starting with X core ∼ 0.72) until the end of the main sequence.Same as in Paper I, the largest impact of the new self-consistent mass loss rates are for the higher initial mass stars: the weaker Ṁsc produces stars that rotate faster.Even though the rotational velocity drops down at TAMS for all models, this braking is slower than the one expected for models with the old mass loss recipe.This weaker braking is more prominent for the low metallicity, where the stars adopting Ṁsc keep their rotational velocity almost constant for a longer time before the final deceleration.This implies that, in our regime of metallicities (between Z = 0.006 and Z = 0.014) mass loss is still a prominent process that impacts the evolution of the rotational velocity.Evolution of rotational velocities is explained by the balance between the transport of angular momentum from the convective core (which is contracting during the main sequence stage) to the stellar surface, and the removal of outer layers due to mass loss (Ekström et al. 2008;de Mink et al. 2013).At extreme low metallicities (Z = 0.0002 = 1/35 Z ), massive stars increase their rot because their mass losses are not strong enough to counteract the transport of angular momentum from the contracting core (Szécsi et al. 2015).Also for Z = 0.006, we see that Ω, which is set initially as = 0.4, increases up to a constant peak of ∼ 0.55 before the final decrease, as a direct consequence of the initial constancy in the equatorial velocity and the decreasing in the evolution of the critical velocity (Eq.5).Weaker winds predicts a nearly exponential increasing in the stellar radii for all initial masses and both metallicities (Fig. A.5), whereas both the stellar mass and Eddington factor show a linear and metallicity dependent difference with respect to the models computed with the Vink's formula (Fig. A.4 and Fig. A.7).Therefore, because they are weaker, self-consistent winds makes our evolution models to approach closer to the break-up limit for rotational velocity in agreement with Ekström et al. (2008), prior to the final drop in rot at the end of main sequence in agreement with Szécsi et al. (2015).Now, let us discuss some more detailed effects of the use of the new mass loss rate for the different initial mass models considered in the present work.

Case M zams = 120 M
This is the range where we find the very massive stars (VMS Yusof et al. 2013), whose properties and evolution have aroused great interest in recent years (Higgins et al. 2022;Sabhahit et al. 2022).Here, the tracks with Z = 0.006 have a special relevance  because the most studied cluster containing VMSs is 30 Dor, located in the LMC (Crowther et al. 2010;Bestenlehner et al. 2014).Same as in Yusof et al. (2013), we consider that the star becomes a Wolf-Rayet (i.e., it changes from optically thin to optically thick wind) when its surface hydrogen abundance goes below 0.3 and T eff > 10 kK.After that point, mass loss rate is assumed to follow the recipe from Gräfener & Hamann (2008) for WR winds, generating the discontinuity in Ṁ after the black cross as observed at Fig. A.3.Notice that prior to the change in the recipe, evolution of the self-consistent mass loss rate proceeds as expected due to the large increment in the He/H ratio (Fig. A.6) as a consequence of Eq. 11.Due to the extreme strength of the winds, the VMS is depleted from the 70% of its surface hydrogen even before H is completely exhausted in its core, as seen from Fig. 6a, generating the stars with spectral type WNh (Hamann et al. 2006;Martins & Palacios 2022).This switch not only in the Ṁ prescription but also in the wind regime (from optically thin to optically thick wind) is also evidenced by Fig. A.5, where the radius needs to the recalculated after the black crosses based on the wind opacity (see, e.g, Meynet & Maeder 2005) due to the absence of a well-defined photosphere (Schaller et al. 1992).
The most remarkable difference between classical and selfconsistent tracks in Fig. 3, is the drift 'redwards' (i.e., towards the cool range of temperatures in the HRD) that the selfconsistent models do since the beginning of the main sequence.For both metallicities, stars with Ṁsc end their MS lifetime as O-type at T eff ∼ 40 kK, before suddenly increasing their temperature and becoming WNh stars.We note that the model using the weaker self-consistent Ṁ shows in general a slightly stronger surface enrichment than the original model.This explains why these models evolve to redder regions of the HRD during the MS phase.It is a well known effect that a star showing a nearly homogeneous internal chemical composition keep blue colours in the HRD diagram (Maeder 1987).A homogeneous chemical composition can result from a very strong internal mixing, from very strong mass losses or from both processes.In contrast an internal non-homogeneous composition will produce redwards evolution on the HRD.The extension to the red is favoured when a larger convective core is present and the radiative envelope is not too efficiently mixed, otherwise we would have a situation approaching the one of a homogeneous internal composition and the track will remain in the blue region of the HRD (Martinet et al. 2021).Here we see that reducing the mass loss favours redwards tracks.In Fig. 5, we show the diffusion coefficients in the new and old mass loss rate models in the middle of the MS phase.For the case of 120 M and Z = 0.014, the model with a lower mass loss rate has a more massive convective core at that stage, even though it represents a slightly minor fraction of the total mass: ∼ 84% for the self-consistent model, against ∼ 87% for the old one.This is a result of the higher difference in the mass retained and also to the possible higher value of the effective diffusion coefficient, D eff , just above the convective core.In The meaning of the different line style is the same as in Fig. 3. Vertical grey lines represent the total mass of both models, which correspond to the final masses at the end of MS tabulated in Table 3.
the rest of the radiative envelope, there is not much difference between the values of D eff .Since the diffusion timescale varies as R 2 /D eff , and since the radius of the self-consistent model is larger, we see that the self-consistent model will take a longer time to reach a given level of surface enrichment than the original one.This imposes redder evolutionary tracks for the model with Ṁsc .Figure 6a compares the chemical composition of the selfconsistent and original models at the end of the MS phase.As expected the model with the lower mass loss rate is showing the greatest contrasts between the surface composition and that of the core.Figure 7a shows the variation of the specific angular momentum j r = r 2 ω inside the 120 M at the beginning and end of the MS phase for the two different mass loss prescriptions.Let us first remind that in the absence of any angular momentum transport, j r remains constant.The evolution that we see between the beginning and end of the core H-burning phase (a global decrease of j r in the whole star) is due to the fact that the transport processes will in general transport angular momentum from the inner to the outer layers of the star.At the surface the angular momentum is removed by stellar winds, thus we see that when the mass loss rates are smaller, j r is larger.Indeed, by decreasing the mass loss rate by about a factor three during the main sequence (compared to Vink's formula), we obtain that the specific angular momentum is shifted by 0.3 dex upwards when Ṁsc is used, whereas the value of the angular velocity of the core is increased by a factor of ∼ 2 (see Fig. 7b).
Finally, we notice that for both metallicities our rotating models reach the WNh stage at T eff 40 kK.Since we know that stars will increase their temperature through the WR stage, this means that rotating models with M zams = 120 M do not cross the temperature regime around ∼ 25 kK during their MS phase, where the so-called bi-stability jump (Vink et al. 1999) takes place.Because of this reason, and contrary to the non-rotating models from Paper I, we do not observe LBV behaviour such as eruptive mass losses at this mass range.The fact that only initially slow rotating models can go through an eruptive mass loss regime at the end of the main sequence supports the idea that the LBV phenomenon is very rare for very massive stars (Gräfener 2021).

Case M zams = 70 M
Same as before, besides the differences in the stellar mass (Fig. A.4) we also observe a divergence bluewards-redwards for the evolution across the HR diagram (Fig. 3).However, unlike the previous case, this time the beginning of the WR winds stage coincides with the end of the H-core burning, which can be appreciated in Fig. 6b where hydrogen is depleted in the core at the end of both (standard and self-consistent) tracks.Again we see the change in the slope of the evolution in mass loss rate (Fig. A.3), this time with a flattening in Ṁ prior to reaching the second black dot, associated with the increasing of the He/H ratio (Fig. A.6).
Concerning rotation, the braking in the equatorial velocity only becomes significant after the star has burnt almost half of the hydrogen in its core: X c ∼ 0.38 for Z = 0.014 and X c ∼ 0.36 for Z = 0.006, which also represents the peak of Ω ∼ 0.42 for Galactic metallicities and Ω ∼ 0.54 for SMC metallicities (Fig. 4).This is interconnected one more time with the retaining of more specific angular momentum for the models adopting Ṁsc , and with the larger stellar radii (Fig. A.5). Indeed, the switch in the wind regime from Ṁsc to ṀVink (marked with an abrupt increase in mass loss rate after the second black dot of Fig . A.3) carries an important drop in the rotational velocity which creates some numerical instability, but does not affect the general trend in the evolution of observed for all the rotation models.
From Fig. 6b, we see that the variation of the element abundances across the stellar structure from core to surface is more prominent for the model with Ṁsc , which is evidence for the less efficient rotational mixing for the weaker wind.However, on the contrary of the case with M zams = 120 M , where the inner structure was almost fully homogeneous, now there are remarkable distinctions between core and surface abundances because the mass loss has been less efficient in reducing the mass of the envelope and revealing the inner layers.
The redward drift of the evolutionary tracks for the reduced mass loss rate models is in line with the results found by Björklund et al. (2022), who implemented a mass loss recipe derived from the wind theoretical prescription from Björklund et al. (2021) where Ṁ is also ∼ 3 times below the values coming from Vink's formula.In their study, they used Mesa to trace the evolutionary track of a star with M zams = 60 M , rot,ini = 350 km s −1 , and solar metallicity; where they found the drift bluewards at T eff ∼ 26 kK (see their Fig. 6).On our side, by making an eye inspection on Fig. 3 and Fig. A.1, we conjecture that the point of turning bluewards for one of our self-consistent models with 60 M and Z = 0.014 must be between ∼ 29−31 kK, i.e., around ∼ 5 kK hotter than found by Björklund et al. (2022).Such a discrepancy is relatively minor considering the differences between Genec and Mesa (check the comparison for the models of these both codes, as done by Agrawal et al. 2022).Therefore, we can conclude that the self-consistent m-CAK prescription for stellar winds predicts evolutionary tracks in fair agreement with the recent study of Björklund et al. (2022), as expected given that both studies are based on mass loss rates of the same order of magnitude.

Cases M zams = 25 and 40 M
Effects of self-consistent mass loss rates over the evolution of stars within these mass ranges are less pronounced, as it is seen in Fig. 3 and in the resulting final masses and radii from Table 3.We also observe a 'redder' evolution for the models adopting Ṁsc , but still 'bluer' compared with non rotating cases (Fig The composition in the radiative envelope shows some significant differences between the two models.This is particularly visible for the 40 M star, where in the low mass loss rate model a convective region is associated to the H-burning shell (between about 17 and 25 M coordinates), while no convective zone is associated to that shell in the high mass loss rate model.We see therefore, that reducing the mass of the envelope disfavours the formation of an intermediate convective shell associated to the H-burning shell.
The specific angular momentum inside our model of 25 M is shown in Fig. 9a.We see along the curves small abrupt variations (as for instance around 9 M at the end of the MS phase).These come from regions where there is a transition between a convective and a radiative zone, where the chemical composition changes, also observed as an abrupt jump in the inner angular velocity (Fig. 9b).Like for the 120 M star, the model computed with the self-consistent mass loss rates ends the MS phase with slightly larger values of the specific angular momentum.The increase is less strong since the mass loss rates are anyway weaker when the initial mass decreases.The changes in the interior angular velocity distribution are also much weaker than for the 120 M model.We note as expected a slightly faster rotating core in the self-consistent stellar model compared to the original one.1. × 10 -5 5. × 10 -5 Fig. 9: Same as Fig. 7, but for our model with M zams = 25 M and Z = 0.014.

Comparison with rotational surveys
Because the self-consistent m-CAK prescription predicts less strong winds than previous studies, stellar models computed with such prescription evolve at higher luminosities and reach larger radii during the MS phase than the models of Ekström et al. (2012) or Eggenberger et al. (2021), as already detailed in Section 4. Besides, evolutionary tracks adopting Ṁsc also predict that the rotational velocity during the main sequence will be higher as seen in Fig. 4, implying a weaker braking in the equatorial velocity of massive stars.This weaker braking is also represented in Fig. 10, where we plot rotational velocities as a function of the effective temperatures for evolutionary tracks adopting old and new winds, and where the timesteps of 0.5 Myr are illustrated with the respective coloured circles.Veloc-ities from evolution models with ṀVink quickly decreases after the first ∼ 1−2 Myr, passing from little variations in temperature for M zams = 120 M to moderate variations for M zams = 25 M .On the contrary, the rot from self-consistent evolutionary tracks keep relatively constant at the beginning of the main sequence for all tracks, at timescales inversely proportional to their initial masses (from ∼ 1.0 Myr for M zams = 120 M to ∼ 6.0 Myr for M zams = 25 M ), whereas T eff decreases ∼ 8 kK per each model prior to the final braking.
In parallel, we plot in Fig. 10 the sample of 285 O-type stars taken from the recent survey of Holgado et al. (2022), whose catalog is available online2 .In their study, they revisited the ro-  tational properties of 285 Galactic massive O-type as part of the IACOB Project (Simón-Díaz et al. 2015;Holgado et al. 2018Holgado et al. , 2020)).This survey covers a range of stellar masses from 15 to 80 solar masses, and a mixture of ages within the main sequence stage.Besides, they compared the rotation of these stars with the state-of-the-art evolution models Brott et al. (2011) and Ekström et al. (2012), finding that neither of the two sets of rotating evolutionary tracks was able to reproduce such features of the survey.
From one side, models from Brott et al. (2011) cannot reproduce the existence of stars with rot 150 km s −1 across the entire domain of O-type stars (see Fig. 9 of Holgado et al. 2022).From the other side, models from Ekström et al. (2012) cannot adequately reproduce the scarcity of stars with rot 75 km s −1 for the left side of Fig. 10, which is also appreciated by the lack of stars matching the old tracks with 70 and 120 M .
In contrast, in Fig. 10 we observe that new self-consistent evolutionary tracks can adequately explain both issues.Tracks from 25 to 70 M match to the set of stars with rotational velocities below 150 km s −1 for the range of temperatures from 27 to 40 kK.Also, the empty region of stars with T eff ≥ 42.5 kK and rot ≥ 75 km s −1 would only be populated by our 120 M model, which is out of the mass range of the sample.Therefore, evolution models adopting weaker winds better interpret the rotational properties of the survey of Holgado et al. (2022), keeping with initial Ω = 0.4 and without the need of decreasing the initial equatorial velocity for our models.Nonetheless, despite this encouraging result, there are still important aspects of the rotational properties of Galactic O-type stars that need to be taken into account.For instance, although our track of M init = 25 M can encircle a larger fraction of stars in the top-right side of Fig. 10, there are still fast rotators ( rot > 300 km s −1 ) outside the scope of any track, which only could be explained by binary interaction effects (de Mink et al. 2013;Wang et al. 2020).Moreover, the slow braking observed for self-consistent evolutionary tracks indicates that stars born with masses between 25 and 70 M should spend the ∼ 75% of the lifetime in the main sequence stage with rot between ∼ 250 km s −1 and ∼ 330 km s −1 .As a consequence, we should expect to find a larger fraction of O-type stars in this range of rotational velocities.However, because of its proximity to the ZAMS, this range also covers an important portion of the region where there is a lack of empirically detected O-type stars, as shown in the Fig. 11 and described by Holgado et al. (2020).For that reason, it is not surprising to find a moderate number of stars in the range from 250 to 330 km s −1 and temperatures between 30 and 42.5 kK.Therefore, even though our evolution models adopting self-consistent winds represent a relevant upgrade, there are still important challenges in the study of O-type stars to deal with, such as the expansion of current samples to more hidden places by the use of infrared observations.

Implications of evolution models with self-consistent winds
The evolutionary tracks across the HRD from our rotating models adopting self-consistent winds exhibit considerable contrasts with the old models adopting Vink's formula, leading into a better description of the rotational properties of the most recent observational diagnostics.In this Section we move one step beyond, and we explore some of the implications that the incorporation of these self-consistent models have at Galactic scale.We remark however that these implications are here discussed at an introductory level, and their respective conclusions are part of forthcoming studies.4.9 × 10 −5 current rotating models allow 26 Al to appear at the surface well before layers enriched by nuclear burning in 26 Al are uncovered by the stellar winds.This is the effect of the rotational mixing that allows diffusion of the 26 Al produced in the core to reach the surface.Thus we see that, like for the case in Paper I, reducing the mass loss rate reduces the maximum value of the abundance of 26 Al reached at the surface.This maximum also occurs at a more advanced age for the model with self-consistent wind, despite the peak of production of 26 Al in the core is found at almost the same age for both (old and new) wind prescriptions.Such time difference between the peak abundance at the centre and at the surface is easily explained by the reduction on mass loss rate, which makes the star to take more time to remove their outer layers and then expose its inner composition, and for the less intense mixing as previously evidenced in Fig. 6a.The total amount of 26 Al released from the stellar surface to the interstellar medium during the MS for each one of our wind prescriptions, together with previous results from Paper I, are tabulated in Table 4.Although the total amount of the aluminium-26 ejected by old winds is almost the same regardless if rotation is considered, for self-consistent winds the rotating model predicts the ejection of ∼ 8 times less fraction of the isotope during the MS, compared to the rotating evolution model adopting Vink's formula.Such a difference is related not only to the mass loss due to the line-driven mechanism.Whereas nonrotating models predict eruptive processes associated to LBVs (Fig. 7 from Paper I), rotating models for M zams = 120 M never reach those magnitudes for Ṁ (Fig. A.3).Therefore, the con-tribution of 26 Al to the interstellar medium predicted by selfconsistent winds is even weaker (compared with models adopting Vink's formula) for rotating models than for non-rotating models.
However, it is important to remark that Fig. 12 covers only the lifetime of the star when the wind is optically thick, then excluding optically thick winds from WNh (H-core burning) and WR (He-core burning) phases.In contrast, the study of Martinet et al. (2022) calculated the evolution of 26 Al for evolution models with mass ranges from 12 to 300 M and metallicities from Z = 0.002 to Z = 0.020, but without making any distinction between thin and thick wind regimes.Instead, they made the distinction between H-core and He-core burning stages, showing that the production of 26 Al abruptly decreases when the star enters into the He-burning.As a consequence, even though here we are selecting only one mass and one metallicity to do a quick comparison with the more complete analysis from Martinet et al. (2022), we infer that the larger contribution of 26 Al from very massive stars to the ISM must stem from the later stages exhibiting optically thick winds, such as H-core WNh stars, to explain the current estimations of 26 Al total abundance in the Milky Way (Knödlseder 1999;Diehl et al. 2006;Wang et al. 2009;Pleintinger et al. 2019).Certainly, we need a extensive study implementing self-consistent evolution models for wider ranges of masses and metallicities, together with updates in the convection criteria (Georgy et al. 2014;Kaiser et al. 2020) and overshooting values (Martinet et al. 2021) in order to have more complete analysis.

Massive stars at the Galactic Centre
The implications of the new evolution models for massive stars with self-consistent mass loss rate are large, not only across the main sequence but also over the subsequent evolutionary stages.One of such consequences is how using models computed with the weaker self-consistent mass loss rates may actually provide estimates of evolutionary masses that are lower than when models with Ṁ from Vink et al. ( 2001) are adopted.
An example of this situation is the study of the Ofpe3 stars located at the Galactic Centre.Massive stars at the Galactic Centre play an important role feeding the supermassive black hole Sgr A* by means of their stellar winds (Cuadra et al. 2015;Ressler et al. 2018;Calderón et al. 2020b).This is a consequence of the remarkable situation of OB and WR stars representing a large fraction of the stars at the Galactic Centre (e.g., Lu et al. 2013;von Fellenberg et al. 2022), even though this is a region thought to be hostile to star formation (e.g., Genzel et al. 2010).The wind properties of these massive stars have been analysed and constrained by Martins et al. (2007) for the O-type and WRs, and by Habibi et al. (2017) for the B-type stars.Martins et al. (2007) calculated the stellar and wind parameters for a set of stars at evolved stages such as Ofpe and WN, which are particularly relevant for feeding Sgr A*.Their analysis was performed by spectral fitting using the CMFGEN code and the evolutionary tracks from Meynet & Maeder (2005).Therefore, it is important to check whether the state-of-the-art evolutionary tracks suggest a revision of the properties of the massive stars at the Galactic Centre, and the consequences these modifications imply for the accretion on to Sgr A*.To illustrate the issue, Fig. 13 shows the hydrogen abundance at the stellar surface as a function of luminosity for standard and self-consistent evolutionary tracks, plus five of the Ofpe stars analysed by Martins et al. (2007, compare with their Fig. 21).From Fig. 13 we can see that for a given hydrogen surface abundance and for a given initial mass, the tracks computed with the weaker self-consistent mass-loss rates are overluminous.These tracks also end the MS phase with higher surface hydrogen abundances.We notice that the lowest luminosity Ofpe star in the sample cannot be reproduced by any of the two families of tracks.The other four stars can reasonably be fitted by both types of models.On average, although the initial masses that will be deduced from comparison with self-consistent mass-loss rate tracks are smaller than the masses deduced from tracks computed with Vink's mass loss rates.
Given that the observed Ofpe stars in the Galactic center would have lower initial masses, we can speculate that their wind parameters (terminal velocity and mass loss rate) will be lower as well.The terminal velocity is an important parameter in this context, since the Galactic centre is a region with large stellar density, and the winds collide with each other, creating a hightemperature medium.However, slower winds ( 600 km/s) produce a plasma of lower temperature ( 4 × 10 6 K) at the collision, which becomes susceptible to hydrodynamical instabilities and ends up forming high-density clumps and streams.These clumps can then be captured by the central black hole, increasing its accretion rate (Cuadra et al. 2005(Cuadra et al. , 2008;;Calderón et al. 2016Calderón et al. , 2020b,a),a).
In a forthcoming paper we will perform a more complete analysis of the Galactic Centre massive star population, extending the tracks to the WR stage and taking into account the nonstandard chemical abundance deduced for this region.From the observational side, we will also include data collected after Martins et al. (2007), which is expected to afford a reduction in the typical error bars for the stellar luminosities from 0.2 down to 0.1 dex (S. von Fellenberg, priv. comm.).The comparison between models and observations will allow us to update the estimates for the initial masses and ages of this population, and also its wind parameters, which have a large influence on the current accretion on to Sgr A*.

Conclusion
We have extended the stellar evolution models performed in Gormaz-Matamala et al. (2022b, Paper I), adopting selfconsistent m-CAK prescription (Gormaz-Matamala et al. 2019, 2022a) for the mass loss rate recipe (Eq.1), by including the rotational effects.Stellar rotation affects the mass loss of a massive star, not only by changing the balance between gravitational, radiative and centrifugal forces (Maeder & Meynet 2000), but also because rotational mixing considerably modifies the internal distribution of the chemical elements and thus impacts the evolutionary tracks in the HRD and hence the mass loss rates.The progressive increase of the helium-to-hydrogen ratio in the wind impacts the line-force parameters (k, α, δ) and henceforth the self-consistent solution of the equation of motion, leading to a decrease of 0.2 dex in the absolute value of the mass loss rate (Eq.11).
The updated mass loss rates are implemented for a set of evolutionary models at different initial stellar masses, for metallicities Z = 0.014 (Galactic) and Z = 0.006 (LMC).New tracks show important differences with respect to the studies of Ekström et al. (2012) and Eggenberger et al. (2021), who adopted the formula from Vink et al. (2001) for the winds of massive stars, but exhibit a fair proximity with the studies of Sabhahit et al. (2022) for VMS, and with Björklund et al. (2022) and their own self-consistent wind prescription.Besides the differences in the tracks across the HRD, already remarked in Paper I, for rotating evolutionary tracks we find as expected that the surface rotation maintains a higher value when the weaker self-consistent mass loss rates are adopted.We observe that for the initial rotation considered here, tracks computed with the self consistent mass loss rates extend more to the red part of the HR diagram during the MS phase Such an effect is more important for masses above about 40 M .
New tracks predict evolution of the surface rotational velocities for O-type stars that, at first sight, appear to be in better agreement with the most recent observational diagnostics.The slow braking in the evolution of rot explains better the rotational properties of the survey of Holgado et al. (2022) for O-type stars, such as the lack of fast rotators for stars with T eff 42.5 kK and the abundance of stars with rot ≤ 150 km s −1 and temperatures between 27 and 40 kK.Besides, the implications of these new evolution models are wide.For example, lower mass loss rate predicts a less important stellar wind contribution to the Aluminium-26 enrichment of the ISM during the main sequence phase, at least whereas the stellar wind is optically thin.Likewise, the fact that self-consistent models are more luminous in the HRD suggests that the initial mass deduced from evolutionary tracks might be lower when the self-consistent tracks are used.This may apply to the mass estimates of Ofpe and WN stars at the Galactic Centre, which leads us to expect that their wind properties (mass loss rate and terminal velocity) might also be overestimated.Nevertheless, a more accurate analysis of the stars surrounding Sgr A* and their stellar winds are deferred to a forthcoming paper.

Fig. 1 :
Fig. 1: Evolutionary tracks across the HRD for models with rotation

Fig. 2 :Fig. 3 :
Fig.2: Differences in the resulting self-consistent mass loss rate ∆ log Ṁsc for the abundances tabulated in Table1, (a) due to the modification of the He/H ratio with respect to the default He/H= 0.085, (b) due to the modification of any of the CNO elements with respect to ξ/ξ = 1.0 (with ξ being carbon, nitrogen or oxygen depending on the respective symbol).

Fig. 4 :
Fig.4: Evolution of surface equatorial velocities and rotation parameter Ω as a function of the mass fraction of hydrogen at the core, which decreases from X core ∼ 0.72 at ZAMS to X core ∼ 0.05 at the end of the H-core burning stage.The symbols (dots and crosses) have the same meaning as indicated in the caption of Fig.3.
Fig. A.1, where the evolutionary tracks with and without rotation are compared.The rotating tracks remain in bluer regions of the HR diagram compared to the non-rotating ones.Similar to Paper I, self-consistent mass loss rates are of a factor about ∼ 3 lower than the mass loss rates from Vink's formula (Fig. A.3), which makes the stars retain more mass (Fig. A.4) and have larger radii (Fig. A.5)

Fig. 5 :
Fig.5: Diffusion factors for convective turbulence (D conv ), shear turbulence (D sh ) and effective turbulence (D eff ) as a function of the Lagrangian mass coordinate, for our evolution models with M zams = 120 M and Z = 0.014, adopting old and new winds at an intermediate point in its main sequence stage (X core = 0.3).Vertical grey lines represent the total mass of both models.Symbology for solid and dashed lines, is the same as in Fig.3.

Fig. 6 :
Fig. 6: Variation in a logarithmic scale of the mass fraction of hydrogen, helium and CNO elements, as a function of the Lagrangian mass coordinate in solar units inside the (a) 120 M model and (b) the 70 M model, for the two prescription of the mass losses studied in the present work at the end of the main sequence stage.The meaning of the different line style is the same as in Fig.3.Vertical grey lines represent the total mass of both models, which correspond to the final masses at the end of MS tabulated in Table3.

Fig. 7 :
Fig. 7: Angular properties for the inner structure of our 120 M model on the ZAMS and at the end of the MS phase for the two prescriptions of the mass losses studied in the present work.(a) Variation in a logarithmic scale of the specific angular momentum as a function of the Lagrangian mass coordinate in solar units.(b) Angular velocity, in radians per second, also as a function of the Lagrangian mass coordinate.
. A.1). Differences in the rotational velocities are also less remarkable, though they still follow the same trend of the more massive mod-els.For the case of Z = 0.006 (LMC), we observe that the end of H-burning is reached inside the region of validity of the selfconsistent tracks (Fig.A.2), and therefore the second dot overlaps again with the cross.Because we observe this situation only with Z = 0.002 (SMC) in the non-rotating cases from the Paper I, it implies that the evolution models adopting m-CAK wind prescription cover a broader range of the main sequence when rotation is incorporated.The impact of the different mass losses on the chemical structures of the 25 and 40 M stellar models can be seen in Fig 8.The convective cores are 8.5 and 18.5 M for the 25 and 40 M stars, respectively, when the high mass loss rates are used (Fig. A.8).

Fig. 10 :
Fig. 10: Rotational velocities as a function of the effective temperature, for evolutionary tracks adopting old (dashed) and new (solid) winds.The coloured dots represent the intervals of age with a step of 0.5 Myr.Grey dots represent the sample of O-type stars taken from the survey of Holgado et al. (2022).

Fig. 11 :
Fig. 11: Spectroscopic Hertzsprung-Russell diagram (Langer & Kudritzki 2014, sHRD, with L := T 4 eff /g), for evolutionary rotating tracks adopting old (dashed) and new (solid) winds.Grey dots represent the sample of O-type stars taken from the surveys of Holgado et al. (2022).The black line represents the ZAMS for all models, whereas dark yellow line represents the region close to the ZAMS where no stars are foundHolgado et al. (2020).

6. 1 .Fig. 12 :
Figure12shows how the new mass loss rates impact the quantity of 26 Al ejected by the stellar winds of our 120 M models computed with the two mass loss rate prescriptions.Compared to the models discussed in Paper I for the non rotating case, the

Fig. 13 :
Fig. 13: Surface hydrogen abundance (in mass fraction) as a function of the stellar luminosity for different evolutionary tracks.Continuous lines are models computed with the self consistent mass-loss rate, while dashed lines are computed with Vink's original recipe.Magenta symbols are the Ofpe stars from the Galactic Centre as plotted by Martins et al. (2007, compare with their Fig.21).

Table 2 :
Summary of the formulae implemented for old and new evo-

Table 3 :
Properties of the stellar models at the end of the Main Sequence.

Table 4 :
Values for the integration of the mass fractions 26 Al surface for our evolution model of M zams = 120 M and Z = 0.014 (see Fig.12), compared with the results from Paper I (non-rotating tracks).