© ESO, 2014

1. Introduction

Very massive stars (VMS) are defined to be initially more massive than 100 M_⊙ (cf. Vink et al. 2013). Whilst a lot of attention has been devoted to the role of VMS at extremely low metallicity in the very early Universe (e.g. Bromm et al. 1999; Abel et al. 2002), the role of VMS in our local Universe has largely been ignored until recently. Hamann et al. (2006) found the existence of a number of highly luminous Galactic objects (up to 10⁷L_⊙). Crowther et al. (2010) confirmed this finding, reporting the existence of VMS up to 320 M_⊙ for the R136 core stars in the Tarantula nebula of the Large Magellanic Cloud (LMC). Given the clustered nature of these objects it was important that Bestenlehner et al. (2011) identified an almost identical twin of R136a3, VFTS 682, in the context of the VFTS. Its isolated nature increased confidence that VMS really exist above the “canonical” upper limit of about 150 M_⊙ (Figer 2005).

Interestingly, the presence of VMS of up to 320 M_⊙ would bring them into the predicted initial mass range of pair-instability supernovae of 140–260 M_⊙ (Heger & Woosley 2002; Langer et al. 2007). As such explosions would disrupt the entire star, all metals would be returned to the interstellar medium (ISM). One such object at the top of the mass function of a massive-star cluster would produce more metals than all the other lower mass stars together (Langer 2009). It could significantly affect chemical evolution modelling of galaxies.

However, the evolution and fate of VMS are highly uncertain. The main culprit is the unknown rate of wind mass loss, given that VMS are thought to evolve close to chemically homogeneously (e.g. Gräfener et al. 2011). Whilst recent evolution models have been published by e.g. Yungelson et al. (2008) and Yusof et al. (2013), the way mass loss has been implemented in these models is entirely different and it is not known a priori which way is more appropriate. The key issue to decide whether pair-instability supernovae are likely to exist in the present day universe is that of stellar mass loss.

Fortunately, there has recently been significant progress in our theoretical understanding of mass loss at the top of the Hertzsprung-Russell diagram (HRD). In particular, Vink et al. (2011) computed mass-loss rates for VMS up to 300 M_⊙ and discovered a “kink” in the mass loss versus Γ dependence, involving a relatively shallow slope for optically thin O star winds, but changing into a steeper slope for optically thick WNh Wolf-Rayet stars. This latter result appears to be in qualitative agreement with theoretical and observational studies of very massive WNh stars by Gräfener & Hamann (2008) and Gräfener et al. (2011). Furthermore, Vink & Gräfener (2012) calibrated currently used mass-loss rates for VMS at the transition point from Of to WNh. Alternative VMS wind models have been proposed by Pauldrach et al. (2012) where it is claimed that VMS winds remain optically thin. In other words, there are not only quantitative, but even qualitative differences in the behaviour of VMS stellar wind models. For these reasons, it is imperative to determine the VMS mass-loss behaviour empirically, thereby testing theoretical concepts.

This has now become possible with the VFTS, providing a sufficiently large dataset at the top end of the HRD. This data-set involves several groups of H-rich Wolf-Rayet stars (WNh), but also transition Of/WN stars, as well as very luminous Of and other O stars. These stars cover a large range in mass-loss rate (which we find to be of more than two orders of magnitude), which enables us to study the transition from optically thin O star winds to optically thick WNh Wolf-Rayet winds.

In Sect. 2, we give an overview of the observational data and describe our target selection. We describe the details of our spectral analysis regarding theoretical model computations and analysing techniques in Sect. 3. On the basis of the results we investigate in Sect. 4 the mass-loss properties, evolutionary stages, and ionising fluxes of our targets. Our conclusions are summarised in Sect. 5.

2. Observational data and target selection

The data used in this work have been obtained in the framework of the VLT-FLAMES Tarantula Survey (VFTS; Evans et al. 2011), a large spectroscopic survey of over 800 O-, B-, and Wolf-Rayet (WR) stars in the 30 Doradus region of the LMC. From this data-set we selected 62 massive and VMS candidates. The data have been augmented by existing archival data and follow-up observations in the infrared (IR) range. In the following we give an overview of the spectroscopic and photometric data used (in Sects. 2.1 and 2.2), and describe the target selection process (Sect. 2.3).

2.1. Spectroscopic observations

Within the VFTS project the entire 30 Doradus region, centred on the dense cluster R136, has been covered with optical fibre spectroscopy using the MEDUSA-GIRAFFE (λ4000–7000 Å) mode of the FLAMES instrument. Roughly 70% of the stars with V< 17 mag have been observed in this programme. For the densely populated region near R 136 additional integral-field observations were obtained using FLAMES/UVES (λ4200–6100 Å) and ARGUS/GIRAFFE (λ4000–4600 Å).

Specifically for the analysis of VMS in the central region of 30 Dor, K-band data were obtained with SINFONI the Spectrograph for INtegral Field Observations in the Near Infrared on the VLT (Program ID 084.D-0980, PI: Gräfener). The spectra cover the central arcsec around R 136 with a seeing constraint ≤0.8 arcsec and a spectral resolution of R = 4000 over the wavelength range of λ2.0 − 2.4 μm. Due to the increased IR wind-emission of VMS compared to lower-mass stars with weaker winds, these data provide a homogeneous set of wind and abundance diagnostics for these objects that is less prone to crowding issues than e.g. optical data. The data were reduced using the standard ESO pipeline. The observed telluric standards were corrected for their intrinsic H and He absorption features with the aid of synthetic telluric spectra from Smette et al. (2010) and Smette et al. (2014).

For some targets close to the centre of 30 Dor, archival data were available from VLT/UVES (Program ID: 70.D-0164, PI: Crowther) and the Hubble Space Telescope (HST) satellite, specifically HST/FOS (Program ID: 6417, PI: Massey), HST/STIS (Program ID: 7739, PI: Massey), HST/GHRS (Program ID: 3030, PI: Ebbets), and HST/COS (Program ID: 11484, PI: Hartig). Some Of/WN and WNh stars had archival data from the International Ultraviolet Explorer (IUE) available. The additional optical archival observations are chosen in a way such that we can use the same diagnostic lines for all our targets, as to allow for a homogeneous spectroscopic analysis. If one diagnostic line was not observed, the star was removed from the discussion (see Sect. 4). However, H_α and He iiλ4686 could be substituted by the near-IR lines He ii at 2.165 μm and Br_γ if they show sufficient emission, or vice versa (see Appendix C). The UV-data are only used to constrain the terminal velocity. An overview of the data used for each target is given in Table B.1. Aliases of the VFTS stars are given in Table B.2.

2.2. Supplementary photometry

To constrain stellar luminosities it is necessary to model their spectral energy distribution (SED) including interstellar extinction. For this purpose, we use optical and near-IR photometry. Optical photometry has been obtained by Evans et al. (2011) using the Wide Field Imager (WFI) at the 2.2 m MPG/ESO telescope at La Silla. Near-IR photometry is provided by the VISTA Magellanic Clouds (VMC) Survey (Cioni et al. 2011). Specifically, we use PSF-fitting photometry from Rubele et al. (2012) and cross-matches to the VFTS sources from Zaggia et al. (in prep.). For objects that are not included in the VMC survey we use the InfraRed Survey Facility (IRSF) Magellanic Clouds catalogue (Kato et al. 2007) as listed in Evans et al. (2011). For stars near the centre of the R 136 cluster a higher spatial resolution is required due to crowding. In this case we utilise optical HST-WFC3 photometry from De Marchi et al. (2011) and near-IR H/K_s-Band VLT/MAD photometry from Campbell et al. (2010). The photometry from Evans et al. has been transformed from the original WFI system to Johnson colours. In the presence of strong emission lines this can lead to large uncertainties. Furthermore uncertainties can arise from crowding in strongly populated fields. Alternative photometry for objects with such uncertainties or objects that are missing in the catalogue of Evans et al. (2011) is given in Table 1.

2.3. Target selection

Our interest in this study is focussed on the most massive stars in 30 Dor sampling the top of the main-sequence, including VMS with masses in excess of ~100 M_⊙. We are interested in objects that are suitable for an emission-line analysis and are brighter than ~10^5.5L_⊙, and we avoided further selection criteria for the O stars to avoid a biassed sample. That said, irrespective of their brightness, we included all of the most promising VMS candidates (namely the Of/WN and WNh stars from Evans et al. 2011).

From the sample of 360 O stars (Evans et al. 2011; Walborn et al. 2014) we pre-selected stars showing emission-lines in He iiλ4686, N iii λ4634/4640, N ivλ4058, or H_αλ6562. Based on our theoretical models (Sect. 3.1) we divided the sample into three temperature classes to estimate their bolometric corrections (class 1: T_eff ≲ 40 kK, class 2: 40kK ≲ T_eff ≲ 45 kK, class 3: T_eff ≳ 45 kK). For each class intrinsic colours B − V and V − K_s are derived from our models. The B, V, and K_s-band magnitudes are determined by applying the filter function for the MPG/ESO/WFI B and V-band filters, and the 2MASS (2 Micron All Sky Survey) K_s-band filter. The obtained intrinsic colours are then used to estimate the extinction parameters E(B − V) and E(V − K_s) from the observed optical and near-IR photometry. As outlined in Bestenlehner et al. (2011) we estimate A_V and R_V for each object adopting the extinction law from Cardelli et al. (1989). For the absolute visual magnitude M_V we adopt a distance modulus of 18.49 mag (Pietrzyński et al. 2013). Stellar luminosities are then obtained using bolometric corrections from Martins & Plez (2006) for each temperature class. Stars with estimated luminosities ≲ 10^5.5L_⊙ are excluded from our sample.

For the stars near R 136 that are only covered by ARGUS integral-field spectroscopy we further added objects based on the presence of emission lines for He ii at 2.165 and 2.189 μm and Br_γ at 2.165 μm in our SINFONI data. We expect that these objects cover a similar mass-loss range as the ones selected on the basis of optical emission lines from the MEDUSA and UVES observations. Figures 1 and 2 show the location of our targets in the Tarantula Nebula.

In total we selected 62 stars combining the VFTS-MEDUSA, ARGUS and UVES, the VLT-SINFONI, and archival HST and IUE observations. The sample is deliberately heterogeneous covering a large range of emission-line strengths, and including a variety of luminosity classes and spectral types. A summary of the selected targets and the available spectroscopy is given in Table B.1. The table also indicates stars that are or may be spectroscopic binaries, higher order multiple systems, or contaminated as a result of crowding.

3. Spectral analyses

In this section we describe the spectral analyses of the 62 objects in our sample. In Sect. 3.1 we describe the model computations performed to produce a large grid of cmfgen non-LTE atmosphere/wind models. In Sect. 3.2 we introduce our method of analysis which combines the classical, subjective approach to fit spectra by eye with a more systematic χ²-approach, followed by a discussion of important diagnostic spectral lines in Sect. 3.3. In Sect. 3.4 we describe our method to fit the observed SED, including the effects of interstellar extinction. The uncertainties arising from our analysis are discussed in Sect. 3.5, and comparisons with previous results are performed in Sect. 3.6.

3.1. Model grid

The atmosphere models used in this work are computed with the cmfgen code by Hillier & Miller (1998), and comprise complex model atoms for H i, He i-ii, C iii-iv, N iii-v, O iii-vi, Si iv, P iv-v, S iv-vi and Fe iv-vii (cf. Table B.3). The models take advantage of the concept of super-levels by Anderson (1989), where the statistical equations of complex atoms are solved by combining many close-by atomic energy levels into few super-levels.

The parameters of our main model grid are described in Table 2. We varied three fundamental parameters, effective temperature T_eff, transformed mass-loss rate Ṁ_t and helium surface mass fraction Y over their relevant parameter range. Furthermore we took wind clumping with a volume-filling factor f_v into account by computing two identical model grids assuming a clumped wind with a volume-filling factor f_v = 0.1 starting at a velocity of 100 km s^-1, and an unclumped wind with f_v = 1¹. We computed one sub-grid for unprocessed material with Y = 0.25 and metal abundances corresponding to half solar metallicity based on solar abundances from Asplund et al. (2005). Furthermore, we computed multiple sub-grids for material that has been processed in the CNO-cycle with a nitrogen-enriched composition and He mass fractions between 0.25 and 0.7 (cf. Table B.4).

In our model grids we keep the luminosity and terminal wind speed fixed to values of L = 10⁶L_⊙ and ν_∞ = 2800 km s^-1. The use of the transformed mass-loss rate Ṁ_t, as defined by Gräfener & Vink (2013), allows us to scale our results to other values of L and ν_∞, and to correct for wind clumping with an arbitrary volume-filling factor f_v. For a synthetic emission-line spectrum with given Ṁ_t the true mass-loss rate Ṁ follows from the relation (1)

We note that this scaling relation is equivalent to the approach using the “transformed radius” by Schmutz et al. (1989) and Hamann & Gräfener (2004).

Furthermore, the model structure is computed assuming hydrostatic equilibrium with a fixed surface gravity of log g = 4 in the inner part of the atmosphere, and a β-type velocity law with a fixed β-parameter of 1.0 in the outer wind. The reference radius (R_ref) and effective temperature (T_eff) of the star are defined at an optical depth of τ = 2/3. The stellar temperature (T_⋆) and radius (R_⋆) of the inner boundary of the atmosphere are defined at τ = 100.

Additional grid computations were necessary to cover the parameter range of WNh stars because of their higher mass-loss rates. The WNh grid covers the range of log  (Ṁ_t/ (M_⊙ yr^-1)) from −4.61 to −4.11 in steps of 1/6 dex with helium mass fractions of 0.625, 0.7, 0.775, 0.85 and 0.925. Furthermore, a small sub-grid was computed for late WN6-8h stars, because their terminal velocities are significantly lower. This sub-grid is computed for a terminal velocity of 1800 km s^-1 and covers the temperature log T [K] from 4.550 to 4.625 in 0.025 dex steps and the transformed mass-loss range log  (Ṁ_t/ (M_⊙ yr^-1)) from −4.61 to −3.61 in 1/6 dex steps with an additional helium abundance of 99.0% compared to the WNh grid. Because of convergence difficulties some of the WNh models required a lower log g of 3.75. As the winds of WNh stars are optically thick this choice has no effect on their emission-line spectra.

As the determination of temperatures and mass-loss rates for O stars is sensitive to the surface gravity we computed sub-grids in the low mass-loss and temperature range with log g = 3.75 and 4.25. We used these models to quantify the impact of log g and estimate the resulting uncertainties. The sub-grids cover a transformed mass-loss range log  (Ṁ_t/ (M_⊙ yr^-1)) from −6.28 to −4.78, temperature range from log T [K] 4.55 to 4.65 and helium mass fractions of 0.25 and 0.325.

In addition, some models with log g = 3.5 and with intermediate helium mass fractions were computed to investigate the influence of the assumed log g on the derived helium abundances. The results of these tests are discussed in Appendix D.

Together with additional test models more than 2500 models were computed.

3.2. Method of analysis

The stellar sample analysed in this work covers a broad parameter range, reaching from O stars with a combination of weak emission and absorption lines to WR stars with their strong emission line spectra. Especially for the latter it is necessary to perform time intensive non-LTE computations with complex numerical codes such as cmfgen. As we are interested in a homogeneous analysis of the transition between the O and WR stars we chose to analyse the complete sample using the complex cmfgen models. A drawback of this approach is that the number of parameters that we are able to vary in our analysis is limited. Our results will thus be less accurate than e.g. the ongoing O star analyses in the VFTS collaboration that are based on the much faster code fastwind by Puls et al. (2005), but they are sufficient for our investigation of the trends of the wind parameters. (The future works will be based on tailored analyses for each object involving more extensive grids and/or advanced fitting methods, however, they are limited to the O star range.)

In this work we use the grid computations described in Sect. 3.1 and try to quantify the quality of the fit by introducing a tailored merit function () that mimics a least-square fit method. This method is based on measuring equivalent widths (EWs) and is thus independent of the adopted terminal wind velocity (ν_∞) in our grid², the projected rotational velocity (ν_esin i), and (micro) turbulent velocity at the stellar surface (see Sect. 3.1). However, due to the variety of spectral types in our sample and their different line diagnostics we found it necessary to perform a final visual assessment of the fit quality thus introducing a subjective component in our approach.

The following continuum subtracted line profiles are taken into account for the analysis: the N iii doublet λ4634/4640, N iv λ4058, the N v doublet λ4604/4620, He i λ4471, He ii λ4686, He ii λ2.189 μm and the Balmer lines H_{α − δ}, and Br_γ. To ensure that the method is independent of the continuum normalisation, we determine the continuum locally and subtract it in exactly the same way from the observations and the synthetic spectra. The EWs are obtained by defining two continuum points in the blue and red direction of the diagnostic line, fitting a straight line through the points, subtracting it as the continuum, and integrating over the spectra (consequently absorption lines have negative, and emission lines positive EWs). Complications can still arise from a contamination of the observed spectra with nebular emission lines. This is usually the case for He i λ4471 and the Balmer lines H_{α − δ}. In these cases we remove the nebular features from the observed spectrum and use a Gaussian fit to the line wings to estimate the line profile shape in between. Again, this is done in the same way for observations and models. Lines that are strongly affected by nebulosity are excluded from our analysis. The fit quality of the grid models can then be quantified for each object by computing the sum (2)where i indicates the different diagnostic lines. Here the division by the square root of the observed EW balances the weight between stronger and weaker lines. The constant (which is of the order 10^-4) prevents divisions by zero and gives a lower relative weight to weak lines that are more affected by noise.

The final decision about the best-fitting model is done by eye, considering the models with the smallest . To this purpose the model spectra are convolved with instrumental and rotational profiles. This allows us to roughly constrain the projected rotational velocities ν_esin i from the observed line profiles. Investigating the parameter range with small also helps to identify possible ambiguities in the parameters that can reproduce the observed spectra.

3.3. Line diagnostics

The determination of mass-loss rates and helium abundances is based on strong H and He lines, mainly He ii λ4686 and H_α in the optical and, if available, He ii λ2.189 μm and Br_γ in the near-IR. For high mass-loss rates other lines of the Balmer series are also useful (H_{β − δ}). As the strength of these lines also depends on T_eff the accuracy of our results depends on how well T_eff can be constrained from other diagnostic lines. The emission line strengths of H_α, Br_γ, and He ii λ2.189 μm roughly show a linear dependence on T_eff. He ii λ4686 shows a much stronger sensitivity at certain temperatures (cf. Fig. C.1) and also depends on gravity (Fig. D.2). In the absence of IR data large uncertainties in the derived He abundance may occur. In Fig. C.2 we show the sensitivity of the diagnostic lines with respect to the helium mass fraction. The mass-loss rates determined from pure emission lines scale with the square root of the adopted volume-filling factor (cf. Eq. (1)). In contrast to this, absorption lines, which form below the region where we assume clumping to start, are affected by the actual mass-loss rate and are thus independent of clumping.

To constrain T_eff, we use He i λ4471, as well as N iiiλ4634/4640, N iv λ4058, N vλ4604/4620. The He i λ4471 absorption line is a classical temperature indicator for O stars, however, it is only available for temperatures up to T_eff ≈ 45 kK. For this reason we use the nitrogen lines as an additional indicator. The ionisation stages of nitrogen cover N iii to N v and can thus be used over the whole temperature range investigated here. However, lines from these ions are usually in emission, which means their analysis is complicated by the effects of mass-loss and abundance variations (cf. Rivero González et al. 2012a). Moreover, their strength depends on detailed atomic physics, such as dielectronic recombination channels (Mihalas 1973). Finally, line overlaps between N iiiλ4634/4640, N vλ4604/4620 and He ii λ4686 can lead to difficulties, particularly for high mass-loss rates.

The terminal wind velocities (ν_∞) are best derived from the P-Cygni profiles of UV resonance lines, in particular the C ivλ1548/1551 doublet. For objects without available UV spectroscopy we estimate ν_∞ from the width of the H_α emission line by a detailed comparison with synthetic line profiles. For VFTS 016, 482, 1022, and MK42 we could use the traditional black absorption edge method (e.g. Prinja et al. 1990; Prinja & Crowther 1998) to test the ν_∞ measurements based on spectral modelling for H_α and the C iv doublet. For WNh stars of late subtype (WN6-8h) the He i λ4471 line has a P-Cygni type profile and can be used to derive ν_∞ with similar accuracy as from the C iv λ1548/1551 doublet. For O and Of stars H_α becomes too weak to estimate ν_∞. In these cases we adopt the values from the relation by Lamers et al. (1995) (cf. Sect. 4.3).

To constrain the surface nitrogen enrichment due to the CNO-cycle we performed a subjective comparison of the overall strength of the observed nitrogen features with our models. We tentatively classified each object as N-enhanced (e), N-normal (n) or partially enhanced (ne). The results are listed in Table 3. For He-enriched stars (Y> 0.25) we found that the strength of the observed nitrogen features is generally in good agreement with an enhanced N abundance (e), as adopted in our models (cf. Table B.4). For stars with normal He (Y = 0.25), we computed N-normal and N-enhanced models (cf. Table B.4) and found that the observed nitrogen features indeed reflect a variety of nitrogen abundances.

3.4. Luminosity and interstellar extinction

To determine the bolometric luminosities (L_⋆) of our sample stars we use the model flux as obtained from our grid computations (which is computed for a luminosity of 10⁶L_⊙) and scale it to match the observed SED in the B, V and K_s bands (cf. Sect. 2.2). Simultaneously we determine the extinction parameters R_V and A_V as described in Sect. 2.3, i.e., we compute the intrinsic brightness in the B, V and K_s bands from the model flux using appropriate filter functions³ and match the resulting values for E(B − V) and E(V − K_s) as described by Bestenlehner et al. (2011). The results are given in Tables 3 and 4.

As the extinction is the smallest in the near-IR range the resulting luminosities rely chiefly on the K_s-band photometry. R_V and A_V mainly follow from B and V in relation to K_s. This method thus leads to very reliable L_⋆ as long as the K_s magnitude is not contaminated by other sources, such as thermal dust emission or nearby stars. To avoid such instances we generally compare the reddened model flux for each object with the complete available photometry, in some cases reaching from the UV to the far infrared (cf. Appendix E). A similar method has previously been used by Hamann et al. (2006) for the analysis of Galactic WR stars. As visual and UV fluxes are strongly affected by deviations from a standard extinction, i.e. R_V = 3.1, this method helps to avoid large errors in the luminosity determination that can arise from using these wavelength ranges alone. A good example for such a case is VFTS 682 for which Bestenlehner et al. (2011) obtained log L/L_⊙ = 6.5 based on the derived R_V = 4.7, while a standard extinction law with R_V = 3.1 would suggest a much lower luminosity of log L/L_⊙ = 5.7.

Uncertainties in our analysis chiefly arise from contributions of unresolved neighbouring stars. Such occurrences are further complicated by differences in the spatial resolution of optical and near-IR photometry. For emission-line stars the use of broad band photometry can cause problems if the line contribution is not correctly taken into account. However, as we know the line contribution from our models and take it into account using appropriate filter functions these uncertainties should be minimised. In some cases we still obtain unrealistic R_V from our analysis (e.g. VFTS 402 and 147). These problems are most likely caused by unresolved companions. Stars for which this happens are excluded from our discussion (cf. Sect. 4.4) and their parameters are given in Table 4. As these problems only affect the optical extinction the derived luminosities may still give reasonable estimates for the main contributor in the IR, where the extinction is low.

Furthermore, many of our sample stars show an apparent UV mismatch in the SED plots in Appendix E. In the majority of cases this is caused by our use of a Galactic extinction law over the full wavelength range from UV to IR. Specific LMC extinction laws (as e.g. from Howarth 1983) are indeed much more appropriate for the UV range and lead to better UV fits. In particular the UV 2200 Å feature appears to depend strongly on metallicity and is much weaker in the LMC. However, deviations from the standard extinction with R_V = 3.1 have only been thoroughly investigated for the Galaxy (e.g. Cardelli et al. 1989). This particularly affects the transition from the IR to the optical range and the overall slope of the extinction law in the UV, which is why we use the extinction law by Cardelli et al. in this work. The UV mismatch does not affect our results as these are only based on optical/IR photometry, and is thus a purely cosmetic issue following from the adopted extinction law.

3.5. Results and error discussion

The results of our analyses are compiled in Tables 3 and 4. To discuss the uncertainties we need to consider that the analysis is performed in two steps. In step one the line spectrum is analysed using the method described in Sect. 3.2. In this step the following parameters are determined: the transformed mass-loss rate (Ṁ_t), the effective temperature (T_eff), the surface helium mass fraction (Y), and the terminal wind velocity (ν_∞). In step two the luminosity (L_⋆) and the extinction parameters A_V and R_V are determined by fitting the observed SED (cf. Sect. 3.4). The actual mass-loss rates (Ṁ) then follow from the scaling relation Eq. (1), which additionally depends on the adopted volume-filling factor f_v. The uncertainties in T_eff, Y and Ṁ_t are dominated by the limited resolution of our model grid. We estimate values of Δlog T_eff = ± 0.02 dex (corresponding to ±2 to 2.5 kK), ΔY = ± 5% and ΔṀ_t = ± 0.1 dex. For the terminal wind speed the measurement error is of the order of Δν_∞ = ± 200 to 400 km s^-1.

Moreover, systematic errors are introduced in our grid analysis because the wind acceleration parameter β and the surface gravity log g are fixed. By changing β from 1.0 to 1.5 (as suggested e.g. by model computations by Vink et al. 2011) the wind acceleration zone is extended and the wind density is increased in the line-forming region. To compensate for this effect Ṁ_t needs to be reduced by 0.1 to 0.2 dex. The slope of the SED is mainly affected in the optical range, impacting R_V but leaving the luminosity almost unchanged (Δlog L ~ + 0.005 dex).

Changes due to variations in log g are investigated in Appendix D. A change in log g can influence the density and ionisation balance in the hydrostatic part of the atmosphere. For stars with optically thin winds this may affect the derived helium and nitrogen abundances. The impact on the latter turns out to be negligible, mainly because we only coarsely estimate the nitrogen enrichment based on the overall strength of the nitrogen lines (see last paragraph in Sect. 3.3). For stars with optically thick winds changes of the hydrostatic structure have practically no effect.

In Fig. D.1, we show the effect of lowering log g by 0.5 dex in a critical temperature range around T_eff = 40 kK. In this range lowering the gravity decreases the strength of the main temperature diagnostic He i λ4471. Consequently a lower temperature by ~2500 K is needed to reproduce this line. While the transformed mass-loss rate remains almost unaffected the derived helium mass fraction is lowered by 2 ... 4%. This systematic uncertainty affects our given errors of ΔY = ± 5% only marginally. As a result of the lower T_eff we derive a lower luminosity. Moreover, the terminal wind speed ν_∞ is affected for the cases where it is derived from the relation by Lamers et al. (1995). The resulting change of the absolute mass-loss rate Ṁ as derived from Eq. (1) is of the order of 0.15 dex and is taken into account in the given error range. The resulting systematic changes are indicated by a vector in Fig. 6 (see Sect. 4). We note that they only affect the O supergiants in our sample, as these objects have significantly lower log g than 4.0.

As we are mainly interested in a homogeneous analysis of the mass-loss behaviour across the transition region between O stars, Of/WN, and WNh stars, and it is not possible to determine log g for the WNh stars, we decided to use the same analysis method for all objects in our sample, i.e., to keep log g fixed. A detailed determination of log g and N-abundances of all O stars in the VFTS sample will be performed in future works (Ramirez-Agudelo et al.; Sabín-Sanjulián et al., in prep.).

The derived luminosities depend on the photometry (mainly in the near-IR) and the model SED, which depends on T_eff and Ṁ_t. The photometric errors are of the order of ±0.05 dex. Considering the resolution of our model grid, the resulting uncertainties are of the order of Δlog L_⋆ ≈ ± 0.1 dex.

In summary, we estimate uncertainties of Δlog Ṁ_t = ± 0.1 dex, Δlog T_eff ≈ ± 0.02 dex, Δν_∞ ≈ ± 200 to 400 km s^-1 (see also Sect. 4.3) and Δlog L_⋆ ≈ ± 0.1 dex. The resulting mass-loss rates Ṁ according to Eq. (1) are uncertain by Δlog Ṁ ≈ ± 0.2...0.3 dex considering the uncertainties in log g and β, and additionally scale with the square root of adopted volume-filling factor .

3.6. Comparison with previous works

In our present work we tried to follow a well-defined method of analysis as described in Sects. 3.2–3.4. In this section we present comparisons with previous works to get a better picture of the systematic uncertainties that can arise from different analysis methods.

A comprehensive study of WR stars in the LMC has recently been published by Hainich et al. (2014), who have analysed 17 stars from our present sample. These authors used a similar grid approach as in our work, based on the Potsdam Wolf-Rayet (PoWR) model atmosphere code (Koesterke et al. 2002; Gräfener et al. 2002; Hamann & Gräfener 2003). An important difference with respect to our work is how the stellar luminosities and interstellar extinction are determined. As described in Sect. 3.4 our results mainly rely on IR K_s-band photometry. The optical flux in B and V is matched simultaneously with the IR by adapting R_V and E(B − V) in the extinction law. Hainich et al. keep R_V fixed and mainly use the optical to UV flux to determine E(B − V) based on the LMC extinction law by Howarth (1983). As discussed in Sect. 3.4 this approach may lead to substantial errors if R_V deviates from the adopted value. The analysis of Hainich et al. largely relies on the UV range where the uncertainties in the extinction are very high. Our results rely on the IR where the extinction is almost negligible, however, there is a danger of contamination by circumstellar dust or other IR sources.

An illustrative case is VFTS 482 (BAT99-99, Mk 39). For this star Hainich et al. obtain a luminosity of log  (L/L_⊙) = 5.9 while our value lies significantly higher (log  (L/L_⊙) = 6.4). As described above Hainich et al. do indeed substantially underestimate the IR flux for this object. However, there is possible evidence for crowding in HST images of this star (cf. Table B.1). Could the excess IR flux thus be due to other sources? Our SINFONI IR data of VFTS 482 suggests that this is not the case as the IR spectrum shows the correct line strength, i.e., it is unlikely that the IR is contaminated by other sources. The same holds for the optical range for which a comparison between HST and UVES spectroscopy shows no sign of contamination. For VFTS 482 we are thus confident that our approach is correct.

Furthermore, the different wavelength coverage and larger signal-to-noise ratio (S/N) ratio of the VFTS data can explain discrepancies arising from differences in the line diagnostics. This partly leads to different He-abundances (VFTS 482 and 545) and temperatures (ΔT ≳ 10 000 K, e.g. for VFTS 545, 617, 1017, and 1025). After correcting for clumping, there is also a systematic offset of 0.1 to 0.2 dex in the mass-loss rates noticeable. This difference lies within our error bars and is most likely caused by different assumptions in the model physics and atomic data.

VFTS  1025 (R 136c) has previously been analysed by Crowther et al. (2010). This star may be affected by crowding as it is located near the core of R 136. Due to the high S/N in the VFTS data we detect a weak photospheric He i λ4471 absorption line in the optical spectrum which is reproduced by our models. Based on this line we determine a lower temperature for this object than Crowther et al. (2010). It is not clear whether this line is intrinsic or stems from a nearby star (see the detailed discussion in Appendix A). We obtain T_⋆ = 42 ± 2 kK, log  (L/L_⊙) = 6.6 ± 0.1 while Crowther et al. obtain T_⋆ = 51 ± 5 kK, log  (L/L_⊙) = 6.75 ± 0.1. The obtained mass-loss rates and He-abundances are very similar. The differences have practically no effect on our results in Sect. 4.

VFTS 016 has been analysed by Evans et al. (2010). Their model has a lower surface gravity which results in a lower T_eff and luminosity compared to our values. The mass loss and He-abundance are in good agreement.

VFTS 072 has been analysed by Rivero González et al. (2012b). They derived a higher T_eff that might be a result of the higher gravity in their models. Mass-loss rate, luminosity and He-abundance are in good agreement.

VFTS 482 has also been analysed by Massey et al. (2005). We derived a ~4000 K higher T_eff. The luminosity is higher by about 0.6 dex. Massey et al. (2005) used a standard R_V = 3.1 and estimated the extinction E(B − V) by averaging the colour excesses in B − V and U − B based on the spectral type. In this work we calculate the R_V (see Sect. 3.4). The mass loss is higher by 0.3 dex while the He-abundance is in agreement.

Mk 42 has been analysed by Puls et al. (1996). The derived T_eff is about 3000 K higher while the He-abundance is a bit lower compared to our values. The luminosity agrees and our mass-loss rate is about 0.1 dex lower.

We conclude that the stellar parameters derived in this work largely agree with previous studies with some important exceptions due to systematics. The main systematic uncertainties arise from the uncertain interstellar extinction (in particular because of its dependence on R_V), and from problems with crowding due to the large distance to the LMC.

4. Discussion

In Sect. 4.1 we start our discussion with a division of our sample into sub-classes based on the spectral classifications in Table 3. The method by which we estimate masses and Eddington factors for each object is discussed in Sect. 4.2. In Sect. 4.3 we compare the observed terminal wind velocities obtained with different methods within this work with commonly adopted relations. In the subsequent sections we discuss the mass-loss properties of our sample stars, namely their wind momenta (Sect. 4.4), the dependence of the mass-loss rate on the Eddington factor (Sect. 4.5), and their wind efficiencies and clumping factors (Sect. 4.6). A comparison with VMS evolutionary models is performed in Sect. 4.7. The discussion closes with a brief section of the ionising fluxes (Sect. 4.8).

4.1. Definition of stellar sub-classes

In the following we divide our sample into sub-classes with different emission line strengths. Based on the spectral classifications in Table 3 we designate stars with spectral types WNh and WN(h) as “WNh stars”, stars with mixed Of/WN spectral types as “Of/WN stars”, O stars in categories f and (f) as “Of stars”, and O stars in categories ((f)) or without category as “O stars”. If we discuss O stars in general we designate them as “O-type stars”. Using these designations the group of O stars is dominated by dwarfs with low Ṁ_t, and the group of Of stars by giants and supergiants with considerably higher Ṁ_t.

4.2. Stellar masses and Eddington factors

To characterise the proximity of our sample stars to the Eddington limit we use a method introduced by Gräfener et al. (2011). The Eddington factor (Γ) is defined as the ratio between radiative acceleration (g_rad) and gravitational acceleration (g), Γ = g_rad/g. In spherical symmetry it can be expressed as (3)Here χ(r) denotes the total flux-mean opacity per gram of material. In general, χ is a complex function of density and temperature. For this reason Γ depends on radius, i.e. it is impossible to find a single value of Γ which characterises the proximity of a star to the Eddington limit. For the hot atmospheres that we consider here H and He are largely ionised and electron scattering is a dominant contributor to the total opacity. As the electron scattering opacity χ_e depends only on the density of free electrons it is almost independent of the radius. For this reason we use the classical Eddington factor (4)in our discussion in Sect. 4.5. Assuming a fully ionised plasma χ_e depends only on the hydrogen mass fraction X_s at the stellar surface, and Γ_e can be expressed by the relation (5)In this equation L and X_s are determined in our spectral analyses, i.e. the only unknown is the stellar mass M.

To estimate stellar masses for given values of L and X_sGräfener et al. (2011) computed mass-luminosity relations for chemically homogeneous stars. In this context it is important to know that hot massive stars have large convective cores. In this case the luminosity of a star with given mass M is chiefly determined by the hydrogen mass fraction X_c in the convective core. As X_c is a priori unknown, Gräfener et al. made two extreme assumptions about X_c. The first assumption is X_c = X_s. In this case the star is chemically homogeneous, and the derived mass provides an upper limit for the real mass of the star. The second assumption is X_c = 0. In this case the hydrogen in the core is exhausted and the star is in the phase of core He-burning. The mass derived in this way provides a lower limit to the real mass of the star.

Chemical homogeneity is expected for stars with strong rotational mixing and/or strong mass-loss, i.e. when the timescales of mixing or mass-loss are shorter than the nuclear burning timescale. E.g., for the extremely massive WNh stars R136a1, a2, a3 and c, Crowther et al. (2010) estimate present day masses of 265, 195, 135 and 175 M_⊙ using evolutionary models. The corresponding homogeneous masses following Gräfener et al. (2011) are 286, 206, 154 and 185 M_⊙. The agreement is better than the error-margins on the stellar masses given by Crowther et al. Mass estimates following Gräfener et al. (2011) for our sample stars are given in Tables 3 and 4. A comparison with evolutionary models for VMS from Köhler et al. (2014) gives values that are on average ~30% lower. This is to be expected, in particular for O stars with H-rich surface compositions. Many of these stars will have evolved away from the zero-age main-sequence (ZAMS), i.e. they are not chemically homogeneous. The Eddington factors estimated here for individual stars are thus likely higher in reality, in particular for relatively low stellar masses.

The use of evolutionary masses for all stars in our sample turned out to be problematic, in particular for stars with helium-enriched surface compositions. In these cases the evolutionary models do not always match the observed abundances, leading to qualitatively wrong estimates of Γ_e following Eq. (5) (the resulting masses are partly higher than the upper limits obtained under the assumption of homogeneity). The strength of our method is therefore to obtain model-independent qualitative information about a Γ-dependence, based on observed luminosities and surface compositions of a large stellar sample. Nevertheless, we test the effect of our assumptions on our results for the H-rich O stars in our sample in Sect. 4.5.

4.3. Terminal wind velocities

The terminal velocities (ν_∞) given in Tables 3 and 4 are determined on the basis of different lines. Dependent on their availability we use P-Cygni line profiles of Civλ1550 or He i λ4471, or the width of H_α (cf. Sect. 3). If no diagnostic lines are available we use the empirical relation for O stars with T_eff> 25 000 K from Lamers et al. (1995), i.e. . Here the effective escape velocity () is defined as (6)with Γ_e from Eq. (5). R_ref is defined as . These values are also given in Tables 3 and 4.

In Fig. 3 we compare our results with this relation. The comparison suggests that the relation by Lamers et al. (1995) can be used for the Of, Of/WN, and WN5h stars in our sample (Fig. 3). For late WN6-8h stars we find lower values than predicted. That said, the relation from Lamers et al. is most suitable for the O and Of stars in our sample for which direct observational values are lacking.

The comparison indicates uncertainties in ν_∞ of ±200 to 300 km s^-1 using the C iv P-Cygni profile. Similar errors are obtained by deriving ν_∞ from the He i line for late WNh stars. A comparison of C iv and He i measurements shows a good agreement. The errors for measurements based on the width of H_α are of the order of ±300 to 400 km s^-1 depending on the quality of the model fit.

4.4. Wind momentum-luminosity relation

The wind momentum of a stellar wind is given by the product of its mass-loss rate and terminal wind velocity (Ṁν_∞). Kudritzki et al. (1999) introduced the modified wind momentum , as this quantity is expected to be almost independent of the stellar mass, i.e. it depends chiefly on the luminosity of the star. Puls et al. (1996) and Repolust et al. (2004) found a tight wind-momentum luminosity relation (WLR) for O stars and supergiants. Martins et al. (2008) studied the most massive stars in the Arches cluster and found two WLRs, for O stars and WNh stars. Mokiem et al. (2007) investigated the WLR for a sample of 28 OB stars in the LMC based on the data of the VLT-FLAMES Survey of Massive Stars (Evans et al. 2005). Their sample included stars in the luminosity range of log L/L_⊙ = 4.5...6.0.

In this work we focus on LMC stars with even higher luminosities in the range between log L/L_⊙ = 5.5 and 6.6, including Of/WN and WNh stars. The mass-loss rates given here are for the case without clumping, to facilitate a comparison with earlier works. This means that even if our modelling indicates a clumped wind structure, e.g. based on the electron scattering wings of some emission lines, we correct them for clumping according to Eq. (1). For our models with f_v = 0.1 this means that Ṁ has to be multiplied with a factor of .

In Fig. 4 we show the WLR for all stars in our sample. Stars whose spectra are affected by multiplicity (i.e. SB2’s and stars in crowded fields) are indicated by squares. Although for the following discussion these objects are removed from our sample, it is remarkable that even though these objects are clearly contaminated by other stars the resulting WLR is hardly affected. Black symbols indicate SB1’s or stars with weak RV variations. We assume that the spectra of these objects are dominated by one star and we treat them as single stars in our analysis. The corresponding errors in the derived luminosities are most likely negligible compared to the impact of the photometric errors.

In Fig. 5 we indicate stars with highly uncertain stellar parameters. Two of these stars have peculiar spectra with possible indications of fast rotation (VFTS 208 and 406). For the rest some essential diagnostic features are missing or are of poor quality (VFTS 562, 1014, 1026, 1028). In total we excluded six objects from the following discussion (see Appendix E for more details).

In Fig. 6 we show the WLR for the remaining sample indicating dwarfs, giants, and supergiants. To investigate the systematic effects of our choice of log  g = 4.0 in our grid computations we highlight the positions of O and Of (super)giants in the WLR. The arrow indicates the estimated shift in the WLR if we adopt a more realistic value of log  (g) = 3.5 for these objects (cf. Sect. 3.5). The resulting shift of –0.1 dex in log  (L) and –0.15 dex in the modified wind momentum is also confirmed by an individual model for the O8 Ib(f) star VFTS 669 with log g = 3.5 (cf. Figs. E.48 and E.49). Therefore the corresponding objects would move along the WLR towards lower luminosities and wind-momenta. Their qualitative position on the WLR is thus unaffected by the change in log  g, i.e. the (super)giants still form a group along the WLR, but with systematically lower wind momenta. For O dwarfs the effect may go exactly in the opposite direction as their gravities may be up to 0.2 ... 0.3 dex higher than our adopted value. Also those stars would remain on the same WLR.

In Fig. 7 we show the WLR for the combination of our sample with the sample of Mokiem et al. (2007). We indicate the four groups of O stars, Of stars, Of/WN stars and WNh stars (see last paragraph of Sect. 2.3) using different symbols. The majority of the O-type stars in both samples follow a clearly defined WLR. We note that the work by Mokiem et al. is based on a different analysis method and different numerical models than ours. However, the good agreement between both samples in the overlapping region between log L/L_⊙ = 5.5 and 6.0 supports the reliability of the comparison. For the stars at the highest luminosities, including the Of/WN, WNh and some Of stars from our sample, the observed wind momenta are clearly enhanced with respect to the overall WLR. In Sect. 4.5 we show that this mass-loss enhancement can be explained by the proximity of these stars to the Eddington limit.

The only object that clearly stands out in Fig. 7 is the WNh star VFTS 108. This star has a high helium surface mass fraction Y = 0.775 and may thus be in the phase of core He-burning. Also the other four WNh stars with Y> 0.75 (427, 695, 758, 1001) have increased modified wind-momenta and may already have reached the core He-burning stage.

Focussing only on the O-type stars we obtain an observed WLR of the form (7)with coefficients m₀ and C₀ given in Table 5. The observed relation is derived in two steps. The first step is to fit a simple linear fit through the O-type stars ignoring the uncertainties in the luminosity and wind momentum. The result is used in the second step as the starting point for a Levenberg-Marquardt fit including error ellipses (Markwardt 2009). The coefficients for the resulting relation are given in Table 5.

In Fig. 7 we compare this relation with the predictions by Vink et al. (2001) in the form given by Mokiem et al. (2007, coefficients are given in Table 5. As previously noted by Mokiem et al. (2007) the observed WLR for O-type stars at low luminosities lies slightly above the theoretical predictions, which can be resolved by adopting moderate wind-clumping to reduce the observed mass-loss rates. For O-type stars with higher luminosities, around log L/L_⊙ ~ 6, the values derived here agree well with the predictions, so that a further downward correction would lead to values below the theoretical expectations. While these effects are relatively small, and will be subject of future dedicated studies within the VFTS collaboration, we find a substantial increase of the wind momenta for the Of/WN and WNh stars at even higher luminosities, similar to the results of Martins et al. (2008) for the luminous WNh stars in the Arches cluster. This increase has been interpreted by Gräfener et al. (2011) as being due to the proximity of these objects to the Eddington limit, and is further discussed in Sect. 4.5.

The average mass-loss relation for the complete sample, including WNh, Of/WN, and O-type stars, would lie 0.2...0.3 dex above the predictions. Moderate volume-filling factors around f_v = 0.25 would thus bring the observed average WLR in agreement with the predictions (cf. also Repolust et al. 2004; Mokiem et al. 2007).

4.5. Mass-loss near the Eddington limit

A primary goal of this work is to investigate the mass-loss properties of VMS at the top of the main sequence, when they are still in the phase of core hydrogen-burning. Gräfener & Hamann (2008) and Vink et al. (2011) predicted a strong dependence on the classical Eddington factor Γ_e in this regime. Gräfener & Hamann (2008) also predicted a dependence on the environment metallicity Z, which mainly decides around which values of Γ_e the strong Γ-dependence sets in. Vink et al. (2011) investigated the transition from the classical WLR for O stars to the Γ-dependent mass-loss of VMS and found an abrupt change in the form of a “kink” between both regimes. Gräfener et al. (2011) confirmed the strong Γ-dependence for the sample of VMS in the Arches cluster near the Galactic Centre. They used the method outlined in Sect. 4.2 to estimate Eddington factors which take advantage of the characteristic dependence of the classical Eddington factor on the hydrogen surface abundance with Γ_e ∝ (1 + X) (cf. Eq. (5)).

In Fig. 8 we show the resulting Ṁ − Γ_e relation for our sample and the sample of Mokiem et al. (2007). To be consistent we also used masses based on the assumption of chemical homogeneity to calculate Γ_e for the sample of Mokiem et al. (2007). As predicted by Vink et al. (2011) we find two branches for high and low Γ_e. Separate fits for both branches are shown using a relation of the form (8)with m₁ and C₁ listed in Table 6. The fit parameters and errors are obtained from a Levenberg-Marquardt analysis.

The low-Γ branch consists only of O-type stars and has a relatively shallow slope with . As the hydrogen surface abundances of the stars on this branch do not vary, their Γ-dependence translates into a dependence on L/M (cf. Eq. (5)). As M is computed from L in our present approach (cf. Sect. 4.2) this translates into a dependence on L only. It is thus not surprising that we find such a relation for the stars on this branch. The stars on the high-Γ branch, on the other hand, are Of, Of/WN, and WNh stars with substantially varying hydrogen surface abundances. The Γ-dependence on this branch is thus a firm indication of a strong intrinsic dependence on the Eddington factor Γ_e. On the high-Γ branch we find a much steeper relation with . The five stars with Y> 0.75 have been excluded from the fit (see paragraph after next for details). The two branches intersect at log Γ_e = −0.58.

We further investigate whether the steep Γ-dependence can be explained within the standard wind theory for OB stars by Castor et al. (1975, CAK). Employing equation (46) from CAK in our fitting algorithm we obtain an (effective) force multiplier parameter α = 0.62 ± 0.08, in good agreement with typical values from the literature (Puls et al. 2008). The resulting relation between Ṁ/M and Γ_e is shown in Fig. 9. The relation matches the O-type stars in our sample fairly well, but the steep increase in Ṁ occurs only for Γ_eextremely close to unity (as a result of the 1/(1 − Γ) term in the CAK relation). In an attempt to match the steep observed increase for high Γ_e with a CAK-type relation we obtained better results for a reduced α = 0.53. However, in this case the mass-loss rates of normal O-type stars are systematically over-estimated. To explain the observations in the framework of the CAK-theory it is thus necessary to adopt a varying α-parameter, or two separate relations for O-type stars and Of/WN, WNh stars.

As outlined in Sect. 4.2 the Γ_e in Fig. 8 are obtained under the assumption of chemical homogeneity, i.e., they are lower limits. For the five He-rich stars with Y> 0.75 we also indicate upper limits for Γ_e based on the assumption that the stars are core He-burning. The most significant change occurs for VFTS 108. The assumption of core He-burning brings this star in agreement with the overall relation. For the other He-rich stars the situation is less clear, i.e., they could be in a core H or core He-burning phase. For VFTS 108 we obtain M = 20.5 M_⊙ and Γ_e = 0.46 under the assumption of core He-burning, and M = 33.1 M_⊙ and Γ_e = 0.29 assuming chemical homogeneity. The lower mass is also supported by the low terminal wind velocity of this star. Under the assumption of chemical homogeneity the escape velocity of this star would be higher by a factor 1.8 and the object would significantly stand out in Fig. 3. We thus conclude that VFTS 108 is most likely a core He-burning object in our sample. The situation is similar for VFTS 427 and VFTS 1001 which also have low terminal wind velocities (~1000 and ~1100 km s^-1). The wind velocities of VFTS 695 and VFTS 758 are higher (~1600 and ~2000 km s^-1) and in agreement with the relation by Lamers et al. (1995) adopting chemical homogeneity.

As noted in Sect. 4.2 the masses derived under the assumption of chemical homogeneity are upper limits, and the true masses of our sample stars are likely smaller. Comparisons with evolutionary models indicate an average difference of ~30% (Schneider et al., in prep.). This holds in particular for the stars with low mass-loss rates and normal He or N abundances for which chemical homogeneity is not expected from single star evolution. The true Eddington factors for these stars may thus be significantly higher than adopted in Fig. 8.

In Sect. 3.5 and Appendix D we discuss that systematically lower log g for the O supergiants in our sample imply slightly lower luminosities for these objects. As this also results in a lower stellar mass the L/M ratio, and thus also Γ_e, remain almost unaffected. However, the absolute mass-loss rates would be affected, and the O supergiants would systematically shift downwards along the ordinate. The effect on the slope of the O-type relation is rather small, as the O supergiants are evenly distributed along the abscissa.

For the O-type relation we obtain under the assumption of chemical homogeneity (Table 6). Using evolutionary masses from Köhler et al. (2014) for the O-type stars with high hydrogen surface abundances the relation changes to . The slope is only slightly shallower and well within the errors bars. The main difference is that the O-type relation is shifted towards higher Γ_e because of the systematically lower evolutionary masses.

Furthermore, clumping could affect our results qualitatively if the clumping properties change systematically, e.g. between O and WR-type stars. However, in Sect. 4.6 we show that extreme volume-filling factors of the order of 0.01 or less (as suggested for O stars e.g. by Bouret et al. 2003; Fullerton et al. 2006) are unlikely for our sample. Moderate volume-filling factors of the order of 0.1 may affect our results moderately if they change systematically, e.g., between different spectral types.

The main result of this work is that the observations qualitatively confirm the existence of a strong Γ-dependence as predicted by Gräfener & Hamann (2008) and Vink et al. (2011). Furthermore we find two branches with different slopes for low and high Γ_e which appear to be connected by a kink as predicted by Vink et al. (2011). The location of the kink in Fig. 8 is at log Γ_e = −0.58 (Γ_e = 0.26). For the reasons outlined in Sect. 4.2 we believe that this is a lower limit for the true value where the transition between normal O star mass-loss and enhanced Γ-dependent mass-loss happens. E.g. the models by Gräfener & Hamann (2008) suggest a higher value around Γ_e ~ 0.5 for LMC metallicity.

The question whether the observed change in the mass-loss properties also marks a change in the underlying wind physics is presently difficult to answer. The models by Gräfener & Hamann (2008) suggest that the wind physics for the Of/WN and WNh stars on the steep branch is different from that for OB stars. The kink found by Vink et al. (2011) also suggests that this is the case. In the present work we confirmed the steep Γ-dependence, but the nature of the transition remains unclear. In this context it is noteworthy that, in terms of mass-loss, the kink in Fig. 8 already appears in the regime of extreme Of stars, well below the mass-loss range of Of/WN stars. As we discuss in Sect. 4.6 the wind performance numbers η in this regime are most likely still below unity (if we assume that O star winds are clumped). The transition to a strong Γ-dependence thus seems to happen for lower wind efficiencies than the value of η = 1 predicted by Vink et al. (2011).

The observed slope with on the high-Γ branch is similar to the prediction by Vink (2006) who gives a relation with for VMS stars approaching the Eddington limit. Gräfener & Hamann (2008) and Vink et al. (2011) also predicted a steep increase but used more complex descriptions than a simple power law. Nevertheless, the relevant result is the qualitative confirmation of a relation between Ṁ and Γ_e.

Finally there is an offset between the majority of the Of stars and the Of/WN stars in Fig. 8. This could indicate systematically lower volume-filling factors (f_v = 0.1) for the Of/WN and WNh stars, compared to more moderate values (f_v = 0.25) for O-type stars (see also Sundqvist & Owocki 2013). Gräfener et al. (2011) did not find a comparable gap in the Arches cluster sample.

4.6. Wind efficiency and clumping

The wind efficiency parameter η = Ṁν_∞/ (L_⋆/c) gives the ratio between mechanical wind momentum Ṁν_∞ and radiative momentum L_⋆/c of a star. η = 1 denotes the single scattering limit, for which each photon is scattered on average once to drive the stellar wind, i.e. η ≈ τ_wind ≈ 1 where τ_wind denotes the flux-mean optical depth at the sonic point (cf. Lamers & Cassinelli 1999, Sect. 7.2). Vink & Gräfener (2012) pointed out that this relation is independent of wind clumping and used it to calibrate the mass-loss rates of VMS at the single scattering limit. Their model computations suggested that, more precisely, at the single scattering limit η = fτ_wind where f depends on the ratio ν_∞/ ν_esc and is typically of the order of f ≈ 0.6. We thus expect stellar winds to become optically thick for η ≈ 0.6. At this point the majority of photons in the flux maximum are absorbed within the wind, with the consequence that the wind recombines and forms emission-line spectra.

Our present stellar sample traces this transition observationally in unprecedented detail. In Fig. 10 we show the wind efficiencies that follow from our present analyses. The dependence of η on the adopted volume-filling factor arises from Eq. (1). The point where the winds become optically thick is indicated for the cases with (f_v = 0.1) and without (f_v = 1.0) clumping. Spectroscopically these points coincide with the domain of Of/WN stars (adopting f_v = 0.1) and the transition between O stars and Of stars (adopting f_v = 1.0), i.e., the regime where stars start to form emission lines. Clumping factors of this order of magnitude are thus plausible. Volume-filling factors of this order of magnitude are commonly found for WR stars (e.g. Hamann & Koesterke 1998) and have been proposed for O stars (e.g. Puls et al. 2006; Šurlan et al. 2013). Extreme volume-filling factors of the order of f_v = 0.01 as proposed e.g. by Bouret et al. (2003); Fullerton et al. (2006) appear very unlikely, as this would move the transition into the region of WNh stars where strong emission line spectra are already fully developed. We thus conclude that volume-filling factors at the transition between optically thin and optically thick winds are likely moderate, i.e., they are likely of the order of f_v = 0.1 or higher.

4.7. Stellar evolution in the upper HRD

The data presented in the previous sections are unique as they provide the largest stellar sample to date with luminosities above ~400 000 L_⊙. They also provide a unique opportunity to investigate the evolution of stars with masses above ~40 M_⊙, for which constraints are otherwise sparse. In the following we draw some straightforward conclusions from comparisons with stellar evolution predictions.

In addition to the effects of mass loss from the stellar wind, internal mixing processes (e.g., convective overshooting, rotationally-induced mixing) also constitute major uncertainties in evolutionary models. Furthermore, the stars considered here are very close to their Eddington limit, which may affect their envelope structure (Köhler et al. 2014) and give rise to instabilities (Gräfener et al. 2012), raising the question whether eruptive mass-loss, as observed e.g. for some LBVs and supernova progenitors, may also affect their evolution. Finally, recent work by Sana et al. (2012), Chini et al. (2012), and Sana et al. (2013) suggests that binary interaction greatly affects the evolution of massive stars, implying that several stars in our sample may have been subject to mass transfer or merger events (see de Mink et al. 2014; Schneider et al. 2014).

The observed HRD positions of our sample stars are indicated in Figs. 11 and 12, distinguishing between different spectral types and surface compositions. They are compared with single star evolutionary tracks from Köhler et al. (2014) with an initial rotational velocity of ν_R,ini = 300 km s^-1. In agreement with our conclusions from Sect. 4.5 the HRD positions of most of our sample stars agree with an evolutionary stage on the main-sequence. The only exception is VFTS 108 which is found to the left of the Zero-Age Main-Sequence (ZAMS), supporting our previous conclusion that this star is most likely in the core He-burning phase.

Considering the surface helium abundances of the main-sequence stars in our sample (Fig. 12), it appears convenient to divide the sample into two luminosity or mass ranges, with the dividing line roughly at 10^6.1L_⊙ or 80...90 M_⊙. The group of 24 stars below this luminosity limit consists of O-type stars of which only three have a somewhat enhanced helium abundance. All other stars still show their original helium abundance. The group of 19 stars above the luminosity limit is dominated by Of/WN and WNh stars. The majority (13 out of 19) of stars in this group have enhanced helium surface abundances. Furthermore, there appears to exist a similar limit for the enrichment in nitrogen. Stars without evidence for nitrogen-enrichment are only found in the low-luminosity group.

Based on their evolutionary models Köhler et al. (2014) identify similar luminosity limits as the ones described above. They find that the minimum luminosity for which the surfaces of main-sequence stars show considerable helium enrichment depends on the initial rotational velocity and varies between log  (L/L_⊙) ≈ 4.8 for fast rotators and log  (L/L_⊙) ≈ 6.3 for slow rotators. Based on the models by Köhler et al. the empirical boundary reported here (log  (L/L_⊙) ≈ 6.1) suggests that the bulk of our sample stars had initial rotational velocities that were slower than ~200 km s^-1.

Projected equatorial rotational velocities ν_esini of the presumed single O-type stars in the VFTS sample have been reported by Ramírez-Agudelo et al. (2013). After correcting for projection effects and macro-turbulence, they find that the rotational velocity ν_esini distribution strongly peaks at 100 km s^-1, and that about 75% of the sample has ν_e ≲ 200 km s^-1.

Based on the synthetic emission line profiles from our models we are further able to constrain ν_esini for 15 Of/WN and WNh stars in our sample. We use narrow emission lines (with an intrinsic width of 50 to 200 km s^-1) that are formed near the stellar surface, and should thus well reflect the rotational speed at the surface. In contrast to the values for O-type stars from Ramírez-Agudelo et al. these are not corrected for additional line broadening due to macro-turbulence. Because of this, and because the intrinsic line profiles are affected by the wind velocity, we are only able to give conservative upper limits for ν_esini of 200 km s^-1. The micro-turbulent velocities adopted in our models are much lower (10 km s^-1) and have thus no effect on the given limits.

The results are presented in Fig. 13 where we plot the observed helium surface mass fractions Y vs. ν_esini for these objects and for the 20 O-type stars from Ramírez-Agudelo et al. (2013) and Sabín-Sanjulián et al. (2014) overlapping with our sample. There is no evidence for a correlation between Y and ν_esini in this figure suggesting that rotationally induced mixing is not the dominant factor responsible for the helium surface enrichment.

In the whole sample there are only two fast rotators with ν_esini> 200 km s^-1. These are VFTS 626 with 278 km s^-1 and VFTS 755 with 286 km s^-1. While VFTS 626 belongs to our group of low-luminosity stars with He-enrichment (Y = 0.4), VFTS 755 does not show any sign of surface enrichment in He or N. The main difference between both stars is their mass-loss rate which is three times higher for VFTS 626.

In Fig. 14 we investigate the relation between surface helium enrichment and mass-loss. There appears to be a very clear correlation between Y and log  (Ṁ/M), indicating that mass loss plays a dominant role in the helium enrichment for these objects. The X-axis in Fig. 14 can also be read as the inverse of the mass-loss timescale τ_Ṁ = M/Ṁ, or more precisely an upper limit for τ_Ṁ as the mass-estimates used are upper limits (cf. Sect. 4.2, some lower limits for potential He-burners are also indicated in Fig. 14). For log  (Ṁ/M) ≳ − 6.5, τ_Ṁ becomes comparable to the main-sequence lifetime⁴ and the correlation becomes very clear.

Our results thus suggest that the observed changes in the surface abundances for luminosities ≳10^6.1L_⊙ are caused by the change in the mass-loss properties discussed in Sects. 4.4 and 4.5. The observed luminosity threshold coincides with the luminosity above which we see evidence for enhanced Γ-dependent mass-loss in our sample (cf. Fig. 7). In the high-luminosity group, with L ≳ 10^6.1L_⊙, we thus find predominantly Of/WN and WNh stars whose hydrogen-rich layers have been stripped off by the enhanced mass-loss during their main-sequence lifetime.

For the stars below the luminosity threshold our previous discussion suggests that mass loss and rotational mixing do not dominate the surface helium enrichment. In this case binary interaction may play a role, and the three He-enriched stars in this group may (predominantly) be binary interaction products. How does this number compare with predictions from binary evolution models by de Mink et al. (2013, 2014)? Assuming continuous star formation and an initial binary fraction of 70% de Mink et al. (2014) predict an incidence rate of 28% of binary interaction products, 50% pre-interaction binaries and 22% single stars. They suggest that these numbers can be compared to the population in 30 Dor where Sana et al. (2013) observed a spectroscopic binary fraction of 35%, and inferred that the intrinsic binary fraction is ≈50% which is comparable to the predicted fraction of pre-interaction binaries from de Mink et al.

The spectroscopic binary fraction in our initial sample is 34% (21 out of 62 stars), i.e., it is comparable to the results of Sana et al. (2013) for 30 Dor. As we removed 8 SB2 systems in our previous discussion we expect a slightly higher incidence rate of binary products than predicted by de Mink et al. (2014). Given the errors of roughly 50% from de Mink et al., and assuming that these objects would show an enhanced He abundance, we thus expect 3...10 stars of this kind in our low-luminosity sample. However, we have to be careful as we also excluded several stars with uncertain parameters from the initial sample. Inspection of our provisional parameters for these objects in Table 4 reveals that we excluded three objects that would belong in the low-luminosity group and show indications of an enhanced helium abundance (VFTS 208, 406, and 562). Including these objects would bring our observed numbers in agreement with the predictions. A thorough investigation of such “problematic” objects may thus be crucial to understand the role of binary interactions in massive star evolution. Furthermore, a higher accuracy in Y than achieved in this work (±5%; cf. Sect. 3.5) may be desirable.

Finally, we note the distribution of spectral types in the upper HR-diagram (Fig. 11). Again, the objects in this diagram can be divided in two groups, above and below 10^6.1L_⊙. The low-luminosity group consists of O-type stars with a clear separation between O stars and Of stars, indicating that Of stars are cooler and more evolved than O stars. In contrast to this the high-luminosity group, which consists of a mixture of O-type stars, Of/WN and WNh stars, shows a much more erratic spectral type distribution.

A comparison with the empirical HRD of the VFTS-MEDUSA region in Doran et al. (2013, their Fig. 12a) suggests that our high-luminosity group covers the complete range of observed effective temperatures in the HRD, while our low-luminosity group appears to be truncated at T_⋆ ≈ 40 kK, omitting cooler temperatures and thus likely the most evolved stars on the main sequence. This selection effect may also affect our previous discussion of the importance of binary interactions in the low-luminosity group.

As the stars in the high-luminosity group cover the majority of the main-sequence phase this group will include relatively young objects evolving towards the red, but also more evolved objects evolving towards the blue. Their observed temperatures are generally higher than ≈40 kK. The evolutionary tracks from Köhler et al. in Figs. 11 and 12 reach much lower temperatures than what is observed. This is the case for all of their models with ν_R,ini ≤ 300 km s^-1. One reason for this is the envelope inflation effect which has originally been predicted for massive stars near the Eddington limit by Ishii et al. (1999). Recent theoretical studies of this effect suggest that it can substantially affect stellar radii and temperatures in the upper HRD and lead to LBV-type instabilities, depending on a variety of physical effects such as the detailed structure and dynamics of the sub-photospheric layers, the outer boundary conditions imposed by a stellar wind, or stellar rotation (Petrovic et al. 2006; Gräfener et al. 2012; Gräfener & Vink 2013; Köhler et al. 2014). The observed HRD positions in Fig. 12 are likely affected by a combination of evolutionary effects and the complex envelope physics in this part of the HRD. Enhanced mass loss in the high-Γ regime, as it has been discussed in the present work, is expected to reduce the inflation effect (cf. Petrovic et al. 2006), and may thus help to resolve the apparent discrepancies.

4.8. Ionising fluxes

Recently Doran et al. (2013) highlighted the pivotal role of VMS in the ionising feedback from massive OB stars in 30 Doradus. In this section we compare the number of ionising photons per second (Q₀) from our work with the values from Doran et al. (2013). In contrast to our grid analysis Doran et al. (2013) used a template method. The ten stars with the highest ionising fluxes in 30 Doradus are listed in Table 10 of Doran et al. (2013). Five stars are in our sample, VFTS 758, 402, 1001, 1025, and Mk42. VFTS 758 (BAT99-122) and Mk42 agree within an uncertainty less than 50%. VFTS 1025 (R136c) has half of the number of ionising photons, but still remains in the top 10. The value for VFTS 1001 (R134) is a factor 4 lower. VFTS 402 (BAT99-095) has a factor 3 lower ionising flux. Both stars would drop out from the 10 stars with the highest ionising flux. However, VFTS 402 is a multiple system and our spectroscopic analysis is not reliable. On the other hand we determined higher values for VFTS 682 and VFTS 695 (BAT99-119) and they would move up into the top 10 with the highest ionising fluxes. By adding up all ionising photons of stars which are in both samples the uncertainties are averaged out and we obtained similar total fluxes for both methods, so that the overall results of Doran et al. (2013) are robust with respect to individual uncertainties.

5. Conclusion

In this work we analysed an unprecedented sample of 62 O, Of, Of/WN and WNh stars with very high luminosities. The sample contains a diversity of spectral types and luminosity classes and covers a range of mass-loss rates of ~2.5 dex. The analysis is based on a large grid of complex stellar atmosphere models computed with the non-LTE radiative transfer code CMFGEN (Hillier & Miller 1998).

We investigated the mass-loss properties of our sample stars based on the classical approach for O stars using the wind-momentum luminosity relation (WLR) from Kudritzki et al. (1999), and the semi-empirical approach by Gräfener et al. (2011) to investigate a dependence on the classical Eddington factor Γ_e. While the majority of O and Of stars in our sample follow the classical WLR, we find that the mass-loss rates of Of/WN, WNh, and some Of stars are enhanced with a steep dependence on Γ_e, in qualitative agreement with theoretical mass-loss predictions for very massive stars (VMS) by Gräfener & Hamann (2008) and Vink et al. (2011). Furthermore, our results suggest that at the transition between both regimes the Ṁ − Γ_e dependence shows a kink as predicted by Vink et al. (2011).

The main uncertainties in our results are due to wind clumping which is poorly constrained, in particular for O-type stars. Assuming that the winds are moderately and uniformly clumped (f_v ≈ 0.1) our results suggest that the mass-loss of VMS is enhanced with respect to standard mass-loss predictions from Vink et al. (2000, 2001). Furthermore, the dependence on Γ_e sets in earlier than expected, before the winds become optically thick. Its quantitative impact on the evolution of VMS may thus be larger than previously thought. Based on our present analysis the transition already occurs for luminous Of stars at Γ_e ≈ 0.3 with a steep mass-loss dependence above this value. We note that these results are based on simplified mass estimates and the true Eddington factors are likely higher than adopted. This may at least form part of the reason why, quantitatively, theoretical models predict the steep Γ-dependence for higher values of Γ_e ≈ 0.5 (Gräfener & Hamann 2008) and Γ_e ≈ 0.7 (Vink et al. 2011).

The observed HRD positions and mass-loss properties of our sample stars suggest that they are predominantly very massive main-sequence stars. Based on the observed stellar properties, the sample can be divided into two luminosity or mass ranges with a dividing line roughly at 10^6.1L_⊙ or 80...90 M_⊙. Above this limit we find predominantly Of/WN and WNh stars whose surfaces are enriched in helium due to their enhanced mass-loss. Below this limit only a few stars are He-enriched. These objects may be binary interaction products. Based on single star evolution models from Köhler et al. (2014) the location of the observed luminosity limit suggests that rotational mixing is relatively unimportant for the He-enrichment in our sample.

Finally, we confirm the important role of VMS in the ionising budget of 30 Dor found by Doran et al. (2013), providing evidence that these objects play a pivotal role in ionising and shaping the interstellar environment in young starburst-like regions.

---

Acknowledgments

We thank the anonymous referee for providing constructive comments, Joachim Puls and Jon Sundqvist for the fruitful discussion about mass-loss, and Danny Lennon and Selma de Mink for providing us with HST/F775W-band images (GO 12499) to check for visual companions. J.M.B., G.G., and J.S.V. thank the Science & Technology Facilities Council (grant No. ST/J001082/1) and the Department of Culture, Arts and Leisure in Northern Ireland for financial support. This work is based on observations with the European Southern Observatory Very Large Telescope, programmes 182.D-0222, and 084.D-0980. We thank John Hillier for CMFGEN and Wolf-Rainer Hamann for WRPLOT. The authors acknowledge the use of the Dell HPC Cluster and associated support services at Queen’s University of Belfast in the completion of this work.

References

Online material

Appendix A: Comments on individual objects

Here we describe the individual properties of the Of/WN and WNh stars in our sample (Sect. A.1). This is followed by a discussion of two peculiarities that occurred during the spectroscopic analysis (Sect. A.2).

Appendix A.1: Of/WN, WNh and WN(h) stars

We now give a short summary of individual stars and their characteristics.

VFTS 108 (BAT99-089/Brey71): a single WN7h-star with no noticeable variation in the peak intensity of the emission lines between the observed epochs. There are no absorption lines suitable for RV measurements to investigate the presence of a companion or the stars’s runaway status. The star is not located in a cluster and lies in the field surrounded by a few fainter objects. The quality of the model is reasonably good and the derived stellar parameters are well constrained (Fig. E.6).

VFTS 147 (BAT99-091/Brey73): the star has the spectral type WN6(h) and is not a known binary or multiple component system. According to Schnurr et al. (2008) it does not show any variability. However, the star is located in the Brey 73 cluster with a bright object nearby. In several epochs absorption features from a possible OB star are visible in the VFTS spectra. The presence/absence of these features is due to slightly different offsets of the MEDUSA fibre position and/or variations in the atmospheric seeing. The mass-loss rate and temperature estimation are uncertain. The low spatial resolution of the available optical photometry unrealistically high flux in the B and V bands. The continuum estimated from the near-IR VISTA photometry is less affected (or unaffected), which results in a negative extinction parameter E(V − K_s). The estimated luminosity is therefore based on the K_s magnitude. photometry (Sect. 3.4), but The luminosity might be up to ~0.2 dex too high depending on the uncertainties in the near-IR photometry and the near-IR flux of the nearby star. The star is excluded from the discussion (Fig. E.8).

VFTS 402 (BAT99-095/Brey80): this WN7(h) star has a probable OB companion. Schnurr et al. (2008) classified it as binary system with a short orbital period of ≈2 days. The model fit quality is very poor as a result of the presence of the secondary. The temperature, mass-loss rate and helium abundance determinations are highly uncertain (Fig. E.18). The impact of the secondary on the total continuum is not known, which leads to an incorrect luminosity. The star is excluded from the discussion.

VFTS 427 (BAT99-096/Brey81): a single star with no noticeable change in the peak intensity of the emission lines between the observed epochs. The WN8(h)-star is surrounded by a few nearby fainter objects, which do not considerably influence the spectra or the continuum. There are no photospheric absorption lines to study RV variations due to any possible secondary. The quality of the model is good and the derived stellar parameters are reasonably accurate (Fig. E.22).

VFTS 457 (BAT99-097): the star is classified as O3.5if*/WN7 and is not a known binary. Small LPVs are observed in the spectra which might be a result of the normalisation process. The fit to the nitrogen lines is poor, in particular for N iii. However, test calculations show that an increase of the nitrogen abundance, moving the starting point of clumping closer to the stellar surface to a wind velocity of 10 km s^-1, and a higher value of the wind parameter β improved the fit quality, but does not fundamentally change the result (see Figs. E.26 and E.27). These different properties might be a result of an extended or asymmetric atmosphere, as a result of rotation. The terminal velocity is between 900 and 1100 km s^-1 lower than the other Of/WN stars. Nevertheless, the derived stellar properties are robust.

VFTS 482 (BAT99-099/Brey78/Mk39): VFTS 482 is classified as O2.5if*/WN6. Massey et al. (2002) identified spectral variability while Massey et al. (2005) detected a composite in the spectra, but provided no classification of the secondary. According to Schnurr et al. (2008) it is a elliptical system with a period of 92.6 days (close to their detection limit). The quality of the fit is good. The temperature, mass-loss rate and He-abundance can be well determined (Fig E.29). The star is located in a relatively crowded field. The luminosity might be over estimated because of its unclear binary nature.

VFTS 509 (BAT99-103/Brey87): this WN5(h) star is part of a binary system with a separated early O star companion (Evans et al. 2011; Doran et al. 2013). Schnurr et al. (2008) give a period of 2.76 days. The fit quality is poor, but gives a fair hint of the actual stellar parameters (Fig. E.31). Still, the star is excluded from the discussion.

VFTS 542 (BAT99-113): this O2if*/WN5 star is a binary system with a period of 4.7 days (Schnurr et al. 2008). It is located in a crowded field. The quality of the model fit is good. The temperature, mass-loss and helium abundance can be determined (Fig. E.37). The luminosity is uncertain, because it is unclear how the continuum is affected by the secondary.

VFTS 545 (BAT99-114): as reported by Hénault-Brunet et al. (2012) there are weak RV variations in the spectra of this O2if*/WN5 star. There is no contamination of the stellar spectra; the fit is good and the results are robust (Fig. E.38).

VFTS 617 (BAT99-117/Brey88): this single star is classified as a WN5ha. The spectra show LPVs particularly in the Balmer series. The quality of the model is good. Ṁ varies as a result of the LPVs, but only within the error bars. (Fig. E.44).

VFTS 682: the WN5h-star is a newly discovered WR-star by the VFTS (Evans et al. 2011; Bestenlehner et al. 2011). The star appears to be single. The model fit is good and the results based on the grid are comparable with the parameters derived from Bestenlehner et al. (2011). The grid fit is shown in Fig. E.50.

VFTS 695 (BAT99-119/Brey90): this is a WR-star with spectral type of WN6h plus a companion. The binary system has a period of 158.8d (Schnurr et al. 2008). The spectra show LPVs in the Balmer series, in particular for H_α. The observations can be reasonably well reproduced by our model, however L_⋆ might be over estimated (Fig. E.51).

VFTS 758 (BAT99-122/Brey91): this is a single WN5h-star. No noticeable line-profile and RV variations were found in the spectra. The fit quality regarding the grid resolution is good and the results are reasonably good (Fig. E.53).

VFTS 1001 (BAT99-100/Brey75): the WN6(h)-star is not a known binary. Hénault-Brunet et al. (2012) found possible LPVs in the He ii λ4542 line. The star is associated with an X-ray source. Portegies Zwart et al. (2002) suggested that these X-rays are potentially due to a wind-wind collision in a binary system. However, our model fits the observation well and we treat the object as if it is a single star (Fig. E.55).

VFTS 1017 (BAT99-104/Brey76): the O2 if*/WN5 spectra of this star show LPVs and weak RV variations (Hénault-Brunet et al. 2012). The possible companion does not impact the fit quality. The results are reliable with slightly larger uncertainties than most well fit stars in this study (Fig. E.57).

VFTS 1022: this star has been reclassified from O4 If+ (Massey & Hunter 1998) to O3.5 if*/WN7 by Crowther & Walborn (2011). The star shows small RV variations. Even though it is likely a binary, the derived stellar parameter are treated as if it were a single star (Fig. E.60).

VFTS 1025 (BAT99-112/Brey82/R136c): this WN5h-star is located in a crowded field and surrounded by multiple sources in HST observations by Massey & Hunter (1998). In their catalogue VFTS 1025 is designated as source No. 10 with M_V = −6.8. One neighbouring HST source at an angular distance ≲0.5 arcsec is bright enough to affect our optical spectrum directly (source No. 57 with M_V = −5.5, classified as O3 III(f*)). Furthermore the bright O star R136b (source No. 9, M_V = −6.9, O4 If+) is relatively close at an angular distance of ≈1.3 arcsec. VFTS 1025 itself shows indications of binarity based on its strong X-ray emission (Guerrero & Chu 2008) and possible low-amplitude RV variations (Schnurr et al. 2009).

In our analysis we determine a significantly lower temperature than Crowther et al. (2010), mainly based on the He iλ4471 Å absorption line in the optical range. Our cooler temperature is further supported by the low terminal wind speed of VFTS 1025 compared to the other WN5(h) stars in R 136 analysed by Crowther et al. Their hotter temperature is based on the N v λ2.10 μm emission line in the IR. It is not clear if the He i line in our optical spectrum is intrinsic to VFTS 1025 or originates from a nearby star. We can exclude star No. 57 from Massey & Hunter as the origin of the He i line, based on the non-detection of strong enough He i in its HST spectrum. R136b has a weak He i absorption line and could possibly affect our spectrum. Hénault-Brunet et al. (2012) found LPVs for the prominent H and He ii emission lines but not for the He i absorption line in the spectrum of VFTS 1025. This is in line with a photospheric nature of the absorption line and wind variability in the emission lines as expected from our cool model, but also with a contribution from a nearby star.

Appendix A.2: Peculiarities

In this section we discuss three peculiarities that arose during our spectroscopic analysis. These peculiarities occurred as discrepancies between the observations and the model predictions. The plots of the model fits are shown in Appendix E.

Firstly, our models predict the C iii lines at λ4647/4650 Å in emission. Martins & Hillier (2012) discussed that these lines may appear in absorption or emission, depending on the detailed model assumptions. In our observations they appear in absorption for VFTS 171, 333, 532, and 669 and in emission for VFTS 216, 404, 422, 603, 608, and 797. Furthermore, VFTS 532 has a cooler companion that could contribute to the C iii absorption. In our analysis these lines are ignored.

Secondly, VFTS 208, 406, and 626 show He ii λ4686 Å emission profiles with a central absorption, which are characteristic for the fast rotating stars in the Onfp sub-class (see Walborn et al. 2010). This type of line profile has been modelled in 2D by Hillier et al. (2012). VFTS 406 has been classified as OVnn by Walborn et al. (2014) who find that its composite spectrum in the VFTS data is likely due to contamination by an adjacent WN spectrum on the detector. For VFTS 208 and 406 we only achieve a bad fit quality and the stars are excluded from our discussion.

Thirdly, for VFTS 259, 457, and 1021 the observed N iii λ4634/4640 Å lines are significantly stronger than in our theoretical models. Furthermore, there is a discrepancy between the diagnostic He i λ4471 Å line, and the N iii λ4634/4640 Å and N iv λ4058 Å lines for the effective temperature. To estimate the uncertainties we carried out test calculations for the O star 259 and Of/WN star 457. We found that the discrepancies can be significantly improved by simultaneously lowering log g by 0.5 dex, increasing the N-abundance by a factor 2, and moving the starting point of clumping closer to the surface (see Figs. E.26 and E.27 for the case of VFTS 457). The resulting parameters agree with those given in Table 3 within the given uncertainties.

Appendix B: Additional tables

Appendix C: Temperature sensitivity and He-abundance

Appendix D: The influence of log g on the stellar parameters

A variation of log g changes the density and ionisation balance in the stellar atmosphere. This can directly affect the derived temperature and helium abundance and with it the luminosity and the mass-loss rate. Here, we investigate the effect of varying surface gravity on the result for the O-type stars. For this purpose, we used our log g sub-grid with some additional smaller steps in helium abundance, and analysed a small sub-sample of O stars in that parameter space. More accurate parameters using a different stellar atmosphere code including a variable log g in the analysis will be provided by Ramirez-Agudelo et al. (in prep.) and Sabín-Sanjulián et al. (in prep.).

As shown in Figs. D.1 and D.2 a lower value of the surface gravity (Δlog g = 0.5 dex) requires a lower model temperature (ΔT_eff = −0.025... − 0.03 dex) to match the main O star temperature diagnostic Hei  λ4471. The transformed mass-loss rate is almost unchanged, within the given error bars. The lower T_eff results in a lower luminosity (~ − 0.1 dex) and lower stellar mass. This gives a lower terminal velocity for the O stars using the relation by Lamers et al. (1995). These changes affect the mass-loss rates derived from Eq. (1) by about ~0.15 dex. The derived He-abundance decreases by 2...4%. The differences are larger for increasing temperature.

A variation of log g to values as low as 3.5 is only possible for the lowest T_eff in our model grid. The reason is that for the hotter models the radiative flux () in the photosphere becomes so high that log g_rad exceeds 3.5, leading to a situation where Γ > 1 (cf. Eq. (3)) and no hydrostatic solution is possible. As Γ is generally high for the hotter models, and the spectroscopically relevant quantity is g_eff = g(1 − Γ), the log g needed to match the H_{β − δ} line profiles of hot supergiants are not significantly lower than the log g of 4.0 adopted in our grid models. E.g., Evans et al. (2010) determined log g = 3.75 for VFTS 016 (O2 III-If*) in agreement with our own test computations.

Appendix E: Spectral modelling for each target star

In the appendix we show model fits for all our targets and describe how we obtained the stellar parameters. The first panel in the figures is the model SED (red solid line) fit to the optical and near-IR photometry (blue boxes). The reddening law and the R_V is given in the bottom left corner. The following panels show the fits of the model spectra (red solid line) to the observations (blue solid line).

All Tables

All Figures