Very Massive Stars Models I. Impact of Rotation and Metallicity and Comparisons with Observations

,


Introduction
Very massive stars (VMS), often considered as stars with initial mass larger than 100 M ⊙ , have been mostly stimulating interest in Population III and early Universe.Indeed, early hydrodynamical simulations predicted preferential formation of very massive stars for Population III stars (Abel et al. 2002;Bromm et al. 2002).More recent simulations predict now a more significant fragmentation process (Stacy et al. 2010;Clark et al. 2011), binaries fraction (Turk et al. 2009) and wide initial mass distributions (Hirano et al. 2014(Hirano et al. , 2015)).Moreover, the IMF might even be redshift-dependent for Population III stars (Hirano et al. 2015) due to the different temperature of formation sites.There is, however, a general agreement on a top-heavy primordial IMF (Greif et al. 2011;Stacy & Bromm 2013;Hirano et al. 2014Hirano et al. , 2015;;Susa et al. 2014;Stacy et al. 2016;Jeřábková et al. 2018;Wollenberg et al. 2020).Population III stars could be a very significant source of ionizing photons and thus be interesting candidates for the reionization of the intergalactic medium (Haehnelt et al. 2001;Faucher-Giguère et al. 2008, 2009;Becker & Bolton 2013;Wise et al. 2014;Sibony et al. 2022).Using the present Pop III models, Murphy et al. (2021) showed that including ini-tial masses up to 500 M ⊙ can increase the total number of ionizing by 30% compared to the case where the upper mass limit is chosen equal to 120 M ⊙ .While VMS have an undeniable impact in the early universe, their contribution across the cosmic time might be underestimated even in higher metallicity environments.
One of the most spectacular observations in the early 80s was the observation in the 30 Doradus Nebulae, in the Large Magellanic Cloud (LMC), of the central object R136, that Cassinelli et al. (1981) showed to potentially be a ∼2500 M ⊙ star according to their spectroscopic analysis.With higher spatial resolution instruments, R136 was then shown later on to be a young cluster, composed of much lower mass stars.The stellar upper mass limit was then found to lie around 150 M ⊙ (Weidner & Kroupa 2004;Figer 2005;Oey & Clarke 2005;Koen 2006).Crowther et al. (2010) used new adaptive optics observations combined with more refined non-LTE atmosphere code to derive the fundamental parameters of these objects.Their results show that for four stars in R136, models suggest masses in the 165-320M ⊙ mass range.This has been confirmed by recent studies using new HST data (Bestenlehner et al. 2020;Brands et al. 2022).
Interestingly, VMS could also be viable candidates to explain the self-enrichment of globular clusters (see Vink 2018).For the last decade, one of the main candidates for this selfenrichment process has been massive AGB stars (D'Ercole et al. 2010) due to their slow winds compare to the fast line-driven winds of massive O-stars.Indeed, such fast winds (up to 2000-3000 km.s −1 ) are not able to deposit the enriched material efficiently into the globular clusters as the winds are too fast to be trapped in their potential well.Alternatives such as rotating massive stars (Decressin et al. 2007), massive binaries (de Mink et al. 2009;Vanbeveren et al. 2012), red supergiants (Szécsi et al. 2018), or supermassive (SMS) stars (Denissenkov & Hartwick 2014) have been proposed and are discussed in Bastian & Lardo (2018).Interestingly, SMS have been proposed as a workaround to the so-called mass-budget problem (Gieles et al. 2018), but the estimated wind velocities of roughly 1000 km s −1 are still too fast compared to the estimated escape velocities (≃ 500 km s −1 ) at the center of young globular clusters (Gieles et al. 2018).Vink (2018) propose then that VMS, after inflating due to their proximity to the Eddington limit (Ishii et al. 1999;Gräfener et al. 2012;Sanyal et al. 2015), lose mass through slower winds, satisfying conditions on both wind velocity for enrichment and massbudget.
Moreover, VMS have been shown to have a potential important impact for the chemical enrichment at high metallicity.Indeed, Martinet et al. (2022) have shown that the galactic production of the short-lived radioactive isotope 26 Al by the winddriven mass loss of massive stars might be increased by up to 150% when including VMS into a Milky Way-like population, and could then account for a large part of the observed 26 Al galactic content.
To explore the evolution of VMS across the cosmic history, and as a follow-up to the work done in Yusof et al. (2013), we computed a grid of both rotating and non-rotating VMS from Population III stars to solar metallicity stars with GENEC (Eggenberger et al. 2008).In this work, we take advantage of the improvements made compared to Yusof et al. (2013) models, such as the improvement of angular momentum conservation (Ekström et al. 2012), the use of another convective criterion (Kaiser et al. 2020), a revised overshooting parameter (Martinet et al. 2021) and the introduction in this work of an equation of state accounting for the pair-creation production.
The paper is organized as follows: in Section 2 we describe the physical ingredients used to compute the present models for very massive stars.The effect of different physics on very massive stars evolution is discussed in Section 3. Comparison between our models and observed VMS in the LMC is discussed in Section 4. In Section 5 we discuss and suggest potential impact of different physics on the final fate of VMS, and we present the main conclusions of this work.

Ingredients of the stellar models
The present models have been computed with the 1D stellar evolution code GENEC (Eggenberger et al. 2008).These models differ mainly in three points with respect to the physics used in the grids by Ekström et al. (2012) and Yusof et al. (2022).First, we adopted the Ledoux criterion for convection instead of Schwarzschild and an overshoot of 0.2 H p instead of 0.1 H p .These changes were made because there are some indications that the Ledoux criterion might be more appropriate for reproducing some observed properties of massive stars (Georgy et al. 2014;Kaiser et al. 2020), and that an increase in overshoot parameter is needed for stars with masses higher than 8 M ⊙ (Martinet et al. 2021;Scott et al. 2021).The third change is that our very massive star models have been computed with an equation of state accounting for electron-positron pair production (Timmes & Swesty 2000).The other ingredients are the same as in Ekström et al. (2012) and Yusof et al. (2022).To make the paper more self-consistent, let us however remind the prescriptions used for rotation and mass loss.We considered here the physics of shellular rotating models as described in Zahn (1992).For the vertical shear diffusion coefficient, we used the expression by Maeder (1997) and the expression of the horizontal shear diffusion coefficient by Zahn (1992).
The radiative mass-loss rate adopted on the MS is from Vink et al. (2001); the domains not covered by this prescription (see Fig. 1 of Eggenberger et al. 2021) use the de Jager et al. (1988) rates.Gräfener & Hamann (2008) prescriptions are used in their domain of application, while Nugis & Lamers (2000) prescriptions are used everywhere else for the Wolf-Rayet phase 1 .We kept here for consistency the same prescriptions as for the overall grids published so far, however more recent prescriptions are now available (see the review of Vink 2021, and Fig. 8).In Sect.6.1, we propose a more in-depth discussion on new mass loss prescriptions and the relevance of the prescriptions used in this work.The radiative mass-loss rate correction factor described in Maeder & Meynet (2000) is applied for rotating models.The dependence on metallicity is taken such that Ṁ(Z) = (Z/Z ⊙ ) 0.7 Ṁ(Z), except during the red supergiant (RSG) phase, for which no dependence on the metallicity is used.This is done accordingly to van Loon et al. (2005) and Groenewegen (2012a,b) showing that the metallicity dependence for the mass loss rates of these stars do appear to be weak.
We computed stellar models with initial masses of 180, 250, and 300 M ⊙ for Z=0.000, Z=10 −5 , Z=0.006 and Z=0.014 metallicities, with no rotation and a rotation rate of V/V c =0.4, where V c is the critical velocity 2 .The nuclear network allows following the abundance variation of 30 isotopes 3 .

The grid of VMS stellar models
Table 1 presents the initial parameters and various physical quantities of the models presented in this paper at different evolutionary stages.We define the beginning of a phase when 1% of the central burning element has been consumed and the end of a phase when the mass fraction of the main fuel at the center is lower than 10 −3 .
Figure 1 presents the Hertzsprung-Russell diagram (HRD) of the higher metallicities models, at Z=0.014 and Z=0.006.The non-rotating models are displayed in solid lines, while the rotating ones are in dashed lines.Starting with the 180 M ⊙ at Z=0.014, we can see that the rotating model evolves more to the blue during the MS due to the rotational mixing.It enters the WR phase at a slightly earlier evolutionary stage.The rotating model enters the core He-burning phase at a much higher effective temperature than the non-rotating ones.The evolution through the advanced phases is then very similar due to the dom-Table 1. Properties of the stellar models at different stages of their evolution.The first columns present the initial mass M ini in M ⊙ , the initial surface velocity ν ini in km.s −1 and the averaged surface velocity over the main-sequence νMS in km.s −1 .For the H-, Heand C-burning phases, the columns present the durations of the phase τ phase in Myrs for H-and He-and in yrs for C-, the current total mass M tot in M ⊙ , the surface velocity ν surf in km.s −1 and the surface mass fraction of helium Y surf at the end of each nuclear burning phases.inant WR mass loss, with the rotating models at slightly lower luminosity.
If we go to a higher initial mass, the evolution is now highly dominated by the larger mass loss, with still the rotating models evolving with slightly lower luminosities.The 300 M ⊙ at Z=0.014 is the perfect example of mass loss dominated evolution, where both the non-rotating and the rotating models undergo so much mass loss that their evolution through the HRD is effectively similar.
If we go to lower metallicity, we now have a lower contribution to the mass loss of the line-driven wind.This means that the evolution of VMS will now be more sensitive to the rotational mixing.Indeed, as we can see for the whole 180-300M ⊙ mass range, the rotating models undergo quasi-chemically homogeneous evolution due to the very strong mixing.Rotating stars reach the WR phase much earlier.This leads the star to undergo large mass loss.The He-burning phase of rotating models occurs at much lower luminosities than for the non-rotating models.The same for the more advanced phases.
To synthesize these results, we can say that at solar metallicity, the tracks of the models for the 250 and 300 M ⊙ nearly always decrease in luminosity as a result of the very strong mass loss.The effects of rotation become less and less strong when the initial mass increases because of the dominating effect of mass loss.At the metallicity Z = 0.006 as well as for the 180 M ⊙ at solar metallicity, i.e. for all the models where the mass losses are not so strong, the effects of rotation are important.As is well known, rotating tracks are bluer on the MS phase.The models enter at an earlier evolutionary phase into the WR phase and thus reach lower luminosities during the post MS phases.
Figure 2 shows the HRD of 180 M ⊙ non-rotating stellar models at different metallicities.One of the main effects of metallicity on the surface properties is, of course, the cooler effective temperature at higher metallicity (see Fig. 2 bottom right panel).Indeed, we can see that the beginning of the H-burning phase, while starting at a roughly similar luminosity, starts at very different T eff .This is mainly due to the fact that the opacity of the outer layers increases with the metallicity.Since the gradient of the radiative pressure scales with the opacity, this means that the outer layers receive a stronger radiative support in high metallicity stars than in low metallicity ones.This produces stars with larger radii and lower effective temperatures.Moreover, a higher metallicity implies more CNO elements in stars, that in turn implies that, for producing a given amount of energy, smaller temperatures are needed.Smaller temperature means that the star needs to contract less to allow nuclear reactions to produce enough energy for compensating the losses at the surface (see Sect. 4.2 in Farrell et al. 2022).This also contributes to making metal-rich star more extended than metalpoor ones.The evolution of the first parts of the MS is similar for the metallicities considered in Fig. 2, until the effect of the linedriven wind, very sensitive to metallicity, truly kicks in.From then, the higher metallicity models start to evolve to the blue due to large quantities of mass loss, quickly forming WR stars.
At low metallicity, the star remains more massive, and this has a strong impact on the evolutionary tracks in the HRD.During the MS phase, tracks cover a larger range both in luminosity and effective temperatures.After the MS phase, since there is still an important H-rich envelope, the star reaches larger radii.
During the post-core He-burning phases, the low metallicity models have larger luminosities than the more metal-rich models.This is a consequence of their larger actual masses.Their effective temperature is not reaching as large values as the stripped core resulting from the higher metallicity models.It is interesting to note that these models are very close to the Eddington  limit and hence undergo large mass loss events for a short time during the advanced phases.
Figure 3 shows the evolution of the surface rotation for a 180 M ⊙ at low metallicities.At these metallicities, the line-driven wind is almost negligible, hence the transport of angular momentum from the core spin-up the surface.In fact, the surface accelerates up to the critical velocity.From this point on, works usually assume that some matter will be lost from the star, forming an equatorial disk (see for example Hirschi 2007).Indeed, at the critical velocity, a small additional force is sufficient to launch the matter in Keplerian orbits around the star.The disk may lose mass (see e.g.Krtička et al. 2011).The net result will be a mechanical loss of mass by the star.
Numerically, this induces large problems of angular momentum conservation at the surface, that are very difficult to overcome when so close to the critical velocity.To manage this problem, we chose to artificially lower the effective velocity needed to trigger this mechanical mass loss to Ω limit /Ω crit =0.8.While the effects of lowering this limit are minimal in this case (≈2% of the total mass is lost by mechanical mass loss) and help tremendously the computation, it can be partly justified also by the uncertainties of the mass loss processes for very fast rotators.If we consider the case of the Be stars, which are stars showing a decretion disk likely resulting from their high rotation, it is not clear at which velocity those stars begin to lose mass (see the review by Rivinius et al. 2013, and references therein).Indeed, in addition to the centrifugal force, some other forces should be involved to kick the matter into an outward-expanding disk.It may be momentum given by radiation, or by pulsations.At least, it is likely that stars may begin to lose mass well before the critical limit is reached.Thus, choosing here a subcritical velocity for triggering the mechanical mass loss is not unrealistic.The star remains at the lower limit of the surface rotation for triggering the mechanical mass loss until the end of the MS phase (see Fig. tation becomes much lower than the limit and the mechanical mass loss stops.Other processes like line-driven winds or mass loss triggered by the continuum when some outer layers have supra Eddington luminosities become dominant.The zero metallicity model has a similar evolution on the MS as the Z=10 −5 since it is not dominated by the line-driven mass loss.However, the Pop III model has not been pushed further in its evolution due to numerical problems linked to mechanical mass loss.

Wolf-Rayet stars from VMS
Wolf-Rayet (WR) stars are classified into different categories from their spectral features.As we do not predict the output spectra of our models, we cannot use here the same spectroscopic criteria as those used by spectroscopists to classify our models among the different WR sub-types.Here we have adopted different theoretical criteria for WR classification.The Wolf-Rayet phase in GENEC is assumed to begin when the model has an effective temperature larger than 10 000 K and a surface mass fraction of hydrogen X 1 H s below 0.3.Assuming the star fulfills these WR conditions: the WNL phase is defined when X the WO phase when it is a WCO with log(T eff ) > 5.25 the WC phase when it is a WCO with log(T eff ) ≤ 5.25 It is important to note that these classification criteria may actually result in differences with the ones obtained from spectroscopy (see Groh et al. 2014).This is why the WO/WC criteria have been defined following the work of Groh et al. (2014).
Figure 4 shows the stacked lifetimes of different spectral phases of massive stars as a percentage of the total lifetime of the concerned star.The lifetimes are displayed for non-rotating (left panels) and rotating at V/V c =0.4 (right panels) models at Z=0.006 on the upper panels and Z=0.014 in the lower panels.The Z=0.006 models below 150 M ⊙ are coming from Eggenberger et al. (2021) and from Ekström et al. (2012) for the Z=0.014 ones.
Most of the life of massive stars is spent as O-type MS stars (note that the y-axis begins at 60%).The post-MS O-type star phase comprises between 10 and 40% depending on the initial mass, metallicity and rotation.In general, this phase covers a larger fraction of the total lifetime when the initial mass, the metallicity, and/or the initial rotation increases.For example, the whole post O-type MS phase lasts about 27% of the total lifetime of a 250M ⊙ non-rotating Z=0.006 model.Increasing the metallicity to Z=0.014 brings that fraction to about 30%.The rotating models at Z=0.006 and 0.014 still increase that fraction to respectively 35 and 40%.Except in the case of the non-rotating Z=0.006 model, where the fractions of the total lifetime spent as a B-type, WNL and WC stars are more or less the same, for the other models, the post O-type MS star phase is only divided in two type of stars, either WNL (the longest ones covering about 2/3 of the remaining time) and the WC phase (one third).This means that statistically, for single stars in a region with a constant star formation rate (for example during the last ten million years or so), there are in these last cases twice more chances to observe a WNL than a WC.
As is well known from previous models of rotating massive stars (see e.g.Fig. 2 in Georgy et al. 2012), rotational mixing increases significantly the mass of the star that has a chemical composition enriched by H-burning products, hence the longer duration of the WN phase in rotating models than in non-rotating ones.

Evolution of the core mass and final fate
The CO core mass is an important quantity that affects the advanced stages of massive stars and especially their final fate.This quantity is highly dependent on the initial mass, but also on rotation and mass loss events.
Figure 5 shows the evolution of the convective core mass of VMS between the beginning and the end of the He-burning phase, for different rotation rates and metallicities.The convective core at the end of He-burning is a good proxy for the CO core mass.The figure also shows the actual total mass of the star at the different evolutionary stages shown on that figure.
We see that for the masses considered in Fig. 5, the total mass significantly decreases during the core He-burning phase and imposes the convective core to also decrease.We note that at the metallicity Z=0.006, the total mass of the non-rotating models at the beginning of the core He-burning phase increases with the initial mass, and the same for the convective core mass.This trend disappears at the end of the core He-burning phase.Typically, the total actual mass of the 300 M ⊙ model is smaller than that of the 250 M ⊙ one.This is a consequence of the fact that the 300 M ⊙ entering into the core He-burning phase with a higher mass, loses more rapidly mass than the 250 M ⊙ , allowing the model to lose a larger fraction of its mass by the end of that phase.This reflects the fact that the mass loss rates increase with the luminosity that itself increases with the actual mass of the star.The rotating models show much fewer differences between the different initial masses.They begin their core He-burning phase with similar actual mass and lose a similar amount of their mass during that phase.Rotational mixing, as shown above, makes the star enter at an early evolutionary stage in the WR phase.From that stage, the mass loss scales with the luminosity and hence the mass, making the models converge towards similar masses.Indeed, as explained above, by starting with more mass, the high luminosity will produce stronger mass losses and thus causes the model to reach a similar mass as models starting with a lower mass.
At solar metallicity, the differences between the non-rotating and rotating models are much smaller than at Z=0.006.The strong winds dominate here, leaving little room for the effects of rotational mixing.
Very often in the literature, the mass of the He core or that of the CO core at the end of He-burning is used to determine whether the star will explode as PISN or PPISNe (see for example Farmer et al. 2019;Costa et al. 2021).As we have seen, the convective core mass at the end of He-burning is a good proxy for the CO core mass.In Fig. 6, the CO core masses limits given by Farmer et al. (2019) are indicated.We see that according to this criterion, the only two models that will explode as a pairinstability supernova would be the 250 and 300 M ⊙ non-rotating models at Z=0.006.The inferior mass limit would be 200 M ⊙ for entering into the domain of PISNe.The result is different from the one by Yusof et al. (2013) (see therein their Fig.18).They find that for the non-rotating Z=0.006, the inferior mass limit for the non-rotating models was found to be around 300 M ⊙ .This is likely due to the absence of models between 150 and 500M ⊙ in the Yusof et al. (2013) study.Indeed, we see that the final mass has a peak for an initial mass around 260M ⊙ in this study, so it is likely that the Yusof et al. (2013) study missed that peak.For rotating models at Z=0.006, the 250 and 300M ⊙ models are expected to undergo PPISNe.
We note also that the models for Z=0.014 are expected to avoid both Pulsational Pair instability phase and Pair instability explosions for all the models computed here, thus at least for models up to an initial mass of 300 M ⊙ .In case those models end their lives as black holes engulfing the whole mass that the star has retained until the final core collapse, black hole masses up to 30 -45 M ⊙ could be formed.to about 60 M ⊙ ) could be formed from the non-rotating Z=0.006 models if the mass loss induced by the pulsational pair instability would be smaller than a few solar masses.
In a paper in preparation, we will present a more in-depth study of the final fate of VMS where we will study both core masses criteria and an average value for the first adiabatic exponent over the whole mass (see Stothers 1999;Marchant et al. 2019;Costa et al. 2021).We will discuss whether indeed those models that should explode as a (P)PISNe according to the core mass criterion show a sufficiently large fraction of their total mass unstable to trigger the explosion mechanism (Martinet et al. in preparation).

Comparison with previous models
Figure 7 presents the evolution of a rotating (V/V c =0.4) 300 M ⊙ at Z=0.014 computed with GENEC, from the computations of Yusof et al. (2013), and the other one from this work.The main differences between these two computations are mainly coming from the improvement to the code, namely the treatment of conservation of angular momentum in the envelope (see Ekström et al. 2012), and the new EOS we have been introducing to take into account the pair-creation effects.As we also have seen in Sect.2, these new VMS models are computed with the Ledoux criterion and an increased overshoot α ov =0.2, while Yusof et al. ( 2013) models are computed with the Schwarzschild criterion and α ov =0.1.The left panel of Fig. 7 shows the HRD of these two models, with the beginning and end of the different burning phases depicted by circles and crosses.Z=0.014 Fig. 5. Evolution of the convective core mass of VMS stars at the beginning and at the end of the He-burning phase, for non-rotating and V/V c =0.4 models at both Z=0.006 and Z=0.014.The light grey region corresponds to the total actual mass of the star at the considered evolutionary stage, and is emphasized by a line of the corresponding color.
The beginning of the MS is very similar for both models, however slight discrepancies appear at the surface during the rest of the MS evolution.This is due to the different treatment of the angular momentum conservation at the surface, combined with the very large mass loss rates occurring for such massive stars at solar metallicity.The larger overshooting parameter does not induce a meaningful increase of the core size during the MS due to VMS having already very large convective cores.The total mass of both models during most of the MS is similar because mass loss dominates its evolution.However, at the end of the Hburning, a large mass loss occurs for both models, leading however to larger quantities lost in the Yusof et al. (2013) models.This leads our models to be 40% larger in mass at the end of the MS.
The He-burning phase occurs at higher luminosity and goes through a slight evolution to the red for the model of this work.The leading consequence is that the VMS will undergo slightly larger mass loss rates during the middle of the He-burning when going to the redder part of the HRD.This reduces the difference in the actual mass between our model and the one by Yusof et al. (2013) from 40% at the end of the MS phase to 35% at the end of the core He-burning phase.Interestingly, the evolution from the C-burning phase is slightly different, with Yusof et al. (2013) model continuing to higher effective temperature while this work model evolves back to the red.To explain this feature, we need to now go to the right panel of Fig. 7 where the central temperature of the star is depicted as a function of the central density.Once again we see that in early stages, the central conditions are very similar, with a slightly higher temperature in the core in the new model (likely due to the larger overshoot parameter).The very interesting features appear at the end of the C-burning phase and are displayed as a zoom in the upper-right part of the diagram.
The oscillations of the central conditions are in fact a direct effect of the pair production inside the star.Indeed, while the center is out of the pair-instability zone, a part of the core slightly off-centered is crossing this instability region.The production of pair inside these parts of the star results in destabilizing this region of the core4 .

Comparison with observed VMS in the LMC
The Tarantula Nebula in the Large Magellanic Cloud hosts one particular young cluster, known as R136, hosting itself very high mass components.Recent work (Bestenlehner et al. 2020;Brands et al. 2022) confirmed these results with new VLT-FLAMES and HST/STIS optic/spectroscopic measurements and provide interesting stellar characteristics to the most massive stars observed.With the new VMS models computed here, we can confront our results to these new observations.We explore here the impact of the physics we have included inside our mod- els on the observable, and compare the results to the ones provided by the most recent observations.As we have seen in the previous section, the evolution of VMS at high metallicity is dominated by mass loss.The prescriptions used in our models are coming from both theoretical and empirical results obtained from lower mass O-type stars, and are highly uncertain, especially when getting close to the Eddington limit (see the review of Vink 2021).We confront in Fig. 8 our models to the observed mass loss rates obtained by Brands et al. (2022) and their corresponding fit using the prescription proposed by Bestenlehner et al. (2020).The mass loss rates are plotted as a function of the Eddington parameter (with Γ edd,e =1 being the Eddington limit).The time-averaged mass loss of our models over the Main-Sequence is shown by colored rectangles.The time-averaged mass loss rates of our nonrotating VMS models with initial masses equal or above 180 M ⊙ are below the fitting line by at most 0.3 dex (a facor two).The contrary for our rotating models that show time-averaged mass loss rates at most about 0.3 dex above the fitting line.At least, this shows that although we did not account for the most recent developments of the mass loss rates in the present models, the fact that the non-rotating and rotating models more or less frame the fitting curve by Brands et al. (2022) with a factor of at most a factor 2 is rather encouraging.Indeed, the scatter of the colored squares is smaller than the scatter of the individual mass loss rate determinations.
Figure 9 presents the HRD of the Tarantula Nebula, a very active star-forming region in the LMC.The compilation of Schneider et al. (2018) is displayed in yellow, while the com-ponents of the R136 cluster, obtained by Brands et al. (2022), are displayed in magenta.The tracks used here are the Z=0.006GENEC tracks from Eggenberger et al. (2021) and the VMS tracks we computed for this work with the new EOS.We can see that this cluster hosts several high-mass components.If we look at the evolution of our VMS models, three stars inside this cluster are above M ini =200M ⊙ (R136a1/2/3).At first glance, one would favor the non-rotating models for these VMS, covering a large range of effective temperatures on the MS, while the rotating models, as we have seen in Sect.3, evolve quasichemically homogeneously and stray directly to the bluer part of the HRD (also true for models rotating at V/V c =0.2 not shown here), hence apparently incompatible with the observed effective temperatures.The problem is however more complicated, the VMS undergo very large mass loss already during the MS.This means, in fact, that the star is embedded inside a "nebula" of very thick wind.When observing such stars, we observe in fact this thick wind, and this blurs the determination of the effective temperature of the surface of the star.Indeed, with the wind being much cooler than the surface of the star, we would obtain a lower effective temperature to be observed if we took this effect into account.We computed corrected effective temperatures and find that both non-rotating and rotating models are compatible with the observations (the corrected effective temperatures go down to 35000K, 6000K below the lowest limit of the observed T eff of VMS see the details of how this correction has been computed in Schaller et al. (1992)).
Figure 10 shows the values of V sin i obtained by Bestenlehner et al. (2020) overplotted over the tracks of VMS models.From the spectroscopic studies, the surface velocities observed do not differentiate the surface velocity from the macro-turbulent velocity.This means that the lower limit of the surface velocity would be when the macro-turbulent is at a maximum.We take the maximum macro-turbulent velocities obtained by Simón-Díaz et al. (2017) for stars close to the Eddington limit to obtain the lower displayed limit.The upper limit is obtained from the angle of view effects.
The theoretical models (see the continuous lines in Fig. 10) shows that the surface velocity decreases rapidly as a function of time as a result of the large amount of angular momentum taken away by stellar winds.
To obtain that the theoretical model goes through the region covered by the error bars given by Bestenlehner et al. (2020) at this stage of the MS, we need quite high rotation rates already at the ZAMS.Actually, since, as recalled above, the V sin i value is a lower limit to the true surface velocity, even larger initial rotations than those shown for the models in Fig. 10 could be  considered.
A 250M ⊙ model with lower rotation rate has been computed to show that the mass loss still dominates the surface velocity evolution and would not be able to reproduce the high velocities observed for these VMS.We can wonder whether assuming a solid body rotation as would do a model accounting for the impact of an internal magnetic field, would allow to start from a lower initial rotation rate.This will be explored in the future.However, we can note that despite the fact that we considered non-magnetic models here, since these stars host very massive convective core they have a solid body rotation anyway in a very large part of the star.We suspect therefore that including magnetic fields in our models would have here a modest effect.(see discussion in Sect.6.1).
Figure 11 shows the observed surface enrichment of nitrogen as a function of the observed surface helium, respectively from Brands et al. (2022) and Bestenlehner et al. (2020).The non-  V/V c =0.4 rotating (in solid lines) and the rotating (in dashed lines) VMS models are overplotted.Interestingly, the evolution of nitrogen at the surface is very similar between the non-rotating and rotating models.This is because the models quickly reach (before the surface Helium fraction reaches 0.4) the maximum quantity of nitrogen produced from the CNO cycle in these stars.This maximum value depends solely on the initial abundances of CNO.
Hence, while the most massive component R136a1 is well reproduced by the models, there is no possibility for these models to explain the very high enrichment value of nitrogen of the two other VMS (R136a2/a3).Understanding these discrepancies is difficult.A first possibility could be that the initial abundances and in particular the initial content of CNO elements could be different in the different stars.We can see that the Z=0.014 mod-  et al. 2017;Sabín-Sanjulián et al. 2014, 2017;McEvoy et al. 2015;Bestenlehner et al. 2014a) and the R136 cluster members from the work of Brands et al. (2022).The GENEC tracks are overplotted in grey, with color emphasis on the VMS models computed here, and the ZAMS is drawn in dark grey.The HRD is plotted for non-rotating (left) and rotating at V/V c =0.4 (right) models.Their effective temperatures are not corrected for the optical thickness of the wind.The three most massive components of the cluster (R136a1/2/3) are pinpointed by arrows on the right panel.els are able to reproduce such values due to the larger initial abundances of CNO.Can we imagine that this region can host populations with different metallicities?This would imply, in a self-enrichment scenario, multiple stellar generations born from material differently enriched in metals by previous generations of stars.There are massive young clusters that show evidence for multiple generations of stars, such as the Orion nebula cluster with potentially generations with an age separation of less than 1 Myr (Beccari et al. 2017).The Tarantula nebula is a much more massive star-forming region than Orion, even more susceptible to having hosted multiple generations.Interestingly, a scenario where the (re)-collapse of the cluster allows the birth of a second generation of stars has been invoked by Domínguez et al. (2022) to explain the observational characteristics of the central region of 30 Doradus.A second possibility to explain the result shown in Fig. 11 is that the error bars have been underestimated.The error bars should be significantly increased to allow tracks at a given fixed metallicity to cross them.This appears not very realistic.Third, we can wonder if the past history of the star, such as binary interaction or early merger, could be invoked to explain the discrepancy in nitrogen surface abundances at so similar values for the surface helium abundance.Of course, some modeling of the evolution of multiple systems is required to answer such a question.

Impact of changing physical ingredients of the models
Changes of some physical ingredients as the criterion for convection, the overshooting, the expression for the diffusion coefficient in rotating models, or the change of the EOS have nonnegligible impacts on the outputs of the stellar models as discussed in Sect.4, where we performed a comparison with the models by Yusof et al. (2013) that used, otherwise also GENEC as the stellar evolution code.
For solar metallicity, models are dominated by mass loss, hence most of the differences come from the changes in angular momentum conservation and the slightly different evolution in the HRD.While non-rotating models are quite similar, rotating ones lead to final masses up to 42% larger than previous models.At Z=0.006, the differences are smaller due to less mass loss occurring.For the rotating models, the final masses are up to 12% larger than in Yusof et al. (2013) models.
We can wonder what would be the impact of changing the angular momentum transport by considering a more efficient process, such as the one proposed in the Taylor-Spruit (TS) dynamo (Spruit 2002;Maeder & Meynet 2004;Eggenberger et al. 2010).The present VMS rotating models at V/V c =0.4 are almost rotating as a solid body at the beginning of the MS due to their very large convective core.It is also the case for slower rotators such as the Z=0.006 180 M ⊙ rotating at V/V c ≃0.2 presented in Fig. 10.The rotation of the envelope and the core stays coupled in fact up until the end of the MS for models at Z=0.006, when V sin(i) [km.s Fig. 10.Evolution of the surface velocity of VMS models.An additional model with slower initial rotation rate has been computed for this plot.The observed velocities of the three most massive components of R136 have been overplotted (R136a2 and R136a3 have similar velocities obtained).The uncertainties on the velocity are taken from the maximum macro-turbulent velocity observed for stars close to the Eddington limit (Simón-Díaz et al. 2017) for the lower limit, and from the maximum velocities obtained from angle of view effect for the upper limit.the mass loss becomes prevalent.A more efficient transport such as the TS dynamo would not change much the rotation profile in that case.For models at Z=0.014, differential rotation between envelope and core sets in from the middle of the MS, due to higher mass loss rates.In this case, a more efficient transport of angular momentum might produce stars with faster surface rotation.
VMS are highly sensitive to mass loss, especially at high metallicities.We chose for this work to keep the same physical ingredients that were used for all the grids published so far (Ekström et al. 2012;Georgy et al. 2013;Groh et al. 2019;Murphy et al. 2021;Eggenberger et al. 2021;Yusof et al. 2022).
In doing so, it allows us to discuss how varying the initial mass, metallicity and rotation impact the outputs.While the mass loss prescription scheme we use here is used in many recent grids for massive stars (see for example Fragos et al. 2023;Simaz Bunzel et al. 2023;Jiang et al. 2023;Brinkman et al. 2023;Song & Liu 2023), recent works presented new predictions for mass loss in massive stars (such as Björklund et al. 2022;Gormaz-Matamala et al. 2023;Yang et al. 2023).Björklund et al. (2022) predict that O-stars mass loss rates are lower by about a factor 3 than the rates typically used in stellar-evolution calculations, however, differences decrease with increasing luminosity and temperature (i.e. less difference for VMS and their high luminosity).Moreover, they find that the remaining key uncertainty regarding these predictions concerns unsteady mass loss for very highluminosity stars close to the Eddington limit, a major mass loss component due to VMS being very close to the Eddington limit (as we have seen in Fig. 8).For RSG mass loss rates, Massey et al. (2023) discussed the relevance of the prescriptions used in this work by comparing the predicted and observed luminosity function of red supergiants, finding a good agreement.Still, Yang et al. (2023) present a new prescription for RSG mass loss rates from a large spectroscopic survey in the Small Magellanic Cloud, finding that it may provide a more accurate relation at the cool and luminous region of the HR diagram at low metallicity compared to previous studies.The mass loss rates prescriptions we use at Z=0.006, although not considering the effects discussed in Bestenlehner et al. (2020), provide a time-averaged mass loss rate during the Main-Sequence that is underestimated by a factor of up to 2 for the non-rotating models and overestimated by a factor up to 2 for the rotating models.The scatter is reasonably small when compared to the scatter of the individual mass loss rate determinations used to establish the fitting mass loss rate prescription.It is however uncertain how our current prescriptions reproduce mass loss rate at solar metallicity, and if we under-or over-estimate the quantity of mass lost by solar metallicity VMS (Bestenlehner et al. 2014b;Vink 2021).For the low metallicity case, while they undergo almost no mass loss during the MS, their mass loss close to the Eddington limit is very uncertain and might have an important impact on their evolution in the later stages (Sander & Vink 2020).
Studying the impact of recent mass loss rate prescriptions on VMS is crucial (see Sabhahit et al. 2022Sabhahit et al. , 2023;;Higgins et al. 2023), but updating important physic such as chemical mixing and transport of angular momentum would be as necessary.Thus, new generations of stellar models should not only update the mass loss rates but also the physics of many other physical ingredients.This is beyond the aim of the present work and will be the topic of future works.

Synthesis of the main results and future perspectives
We computed a grid of VMS from Population III to solar metallicity with and without rotation.We showed that the low metallicity models are highly sensitive to rotation, while the evolution of higher metallicity models is dominated by mass loss effects.The mass loss affects strongly their surface velocity evolution, quickly at high metallicity while reaching the critical velocity for low metallicity models.The comparison to observed VMS in the LMC showed the mass loss prescriptions used for these non-rotating/rotating models are slightly under-/over-estimating compared to the observed mass loss rates.According to the treatment of rotation in our models, observed VMS need a high initial velocity to be able to account for their observed surface velocity.The surface enrichment of these observed VMS is difficult to explain with only one initial composition and could suggest multiple populations in the R136 cluster or maybe a more complex scenario involving multiple star interactions.
The tremendous dependence of the evolution of VMS upon the physic input will have of course a crucial impact on the final fate and the nature of the remnants that such stars can produce.As shown in multiple studies (Woosley & Heger 2021;Costa et al. 2021;Marchant & Moriya 2020;Farmer et al. 2019), the final fate of VMS is strongly linked to the CO core mass at the end of He-burning.This work provides a new illustration of the fact that rotation, mass loss, and other physic inputs play an important role in the CO core mass value, and could be essential to understand the final fate and the maximum mass of the remnants produced by VMS.A first discussion shows that at a metallicity typical of R136, only the non-or slowly rotating VMS models may produce Pair Instability supernovae.The most massive black holes that could be formed then are expected to be less massive than about 60 M ⊙ .In a dedicated follow-up paper, we will discuss the impact of the different physics introduced here on the final fate of the stars.Thanks to the addition of an equation of state following the pair production inside the star, we can now track the evolution of the pair-instability inside VMS.This will allow us to determine for the whole grid of models which ones are either undergoing core collapse, pulsational-, pair-instability supernovae or direct collapse through pair-instability.

Fig. 1 .
Fig.1.HRD of the non-rotating (solid lines) and rotating at V/V c =0.4 (dashed lines) models at Z=0.014 (solar) and Z=0.006 (LMC).The top/middle/bottom row shows model with initial masses of 180/250/300M ⊙ .The tracks are color-coded according to their mass loss rates.The beginning and the end of the burning phases are indicated respectively by circles and crosses.The beginning of the Wolf-Rayet phase is indicated by light blue circles.

Fig. 2 .
Fig.2.HRD for 180M ⊙ non-rotating models at Z=0.014, Z=0.006, Z=10 −5 and Z=0.The tracks are color-coded according to their current total mass.The MS of the models at Z=0.014, Z=0.006, Z=10 −5 are plotted in dashed grey in the lower right panel for comparison.

Fig. 3 .
Fig.3.Evolution of the surface angular velocity divided by the critical angular velocity for 180 M ⊙ rotating models initially at V/V c =0.4 (Ω/Ω crit ≃0.58), for Z=10 −5 and Z=0.In grey dashed line is shown the numerical limit used to trigger the mechanical mass loss.
1 H s > 10 −5 the WNE phase when both X 1 H s ≤ 10 −5 and the surface mass fraction of carbon X 12 C s is smaller or equal to the one of ni-

Fig. 6 .
Fig.6.Convective core mass at the end of He-burning as a function of the initial mass for non-rotating and rotating stars at Z=0.006 and Z=0.014.The models with initial masses lower than 180M ⊙ are fromEkström et al. (2012) for Z=0.014 andEggenberger et al. (2021) for Z=0.006, and share the same physic, albeit the differences described in Sect. 2. The range of CO core masses fromFarmer et al. (2019) where (P)PISNe can occur is shown as an example to underline the very different potential final fate of VMS when considering rotation and metallicity effects.

Fig. 7 .
Fig.7.Left: HRD of a rotating 300M ⊙ at Z=0.014 from the work ofYusof et al. (2013) and from this work.The beginning of each burning phase (defined when the central abundance of the concerned element decreased by 1%) and the end of it (defined by the central abundance being lower than 10 −5 ) are displayed respectively as open circles and crosses.Right: Central temperature as a function of the central density for the same models.A zoom can be seen in the upper left displaying the region where the model starts undergoing pair-creation inside the core, leading the model with the new EOS taking this effect into account to contract and expand rapidly as the pair-creation reduces the radiative pressure.

Fig. 8 .
Fig. 8. Time-averaged mass loss rates over the Main-Sequence as a function of the Eddington parameter Γ edd,e .The models of 85-120M ⊙ are from Eggenberger et al. (2021), and the 180-300M ⊙ from this work.The observations of Brands et al. (2022) of the R136 components are displayed in gold.The fit they obtained is displayed in purple.The three most massive components of the cluster (R136a1/2/3) are pinpointed by arrows.

Fig. 9 .
Fig. 9. HRD of the Tarantula Nebula from the compilation of Schneider et al. (2018) (regrouping results from Ramírez-Agudeloet al. 2017;Sabín-Sanjulián et al. 2014, 2017;McEvoy et al. 2015;Bestenlehner et al. 2014a) and the R136 cluster members from the work ofBrands et al. (2022).The GENEC tracks are overplotted in grey, with color emphasis on the VMS models computed here, and the ZAMS is drawn in dark grey.The HRD is plotted for non-rotating (left) and rotating at V/V c =0.4 (right) models.Their effective temperatures are not corrected for the optical thickness of the wind.The three most massive components of the cluster (R136a1/2/3) are pinpointed by arrows on the right panel.

Fig. 11 .
Fig. 11.Surface enrichment of Nitrogen over Hydrogen as a function of the surface Helium.The observed values of [N/H] are from Brands et al. (2022), and the surface Helium fraction from Bestenlehner et al. (2020).
Still higher mass black holes (up