© ESO, 2014

1. Introduction

Massive stars, with initial masses above ~8 M_⊙, are powerful cosmic engines (Bresolin et al. 2008). They produce copious amounts of ionizing photons, strong stellar winds, energetic final explosions, and most of the heavy elements in the Universe. The most massive amongst them conduct a life very close to their Eddington limit, and may thus be prone to become unstable. They are thought to be able to produce the most spectacular stellar explosions, like hypernovae, pair-instability supernovae, and long-duration gamma-ray bursts (Langer 2012).

While the value of the upper mass limit of stars is presently uncertain (Schneider et al. 2014), there is a great deal of evidence of stars with initial masses well above 100 M_⊙ in the local Universe. A number of close binary stars have been found with component initial masses above 100 M_⊙ (Schnurr et al. 2008, 2009; Sana et al. 2013b). Crowther et al. (2010) proposed initial masses of up to 300 M_⊙ for several stars in the Large Magellanic Cloud (LMC) based on their luminosities.

We present stellar evolution models of very massive rotating core-hydrogen burning stars, with initial masses up to 500  M_⊙. These models will be used for comparison with the VLT FLAMES Tarantula survey (Evans et al. 2011), which observed more than 800 O and early-B stars in the 30 Doradus region located in the LMC, to study the effects of rotational mixing and mass loss on the evolution of very massive stars.

Stellar evolution models for very massive stars have already been presented (see Table 1 in Maeder & Meynet 2011 and references therein). In particular, we refer to the rotating stellar models of Crowther et al. (2010); Yusof et al. (2013), who calculated models with initial masses up to 500  M_⊙. It is the aim of this paper to present a grid of stellar evolution models that has a dense spacing in mass and rotation rates, such that it is suitable for forthcoming population synthesis calculations. In this sense, our models are an extension of the LMC models of Brott et al. (2011a) to higher masses. In the mass range considered here, stellar wind mass loss and the proximity to the Eddington limit play prominent roles.

In Sect. 2 we briefly describe the employed stellar evolution code, including the input parameters for our calculations. Following that, we present and discuss the stellar evolution models in Sect. 3, including the evolution of their mass and their surface abundances, with emphasis on chemically homogeneous evolution, and we compare our models with previous work in Sect. 4. Section 5 contains our short summary. In Appendix A, we present isochrones derived from our models, as functions of age and initial stellar rotation rate.

2. Input physics and assumptions

Brott et al. (2011a) computed three grids of stellar evolution models for different metallicities (Galactic, SMC, LMC) for initial masses up to 60  M_⊙, which were compared with the results of the VLT-FLAMES Survey of Massive Stars (Evans et al. 2005, 2006). In particular, the convective core overshooting parameter and the efficiency parameters for rotationally induced mixing used by Brott et al. (2011a) were calibrated to reproduce the results of this survey (Hunter et al. 2008). The dense spacing in initial mass and rotational velocity used for the grid of Brott et al. (2011a) opened the door for statistical tests of the theory of massive star evolution (Brott et al. 2011b; Schneider et al., in prep.). Here, we extend the LMC grid of Brott et al. up to initial masses of 500  M_⊙, using the same numerical code, physical assumptions, and a similarly dense grid spacing.

To calculate rotating stellar evolution models, we use our one-dimensional hydrodynamic binary stellar evolution code (BEC) which is described in Heger et al. (2000); Petrovic et al. (2005), and Yoon et al. (2012). It contains a detailed treatment of rotation, angular momentum transport due to internal magnetic fields, and stellar wind mass loss. The code solves all five stellar structure equations throughout the stellar interior, including the stellar envelope up to a Rosseland optical depth of τ = 2 / 3. Convection is considered throughout the star using the non-adiabatic mixing-length theory. Our code is suited to treating stars close to the Eddington limit, and to describe the effects of envelope inflation which occur in this situation (Ishii et al. 1999; Petrovic et al. 2006).

2.1. Chemical composition

The initial chemical composition for our models is chosen according to corresponding observations of young massive stars and of H ii-regions in the LMC (Brott et al. 2011a). We thus adopt initial mass fractions for hydrogen, helium, and the sum of all metals of X = 0.7391, Y = 0.2562, and Z = 0.0047, respectively, with a non-solar metal abundance pattern (see Tables 1 and 2 in Brott et al.). We note that the applied opacities and mass loss rates (see below) are scaled with the LMC iron abundance, not with the total metallicity, which is reduced by 0.35 dex with respect to that of the Sun.

2.2. Convection and rotational mixing

Convection with a mixing-length parameter of α_MLT = 1.5 (Böhm-Vitense 1958; Langer 1991) is applied as well as semi-convection with an efficiency parameter of α_SEM = 1 (Langer et al. 1983; Langer 1991). In addition to convection, convective core overshooting is included with α_over = 0.335 local pressure scale heights, as calibrated in Brott et al. (2011a) with the rotational properties of B-type stars (Hunter et al. 2008; Vink et al. 2010) from the VLT-FLAMES survey. While no observational calibration of the overshooting parameter exists for stars of the considered mass range, we point out that the role of overshooting in our models is minor because of the large convective core mass fractions of very massive stars.

Rotational mixing (Heger et al. 2000) is considered with the efficiency parameters of f_c = 0.0228 and f_μ = 0.1 (Brott et al. 2011a). The most significant process causing rotationally induced mixing in our models is the Eddington Sweet circulation. Furthermore, the transport of angular momentum by magnetic fields due to the Spruit-Taylor dynamo (Spruit 2002) is applied, which is assumed here not to lead to an additional transport of chemical elements (Spruit 2006). Since the magnetic torques lead to a shallow angular velocity profile in our models, the effects of the shear instability, although included, are quite limited during the main sequence evolution.

2.3. Mass loss

The evolution of very massive stars is intimately connected to their mass-loss behaviour. Mass loss of very massive stars, for which few observational constraints exist, is a very active field of research. Here, we provide a brief overview of the mass loss prescriptions adopted in our calculations.

From the zero-age main sequence (ZAMS) up to a surface helium fraction Y_s = 0.4, we use the mass-loss predictions by Vink et al. (2000, 2001) for O and B stars. These rates are valid for stars of one million solar luminosities or less, i.e. stars below 80 M_⊙ (Mokiem et al. 2007), which are not very close to their Eddington limit (i.e. have Γ ≲ 0.3; see Eq. (3)). Empirical tests of these predictions depend critically on the presence and properties of small scale structures in the outflows, known as clumping. Mokiem et al. (2007) showed on the basis of Hα and He ii λ4686 analysis that the predicted rates agree with observations if the material is concentrated in clumps that have a 3–4 times higher density than in a smooth outflow. Other analyses, which include modelling of ultraviolet resonance lines, derive clumping factors that may reach values of 10 (e.g. Bouret 2004; Bouret et al. 2005, 2012; Fullerton et al. 2006). In this last case, the Vink et al. prescription adopted here may overestimate Ṁ by about a factor of 2. An improved hydrodynamical treatment shows that for normal O stars the mass loss rates may be somewhat lower (Müller & Vink 2008; Muijres et al. 2012).

Predictions for stars in the 40–300 M_⊙ range have been presented by Vink et al. (2011). Objects at the upper mass end of this range may have a very high Eddington factor Γ. It has been found that, at solar metallicity, the wind strengths agree with the standard Vink et al. recipe for objects that have an Eddington factor for Thomson scattering of Γ_e ≲ 0.7. For higher Γ_e values these new predictions show an upturn in the mass-loss rate, leading to rates that are higher by up to about a factor of 2 compared to the Vink et al. values used here. The authors associate the upturn with the stellar wind becoming optically thick, i.e. leading to spectral morphology which is typical for Wolf-Rayet stars of nitrogen subclass (WN) showing hydrogen emission lines. Interestingly, this predicted upturn may have been confirmed observationally (Bestenlehner et al. 2014). The relation between mass loss and Eddington factor has also been explored for late-WN stars by Gräfener & Hamann (2008) and Gräfener et al. (2011). They find a behaviour that is similar to the results of Vink et al. (2011), though the onset of Wolf-Rayet type outflows occurs at lower Γ_e values. Gräfener & Hamann however report a temperature dependence of their mass-loss rates that is steeper than that of Vink et al. (2011).

Since the mass loss predictions for large Eddington factors cannot yet be implemented unambiguously into stellar evolution calculations, we extrapolate the Vink et al. (2001) rates for stars above 10⁶ L_⊙. We note that Crowther et al. (2010) found the Vink et al. rates to agree within error bars with those observed in stars of up to ~10⁷ L_⊙ found in 30 Doradus. For objects in the range 60...100 M_⊙ the objects are close to the model-independent mass-loss transition point between optically thin and thick winds. In this range mass-loss rates have recently been calibrated with an uncertainty of only ~30% (Vink & Gräfener 2012).

The Vink et al. mass-loss prescription shows a bi-stability jump at about 25 000 K, leading to an increase of the mass-loss rate by a factor of 5 for stars of spectral type B1.5 or later. We include this bi-stability jump in our calculations (cf. Brott et al. 2011a). Additionally the Nieuwenhuijzen & de Jager (1990) empirical mass-loss rate is applied to cope with an increase in mass loss when approaching the Humphreys-Davidson limit (HD limit). The transition from Vink et al. (2000, 2001) to Nieuwenhuijzen & de Jager (1990) occurs at any effective temperature smaller than 22 000 K where the Nieuwenhuijzen & de Jager (1990) mass-loss rate exceeds the Vink et al. (2000, 2001) mass-loss rate.

To account for Wolf-Rayet mass loss, the Hamann et al. (1995) mass-loss rate divided by a factor of 10 is applied for the surface helium fraction Y_s ≥ 0.7. Figure 1 in Yoon et al. (2010) shows that this corresponds well to the Wolf-Rayet mass-loss rate proposed by Nugis & Lamers (2000) for Wolf-Rayet masses in the range 5  M_⊙ to 20 M_⊙. For surface helium mass fractions from 0.4 ≤ Y_s ≤ 0.7 the mass-loss rate is linearly interpolated between either the Vink et al. (2000, 2001) or the Nieuwenhuijzen & de Jager (1990) mass-loss rate and that of Hamann et al. (1995).

We apply a mass-loss enhancement for stars near critical rotation as in Yoon & Langer (2005) which considers a reduction of the critical rotational velocity for stars near their Eddington limit. It is still unclear whether rapid rotation per se leads to an enhanced mass loss (Müller & Vink 2014), but it appears reasonable to consider that the mass-loss rate increases close to the Eddington limit, which is indeed reached sooner for rotating objects (Langer 1997).

2.4. Model grid

Figure 1 gives an overview of our grid of evolution models by indicating the initial masses and initial surface rotational velocities of all computed model sequences. Because of the increase in the Eddington factor with mass, we decrease the maximum initial rotational velocity for higher masses in order to avoid strong rotationally induced mass loss already on the ZAMS (Langer 1998). Whereas most sequences are computed to core-hydrogen exhaustion, some of the most massive and most rapidly rotating models were stopped shortly before their proximity to the Eddington limit caused numerical difficulties.

3. Results

3.1. Evolutionary tracks in the Hertzsprung-Russell diagram

A selection of stellar evolution tracks is presented in Fig. 2. The luminosity of the stellar models for a given initial mass and initial surface rotational velocity is shown as a function of their effective temperature. Tracks are shown for nine different initial masses from 60  M_⊙ to 500  M_⊙, with initial surface rotational velocities of 0 km s^-1, 400 km s^-1, and 500 km s^-1.

Up to a luminosity of log  L/L_⊙ ≃ 6.5, at about 190  M_⊙, the effective temperature of the ZAMS increases with increasing luminosity. For higher initial masses this behaviour changes. The ZAMS moves towards lower effective temperature with increasing luminosity as a result of stars having significantly increased radii and as a consequence of their proximity to the Eddington limit. This effect is discussed in more detail in Sect. 3.6.

Figure 2 shows only the tracks of the slowest and fastest rotators in our grid. However, we point out that in the Hertzsprung-Russell (HR) diagram, the evolutionary tracks do not change much below an initial rotational velocity of ~250 km s^-1. This is demonstrated by the example of 80 M_⊙ models in Fig. 3, for which only the tracks with initial rotational velocities of ~300 km s^-1 and higher deviate significantly from the track of the model without rotation. As shown in Sect. 3.3, this is due to the absence of significant rotationally induced mixing of helium below the threshold rotational velocity for quasi-chemically homogeneous evolution (cf. Brott et al. 2011a). We note that although the evolutionary tracks of the sequences computed with initial rotational velocities below ~250 km s^-1 are almost identical, the corresponding models show significant differences concerning the evolution of the surface abundances of trace elements, e.g. boron and nitrogen (cf. Sect. 3.3).

The models without significant mixing of helium, in the mass range between 60 M_⊙ and 80 M_⊙, expand during core hydrogen burning to surface temperatures as low as 5000 K. This occurs partly because of the relatively large amount of convective core overshooting in our models. A second reason is that when these models evolve through core hydrogen burning, they approach the Eddington limit because their L/M ratio is increasing, and at the same time their envelope opacity is increasing as cooler surface temperatures are achieved (cf. Sect. 3.6). For higher initial masses, the redward evolution is truncated at T_eff ≃ 25 000 K because of the bi-stability mass-loss enhancement in the Vink et al. (2000) mass-loss recipe and the assumed mass-loss enhancement for stars near their Eddington limit (Sect. 2).

Figure 2 shows that the evolutionary tracks of even our slowly rotating models avoid the upper-right corner of the HR diagram. This is remarkable, since our very massive star models evolve very close to their Eddington limit, which leads to an inflation of the envelope (cf. Sect. 3.6). We find that the high mass-loss rates at temperatures below ~30 000 K lead to significant helium enrichments for all stars above ~60 M_⊙ (cf. Sect. 3.3) such that, as core hydrogen burning continues, they evolve towards hotter rather than cooler surface temperatures.

Figure 2 also contains the empirical upper luminosity boundary of stars in the Milky Way, as derived by Humphreys & Davidson (1994). As we see, our slowly rotating models do penetrate the Humphreys-Davidson (HD) limit and spend a significant amount of time at cooler temperatures. Whether this prediction is in contradiction with observations for the LMC is currently unclear. The stellar statistics near the HD limit for LMC stars from published work is not very good (Fitzpatrick & Garmany 1990). And even in the Milky Way, stars above the HD limit are observed, the most prominent example being ηCarinae with a luminosity of log  L/L_⊙ = 6.7 and T_eff ≃ 30 000 K (Smith 2013).

In any case, our models predict a short lived (~10⁵ yr) yellow or red supergiant phase of stars in the mass range 60 M_⊙ to 80 M_⊙ during which the core is still burning hydrogen. Our models obtain a stellar wind mass-loss rate of the order of 10^-4 M_⊙ yr^-1 during this stage. However, in this part of the HR diagram, mass-loss rates are very uncertain. We note that in particular higher mass-loss rates (which may be due to the pulsational instability of these models; Sanyal et al., in prep.) would lead to shorter life times in this evolutionary stage.

The higher the initial mass of our model, the higher is their convective core mass fraction. While it exceeds 80% in ZAMS stars above 100 M_⊙, it reaches 90% at 300 M_⊙, and 500 M_⊙ stars with a convective core mass fraction of 95% are almost fully convective. This puts the convective core boundary far out into the stellar envelope, where the pressure scale height is rather small. Consequently, convective core overshooting does not have the same importance as at smaller stellar mass and plays only a minor role in most of the models presented here.

The effect of rotation on the evolution of massive stars is discussed previously, for example in Maeder & Meynet (2000), Heger & Langer (2000), Brott et al. (2011a), Chieffi & Limongi (2013). In our models, there are two main changes in the stellar evolution tracks in the HR diagram. First, the effective gravity is reduced as a result of the centrifugal acceleration, which leads to a decrease in the effective temperature and luminosity of a star compared to a non-rotating model. Second, above a threshold rotational velocity for which the timescale of rotational mixing becomes smaller than the nuclear timescale, the stars evolve quasi-chemically homogeneously (Yoon & Langer 2005; Woosley & Heger 2006), and the corresponding models evolve directly towards the helium main sequence in the HR diagram (Brott et al. 2011a).

Both effects are clearly visible in our models. For stellar evolution tracks with initial rotational velocities below 250 km s^-1, the first effect mentioned is dominant. It can be recognized most easily from the ZAMS position of our stellar models (Figs. 2 and 3), where the models are dimmer and cooler the faster they rotate.

The most rapidly rotating stellar models undergo quasi-chemically homogeneous evolution. As displayed in Fig. 2 for ν_ZAMS = 400,500 km s^-1, those with initial masses below ~125 M_⊙ show a strong increase in luminosity and effective temperature. For more massive models, however, homogeneous evolution also leads to higher luminosities, but the surface temperature is decreased. Towards core helium exhaustion, strong Wolf-Rayet mass loss leads these stars to lower luminosities and larger surface temperatures (cf. Sect. 3.2). We return to the discussion of quasi-chemically homogeneous evolution in Sect. 3.4. Isochrones in the HR diagram based on tracks discussed here are presented in Appendix A.

3.2. Mass loss and surface rotational velocity

The evolution of the models for the most massive stars is strongly affected by mass loss. According to our mass-loss prescription (Sect. 2), our models undergo mass loss of three different strengths. Initially, the mass-loss rate proposed by Vink et al. (2000, 2001) is used. For our adopted composition, the mass-loss rate is initially of the order of 3 × 10^-6 M_⊙ yr^-1 for the 100 M_⊙ models, while for the 500 M_⊙ models it is ~6 × 10^-5 M_⊙ yr^-1. At this rate, only 10–20% of the initial mass will be lost over the lifetime of the stars (Fig. 4, right panel).

The models that do not experience quasi-homogeneous evolution undergo an increase in their mass-loss rate at effective temperatures of ~25 000 K, which is referred to as the bi-stability jump (Vink et al. 1999), according to the Vink et al. prescription (Fig. 4). The corresponding mass loss is so intense that the models above 100 M_⊙ stop evolving redward at this stage, since their surfaces become strongly helium-enriched. For our most massive models, the 500 M_⊙ sequences, the maximum mass-loss rate at this stage is about 2 × 10^-4 M_⊙ yr^-1.

Once the surface helium mass fraction has increased owing to mass loss to more than 40%, the Wolf-Rayet mass-loss rate is phased in, reaching its full strength at a surface helium mass fraction of 70%. As a consequence, our 100–500 M_⊙ models lose 50–80% of their initial mass before core hydrogen exhaustion. During the phase of Wolf-Rayet winds, the models evolve at decreasing luminosity and all accumulate in a narrow region in the HR diagram (Fig. 2).

The three different phases of mass loss can be clearly seen in both panels of Fig. 4. Again, for our most extreme models, the slowly rotating stars with an initial mass of ~500  M_⊙, the largest obtained mass-loss rate is ~5 × 10^-4 M_⊙ yr^-1, according to our Wolf-Rayet mass-loss prescription. At this moment, the stars have a luminosity of log  L/L_⊙ ≃ 7.1. Assuming a terminal wind speed of 1000 km s^-1, we compute a wind-momentum-to-photon-momentum ratio of η ≃ 2 for this situation. With the same numbers, we obtain a wind-kinetic-energy-to-luminosity ratio of 0.003. The corresponding wind darkening is therefore negligible (Heger & Langer 1996) and significantly smaller than in Galactic Wolf-Rayet stars (Lucy & Abbott 1993).

The models that evolve quasi-homogeneously first increase their luminosity more than the models described above. Consequently, their mass-loss rate becomes larger than that of the inhomogeneous models. The homogeneous models do not reach the bi-stability limit as the increased helium surface abundance due to rotational mixing keeps their surface temperature above 25 000 K. However, they can reach the Wolf-Rayet stage earlier, and in the end lose similar amounts of mass as the inhomogeneous models. During the Wolf-Rayet stage, there are no clear characteristics that can distinguish between the two types of evolution.

Figure 5 compares the mass evolution of the homogeneously and inhomogeneously evolving 60 M_⊙ to 80 M_⊙ models. It shows that in this mass range, the homogeneously evolving models lose more mass, with more than half of the mass being lost during the Wolf-Rayet phase. The inhomogeneously evolving models lose almost equal amounts during their evolution through the hot part of the HR diagram as at cool surface temperatures, where the latter mass loss occurs on a timescale of only 10⁵ yr.

As a consequence of the initially rather moderate stellar wind mass loss, most of our models evolve at a constant rotational velocity for most of the core hydrogen burning. In this case the loss of angular momentum through the stellar wind (Langer 1998) and the increased momentum of inertia due to expansion is compensated by the transport of angular momentum from the contracting core to the envelope. As shown in Fig. 6, this holds even for the most massive stars in our grid. Only the models which undergo chemically homogeneous evolution, as a result of their increased luminosity and mass-loss rate, show a moderate decline of their surface rotational velocity before entering the phase of Wolf-Rayet winds.

The strong mass loss for T_eff< 25 000 K for the inhomogeneously evolving models, and the Wolf-Rayet type mass loss for the homogeneously evolving ones both lead to a strong spin-down of the stars towards the end of their core-hydrogen burning evolution. In fact, this is only avoided for stars initially below ~30 M_⊙ (cf. Fig. 3 of Vink et al. 2010). Since the magnetic coupling in our models ensures close-to-rigid rotation during core hydrogen burning, this implies that all our LMC models above ~80 M_⊙ lose so much angular momentum that they cannot be considered to produce candidates for long-duration gamma-ray bursts. In that respect, models that undergo chemically homogeneous evolution of lower mass are better suited (cf. Sect. 3.4).

Effective temperatures and radii of observed massive stars may be affected by optically thick stellar winds, making the star appear larger and cooler. We use Eq. (14) from Langer (1989)(1)to estimate the optical depth of the stellar winds of our stellar models. Here, R designates the radius of the stellar model without taking the wind into account. This equation is derived from a β-velocity law with β = 1. In this case, we use the electron scattering opacity (κ = σ⁽1 + X⁾, where σ is the Thomson scattering cross section), an expansion velocity ν₀ = 20 km s^-1 at the surface of the stellar model, and a terminal wind velocity of ν_∞ = 2000 km s^-1.

Figure 7 shows the estimated optical depth of stellar winds as a function of time for several model sequences. The behaviour of the optical depths seen in the figure is mostly related to the change in the mass-loss rate. The optical depth increases for higher initial mass and higher rotation rate as a result of a corresponding increase in the mass-loss rate. While these numbers are only approximations, it shows that the winds of the most massive stars might already be optically thick (τ> 1) on the ZAMS, which implies that spectroscopically, these stars may already show Wolf-Rayet characteristics at this time (see also Gräfener & Hamann 2008; Crowther et al. 2010). Furthermore, while the stars of 100  M_⊙ and below are expected to show optically thin winds for most of the core hydrogen burning evolution, the optical depths of the winds for models for which a Wolf-Rayet type wind has been assumed may be quite large.

3.3. Surface abundance evolution

Deriving masses or ages of stars using stellar evolution tracks or isochrones in the HR diagram is no longer straightforward when rotational mixing is important. A given pair of effective temperature and luminosity can be explained by several combinations of masses, rotational velocities, and ages. To uniquely determine the initial parameters and the age of a star, it is necessary to include additional observables. The mass fractions or abundances of elements at the stellar surface are tracers of the rotational mixing. They can be used to uniquely identify the initial parameters and the age of a star in a three (or more) dimensional space of observables (cf. Köhler et al. 2012). We therefore discuss the influence of time, initial rotational velocity and initial mass on the surface composition of our models.

The rotational mixing in our models is sensitive to the gradient of the mean molecular weight μ. Essentially, any significant μ-barrier prevents rotational mixing. Therefore, mixing of helium can only occur in our models as long as they are quasi-chemically homogeneous. This, in turn, requires the timescale for rotational mixing to be shorter than the nuclear timescale, as the star attempts to establish a μ-barrier on the nuclear timescale.

In the fastest rotators, e.g. in the 100 M_⊙ and 300 M_⊙ models initially rotating with ~400 and 500 km s^-1, the surface helium abundance goes almost all the way to Y_s = 1 towards core hydrogen exhaustion (cf. Fig. 8). However, the same figure shows that this evolution is truncated for the rapidly rotating 60  M_⊙ models, at Y_s = 0.68 and Y_s = 0.80 for initial rotational velocities of ~400 and 500 km s^-1, respectively. These models spin down during core hydrogen burning (Fig. 6) such that the timescale for rotational mixing eventually becomes larger than the nuclear timescale.

Figure 8 shows that the non-rotating 300 M_⊙ sequence also evolves to Y_s = 1. The reason is the large convective core fraction of this model, and its large mass-loss rate. Whereas the model always keeps a small radiative envelope mass, the equivalent amount of mass is lost on a short timescale such that core and surface abundances become almost equal (cf. Sect. 3.4; Eq. (2)). For the rapidly rotating 100 M_⊙ and 300 M_⊙ models, mass loss rather than rotational mixing must take over to enforce chemical homogeneity, since these models also spin down quite dramatically (Fig. 6).

Figure 9 (left panel) demonstrates the role of rotational mixing of helium in greater detail using the example of 80  M_⊙ sequences. In the models initially rotating with 0, 50, 100, 150, and 200 km s^-1, such mixing is essentially negligible, and only mass loss can increase the helium surface abundance close to core hydrogen exhaustion, which is why all of these models evolve along the same track in the HR diagram (Fig. 3). The models between 250 and 400 km s^-1 start out with chemically homogeneous evolution (though the one with ν_rot,i = 250 km s^-1 only for a short amount of time), and truncate the homogeneous evolution later in time for higher initial rotation. The models which initially rotate with about 450 and 500 km s^-1 undergo chemically homogeneous evolution all the way (though with the help of mass loss in the end), and their evolution in the HR diagram is also practically identical (Fig. 3).

The evolution of the surface nitrogen abundance in our models of rotating and mass losing stars follows different rules. During a phase of chemically homogeneous evolution, nitrogen is quickly mixed from the core to the surface establishing the CNO-equilibrium value in atmospheric layers (Figs. 8 and 9, left panels). Mass loss can also lead to nitrogen enhancements, even without rotational mixing (Fig. 8). Figure 9 shows that a strong nitrogen surface enhancement can be obtained even in models with initial rotational velocities well below the threshold value required for chemically homogeneous evolution. For example, the 80  M_⊙ models with ν_rot,i = 150 and 200 km s^-1 enrich nitrogen at the stellar surface by factors of 4 and 12 by rotation, i.e. before mass loss kicks in. This order of magnitude of nitrogen enrichment is as expected, since our models were calibrated to increases the surface nitrogen abundance by a factor of about 3 for stars of 13–15  M_⊙ for an initial rotational velocity of 150 km s^-1, and rotational mixing is stronger in more massive stars.

The reason is that nitrogen, carbon, and oxygen, with a maximum mass fraction of less than one per cent, are just trace elements, such that even significant internal gradients do not lead to strong μ-barriers. Therefore, at the beginning of core hydrogen burning, before significant amounts of hydrogen have been converted to helium, but after the CNO-cycle has already strongly enhanced the nitrogen abundance in the convective core, rotational mixing can bring nitrogen out of the core into the radiative envelope and later on also to the surface.

The surface enrichment of helium and nitrogen explained above allows a simple understanding of the occurrence of surface enrichments in the HR diagram. As shown in Fig. 10, stars initially rotating with rotational velocities of 200 km s^-1 or less that are less luminous than log  L/L_⊙ ≃ 6.2 are not expected to show any helium surface enrichment (i.e. Y_s< 0.28) as long as their surface temperature is higher than ~20 000 K. The same is true for stars initially rotating slower than 300 km s^-1 below log  L/L_⊙ ≃ 5.6, and for those initially rotating slower than 400 km s^-1 below log  L/L_⊙ ≃ 4.8. Similar higher thresholds can be read off from Fig. 10 for larger surface helium mass fractions. Given that Ramírez-Agudelo et al. (2013) found that 75% of all O stars rotate slower than 200 km s^-1 and that amongst the 31 O2 to O5 stars none was found to rotate faster than ~300 km s^-1, the helium enrichment in LMC stars below ~100 M_⊙ is expected to be quite small during the main-sequence phase.

In contrast, as more than half of all O stars were found to rotate faster than 100 km s^-1, nitrogen enrichment by at least a factor of 2 is expected to be almost ubiquitous above log  L/L_⊙ ≃ 6.0, and quite frequent in main-sequence O stars in general (Fig. 11).

Changes in the amount of neon, sodium, aluminum, and magnesium at the surface also occur, indicating that in very massive stars the MgAl- and the NeNa-cycles are active. In our models, the amount of magnesium is reduced while aluminum is produced. The mass fraction of neon increases in our models, while sodium decreases. For more details, we refer to the electronic data published with this paper.

Lithium, beryllium, and boron are destroyed in the stellar interior (McWilliam & Rauch 2004) at temperatures above ~2.5 × 10⁶ K for lithium, ~3.5 × 10⁶ K for beryllium, and ~5.0 × 10⁶ K for boron. In our models, temperatures below 5 × 10⁶ K are only found in the outer envelope. Therefore lithium, beryllium, and boron can exist in the outer envelope, while they are destroyed in deeper layers. For the slowest rotating models, the surface abundances of these three elements remain constant until layers which were exposed to higher temperatures earlier in the evolution are exposed to the surface due to mass loss. When that happens, the mass fractions of beryllium, boron, and lithium are quickly reduced. For fast rotating models, rotational mixing and mass loss lead to a gradual decrease in lithium, beryllium, and boron over time.

3.4. Chemically homogeneous evolution

As mentioned in Sect. 3.1, we can divide our models into three classes. The first one (which we call Class O) likely corresponds to most stars in observed stellar samples and contains the models that evolve in the normal way, i.e. in which the rotationally induced mixing of helium is negligible. The second one describes the models that undergo quasi-chemically homogeneous evolution (Class H), which correspond to the initially fastest rotators. As the third class, we have the models that start out evolving homogeneously, but which spin down such that the rotational mixing of helium stops (Class HO), after which they evolve redward in the HR diagram as the ordinary models. These models retain a memory of their past homogeneous evolution, in that they will keep an enhanced helium surface abundance and a higher luminosity-to-mass ratio (Langer 1992) compared to ordinary stars throughout the rest of their core hydrogen burning evolution.

At the highest masses considered here (M ≳ 150 M_⊙), our models may also undergo quasi-chemically homogeneous evolution without rotationally induced mixing, which is due to a combination of an extremely high fraction of the convective core to the total stellar mass and a very high Wolf-Rayet type mass-loss rate. We illustrate this by the example of our 200  M_⊙ sequence with v_ZAMS ≃ 300 km s^-1 in Fig. 12. For such very massive stars, the convective core comprises a major fraction of the stellar mass.

We can divide the core hydrogen burning evolution of this model into three parts according to the mass evolution of the convective core. During the first part, the mass of the convective core decreases as a function of time, but more slowly than the mass of the star, such that the mass of the radiative envelope decreases and the convective core mass fraction increases. During this phase, the model undergoes quasi-chemically homogeneous evolution (cf. Fig. 13). At t ≃ 1.2 Myr, the model transitions to ordinary evolution, from which time on the convective core mass decreases somewhat faster than the total mass. During this second part of its evolution, the helium surface abundance increases, even though at a much lower rate than the central helium abundance (Fig. 13). Finally, in part three, the Wolf-Rayet mass loss kicks in, which leads to a very small mass of the non-convective stellar envelope (M_env ≃ 5 M_⊙). The amount of mass corresponding to the non-convective envelope is lost (at a rate of ~2 × 10^-4 M_⊙ yr^-1) on a timescale of ~20 000 yr which corresponds to only about 1% of the core hydrogen burning time, τ_H. Consequently, the surface helium mass fraction during this stage is roughly equal to (within 1%) the central helium mass fraction. Thus, regardless of rotation, the model is extremely chemically homogeneous during this phase.

We define a fourth type of evolution (Type M) by (2)where τ_H is the hydrogen burning timescale. Since this condition ensures that the surface-to-central helium abundance ratio is above ~0.9 – i.e. similar to or smaller than in the case of rotationally induced quasi-chemically homogeneous evolution – we can subdivide our models into further classes. For example, the 200  M_⊙ sequence discussed here can be considered of Class HOM, as it first undergoes chemically homogeneous evolution due to rotational mixing, then ordinarily, and finally chemically homogeneous evolution due to mass loss.

It is instructive to consider Fig. 14 to understand which evolutionary classes are realized by our models. The non-rotating sequences with initial masses of 300  M_⊙ and 500  M_⊙ achieve chemical homogeneity at a central helium mass fraction of 0.85 and 0.75, respectively, and thus belong to Class OM. The non-rotating 100  M_⊙ sequence does not reach homogeneity. Consequently, although it achieves a quite high final surface helium mass fraction of ~0.65, it belongs to Class O. The most rapidly rotating models depicted in Fig. 14 never show a significant discrepancy between their central and surface helium abundances, which implies that they evolve from H-type to M-type evolution and belong to Class HM. The stars rotating initially with ~300 km s^-1 all undergo an H→O transition, where the two more massive ones move further on to Type M evolution. As stars in the Classes HOM and HM both start out and end their chemically homogeneous evolution, we will not distinguish them further and designate them both as Class HM.

Thus, we have five different classes of evolution, H, HM, HO, OM, and O, where all but the last involve quasi-chemically homogeneous evolution. Class O can be identified in Fig. 15 (right panel), which shows the close-to-final surface helium mass fraction for all our model sequences. Looking at the slowly rotating models, this figure shows that mass loss starts affecting the surface helium abundance above an initial mass of ~65 M_⊙ (log M/M_⊙ ≃ 1.8), and that M-type chemically homogeneous evolution is obtained above ~160 M_⊙ (log M/M_⊙ ≃ 2.2). According to our definition, the slow rotators below ~160 M_⊙ are thus in Class O, the more massive ones in Class OM. For the fast rotators, the dividing line between homogeneous and inhomogeneous evolution depends on stellar mass. At 20  M_⊙ (log M/M_⊙ ≃ 1.3), stars rotating initially faster than ~500 km s^-1 evolve homogeneously (Class H), while at ~160 M_⊙ the critical velocity is at ~350 km s^-1, and more massive stars above this velocity are in Class HM. Below the dividing line defining Class H, between 20 M_⊙ and 160 M_⊙ is a stretch of Class HO models, which comprises a larger initial mass range for higher initial mass.

As mentioned in Sect. 3.1, rapidly rotating stellar models with initial masses greater than 125  M_⊙ show significantly different behaviour in the HR diagram (see Fig. 2) to the stellar models with M ≤ 125 M_⊙ (see also Sect. 3.6). Nevertheless, they behave in a similar way in the Kippenhahn diagram and when comparing the helium mass fraction in the core and at the surface.

Stellar evolution models above 80  M_⊙ have strong mass loss. Even the fastest rotators are slowed down significantly when the Hamann et al. (1995) mass-loss rate is applied. For models above M_ZAMS ≃ 100 M_⊙, the rotation rate at 95% helium mass fraction in the core is below 100 km s^-1. The surface rotational velocity goes down to less than 50 km s^-1 for models above 150  M_⊙, independent of the initial surface rotational velocity (Fig. 15; right panel).

From the above, we can come to several conclusions which are relevant in comparison to observed stars. First, helium-enriched single stars below ~65 M_⊙ are prime suspects of H- or HO-type chemically homogeneous evolution, whereas M-type evolution can be excluded for them. These stars are expected to preserve their rapid rotation throughout their core hydrogen burning evolution. As a consequence, we may expect some correlation of the helium surface abundance with the stellar rotation rate in the considered mass regime.

For very massive stars (above ~200 M_⊙) rotation, whether initially fast or not, makes little difference, i.e. the evolution becomes almost independent of the initial rotational velocity. These stars are efficiently spun down by mass loss, and all undergo chemically homogeneous evolution during their advanced hydrogen-burning evolution. Their surface helium abundance is always very close to their central helium abundance because of a rapid loss of the thin outer radiative envelope through mass loss. As a consequence, these sequences can produce models that have a very high surface helium mass fraction (Y_s> 0.9) and are still undergoing core-hydrogen burning.

3.5. Mass-luminosity relation

Stars of higher mass are increasingly luminous. The most massive models in our grids, at 500  M_⊙, radiate at more than 10⁷ L_⊙. Figure 16 shows the evolution of selected non-rotating and of rapidly rotating model sequences in the mass-luminosity plane, from core hydrogen ignition to core hydrogen exhaustion. At the lowest considered masses (15 M_⊙), the models evolve vertically upward since their mass loss is negligible. The more massive, rapidly rotating models evolve to higher luminosity, and they turn left toward lower mass.

A comparison with the mass-luminosity relations for chemically homogeneous stars of Gräfener et al. (2011) reveals that, except for the depicted 15 M_⊙ model, the rapidly rotating sequences shown in Fig. 16 undergo quasi-chemically homogeneous evolution. They start at the mass-luminosity (ML) relation for ZAMS stars and finish close to the ML relation for chemically homogeneous helium stars, indicating their nearly homogeneous chemical structure. This also demonstrates that the agreement of the ML relations of Gräfener et al. (2011, which represent fits to stellar models computed with a different code) with our models is very good.

We consider a mass-luminosity relation for chemically homogeneous stars as L ~ M^αμ^β, where μ is the mean molecular weight of the stellar gas. The power-law exponent α thus describes the slope of the ML relation in the log M − log L-plane. According to Kippenhahn & Weigert (1990), we have α> 1, and α → 1 for M → ∞, due to the increasing radiation pressure for higher masses (P_gas/ (P_gas + P_rad) → 0). The inequality α> 1 implies that, in the frame of simple opacity laws, massive main sequence stars will never exceed the Eddington limit. The limit of α → 1 implies that for sufficiently massive stars, their hydrogen burning life time τ_H = E/Ė becomes independent of mass, since their nuclear energy reservoir E is proportional to their mass, just as their energy loss rate Ė = L. Figure 16 shows that our most massive models reach α-values very close to one. We find the limiting stellar life time to be close to 1.9 Myr (cf. Table B1).

Figure 17 also shows that, at a given mass, the slope of the mass-luminosity relation is smaller for the homogeneous hydrogen-free stars, compared to that for the ZAMS stars because the helium stars are much more luminous and thus more radiation pressure dominated. Furthermore, the power-law exponent β also approaches one for the highest considered masses. While a 15 M_⊙ helium star is about 20 times as luminous as a 15 M_⊙ ZAMS star (, and μ_tams/μ_zams ≃ 2.23), the corresponding factor is only 2.4 at 500 M_⊙.

3.6. Near the Eddington limit

In its general form, the Eddington limit is complex as it involves an appropriate mean of the total of all opacity sources. This total opacity κ may be a function of depth in the stellar atmosphere; therefore, it is non-trivial to uniquely define the circumstances in which the star encounters its Eddington limit, nor to assess the consequences such an encounter may have. The proximity of a star of mass M and luminosity L to the Eddington limit is usually expressed in terms of the Eddington factor (3)where the constants have their usual meaning and T_eff is the effective temperature and g the surface gravity.

Considering only the dominant contributor to the opacity, i.e. photon scattering on free electrons, greatly simplifies the concept of Eddington limit as electron scattering is a grey process and, for the hot stars considered here, independent of depth in the atmosphere. For a star of mass M and luminosity L we obtain Γ_e = κ_eL/ (4πcGM), where κ_e = σ_e(1 + X) is the opacity due to Thomson scattering. It holds that Γ > Γ_e as κ>κ_e. Even in the most massive stars, the Eddington factor Γ_e will not exceed unity (see Sect. 3.5). Even so, since the mass-luminosity exponent α> 1, main sequence stars will get ever closer to the Eddington limit for electron scattering the higher their mass.

Considering the true Eddington factor, we must expect that the situation Γ = 1 is actually achieved in very massive stars. Indeed, while κ_e ≃ 0.34 in our zero-age models, their true surface opacity is typically κ_surface ≃ 0.5. We can thus expect that Γ = 1 is achieved at Γ_e ≃ 0.7.

In Fig. 18, we plot the quantity as a function of the effective temperature of selected evolutionary sequences. We note that L is normalized to the solar value, log L_⊙ ≃ 10.61 for convenience. For instance, for the initial composition of our stars we have log L/ L_⊙ ≃ 4.6 + log Γ_e. The quantity L is proportional to Γ_e and can be derived from spectroscopic observation of stars without knowledge of their distance (cf. Langer & Kudritzki 2014).

Since for stars of constant mass L~ L, the quantity L behaves in a similar way to the stellar luminosity, and the evolutionary tracks in the L− T_eff-diagram partly resemble those in the HR diagram. However, unlike the latter, the L− T_eff-diagram has an impenetrable upper limit: the Eddington limit.

Comparing Fig. 18 with Fig. 2, we see that while the tracks of the slow rotators in Fig. 2 are almost horizontal, the corresponding tracks in Fig. 18 are significantly steeper. Their luminosities are nearly constant despite strong mass loss, because their increasing mean molecular weight μ compensates for their decreasing mass in the mass-luminosity relation (Sect. 3.5). While it remains hidden in Figs. 2, and 18 reveals nicely that the mass loss drives these stars towards the Eddington limit. The same is true for the rapidly rotating stars; their tracks move steeply upward in both diagrams.

In Fig. 18, we have drawn the horizontal line indicating Γ_e = 0.7 for our initial chemical composition – which corresponds to log L/ L_⊙ ≃ 4.445. From the discussion above, we should expect that stellar models with unchanged surface composition should not be found above this line (as long as their hydrogen and helium can be considered fully ionized). This is confirmed by the triangles placed on the evolutionary tracks, indicating the position at which a surface helium mass fraction of Y = 0.3 is achieved, and which are all found below Γ_e = 0.7. Similarly, even the models with the most helium-enriched surfaces stay below the upper horizontal line for Γ_e(X = 0) = 0.7; in other words, the diagram shows that it is not the electron-scattering Eddington limit which constrains the evolution of the most massive stars, but the true Eddington limit. Furthermore, Fig. 18 indicates that the limiting electron-scattering Eddington factor – which can be read off from Fig. 18 to be a little bit below Γ_e = 0.7 – is independent of the hydrogen/helium surface abundances.

We find our models near the Eddington limit to undergo a significant envelope inflation. This phenomenon has been described before by Kato (1986) and Ishii et al. (1999) for ZAMS stars of different metallicities, and by Ishii et al. (1999), Petrovic et al. (2006), Gräfener et al. (2012) for Wolf-Rayet stars. Figure 19 shows the location of our ZAMS stars, and of zero-age helium stars with the same metallicity, in the HR diagram. Both curves show a maximum effective temperature, which is reached at about 30 M_⊙ (140 000 K) for the helium stars, and at 200 M_⊙ (54 000 K) for the hydrogen-rich stars. Above these masses, the effective temperature decreases again owing to the envelope inflation.

As seen in Fig. 19, the line describing the helium stars bends sharply to cooler temperatures at about 1 000 000  L_⊙, and is even expected to cross the hydrogen-ZAMS (see also Ishii et al. 1999), implying that more luminous helium stars are cooler than equally bright ZAMS stars. This has the consequence that stars undergoing chemically homogeneous evolution above ~125 M_⊙ evolve to cooler surface temperature (see Fig. 2).

The envelope inflation occurs because layers inside the stellar envelope reach or exceed the true Eddington limit. The inflation begins before the critical value of Γ_e ≃ 0.7 is reached, i.e. the electron-scattering Eddington factors are only Γ_e ≃ 0.31 and Γ_e ≃ 0.41 for the 15 M_⊙ helium star and 190 M_⊙ ZAMS star, respectively, because of the opacity peaks of iron and helium which can reach values of 1 cm²/g or more. While the subsurface layers with such opacities always turn convectively unstable, convection, according to the mixing-length theory, is mostly very inefficient in transporting energy in these layers. Consequently, the high radiation pressure pushes the overlying layers outwards, until the opacity drops such that the true Eddington factor falls short of Γ = 1.

Many of these models also develop density inversions near the layer with the maximum Eddington factor (Joss et al. 1973). The corresponding inward directed gas pressure gradient then allows the star to retain layers whose Eddington factor exceed the critical value of Γ = 1. While the stability of these structures remains to be investigated, any instability which would tend to iron out the density inversion will likely lead to a further inflation of the overlying layers. A further degree of complexity is added by the fact that we find many models close to the Eddington limit to be unstable at least to radial pulsations (the only pulsation mode which we can see with our one-dimensional hydrodynamic stellar evolution code), and that convective velocities close to the sound speed may imply significant acoustic fluxes (Goldreich & Kumar 1990). A deeper investigation of these phenomena exceeds the scope of the present paper and will be pursued elsewhere.

4. Comparison with previous results

Recently, Yusof et al. (2013) published a grid of stellar evolution models for very massive stars. Eight of their evolutionary sequences were modelled using a similar initial composition to the one we used here, except that Yusof et al. used a scaled solar metal mix, while we take the measured LMC composition where possible (cf. Brott et al. 2011a). Yusof et al. present non-rotating models with initial masses of 120  M_⊙, 150  M_⊙, and 500  M_⊙ as well as rotating models at 120  M_⊙, 150  M_⊙, 200  M_⊙, 300  M_⊙, and 500  M_⊙. The rotating models of Yusof et al. had initial rotational velocities of 400 km s^-1 or more (cf. their Table 2).

4.1. Subsurface convection and envelope inflation

A major difference between the models of Yusof et al. and our models occurs because of the different treatment of convection in the stellar envelopes. Starting at zero age, all our models contain subsurface convection zones which occur because of the opacity peaks of iron and helium (cf. Cantiello et al. 2009). When using the standard mixing-length theory, where the mixing length is assumed to be proportional to the pressure scale height, we find the convective energy transport in these zones to be mostly inefficient because of the small local thermal timescales imposed by the low densities. As a consequence, we obtain a rich envelope phenomenology as described in Sect. 3.6, including envelope inflation and density inversions.

Yusof et al. describe this situation as unphysical, and avoid it by assuming the convective mixing length to be proportional to the density scale height. As a consequence, when the star attempts to establish a density inversion, the density scale height tends to infinity, which imposes the convective energy transport to be extremely efficient, and opacity peaks will have few consequences. While the existence of inflated envelopes in nature remains to be shown (cf. Gräfener et al. 2012; Gräfener & Vink 2013), we consider it likely, since the short thermal timescales of the envelopes of very massive stars cannot be in accordance with efficient convective energy transport. While this will surely be studied in more detail in the near future, we restrict ourselves here to pointing out the consequences of the different assumptions on the convective efficiency.

Strong differences between the models of Yusof et al. and our models already occur at the ZAMS. As described in Sect. 3.1, our zero-age models (with our metallicity of Z = 0.0047) reach a maximum effective temperature of about 56 000 K at ~190 M_⊙, which agrees very well with the results of Ishii et al. (1999) for Z = 0.004. In Fig. 20, we compare the location of the rotating zero-age models of Yusof et al. (according to their Table 2) with the early hydrogen burning evolution of our models with an initial rotational velocity of ~300 km s^-1. Above ~120 M_⊙, the models of Yusof et al. are hotter than ours, and the more so the larger the mass.

This discrepancy is likely explained by the difference in the treatment of the envelope convection, since our treatment yields a larger radius inflation the closer the star is to its Eddington limit (Sect. 3.6). At a given initial mass, the discrepancy between the models of Yusof et al. and our models is therefore expected to increase with time at least initially, since the models increase their luminosity and decrease their mass (cf. Fig. 16), and thus their Eddington factors grow. The biggest consequence of this is that – although we use a very similar mass-loss prescription to that of Yusof et al. – our mass-loss rates are higher because our stars have cooler surfaces.

4.2. M_ZAMS versus M_TAMS relation

Figure 21 shows that as a result of the different treatment of envelope convection in our models and those of Yusof et al., our non-rotating models above ~150 M_⊙ undergo a dramatic mass loss during core hydrogen burning, leading to masses at core hydrogen exhaustion of the order of 80 M_⊙ even up to the highest initial mass of 500 M_⊙. Yusof et al. did not consider non-rotating models between 150 M_⊙ and 500 M_⊙. However, their non-rotating models below 150 M_⊙ and at 500 M_⊙ remain close to the diagonal in Fig. 21. While our 500 M_⊙ model ends hydrogen burning with a mass of 87 M_⊙, theirs does so with 239 M_⊙. We discuss the consequences for the final fate of the stars in Sect. 4. Remarkably, because of the low mass-luminosity exponent of the most massive stars considered here, the core hydrogen burning life time of the 500 M_⊙ sequence of Yusof et al. (1.9 Myr) resembles ours, indicating that the life time of the most massive stars is not strongly affected by the differences discussed above.

The non-rotating 500 M_⊙ sequence of Yusof et al. ends core hydrogen burning as a nearly chemically homogeneous helium star with a surface helium mass fraction of 0.97. Our Fig. 19 above shows that at the considered metallicity, our helium star models start to have inflated envelopes above ~30 M_⊙ (see also Ishii et al. 1999; Petrovic et al. 2006), and that a 239 M_⊙ helium star model would be so enormously inflated that it would appear as a red supergiant, if a stable envelope structure existed at all. This demonstrates that while the envelope physics used in our models leads to pushing the outer layers of stars at the Eddington limit to large radii, this seems to occur to a much lesser degree in the models of Yusof et al., with the consequence of a smaller mass-loss rate.

When looking at the Eddington factor for electron scattering Γ_e, the models of Yusof et al. behave as our models do, i.e. they remain mostly lower than Γ_e = 0.7 (cf. Sect. 3.6). However, their non-rotating 500 M_⊙ sequence is an exception, where Γ_e = 0.82 is reached at core hydrogen exhaustion. In comparison to the stellar wind calculations of Vink et al. (2011), who found a strong mass-loss enhancement for stars above Γ_e = 0.7, we would again expect that the mass at core hydrogen exhaustion of this sequence could be considerably lower than that found by Yusof et al.

Our models show the largest mass-loss rates during their Wolf-Rayet phase (cf. Sect. 3.2). Nevertheless, it is not the choice of the Wolf-Rayet mass-loss prescription which makes our non-rotating models lose more mass than the comparable models of Yusof et al. We find, for example, that for our slowly rotating 500 M_⊙ models at the time of the maximum mass-loss rate, the prescription by Nugis & Lamers (2000) leads to a mass-loss rate for our stars which is two times larger than the one we use. We conclude that the high pre-Wolf-Rayet mass loss imposed by the inflated envelopes of our models plays a crucial role in explaining the smaller TAMS-mass of our very massive slow rotators compared to those of Yusof et al.

4.3. Rotational mixing of helium

Another difference between the models of Yusof et al. and our models concerns the rotational mixing of helium. As shown in Sect. 3.3, the mixing of helium in our models is practically absent below a threshold rotational velocity, and nearly complete above the threshold. The mixed models suffer more mass loss since they become Wolf-Rayet stars earlier. As the threshold velocity decreases for higher masses, this results in a bimodal behaviour of the total mass lost during core hydrogen burning (Fig. 21). Our models with an initial rotational velocity of ~300 km s^-1 follow the non-rotating models below ~150 M_⊙, while they behave in a similar way to our fastest rotators for higher mass. In the models of Yusof et al., the mixing of helium is gradually increased for larger rotational velocities. Therefore, their rotating models end core hydrogen burning with masses that exceed those of our homogeneously evolving models.

A more detailed comparison with the models of Yusof et al. is difficult because they considered only two rotational velocities, i.e. ν_rot,i = 0 km s^-1 and ν_rot,i ≃ 400 km s^-1. This may also prevent a clear discrimination of their models from ours through observations, since Ramírez-Agudelo et al. (2013) showed that only very few O stars are expected with ν_rot,i = 0 km s^-1 and ν_rot,i ≃ 400 km s^-1, while the majority rotate with ν_rot,i = 100...200 km s^-1.

5. Summary

We present a detailed grid of stellar evolution models for single stars with initial masses from 70 to 500  M_⊙ and rotation rates up to 550 km s^-1. We used the same physics and assumptions as Brott et al. (2011a) did for stars in the mass range 5–60 M_⊙, and this new grid therefore is an extension of their work. The initial composition of our models corresponds to abundance measurements for massive stars in the LMC. We follow the evolution of the stellar models through their core hydrogen burning phase, with some of them computed well beyond this stage.

Given the high fraction of close binaries in massive stars (Sana et al. 2012), we cannot hope to obtain a complete picture of massive star evolution from our models (de Mink et al. 2014), which is true in particular for the most massive stars (Schneider et al. 2014). On the other hand, our rotating models may constitute a fair approximation for those mass gainers and mergers in close binaries which rejuvenate after the binary interaction (Braun & Langer 1995), even though their age, their detailed surface abundances, and their spin may still be peculiar (de Mink et al. 2013). In addition, since most mass donors in post-interaction close binaries either remain unobserved (de Mink et al. 2014) or show extreme surface enrichments (Langer 2012), and as Sana et al. (2012, 2013a) suggest that 30–50% of the massive stars in 30 Dor appear not to have a close companion, a comparison of our single-star models with observed very massive stars will still be meaningful.

We find that stellar rotation can influence the evolution of our models significantly, mostly through mixing of helium. We find two threshold initial rotational velocities, which both decrease with increasing initial mass, one below which the rotational mixing of helium is negligible, and a second one above which the models undergo quasi-chemically homogeneous evolution. For initial rotational velocities in between the threshold values, the models start chemically homogeneously and develop helium-enriched surfaces, but transit to normal evolution as a result of spin-down.

Above an initial mass of ~160 M_⊙, we find that quasi-chemically homogeneous evolution can also be achieved through mass loss, which in fact ensures that all our more massive models end core-hydrogen burning as nearly pure helium stars, independent of their initial rotation rate. We find that single stars with initial rotational velocities below 300 km s^-1 need to be more massive than ~100 M_⊙ in order to achieve a surface helium mass fraction above 35% (cf. Fig. 10).

Significant mixing of trace elements, including nitrogen and boron, is found well below the threshold velocities mentioned above. Because of the high core temperatures in our very massive models, the NeNa- and MgAl-cycles are also activated and lead to surface abundance changes of the involved isotopes. A tripling of the surface nitrogen abundance occurs in our models above an initial rotational velocity of 150 km s^-1 and a mass of ~25 M_⊙, while a doubling at this rotation rate occurs even down to 6 M_⊙ (Fig. 11). Because of the strong mass loss and large convective core masses, the surfaces of our models above ~100 M_⊙ become nitrogen-rich irrespective of their rotation rate.

Our zero-age models above ~150 M_⊙, as well as more evolved models down to ~40 M_⊙ show significant envelope inflation because the true Eddington factor inside the stellar envelope approaches or even exceeds the value of unity, while the electron-scattering Eddington factor remains below or near a value of Γ_e = 0.7 (cf. Sect. 3.6). Consequently, our ZAMS models show a maximum surface temperature of 56 000 K at about 180 M_⊙, with lower values for higher masses (Fig. 2).

During core hydrogen burning, the envelope inflation drives our slow rotators in the mass range 50–100 M_⊙ to effective temperatures near or below 10 000 K. At these cool temperatures, partial hydrogen recombination leads to high opacities, which in turn produces true Eddington factors well above one in the envelopes of these models. While a detailed study of this phenomenon remains to be done, we speculate that this feature is related to the outbursts of luminous blue variables.

The inflation and the corresponding decrease of the surface temperature lead to enhanced mass loss compared to the models of Yusof et al. (2013). The enhanced mass loss has strong consequences for the final fate of very massive stars. In this context it is important to emphasize the metallicity dependence of inflation (Ishii et al. 1999). The reduced masses (Fig. 21) exclude very massive stars of the considered metallicity as progenitors of pair-instability supernovae. This is in agreement with Langer et al. (2007), who demonstrated that models with the same envelope physics as used here may produce pair instability supernovae at a metallicity of Z_⊙/ 20 (Kozyreva et al. 2014). This is in contrast to the finding of Yusof et al., who suggest that slow rotators above ~300 M_⊙ end up with CO-cores above 60 M_⊙ at a metallicity of that of the LMC.

A further consequence of the strong mass loss of our inflated models is their distinct spin-down during core hydrogen burning (Fig. 6). This prevents all our models above ~60 M_⊙ from producing a long-duration gamma-ray burst upon collapse. From the models considered here and including those from Brott et al. (2011a), the highest chance to produce long-duration gamma-ray burst at LMC metallicity from the chemically homogeneous evolution scenario (Yoon & Langer 2005; Woosley & Heger 2006) may occur in the mass range 20–30 M_⊙ (Langer 2012).

The VLT-FLAMES Tarantula Survey (Evans et al. 2011) provides one of the major motivations for this study. Our models will be extensively compared with these and other observations in the near future. To this end, we provide the main outputs of our stellar evolution models as electronic tables (stored at CDS) in the same format as those of Brott et al. (2011a). Additionally, we provide sets of isochrones in Appendix A.

Acknowledgments

We want to thank Durand D’Souza for fruitful discussions, and Jesús Maíz Apellániz, Matteo Cantiello, Joachim Puls, Jon Sundqvist and Paco Najarro for useful comments. S.d.M. acknowledges support for this work by NASA through an Einstein Fellowship grant, PF3-140105. We are grateful to the referee for improving this paper.

References

Online material

Appendix A: Isochrones

The evolution of our rapidly rotating stellar models is influenced significantly by rotation (see Sect. 3.1). Therefore, isochrones of rotating models calculated for the same age differ from isochrones of non-rotating models.

Figure A.1 depicts isochrones of non-rotating models for ages up to 5.4 Myr in the HR diagram. Eight different isochrones are shown with age steps of 0.8 Myr. The core hydrogen burning stellar evolution models presented in Sect. 3 are used to generate the isochrones. Because more massive stellar models have shorter lifetimes, older isochrones terminate at the less massive model at the upper end of the track.

For a given initial composition and age, isochrones span an area in the HR diagram when different initial surface rotational velocities are considered simultaneously (Brott et al. 2011a). Isochrones of 16 different ages from 0.2 to 6.2 Myr for rotating stellar models are shown in Figs. A.2 and A.3 for different ages. Switching from inhomogeneous to chemically homogeneous evolution, the isochrone is located at higher effective temperatures and luminosities than for the non-chemically homogeneous evolution.

The more massive a stellar model is, the earlier it reaches the point in the HR diagram where the stellar evolution track starts to evolve blueward. The related turn in the isochrones in Fig. A.2 is first visible for 1.4 Myr (bottom left panel). The minimum effective temperature depends on the surface rotational velocity at the ZAMS, which is reflected by the isochrones of different rotation rates. The 2.2 Myr and 2.6 Myr isochrones show that the more massive a stellar model is, the earlier the Hamann et al. (1995) mass-loss rate is applied to the stellar evolution calculation. The isochrones therefore show a decrease in luminosity for the most massive models.

Figure A.4 depicts a population synthesis of stars with the age and distributions of mass and surface rotational velocity given in Table A.1. The calculation was done using the code Starmaker (Brott et al. 2011b), using the parameters listed in Table A1.

It can be seen that the randomly drawn stellar models with the same age do not lie on one line, but instead spread over a certain area in the HR diagram. The initial mass distribution determines the model density along the line corresponding to the isochrone of non-rotating stellar models of this age. We choose a flat mass distribution for the simulation to have a better view of the behaviour of the most massive stars. The velocity distribution on the other hand shapes the stellar model density as a function of the effective temperature for (roughly) constant luminosity. The colour coding indicates the number of stars within one pixel of 500 K and log  (L/L_⊙) = 0.05.

Figure A.4 shows similar information as discussed in Fig. A.2. Additionally, it gives the probability of observing a star for the conditions given in Table A.1 depending on the surface rotational velocity distribution.

Isochrones for slow rotating stellar evolution models that undergo non-chemically homogeneous evolution do not differ significantly in the HR diagram. The probability of observing a star along these isochrones is highest. The change from non-chemically homogeneous to quasi-chemically homogeneous evolution occurs at a small velocity range related to a significant change of the position in the HR diagram.Therefore, the probability of observing a star between isochrones of slow-rotating and rapidly rotating models is small. Again, isochrones of rapidly rotating stellar evolution models are close in the HR diagram, leading to a higher probability of observation.

Appendix B: Summary of data

All Tables

All Figures