© ESO 2021

1. Introduction

The recent identification of massive black-hole binaries with masses up to 80 M_⊙ as powerful sources of gravitational waves (Abbott et al. 2016a,b) opened the first direct view onto the end products of the evolution of the most massive stars. To be able to exploit this new information, it is necessary to model stellar populations that are young enough to contain very massive stars (VMS) with masses ≳100 M_⊙. Because of the steep increase in stellar luminosities with increasing mass, VMS in this regime have luminosities well in excess of 10⁶ L_⊙. The resulting high L/M ratios ≳10⁴ L_⊙/M_⊙ bring such VMS close to the Eddington limit, where the radiative acceleration a_rad equals the gravitational acceleration g.

Due to the resulting low effective surface gravities, VMS are subject to physical effects that impact their outer envelope structure and stellar winds. The envelopes of VMS are predicted to be subject to the envelope inflation effect, which leads to the formation of radially extended low-density envelopes (Ishii et al. 1999; Petrovic et al. 2006). In the limit of low convective efficiency, inflated envelopes become unstable above a specific L/M ratio because of a lack of hydrostatic envelope solutions (Gräfener et al. 2012). This leads to numerical problems in stellar evolution models, which are often tackled by an artificial increase in the convective efficiency (cf. Paxton et al. 2013). As a result, the inflation effect is suppressed in most evolutionary models. Relaxing this assumption leads to the formation of relatively cool stars with inflated envelopes at the top of the main sequence (Gräfener et al. 2012; Köhler et al. 2015; Sanyal et al. 2015).

One key question addressed in this work is whether the existence of the inflation effect can be corroborated by observations. This is particularly interesting because the inflation effect may be related to the still enigmatic group of luminous blue variables (LBVs) and their S-Dor type (radius) variability on timescales of months to years (cf. Gräfener et al. 2012; Grassitelli et al. 2020).

The enhanced mass loss of VMS and its transition from optically thin O-star winds to optically thick Wolf-Rayet (WR) winds have been investigated by numerical wind models (Gräfener & Hamann 2008; Vink et al. 2011) that are supported by systematic empirical studies of samples of VMS (Gräfener et al. 2011; Bestenlehner et al. 2014). Gräfener et al. (2011) pointed out that the underlying reason for the enhanced mass loss of WR stars is their proximity to the Eddington limit, and not their surface chemistry as is usually assumed in stellar evolution models. The enhanced mass loss is the reason why VMS often appear as emission-line stars, namely as hydrogen-rich Wolf-Rayet (WNh) stars, even during their main-sequence evolution (Hamann et al. 2006; Schnurr et al. 2008, 2009; Crowther et al. 2010). Due to the extreme luminosities of such very massive WNh stars, their spectral signatures are detectable in young starbursts even at cosmological distances (Wofford et al. 2014; Gräfener & Vink 2015).

The winds of WNh stars are often confused with those of classical WR stars. Classical WR stars are massive stars near the end of their evolution, predominantly in the phase of core He-burning. These stars lost their hydrogen-rich envelopes during their previous evolution and have much lower masses. They reach high Eddington factors due to the high luminosity of their He-burning cores. Classical WR stars have different physical and wind properties than the more massive core hydrogen-burning WNh stars, but they are usually treated with the same mass-loss prescriptions in evolutionary models.

In this work, we construct evolutionary models that take the enhanced mass-loss rates and the envelope inflation of VMS into account. We compute model grids for massive main-sequence stars up to 500 M_⊙ and compare the results with the most massive stars in the VLT-Flames Tarantula Survey (VFTS), the largest and best studied sample of VMS obtained to date (Evans et al. 2011). The VFTS is a near-complete spectroscopic survey of ≈1000 hot luminous stars in 30 Dor in the Large Magellanic Cloud (LMC), the closest analogue of a young starburst in our neighbourhood, with a well-controlled brightness cut-off. All single stars in the VFTS sample have been analysed individually, and their star-formation histories (SFHs) and rotational-velocity distributions have been investigated in detail (Schneider et al. 2018a,b; Ramírez-Agudelo et al. 2013; Dufton et al. 2013). This makes the sample an ideal target for detailed population synthesis modelling, enabling quantitative predictions of the expected numbers of stars in specific evolutionary phases.

We describe our evolutionary models in Sect. 2 and the population synthesis in Sect. 3. In Sect. 4 we compare our models with the observed population of VMS in 30 Dor, and in Sect. 5 we discuss the results. The most important conclusions are summarised in Sect. 6.

2. Stellar evolution models

We use the Modules for Experiments in Stellar Astrophysics (MESA) version 11701 (Paxton et al. 2011, 2013, 2015, 2018, 2019) to compute single-star evolutionary models covering the main-sequence phase from zero to 99% hydrogen exhaustion in the stellar core. The same input physics as in Schootemeijer & Langer (2018), Schootemeijer et al. (2019) are used, except for the treatment of convection in the outer stellar envelope and stellar-wind mass loss (see below). In particular, we adopt a wind-scaling as a function of X_Fe/X_Fe, ⊙, as described by Schootemeijer et al. (2019) We use the MESA default chemical composition from Grevesse & Sauval (1998) scaled to an initial metallicity of Z = 0.008 (compared to the corresponding solar metallicity of Z_⊙ = 0.0169) to match the metallicity, and in particular X_Fe, in the LMC.

Because of our interest in the envelope inflation effect we employ standard mixing-length theory (MLT) instead of the modified treatment of superadiabatic convection in radiation-dominated regions (MLT++, Paxton et al. 2013) that is used by default in MESA models. Furthermore, we introduce a new treatment of stellar-wind mass loss to meet the requirements for VMS near the Eddington limit. In the following sections these new mass-loss relations are described.

2.1. The winds of OB stars

The most widely used mass-loss prescription for OB stars is the relation by Vink et al. (2000, 2001). These theoretical mass-loss estimates are based on a Monte-Carlo approach that makes use of the first-order Sobolev approximation. In the Sobolev approximation spectral lines are assumed to be infinitely narrow, but their flux-mean opacity is high due to the effective line broadening through Doppler shifts in the wind. Consequently, the radiative acceleration depends predominantly on the velocity gradient ∂v/∂r, leading to the same wind physics as in Castor et al. (1975, CAK) who described the first closed theory for the formation of optically thin OB-star winds.

Notably, the Sobolev approximation does not apply in the limit of large wind optical depths where the photon mean-free path becomes small, and Doppler shifts become inefficient (cf. the discussion in Gräfener et al. 2017). In this work we therefore only apply the Vink et al. (2001) relation in the optically thin regime. Furthermore, we change the absolute value and the slope of Ṁ as a function of luminosity L, to match the average observed mass-loss rates of O stars in 30 Dor as in Bestenlehner et al. (2014).

For the absolute calibration of our mass-loss relations with empirical results in Sect. 4.2 we adopt a wind clumping factor of D = 10. This means that we assume that the wind consists of dense clumps with a density that is a factor D higher than the mean density, while the inter-clump medium is void. Clumping constitutes the main uncertainty in the empirical mass-loss rates of OB and WR stars. It increases the line strengths of recombination lines, and therefore reduces the empirical mass-loss rates derived from this type of lines by a factor $\sqrt{D}$ (Hillier 1991; Hamann & Koesterke 1998).

For the most massive O and WNh stars in 30 Dor Bestenlehner et al. (2014) found clumping factors of the order of D = 10, based on the relative strengths of the electron-scattering wings of strong recombination lines (cf. also Vink & Gräfener 2012). For O stars with weaker winds even higher clumping factors up to D ≈ 100 have been discussed, based on the strength of non-saturated UV resonance lines (Fullerton et al. 2006). However, considering the influence of optically thick clumps on the UV line profiles Oskinova et al. (2007), Sundqvist et al. (2014) derived more moderate values of the order of 10. Here we adopt D = 10 for all stars in our sample, resulting in a reduction of the empirical mass-loss rates by a factor of 3.3 compared to un-clumped winds.

Our comparison with observations in Sect. 4.2 leads to a considerable reduction of O-star mass-loss rates, and a shallower slope as a function of L. For the slope we adopt the value of 1.45 derived by Bestenlehner et al. (2014) and replace it in Eq. (24) of Vink et al. (2001) instead of the original value of 2.194. This substitution already lowers the mass-loss rates in the VFTS O-star range considerably. To match the observations, we further apply a constant factor of 0.8 to the altered mass-loss relation. For the dependence of Ṁ on the rotational velocity, we use the standard relation implemented in MESA (Eq. (26) in Paxton et al. 2013).

2.2. The winds of WNh stars

The most massive and luminous main-sequence stars are thought to appear spectroscopically as WNh stars (Gräfener & Hamann 2008; Crowther et al. 2010). WNh stars are therefore the extension of the O-star regime towards higher masses and luminosities. Previous implementations of wind mass loss in stellar evolution models do not take the transition from the optically thin winds of OB stars to the optically thick winds of WNh stars into account. Instead, they extrapolate the Vink et al. (2000) relations to higher luminosities. Only when the surface composition becomes substantially deficient in hydrogen, the mass-loss relations switch to WR type mass loss, where no distinction between the mass loss of very massive WNh stars, and the mass loss of classical core He-burning WR stars is made. For instance, in the default ‘DUTCH’ MESA wind-scheme the empirical mass-loss relation from Nugis & Lamers (2000) is used for hot stars (T_eff > 10 kK) with hydrogen surface mass fractions X_s < 0.4, otherwise the Vink et al. (2000) relations are used.

In this work we follow previous theoretical (Gräfener & Hamann 2008; Vink et al. 2011) and empirical (Gräfener et al. 2011; Bestenlehner et al. 2014) studies of the mass-loss properties of VMS that suggest that the proximity to the Eddington limit determines the formation of optically thick WR-type winds. According to these studies WR-type winds show a distinct, strong dependence on the proximity to the Eddington limit, which is most easily expressed via a dependence on the classical Eddington factor

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{e}}=\frac{{\chi }_{\mathrm{e}}L}{4\pi cGM}={10}^{-4.813}\times (1+X_{\mathrm{s}})\frac{L/L_{\odot }}{M/M_{\odot }}·\end{array}$$(1)

Here we define Γ_e as the ratio between the radiative acceleration on free electrons (with the electron-scattering opacity χ_e for a fully ionised plasma) and the gravitational acceleration of the star. Because χ_e depends only on the number of free electrons per gram, Γ_e depends only on fundamental stellar parameters, namely the hydrogen surface abundance X_s and the L/M ratio. In this picture the mass loss of VMS, which are still in the core hydrogen-burning phase, differs from the mass loss of classical core He-burning WR stars because of their different L/M ratios, effective temperatures, and chemical compositions.

In the following we concentrate on VMS at LMC metallicity. We employ the empirical mass-loss relation for VMS in 30 Dor from Bestenlehner et al. (2014) who derived a relation of the form

$$\begin{array}{c}\hfill log({\dot{M}}_{\mathrm{thick}})=5.22\times {\mathrm{\Gamma }}_{\mathrm{eff}}-0.5\times log(D)-2.6\end{array}$$(2)

for the optically thick winds of WNh stars. Here D denotes the wind clumping factor (cf. Sect. 4.2), and Γ_eff the effective Eddington factor including centrifugal acceleration (cf. the discussions in Gräfener & Hamann 2008; Gräfener & Vink 2015). As Γ_eff varies over the stellar surface, we use an intermediate value that lies between the value at the poles, where the centrifugal acceleration is zero, and the equator, where g_rot = Ω²R:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{eff}}={\mathrm{\Gamma }}_{\mathrm{e}}+0.5\times {\mathrm{\Gamma }}_{\mathrm{rot}}\equiv {\mathrm{\Gamma }}_{\mathrm{e}}+0.5\times \frac{{\mathrm{\Omega }}^2{R}^3}{\mathit{GM}}·\end{array}$$(3)

As a result, our mass-loss prescription is sensitive to a combination of the proximity to the Eddington limit and the critical rotation rate similar to the ΓΩ-limit discussed by Maeder & Meynet (2000), Langer (1997).

As a criterion to switch between the regimes of WR-type winds and optically thin OB-star winds, we evaluate their sonic-point optical depth τ_s. This approach is motivated by the theory of optically thick winds, where the properties of winds with large τ_s are largely constrained by the conditions at the sonic point and τ_s itself (Nugis & Lamers 2002; Gräfener et al. 2017). As described in Sect. 2.1, the different wind physics for the optically thin winds of OB stars (with small values of τ_s) leads to lower mass-loss rates.

To decide whether a stellar model with given surface values of M, L, R, Ω, and X_s can develop a WR-type wind, we estimate the value of τ_s that its wind would have if it had an optically thick wind. To this purpose we compute Ṁ_thick from Eq. (2) and estimate the terminal wind velocity v_∞ using the empirical relation from Lamers et al. (1995) with $v_{\infty }/v_{\mathrm{esc}}^{\mathrm{eff}}=2.51$ (cf. also Bestenlehner et al. 2014). For given Ṁ and v_∞ we estimate τ_s using a relation between τ_s and the wind efficiency factor η ≡ Ṁv_∞/(L/c) from Gräfener et al. (2017). η denotes the ratio between the mechanical wind momentum and the momentum of the radiation field driving the wind. It is related to τ_s via the number of absorption and re-emission processes needed to drive the wind. Gräfener et al. (2017) derived an approximate relation of the form

$$\begin{array}{c}\hfill {\tau }_{\mathrm{s}}\approx \frac{\dot{M}{v}_{\infty }}{L/c}(1+\frac{v_{\mathrm{esc}}^2}{v_{\infty }^2})·\end{array}$$(4)

In our models we use τ_s from Eq. (4) as a criterion for the onset of WR-type mass loss. For cases where τ_s < τ₁ we use the mass-loss prescriptions for optically thin winds from Sect. 2.1, and for τ_s > τ₂ we use the Γ-dependent mass-loss relation Eq. (2). For τ₁ < τ_s < τ₂ we interpolate between both cases. As demonstrated in Sect. 3 we can reproduce the observed mass-loss rates for VMS in 30 Dor from Bestenlehner et al. (2014) very well with values of τ₁ = 2/3 and τ₂ = 1.

3. Population synthesis models

In this section we describe the population synthesis models that we use to compare the models from Sect. 2 with observations.

3.1. Stellar evolution models

Each of our population synthesis models is based on a set of stellar evolution models for massive stars with initial masses in the range of 9−500 M_⊙ and initial rotational velocities between 0 and 90% of the critical rotation rate. Typically, each model grid consists of models with 18 different masses and 11 different rotational velocities, that is, a total of 198 models. As start models we use rotating, chemically homogeneous models in thermal equilibrium, with a default initial composition with Z = 0.008 (cf. Sect. 2). The models cover the main-sequence phase up to 99% hydrogen exhaustion.

3.2. Population synthesis

The goal of our population synthesis models is to compute the fraction of stars in a specific stellar population that show one or more observables o in a given parameter range o_i. The population is defined through a given initial mass function $\xi (M)=\frac{\mathrm{d}N}{\mathrm{d}M}$, initial rotational velocity distribution $P(v_{\mathrm{ini}})=\frac{\mathrm{d}N}{\mathrm{d}{v}_{\mathrm{ini}}}$, the SFH $S(t)=\frac{\mathrm{d}N}{\mathrm{d}t}$, and the total number of stars N.

Our models provide stellar parameters p as a function of time t for each evolutionary track with a given initial mass M and initial fraction of the critical rotation rate ω ≡ Ω/Ω_crit. We further determine the main-sequence lifetime T_MS for each evolutionary track and define the relative age t_r ≡ t/T_MS, so that the p are given as a function of the form p(t_r, M, ω) on the model grid. To compute the o_i we transform the p(t_r, M, ω) onto a homogeneous grid that is logarithmic in M and linear in ω and t_r, through linear interpolation between evolutionary tracks. In this step we typically use 100 grid points in M and t_r, and 64 grid points in ω, where t_r lies between 0 and 1, M between 9 and 500 M_⊙, and ω between 0 and 0.9. Based on the p(t_r, M, ω) we then compute any observable o(t_r, M, ω) of interest on the same grid.

To compute the probability distribution of the different o_i for the adopted stellar population we need to compute a weight function W(t_r, M, ω) that takes S(t) and the initial distributions of M and ω into account. The probability P_i to find an observable o in a given parameter range o_i is then

$$\begin{array}{c}\hfill P_i=\frac{{\sum }_{t_{\mathrm{r}},M,\omega }W(t_{\mathrm{r}},M,\omega )|{}_{o(t_{\mathrm{r}},M,\omega )\in o_i}}{{\sum }_{t_{\mathrm{r}},M,\omega }W(t_{\mathrm{r}},M,\omega )},\end{array}$$

where the summation is a sum over all grid points, and the weight W(t_r, M, ω) takes the initial distribution of M and ω and the sampling of t_r, M and ω into account. We write W(t_r, M, ω) in the form

$$\begin{array}{c}\hfill W(t_{\mathrm{r}},M,\omega )=W_t(t_{\mathrm{r}})\times W_M(M)\times W_{\omega }(\omega )\end{array}$$

and compute W_t, W_M and W_ω in the following way.

As we use linear time steps Δt_r, W_t is given by the product of S(t) and the derivative dt/dt_r = T_MS

$$\begin{array}{c}\hfill W_t=S(t)\times T_{\mathrm{MS}}(M,\omega ),\end{array}$$

with t = t_r × T_MS(M, ω).

For M we use a logarithmic grid with dM/dlog(M) = M, and adopt an initial mass function ξ(M) ∝ M^−γ, so that

$$\begin{array}{c}\hfill W_M=\xi (M)\times \frac{\mathrm{d}M}{\mathrm{d}log(M)}\propto M^{1-\gamma }.\end{array}$$

For ω we use a linear grid, but the initial rotational distribution P(v_ini) is given as a function of the rotational velocity, so that

$$\begin{array}{c}\hfill W_{\omega }=P(v_{\mathrm{ini}})\times \frac{\mathrm{d}{v}_{\mathrm{ini}}}{\mathrm{d}\omega }·\end{array}$$

Here we compute dv_ini/dω numerically, where we define v_ini(M, ω) ≡ v_rot(t_r = 0.01, M, ω) to avoid problems due to the rapid adjustment phase of the internal rotation profile at the start of each evolutionary sequence.

4. Very massive stars in 30 Dor

In this section we compare the results of our population synthesis with the empirically derived properties of the stellar population in 30 Dor. The empirical data consist of the sub-sample of single main-sequence stars from the VFTS as compiled by Schneider et al. (2018a). In particular, the empirical data are based on detailed spectroscopic analyses of the most massive stars in 30 Dor from Bestenlehner et al. (2014), O-type giants and supergiants from Ramírez-Agudelo et al. (2017), O-type dwarfs from Sabín-Sanjulián et al. (2017), B-type supergiants from McEvoy et al. (2015) and B dwarfs from Schneider et al. (2018a).

For the population synthesis we used the SFH S(t) and initial mass function (IMF) ξ(M) for 30 Dor from Schneider et al. (2018a), together with the rotational velocity distribution P(v) for single O stars from Ramírez-Agudelo et al. (2013). Analogous to the work of Schneider et al. (2018a) the present work thus focuses on single stars with a luminosity ≳10^4.5 L_⊙, well above the VFTS brightness cut-off. This means that only stars with masses ≳15 M_⊙ and ages ≲12 Myr are considered in the present analysis.

We further employed evolutionary models for single stars as described in Sect. 2, tailored to match the empirical properties of the massive single-star population in 30 Dor. Because of our focus on the physics of VMS these models are designated as VMS models. We compare these models with empirical data from VFTS and models based on alternative input physics and/or common standard assumptions.

4.1. Mass distribution

An important input parameter for our population synthesis models is the exponent γ of the IMF with ξ(M) ∝ M^−γ. In Fig. 1 we compare the present-day luminosity distribution resulting from our population synthesis with the observed distribution in the VFTS sample. In this plot we normalised the synthetic distribution to match the total number of 288 observed, putatively single stars above log(L/L_⊙) = 4.5. We further adjusted γ to match the number of 20 observed VMS in the range 6.0 < log(L/L_⊙) < 6.6 (roughly corresponding to a mass range of 70−200 M_⊙). For our VMS models this was achieved using a value of γ = 1.977, slightly higher than the value of $\gamma =1.{90}_{-0.26}^{+0.37}$ derived by Schneider et al. (2018a). Using alternative models employing the DUTCH MESA wind scheme (cf. Sect. 2.2) we obtain γ = 1.935, in agreement with Schneider et al. (2018a). Altogether our models confirm a top-heavy IMF in comparison to the standard Salpeter (1955) IMF with γ = 2.35, but with a moderate impact of the adopted mass-loss prescription on the present stellar masses and therefore on γ. In the remainder of this work, we use γ = 1.977 to achieve consistency between our models and observations.

Fig. 1. Observed luminosity distribution of single stars in 30 Dor compared to theoretical predictions. The histogram shows the observed numbers of single stars in the VFTS sample as function of luminosity (blue) compared to the theoretical distribution resulting from our VMS population synthesis models described in Sect. 4 (green). Both distributions are plotted transparent so that overlaps appear in dark green.

Notably, there are no stars with luminosities higher than log(L/L_⊙) = 6.6 in the VFTS sample while our VMS models predict 4.7 stars above this luminosity. This upper luminosity limit corresponds to an upper initial-mass limit of M_ini ≈ 200 M_⊙. We will discuss the relevance of this upper mass limit in Sect. 5.2. For the investigation of the physical properties of VMS in the remainder of this work, we consider the consistency between observed and predicted stellar numbers below 200 M_⊙ as satisfactory.

4.2. Mass-loss rates

To test the validity of the mass-loss prescriptions employed in our models we compare the mass-loss relations resulting from our population synthesis models directly with empirical data. In Figs. 2 and 3 the resulting modified wind momenta $D_{\mathrm{mom}}=\dot{M}{v}_{\infty }\sqrt{R/R_{\odot }}$ are shown as a function of luminosity L. According to the CAK theory for optically thin winds, OB stars are expected to follow a tight wind-momentum luminosity relation (WLR; Kudritzki et al. 1999), largely independent of secondary parameters, such as mass or radius.

Fig. 2. Wind-momentum – luminosity relation for O and WNh stars in 30 Dor. Different symbols indicate empirical results for the most massive stars in 30 Dor from Bestenlehner et al. (2014) in blue, O-type giants and supergiants from Ramírez-Agudelo et al. (2017) Tables C.4 (dark green) and C.5 (light green), and O-type dwarfs from Sabín-Sanjulián et al. (2017) Tables A.1 (red) and A.2 (pink). Upper limits are indicated in grey. Evolutionary tracks are indicated by coloured lines analogous to Fig. 5 but restricted to the O star range. Predicted absolute numbers of stars according to our population synthesis models are mapped and colour coded as indicated by the colour bar on the right.

Fig. 3. Wind-momentum – luminosity relation as in Fig. 2 compared with population synthesis models based on standard MESA models using the DUTCH wind scheme.

In Figs. 2 and 3, we indicate the expected numbers of stars in this diagram according to our population synthesis models by colour maps. The synthetic D_mom were computed based on mass-loss rates Ṁ, radii R and stellar masses M from the underlying stellar evolution models. The v_∞ were computed from the resulting effective escape velocities using a ratio of $v_{\infty }/v_{\mathrm{esc}}^{\mathrm{eff}}=2.51$ as described in Sect. 2.2. The empirical wind momenta were compiled from Ramírez-Agudelo et al. (2017), Sabín-Sanjulián et al. (2017), Bestenlehner et al. (2014) and cover the complete sample of single O and WNh stars in the VFTS, meaning that they include all apparently single main-sequence stars with effective temperatures > 29 kK. The same temperature cut-off has been applied to the theoretical results.

Our VMS models in Fig. 2 show a good agreement with observations. In particular, the location of the transition from O-star winds to the enhanced mass loss in the WNh phase is reproduced well. On the other hand, the WLR resulting from the oft-used DUTCH wind scheme in Fig. 3 matches the WNh stars well, but over-estimates the mass-loss rates of the O stars in the VFTS sample. Also, the slope of the WLR appears to be too steep in the O-star range. For both wind schemes it is noticeable that the synthetic WLRs are much narrower than the observed one. In particular for O stars the observed WLR shows a large spread in D_mom, where dwarfs tend to have lower wind momenta than giants and supergiants by up to a factor of ten. As we will show later, this may hint on substantially reduced mass-loss rates for younger O stars.

A novelty in our VMS models is the treatment of the transition between the regimes of optically thin O-star winds and optically thick WNh winds. This transition occurs qualitatively differently in the VMS and DUTCH approaches. In the VMS wind scheme in Fig. 2 stars with luminosities above log(L/L_⊙) ≈ 5.5 reach high enough Eddington factors at some point during their main-sequence evolution to perform a quick transition to enhanced WR-type mass-loss rates. In the WNh regime the evolution then proceeds at nearly constant luminosity (cf. Fig. 5). On the other hand, the DUTCH models in Fig. 3 show no marked transition between the two regimes but change their behaviour at some point during the WNh stage, when they reach a helium surface mass fraction of Y_s = 0.4. Before this point they evolve continuously with increasing luminosity from the O-star regime into the WNh stage. After they reached Y_s = 0.4 their mass-loss rates increase, and their evolution continues with decreasing luminosity.

As noted earlier, the absolute mass-loss rates in the WNh range are reproduced well by both mass-loss prescriptions. One notable outlier in Figs. 2 and 3 is VFTS 108, which has a WR-type mass-loss rate despite of its low luminosity of log(L/L_⊙) = 5.7. This star is highly enriched in helium (Y_s ≈ 0.8) and most likely already in the core helium-burning phase. Bestenlehner et al. (2014) demonstrated that this star matches the adopted Ṁ(Γ_eff) relation (Eq. (2)) precisely, if a lower mass, corresponding to the core helium-burning case, is adopted.

In Fig. 4 we put the Γ_e-dependent mass-loss prescription (Eq. (2)) in our VMS models to a direct test. As Γ_e is not directly observable, we estimate the Eddington factor ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}(L,X_{\mathrm{s}})$ for the chemically homogeneous case instead, where we compute ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}(L,X_{\mathrm{s}})$ from Eq. (1) using the mass estimate M_hom(L, X_s) for a chemically homogeneous star from Gräfener et al. (2011). Notably, ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}(L,X_{\mathrm{s}})$ depends only on observable quantities. In reality, most stars are of course not chemically homogeneous, but have a chemical profile with an increasing mean molecular weight towards their centre. For fixed values of M and X_s this results in a higher luminosity than for a chemically homogeneous star. The other way round, for given L and X_s, M_hom(L, X_s) constitutes a strict upper limit to the real mass of a star, and ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}$ a lower limit to Γ_e.

Fig. 4. Mass loss as a function of Eddington factor. Mass-loss rates Ṁ as a function of the estimated classical Eddington factor for chemically homogeneous stars ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}$ (see text) are indicated analogous to Fig. 2. The Ṁ(Γe) relation for R 136 from Bestenlehner (2020) is indicated by the dashed curve.

In Fig. 4 we plot the synthetic and observed mass-loss rates Ṁ for our O/WNh sample as a function of ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}$. To allow for a direct comparison between theory and observations, the ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}(L,X_{\mathrm{s}})$ are computed in an identical manner for models and observations. Analogous to the WLR in Fig. 2 the theoretical mass-loss rates increase rapidly as soon as the formation of optically thick winds becomes possible. The transition between O and WR-type winds happens at $log({\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}})\approx -0.5$ and continues with a steeper slope for high $log({\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}})$, in good agreement with the observations. As our Γ_e-dependent mass-loss prescription (Eq. (2)) was initially derived from the same data under the assumption that ${\mathrm{\Gamma }}_{\mathrm{e}}={\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}$ (cf. Bestenlehner et al. 2014) this can be seen as a confirmation that most WNh stars in our sample are in fact evolving chemically homogeneously.

Again, the overall agreement between models and observations in the WNh stage (with log(Ṁ/M_⊙ yr⁻¹) ≳ − 5.2) is very good, except for a short phase directly after the transition. In this short phase the surface is still more hydrogen-enriched than the stellar core, and thus ${\mathrm{\Gamma }}_{\mathrm{e}}>{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}$. However, as soon as the enhanced mass loss in the WNh phase sets in, the H-rich surface layers are quickly removed, and the stars become nearly chemically homogeneous.

While the O-star mass-loss rates in Fig. 4 (with log(Ṁ/M_⊙ yr⁻¹) ≲ − 5.2) show a good agreement with observations for ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}\gtrsim -1$, the agreement gets worse for lower values of ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}$. For the O dwarfs in this regime Sabín-Sanjulián et al. (2017) derived upper limits for Ṁ which are indicated in grey. For O (super)giants Ramírez-Agudelo et al. (2017) provided explicit values for Ṁ that tend to be higher than those in our models. Notably, Ramírez-Agudelo et al. (2017) did not provide absolute values for D_mom for the same objects (cf. Figs. 2 and 3) because they considered the results as uncertain (H. Sana and A. de Koter, priv. comm.). Depending on whether the mass-loss rates in this parameter range can be believed or not the differences between O dwarfs and (super)giants, which were already apparent in Fig. 2, may be even more pronounced for these objects. This is further supported by recent results of Bestenlehner (2020, 2020) who investigated the mass-loss properties of the most massive stars in the R 136 cluster in the centre of 30 Dor. These authors determined a very young age of 1−2 Myr for R 136 and derived an Ṁ(Γ_e) relation that lies up to two orders of magnitude below the values derived for O stars the VFTS sample as indicated by the dashed curve in Fig. 4.

4.3. Location in the Hertzsprung-Russell diagram

In Fig. 5 we show the theoretical HRD for 30 Dor including a colour-coded map of the expected numbers of stars according to our population synthesis models. For comparison, we show the observed HRD in Fig. 6 with the observed number of stars colour-coded in the same way. Generally, the observed HRD positions of stars above log(L/L_⊙) ≈ 5 are well reproduced by our VMS models. The lack of evolved stars below this luminosity in Fig. 5 is caused by the maximum age of 12 Myr adopted in our population synthesis models. Beyond this, our VMS models reproduce the observed bending of the zero-age main sequence (ZAMS) for the highest masses as well as the lack of stars with cooler temperatures than log(T_eff/K) ≈ 4.5 above log(L/L_⊙) ≈ 6 very well.

Fig. 5. Theoretical HRD for VMS models with standard MLT. Coloured areas indicate the expected number of stars according to our population synthesis models for 30 Dor (cf. Sect. 4). Evolutionary tracks for Ω/Ωcrit = 0.3 are indicated, labelled with ZAMS masses in M⊙.

Fig. 6. Observed HRD for 30 Dor. Coloured areas indicate the observed number of stars, evolutionary tracks are the same as in Fig. 5.

The bending of the ZAMS towards cooler temperatures at the top of the main sequence is a consequence of the envelope inflation effect. This is demonstrated in Fig. 7, where the inflation effect has been suppressed by artificially increasing the convective efficiency in radiation-dominated stellar envelopes (MESA option okay_to_reduce_gradT_excess = .true.). As a consequence, the most massive stars above ≈80 M_⊙ become much hotter and evolve towards even hotter temperatures during their main-sequence evolution. The reason is that stars in this regime become chemically homogeneous and therefore more and more compact throughout their evolution. Only when standard MLT is used (Fig. 5) does the inflation effect compensates this effect and brings the effective temperatures of stars in this mass range into agreement with observations.

Fig. 7. Influence of envelope inflation. The theoretical HRD for models with reduced envelope inflation effect through an artificially enhanced convective efficiency (MLT++; MESA option okay_to_reduce_gradT_excess = .true.) is indicated analogous to Fig. 5.

In Figs. 8 and 9 we focus on the distribution of stars in the upper HRD, relative to the ZAMS. To demonstrate how the predicted HRD positions depend on mass loss, we show models with alternative mass-loss prescriptions in Figs. 10 and 11. In Fig. 10 the onset of the enhanced mass loss of WNh stars has been completely omitted while otherwise the same wind scheme as in our VMS models was used. Consequently, the mass loss of the most massive stars is reduced, and the outer H-rich layers are removed much more slowly, leading to higher hydrogen surface mass fractions X_s and therefore higher values of Γ_e (cf. Eq. (1)). Due to the strong dependence of the inflation effect on Γ_e (cf. Gräfener et al. 2012, Eqs. (20), (28), and (29)), this leads to an enhanced envelope inflation and correspondingly cooler stellar temperatures.

Fig. 8. Effective temperature distribution. The theoretical HRD for the upper main sequence is indicated analogous to Fig. 5 but as a function of the effective temperature relative to the ZAMS temperature.

Fig. 9. Observed HRD analogous to Fig. 8.

Fig. 10. Influence of mass loss. The theoretical HRD is indicated analogous to Fig. 8 but adopting a pure OB-type mass-loss prescription without enhanced Γ-dependent mass loss.

The models in Fig. 11 are computed using the standard DUTCH wind scheme with much higher mass-loss rates than our VMS models. Because of the higher mass loss the most massive stars become more compact and less massive, and therefore hotter and less luminous. The qualitatively correct reproduction of the observed HRD positions of the most massive stars with log(L/L_⊙) ≳ 6 by our VMS models can therefore be seen as an indirect confirmation of our new mass-loss prescription, and in particular of the adopted clumping factor D in Eq. (2).

Fig. 11. Influence of mass loss. The theoretical HRD is indicated analogous to Fig. 8 but adopting the DUTCH MESA wind scheme.

For stars with log(L/L_⊙) ≲ 6, our models predict very extended envelopes and correspondingly cool stellar temperatures towards the end of the main-sequence evolution, irrespective of the mixing approach (cf. Figs. 5 and 7). In the observed HRD (Fig. 9) this region is populated by B supergiants whose evolutionary origin is still unclear (e.g., Vink et al. 2010; McEvoy et al. 2015). In Fig. 9 there are only a few objects of this type that are clearly located in the Hertzsprung-gap, that is, outside the expected parameter range for main-sequence stars.

To quantify the contribution of the envelope inflation effect, we compare in Fig. 12 the distribution of stars between 5 < log(L/L_⊙) < 6 as a function of their distance from the ZAMS (log(T_eff/T_ZAMS)) with observations. Our models predict 22.4 stars with log(T_eff/T_ZAMS) <  − 0.2 compared to 30 stars that are observed. This means that we expect the majority of stars in this regime to be evolved main-sequence stars with inflated envelopes. Setting the maximum effective temperature for stars to be classified as B-type stars to 29 kK, we predict 22.4 stars in this phase to appear as B supergiants, while 34 stars are observed. Again, according to our models the majority of B supergiants between 5 < log(L/L_⊙) < 6 are evolved main-sequence stars.

Fig. 12. Effective temperature distribution. Histogram of the observed numbers of single stars in the VFTS sample with 5 < log(L/L⊙) < 6 as a function of their distance from the ZAMS (log(Teff/TZAMS)) (transparent blue) compared with the theoretical distribution from our VMS population synthesis models (transparent green).

While the overall agreement between theory and observation in Fig. 12 appears satisfactory, it is notable that there is a small systematic shift in T_eff/T_ZAMS between the observed and predicted median locations of the main sequence, by 0.025 dex. In principle, this could be an indication that the predicted core masses in our models are too small, similar to what has been claimed, for instance, by Tkachenko et al. (2020). For this reason, we varied the overshooting parameter analogous to Tkachenko et al. (2020) as well as the parameters for rotational mixing in our models but did not achieve a satisfactory agreement between theory and observations. Also, the relatively high oxygen abundance in our models (based on the solar composition from Grevesse & Sauval 1998) does already bring our models towards slightly lower effective temperatures than what would result from more recent estimates (such as Asplund et al. 2009). Therefore, it remains unclear whether the observed shift between theory and observations is due to deficiencies in our models or systematic errors in the spectroscopic analyses of the VFTS sample.

4.4. Helium surface enrichment

In Fig. 13, we show a comparison between the observed and predicted helium surface mass fractions Y_s for stars with T_eff > 29 kK, that is, for O stars and WNh stars in the VFTS, as a function of their estimated inverse mass-loss timescale Ṁ/M_hom. Again, for a meaningful comparison between models and observations, the mass estimate M_hom(L, X) from Gräfener et al. (2011) has been computed in an identical manner for models and observations.

Fig. 13. Helium surface enrichment through mass loss. The predicted helium surface mass fraction Ys, as a function of the estimated inverse mass-loss timescale Ṁ/Mhom for stars with Teff > 30 kK and log(L/L⊙) > 5, is compared with VFTS O and WNh stars labelled as in Fig. 2. The colour scale indicates the number of objects predicted by our population synthesis. Evolutionary tracks from Fig. 5 are indicated by coloured lines.

Our population synthesis models show two regions with substantially enhanced values of Y_s in Fig. 13. The WNh stars with log Ṁ/M_hom ≳ − 7.5 follow a clear trend of increasing Y_s with increasing Ṁ/M_hom (cf. also Bestenlehner et al. 2014, 2020). The reason is that the enhanced mass loss in the WNh regime efficiently removes the hydrogen-rich outer layers, leading to quasi chemically homogeneous evolution with a strong dependence of Ṁ on Γ_e, and therefore on Y_s.

In the O-star range, with lower values of Ṁ/M_hom, our models predict substantial He enrichment for initially fast-rotating stars. This is demonstrated in Fig. 14 where we show the He enrichment for O stars with 5 < log(L/L_⊙) < 5.9 (that is, excluding WNh stars) as a function of rotational velocity. According to our models, fast-rotating stars in this regime develop strong winds due to their proximity to the ΓΩ-limit. This leads to a phase of enhanced mass loss that removes angular momentum until the rotational velocity falls below a value of typically 300 km s⁻¹. As a consequence, there is a peak in the rotational velocity distribution for O stars at v_rot ≈ 300 km s⁻¹ (cf. the orange step function in Fig. 14). Notably, Ramírez-Agudelo et al. (2013) found a similar feature in the observed rotational velocity distribution of the VFTS O star sample, however, at a higher velocity of ≈400 km s⁻¹.

Fig. 14. Helium surface enrichment through rotation. The predicted helium surface mass fraction Ys, as a function of rotational velocity for O stars with Teff > 29 kK and 5.0 < log(L/L⊙) < 5.9, is compared with observed Ys and v sin i of He-enriched VFTS O stars (with Ys > 0.35) labelled as in Fig. 2; stars with log(L/L⊙) < 5.0 are plotted in grey. The colour scale indicates the number of objects predicted by our population synthesis models. Evolutionary tracks for an initial mass of 30 M⊙ are indicated by coloured lines. Step functions indicate the initial (blue) and present (orange) total number distributions in arbitrary units, as well as the present distribution of He-enriched stars with Ys > 0.35 (red).

According to our models, the rotational braking associated with this process induces substantial mixing of helium-rich material to the surface, with a maximum in the number distribution of He-enriched stars at ≈300 km s⁻¹ (cf. the red step function in Fig. 14). The same objects constitute the group of helium-enriched stars with log Ṁ/M_hom ≲ − 7.5 in Fig. 13. According to our models, we expect 14.6 O stars with Y_s > 0.35 in the luminosity range 5.0 < log(L/L_⊙) < 5.9 in the VFTS sample that have been enriched by rotational mixing. This compares to seven observed objects in the same regime, which are indicated by coloured symbols in Fig. 14. For lower luminosities, we expect 0.7 objects, compared to five observed objects, which are indicated by grey symbols in Fig. 14.

Notably, VFTS 285, the star with the highest value of v sin i = 600 km s⁻¹ in our sample and Y_s = 0.36 belongs to the group of low-luminosity stars for which no He-enhancement is predicted (VFTS 285 is not visible in Fig. 14). In both, the high- and low-luminosity regime, we expect the vast majority of He-enriched objects near a rotational velocity of ≈300 km s⁻¹. Due to the effect of random inclination angles the expected average value of v sin i is π/4 (cf. Brott et al. 2011), with a decreasing distribution function for small values of sin i. We therefore expect the distribution of He-enriched stars to peak near 235 km s⁻¹, while the observed v sin i distribution appears evenly distributed, with a potential peak near 100 km s⁻¹. This means that neither the predicted numbers of rotationally enriched stars, nor their distribution agrees with observations. The implementation of rotational mixing and/or rotationally induced mass loss in our models therefore seems to leave space for improvements (cf. the discussion in Sect. 5.4).

5. Discussion

In the following we discuss the implications of our results for the physics and evolution of the most massive stars (Sect. 5.1), their upper mass limit (Sect. 5.2), the evolution of OB stars (Sect. 5.3), and the mechanisms for surface He-enrichment in the upper HRD (Sect. 5.4).

5.1. WNh stars: The most massive stars

According to our analysis, at LMC metallicity, stars in excess of ≈80 M_⊙ enter the WNh stage at some point during their main-sequence evolution. This means that they reach high enough Eddington factors to initiate optically thick winds, so that they appear as WR-type emission-line stars. Because they still have hydrogen on their surface, they are traditionally classified as H-rich WN stars, or WNh stars. As shown in Sect. 4.2, the absolute mass-loss rates, as well as the onset of WR-type mass loss in this regime, are well described by our VMS models employing a Γ_e-dependent mass-loss relation.

Because of their enhanced mass loss, the outer radiative layers of WNh stars are quickly removed so that they develop a quasi chemically homogeneous structure, with hydrogen surface abundances starting at near-solar values and almost reaching hydrogen depletion at the end of their main-sequence evolution. This is in stark contrast to the traditional treatment of WR-type mass loss in stellar evolution models, where enhanced mass loss only sets in beyond a limit of Y_s > 0.4. In Sect. 4.4, we showed that our VMS models describe the resulting He-enrichment through the enhanced mass loss in the WNh stage well, in good agreement with observations.

Because of their proximity to the Eddington limit, WNh stars are subject to the envelope inflation effect. Therefore, their HRD positions are shifted towards cooler temperatures. For masses ≳80 M_⊙ the ZAMS is bent towards cooler temperatures due to this effect. As we showed in Sect. 4.3, our VMS models with standard MLT reproduce this effect very well, in agreement with the observed HRD positions. In particular, there appears to be no need to invoke a clumped envelope structure to explain their radii as it has been proposed for the inflated envelopes of classical WR stars (Gräfener et al. 2012; Gräfener & Vink 2013).

A further consequence of the inflation effect is that WNh stars evolve towards cooler temperatures, that is, towards the right side in the HRD, although they are close to chemical homogeneity. Again, the HRD positions following from our VMS models agree with observations. Models that employ the oft-used modified treatment of super-adiabatic convection in radiation-dominated regions (MLT++, Paxton et al. 2013) clearly fail to reproduce the observed HRD positions of WNh stars. They show much hotter temperatures than observed and evolve towards the left side in the HRD (cf., Figs. 6 and 7).

The general agreement of our models with the observed population of WNh stars in 30 Dor in terms of numbers, HRD positions, mass loss, and surface enrichment strongly support their evolutionary status as VMS in the phase of core H-burning.

5.2. A variable upper stellar mass limit

The first estimate of an upper limit for the masses of stars of ≈150 M_⊙ was determined by Figer (2005), based on photometric data of the Arches cluster in the Galactic Centre, one of the youngest and densest clusters in the Galaxy. Later Crowther et al. (2010) analysed the brightest WNh stars in R 136, the central cluster of 30 Dor, and found luminosities up to log(L/L_⊙) = 6.94 ± 0.1, corresponding to initial masses up to 300 M_⊙. Because of its high stellar density R136 has not been included in the VFTS sample but Bestenlehner et al. (2020) analysed the 55 brightest stars in R 136 based on new HST spectroscopy of Crowther et al. (2016) and revised the upper luminosity limit to log(L/L_⊙) = 6.8 ± 0.1, corresponding to an upper mass limit of ≈250 M_⊙.

According to our present models the brightest WNh stars evolve with near-constant luminosity, implying slightly higher initial masses than the ones derived by Bestenlehner et al. (2020). As a consequence, a value of log(L/L_⊙) = 6.8 corresponds almost precisely to M_ini = 300 M_⊙. As we showed in Sect. 4.1, the VFTS sample displays an upper luminosity cut-off at a slightly lower value of log(L/L_⊙) ≈ 6.6 (cf. Fig. 1), corresponding to an initial mass of ≈200 M_⊙. According to our population synthesis models we would expect 4.7 stars above log(L/L_⊙) = 6.6 and 1.8 stars above log(L/L_⊙) = 6.8 without assuming an upper mass limit. If the R 136 upper mass limit would be universal, we would therefore miss 3.1 stars with 6.6 < log(L/L_⊙) < 6.8 in the VFTS sample. Given that there are 20 stars with log(L/L_⊙) > 6.0, the probability to miss 3 stars with 6.6 < log(L/L_⊙) < 6.8 by chance would be (17/20)²⁰ = 0.039. If there were no uncertainties in the individual luminosity measurements this would indicate a significant difference of the upper mass limits in R 136 and the remainder of 30 Dor.

Now the relevant most massive stars in both samples are spectroscopically remarkably similar objects of subtype WN5h, which have been analysed by Bestenlehner et al. (2014) for the VFTS sample and Bestenlehner et al. (2020) for R 136. In both works individual uncertainties of ±0.1 dex in log(L) were estimated which are of course significant compared to the difference in the maximum luminosities of 0.2 dex. However, we think that the relative uncertainties in the derived log(L) may be smaller than the individual errors because of the similarity of the objects, and because they have been analysed by the same person using the same method of analysis. Moreover, the brightest star in the VFTS sample, VFTS 1025 alias R 136c is in fact a member of R 136 with a luminosity of log(L/L_⊙) = 6.58. The second brightest star, VFTS 682 has a much lower luminosity of log(L/L_⊙) = 6.51, and has been discussed as a slow runaway from R 136 (Renzo et al. 2019). In any case the VFTS population excluding R 136c clearly displays lower luminosities than log(L/L_⊙) = 6.6 compared to R 136 with its brightest member R 136a1 with log(L/L_⊙) = 6.78.

As Bestenlehner et al. (2020) report 3 stars with 6.6 < log(L/L_⊙) < 6.8 in R 136, in agreement with expectations for a top-heavy IMF, dynamical mass segregation cannot account for the lack of the same number of stars in the VFTS sample. Instead, the higher mass limit for R 136 could imply that R 136 initially formed stars with higher masses, depending on the local conditions in the star-forming region. This could mean that the stars in 30 Dor and R 136 already formed in a mass-segregated way (cf. Pavlík et al. 2019). Alternatively, the most massive stars in R 136 could be the result of successive stellar mergers due to the high stellar density in the cluster (e.g., Portegies Zwart et al. 2004).

5.3. OB stars

5.3.1. Mass-loss rates

For the mass-loss rates of OB stars we employed a modified form of the Vink et al. (2001) recipe adjusted to match the empirical mass-loss rates in the VFTS sample with an adopted clumping factor of D = 10. To this purpose we adopted a shallower dependence on L from Bestenlehner et al. (2014) and used overall reduced mass-loss rates. In Sect. 4.2, we showed that this mass-loss relation does indeed match the average O star mass-loss rates in the VFTS sample, however, it fails to reproduce their observed wide spread. In particular, the O dwarfs analysed by Sabín-Sanjulián et al. (2017) show systematically lower wind momenta and mass-loss rates than the O (super)giants analysed by Ramírez-Agudelo et al. (2017). This trend of reduced mass-loss rates for younger objects becomes even more prominent when the VFTS results are compared with the mass-loss relation from Bestenlehner (2020) for O and WNh stars in R 136. The stars in R 136 are only about 2 Myr old, and therefore even younger than most of the stars in the VFTS sample. Notably, the mass-loss rates for O stars in R 136 are up to two orders of magnitudes lower than those found in the VFTS sample (cf. Fig. 4), possibly indicating a very strong trend of reduced mass-loss rates for young O dwarfs.

A reduction in O-star mass-loss rates is also reported in recent theoretical works that employ more elaborate techniques to compute the radiative force than the oft-used first-order Sobolev approximation as in Vink et al. (2001) Employing higher-order Sobolev models, Owocki & Puls (1999) reported a drop in the line source function near the sonic point due to rapid changes in the velocity gradient. Gräfener (2003) reported a similar effect in non-LTE models employing a co-moving frame (CMF) radiative transfer causing a local drop in the radiative force near the sonic point. More recent mass-loss predictions for OB stars employing higher-order Sobolev or CMF techniques by Lucy (2007), Krtička & Kubát (2017), Sundqvist et al. (2019) report similar effects and predict much lower mass-loss rates than those of Vink et al. (2000). One could speculate that young stars with high log g are more susceptible to this effect because a high log g implies stronger changes in the velocity gradient. Physical effects occurring in more evolved stars that are not included in current wind models, such as turbulence caused by sub-surface convection (Cantiello et al. 2009), could help to overcome the associated local drop in the radiative force for these objects, and help to bring their mass-loss rates closer to the values of Vink et al. (2000) as it is observed for the O supergiants in the VFTS sample (cf. Ramírez-Agudelo et al. 2017).

5.3.2. HRD positions

In agreement with previous studies by Köhler et al. (2015), Sanyal et al. (2015) we find that luminous OB stars are subject to the inflation effect, albeit less dependent on the detailed treatment of convection than for the more luminous WNh stars. As a consequence, such objects appear as B supergiants during their main-sequence evolution, in regions of the HRD that previously may have been attributed to the Hertzsprung-gap. In Sect. 4.3, we showed that, according to our population synthesis models, 22 out of the 34 bright B supergiants in the VFTS sample are expected to be inflated main-sequence stars. However, there are some (at least four) objects in the VFTS sample whose HRD positions cannot be explained by our models, so that alternative explanations, such as blue loops or binary interaction, may still be needed to explain the complete B supergiant population.

Notably, our population synthesis models are also able to reproduce the lack of B supergiants with high luminosities (log(L/L_⊙) ≳ 6) in the VFTS sample. In Sect. 4.3 we identified the removal of H-rich layers due to the onset of enhanced WR-type mass loss as the reason for this effect. The successful reproduction of the VFTS population in the HRD is therefore also an indirect confirmation of the absolute mass-loss rates and clumping factors adopted in our models.

Previous works, such as Gräfener et al. (2012) and Grassitelli et al. (2018), associated stars with inflated envelopes with the S Doradus type variability of LBVs. The fact that we expect about 20 objects in the VFTS sample to be in this phase, but none of the observed stars shows LBV-type behaviour, implies that envelope inflation is not generally associated with the LBV phenomenon, in particular not in the H-burning phase. This could still mean that the inflation effect is associated with LBVs, but under rarer conditions, such as fast rotation, or in very late (pre-SN) evolutionary stages (cf. Groh et al. 2006, 2013; Gräfener et al. 2012).

Finally, a clumped density structure inside inflated stellar envelopes has been proposed by Gräfener et al. (2012), Gräfener & Vink (2013) to explain the observed radii of classical core He-burning WR stars. Our present analysis shows that there is no need to invoke similar clumped sub-surface regions to explain the radii of B supergiants through the inflation effect.

5.4. Surface helium enrichment through mass loss and rotation

As outlined in Sect. 5.1 the surface He-enrichment through the strong mass loss of WNh stars is well reproduced by our models. On the other hand, the He-enrichment through rotational mixing appears to be over-predicted by our models. According to our models we expect to find 15 O stars with Y_s > 0.35 in the VFTS sample that are He-enriched through rotational mixing. These objects are expected to have high rotational velocities (≈300 km s⁻¹) and high luminosities 5.0 < log(L/L_⊙) < 5.9. In reality, we find 12 He-enriched objects that are concentrated at low values of v sin i of ≈100 km s⁻¹, five of which are located at luminosities of log(L/L_⊙) < 5.0. Therefore, the observed sample appears to be largely incompatible with our model predictions, implying a much lower number of truly rotationally enriched O stars.

This result implies that either the internal rotational mixing is over-predicted in our models, or that the rotationally enhanced mass loss, which induces mixing, is too high. A large part of the 15 He-enriched stars in the VFTS sample could therefore be false detections (such as unidentified multiple systems) or products of alternative processes, such as binary interaction or mergers (cf. de Mink et al. 2014). Whether the observed discrepancy can be explained through changes in internal mixing parameters, or the previously discussed lower mass-loss rates for young O stars, remains to be shown in future works.

6. Conclusions

The implementation of new mass-loss relations in combination with standard-MLT envelope physics in our present evolutionary models leads to a successful reproduction of the properties of the most massive stars (≳80−100 M_⊙) in 30 Dor, including their HRD positions, mass-loss rates, and surface He-enrichment. In particular the transition between optically thin O-star winds and the optically thick winds of WNh stars is well reproduced.

For the less massive O stars (with ≈25−100 M_⊙) we find much more diverse mass-loss properties than previously thought, with indications for substantially reduced mass-loss rates for young O dwarfs. Also, the surface helium enrichment in the O-star range does not appear to be well understood as we find no clear indications for rotational mixing of helium in contrast to our model predictions. Both phenomena may be related as reduced mass-loss rates for young fast-rotating stars could lead to less rotational braking and mixing.

According to our models the distribution of stellar luminosities in 30 Dor is in good agreement with the top-heavy IMF from Schneider et al. (2018a), with a very modest dependence on the adopted mass-loss prescription. However, the upper luminosity cut-off in the greater 30 Dor region is lower than for its central cluster R 136, implying an upper initial mass limit of ≈200 M_⊙ compared to ≈300 M_⊙ in R 136. This could be an indication of a general dependence of the upper stellar mass limit on the conditions in the star-forming environment.

Finally, the observed number distribution of stars as a function of effective temperature supports the existence of inflated stellar envelopes that are theoretically predicted to form when main-sequence stars reach high L/M ratios. This puts the majority of bright B supergiants in the end phase of core hydrogen-burning, where they display cool temperatures due to their large, inflated radii.

Acknowledgments

GG thanks the Deutsche Forschunsgemeinschaft (DFG) for financial support under grant No. GR 1717/5-1, the Stellar Evolution group at the Argelander Institute for Astronomy in Bonn for sharing their expertise in stellar modelling and providing computational resources, and the anonymous referee for reviewing this work.

References

All Figures

Fig. 1. Observed luminosity distribution of single stars in 30 Dor compared to theoretical predictions. The histogram shows the observed numbers of single stars in the VFTS sample as function of luminosity (blue) compared to the theoretical distribution resulting from our VMS population synthesis models described in Sect. 4 (green). Both distributions are plotted transparent so that overlaps appear in dark green.

In the text

Fig. 2. Wind-momentum – luminosity relation for O and WNh stars in 30 Dor. Different symbols indicate empirical results for the most massive stars in 30 Dor from Bestenlehner et al. (2014) in blue, O-type giants and supergiants from Ramírez-Agudelo et al. (2017) Tables C.4 (dark green) and C.5 (light green), and O-type dwarfs from Sabín-Sanjulián et al. (2017) Tables A.1 (red) and A.2 (pink). Upper limits are indicated in grey. Evolutionary tracks are indicated by coloured lines analogous to Fig. 5 but restricted to the O star range. Predicted absolute numbers of stars according to our population synthesis models are mapped and colour coded as indicated by the colour bar on the right.

In the text

	Fig. 3. Wind-momentum – luminosity relation as in Fig. 2 compared with population synthesis models based on standard MESA models using the DUTCH wind scheme.
In the text

	Fig. 4. Mass loss as a function of Eddington factor. Mass-loss rates Ṁ as a function of the estimated classical Eddington factor for chemically homogeneous stars ${\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{hom}}$ (see text) are indicated analogous to Fig. 2. The Ṁ(Γe) relation for R 136 from Bestenlehner (2020) is indicated by the dashed curve.
In the text

	Fig. 5. Theoretical HRD for VMS models with standard MLT. Coloured areas indicate the expected number of stars according to our population synthesis models for 30 Dor (cf. Sect. 4). Evolutionary tracks for Ω/Ωcrit = 0.3 are indicated, labelled with ZAMS masses in M⊙.
In the text

	Fig. 6. Observed HRD for 30 Dor. Coloured areas indicate the observed number of stars, evolutionary tracks are the same as in Fig. 5.
In the text

	Fig. 7. Influence of envelope inflation. The theoretical HRD for models with reduced envelope inflation effect through an artificially enhanced convective efficiency (MLT++; MESA option okay_to_reduce_gradT_excess = .true.) is indicated analogous to Fig. 5.
In the text

	Fig. 8. Effective temperature distribution. The theoretical HRD for the upper main sequence is indicated analogous to Fig. 5 but as a function of the effective temperature relative to the ZAMS temperature.
In the text

	Fig. 9. Observed HRD analogous to Fig. 8.
In the text

	Fig. 10. Influence of mass loss. The theoretical HRD is indicated analogous to Fig. 8 but adopting a pure OB-type mass-loss prescription without enhanced Γ-dependent mass loss.
In the text

	Fig. 11. Influence of mass loss. The theoretical HRD is indicated analogous to Fig. 8 but adopting the DUTCH MESA wind scheme.
In the text

	Fig. 12. Effective temperature distribution. Histogram of the observed numbers of single stars in the VFTS sample with 5 < log(L/L⊙) < 6 as a function of their distance from the ZAMS (log(Teff/TZAMS)) (transparent blue) compared with the theoretical distribution from our VMS population synthesis models (transparent green).
In the text

Fig. 13. Helium surface enrichment through mass loss. The predicted helium surface mass fraction Ys, as a function of the estimated inverse mass-loss timescale Ṁ/Mhom for stars with Teff > 30 kK and log(L/L⊙) > 5, is compared with VFTS O and WNh stars labelled as in Fig. 2. The colour scale indicates the number of objects predicted by our population synthesis. Evolutionary tracks from Fig. 5 are indicated by coloured lines.

In the text

Fig. 14. Helium surface enrichment through rotation. The predicted helium surface mass fraction Ys, as a function of rotational velocity for O stars with Teff > 29 kK and 5.0 < log(L/L⊙) < 5.9, is compared with observed Ys and v sin i of He-enriched VFTS O stars (with Ys > 0.35) labelled as in Fig. 2; stars with log(L/L⊙) < 5.0 are plotted in grey. The colour scale indicates the number of objects predicted by our population synthesis models. Evolutionary tracks for an initial mass of 30 M⊙ are indicated by coloured lines. Step functions indicate the initial (blue) and present (orange) total number distributions in arbitrary units, as well as the present distribution of He-enriched stars with Ys > 0.35 (red).

In the text