© ESO 2022

1 Introduction

The star cluster R136 inside the Large Magellanic Cloud (LMC) hosts some of the richest populations of high-mass stars in the local Universe. Nine stars within this cluster have masses around or exceeding 100 M_⊙ and a few even surpass ~150 M_⊙ (Crowther et al. 2010, 2016; Bestenlehner et al. 2020). These massive, very luminous, hot stars play an important role in the Universe. They strongly influence their surroundings through the direct injection of mass, momentum, and energy from winds (Weaver et al. 1977), radiation pressure (Mathews 1967), and thermal expansion caused by photoionisation from extreme ultraviolet (UV) photons (Kahn 1954; Spitzer 1978). This feedback can act to both disperse star-forming clouds (e.g. Dale et al. 2014; Kim et al. 2018) and cause compressive flows that lead to the formation of new stars (e.g. Inutsuka et al. 2015; Rahner et al. 2018; Fujii et al. 2021). The mass-loss rates and terminal velocities of winds have a direct impact on the ability for stars to regulate their environment. The deposition rates of stellar wind energy and ionising photons from massive stars are important quantities in driving the multi-phase structure of the interstellar medium and the regulation of future star formation. Furthermore, massive stars end their lives in supernovae, thereby enriching the interstellar medium with newly formed chemical elements and leaving behind compact remnants such as black holes (see, e.g. Smartt 2009; Langer 2012, for a review). Moreover, with their high masses and at half solar metallicity, the stars in R136 are close observable counterparts to the first stars, which are estimated to have a characteristic mass of tens to hundreds of solar masses (see, e.g. Hirano et al. 2014, 2015; Sugimura et al. 2020; Chon et al. 2021; Park et al. 2021, and references therein).

R136 resides inside the Tarantula Nebula or 30 Doradus. This nearby, unobscured, large and intrinsically very bright star forming region (Kennicutt 1984; Doran et al. 2013), hosting many hundreds of massive stars (M ≳ 8 M_⊙), resembles giant starbursts observed in distant galaxies (Cardamone et al. 2009; Crowther et al. 2017). The massive star content of the Tarantula Nebula was studied in great detail in the VLT Flames Tarantula Survey (VFTS, e.g. Evans et al. 2011, Evans et al. 2020; Bestenlehner et al. 2014; Schneider et al. 2018b); however, due to severe crowding, the central core of the R136 cluster was not part of the observing campaign. With the spatial resolution of the Hubble Space Telescope (HST), individual stars in the R136 core can be resolved. This was employed by Crowther et al. (2016), who used HST to collect optical and UV spectroscopy of the cluster, hereby complementing the VFTS survey and extending the coverage to the most massive stars. Focussing on the UV spec-troscopy of this dataset, Crowther et al. (2016) derived a cluster age of $1.5{\pm }_{0.7}^{0.3}$ Myr and found that the He ii λ1640 emission is completely dominated by stars with initial masses ≳100 M_⊙. Bestenlehner et al. (2020) focussed on the optical spectra and derived detailed spectral parameters for all sources. Their findings include a top-heavy initial mass function (IMF) for massive stars in the cluster, and a strong helium enrichment for the most luminous stars.

In order to understand massive star evolution, it is key to know the mass-loss rates of these very massive stars. Moreover, by calibrating theoretical models with observed mass-loss rates and stellar properties, we can improve future large-scale studies of stellar feedback, and hence obtain a more complete picture of how stars shape our Universe. For very massive stars, the effects of mass loss become especially important, as the rates generally increase with luminosity and thus with mass (e.g. Puls et al. 2008; Vink et al. 2011; Vink 2015). Moreover, their large convective cores ensure that they evolve close to homogeneously, diminishing the relative effect of other processes such as rotational mixing and magnetic fields (Yusof et al. 2013; Köhler et al. 2015; Ramachandran et al. 2019). Unfortunately, due to a lack of empirical constraints and the proximity to the Eddington limit, the mass-loss rates in this regime are uncertain (e.g. Langer 2012). Furthermore, obtaining accurate, empirical mass-loss rates is hampered by the presence of small-scale inho-mogeneities in the wind, also called ‘clumps’ (see e.g. Puls et al. 2008, for a review). The origin of this wind structure, or so-called wind clumping, is theoretically attributed to the line-deshadowing instability (LDI), an inherent property of the line force that drives the winds of these stars (e.g. Owocki & Rybicki 1984, 1985, and references therein).

Since wind clumping determines how our diagnostics respond to mass-loss rates, it is imperative to take it into account properly when studying massive star winds. The simplest approach to account for wind clumping in diagnostic models is to assume that the outflowing gas is concentrated in clumps that are small and rarefied enough so that they stay optically thin, and that the interclump medium is void (e.g. Hamann & Koesterke 1998; Hillier & Miller 1999; Puls et al. 2006). For O-stars, this ‘optically thin clumping’ or ‘micro-clumping’ approach leads to a downward revision of empirical mass-loss rates compared to the assumption of a smooth wind for processes that are dependent on the square of the wind density, such as the formation of the Hα line, but it can also affect lines indirectly due to changes in the ionisation and excitation equilibrium (see e.g. Puls et al. 2008, and references therein). If clumps become optically thick (for the considered process), the clumping affects diagnostics in a more complicated way. In this case, light can be blocked by clumps, but can also leak through porous channels in the wind and in this way escape without interacting (Shaviv 1998, Shaviv 2000; Owocki et al. 2004). These velocity-porosity or ‘vorosity’ effects mostly impact resonance lines; neglecting these phenomena can lead to an underestimation of mass-loss rates (e.g. Fullerton et al. 2006; Oskinova et al. 2007; Sundqvist et al. 2011; Šurlan et al. 2013). In order to obtain reliable mass-loss rate measurements, it is thus essential to consider all the aforementioned effects. To date, only one sample of O4–O7.5 supergiants has been studied using an optically thick clumping description in a model atmosphere code (Hawcroft et al. 2021).

In this paper, we reanalyse the R136 sample of Crowther et al. (2016) and Bestenlehner et al. (2020), but now combine the optical and UV spectroscopy, allowing us to study the mass-loss rates and wind structure in detail. For the wind structure, we assume the two-component formalism of Sundqvist et al. (2014) implemented in Fastwind (Sundqvist & Puls 2018), allowing for optically thick clumps and thus including the effects of porosity, velocity-porosity, and a non-void interclump medium. This yields the most accurate mass-loss rate measurements possible with current model atmosphere codes, and furthermore will allow us for the first time to investigate the wind structure for a wide range of stellar properties.

The remainder of this paper is structured as follows. We start by presenting the R136 sample and our dataset in Sect. 2. In Sect. 3 we lay out our methodology. Here, we introduce our fitting algorithm Kiwi-GA and describe the model atmosphere code Fastwind. In particular, we emphasise the parameterisa-tion of the wind structure parameters (Sect. 3.1.1). The results of our analysis are presented in Sect. 4. This section is concluded with several tests of robustness (Sect. 4.7). We discuss our results in the context of theoretical predictions and evolutionary models in Sect. 5; Sect. 5.1 is dedicated to mass-loss rates and wind momentum; Sect. 5.2 pertains to the potential trends that we observed in the wind structure parameters; and in Sect. 5.3 we consider our findings in the context of stellar evolution. Two methods for measuring terminal velocities are compared in Sect. 5.4. We conclude with a summary and outlook (Sect. 6).

2 Sample and Data

Our sample consists of 56 stars residing in the core of the R136 cluster. Spectral types range from late to early O-type, plus three hydrogen-rich Wolf-Rayet (WNh) stars of subtype WN5h. Of the O-type stars, four are supergiants, five are giants, and the rest are dwarfs (Bestenlehner et al. 2020, Caballero-Nieves, in prep.). Figure 1 shows their positions with their Hunter et al. (1995) or Weigelt & Baier (1985) identification, projected onto an HST/WFC3 image (O'Connell 2010). The figure shows that the area is very crowded; the high spatial resolution of HST is thus a necessity to resolve individual stars in the core of the cluster. Recent advances in adaptive optics have made it possible to obtain such high-resolution observations of the core of R136 also with ground based instruments, as is done by Khorrami et al. (2017, 2021, imaging with SPHERE) and Castro et al. (2021, optical spectroscopy with MUSE-NFM).

Previous spectroscopic analyses of several sample stars have been carried out. A sub sample consisting of bright members of the cluster core has been studied by de Koter et al. (1997, 1998), who measured stellar parameters as well as mass-loss rates of 14 sources using HST/GHRS/FOS optical and UV spectroscopy (overlap with our sample is indicated in Fig. 1). Massey & Hunter (1998) analyse optical HST/GHRS/FOS spectra, focussing on stars in the outskirts and surroundings of the cluster. Schnurr et al. (2009) obtained time-resolved NIR spectra of five stars in the core, searching for binarity and reporting a dearth of short period binaries in their sample. We assume in this study that the sources we observe are either single or that the light is dominated by the brightest component; however, the multiplicity properties of the sample remain an open question and require further investigation (Shenar et al. 2019). Combining aforementioned UV, optical and NIR spectroscopy, Crowther et al. (2010) re-derived physical properties of the WNh stars, finding a present day mass of 265 M_⊙ for the most massive star, R136a1.

A comprehensive view of the cluster core was first given by Crowther et al. (2016), who secured optical and UV HST/STIS spectroscopy of the central cluster and obtained temperatures, wind velocities and spectral types from the UV spectra. The 55 optical spectra of this dataset¹ were analysed by Bestenlehner et al. (2020), who determined detailed stellar parameters for all stars and found that at least seven stars have current masses of 100 M_⊙ or more, reporting 215 M_⊙ for R136a1.

2.1 HST/STIS data

For our analysis we use blue-optical, Hα, and far-UV HST-STIS spectroscopy (PI: Crowther, Crowther et al. 2016). For this, 17 HST-STIS long-slit (52′′ × 0.2′′) contiguous pointings were done for six different gratings; technical details are summarised in Table 1. The setup is depicted in Fig. 1. The orientation of the slits (at position angle 64° or 244°) was chosen to align with R136a1 and R136a2, that lie only 0.1′′ apart. The image shows that crowding can occur elsewhere, which can cause contamination of spectra of stars that lie close together. While the spectral extraction process was designed to avoid this, several spectra might still be affected (see Table I.1). The contamination is severe only in the case of R136a6. This source can be resolved into two sources, H19 and H26, having a flux ratio of 0.78 in the V-band and a separation of only 70 mas² (Hunter et al. 1995; Khorrami et al. 2017). The brighter component was likely located partially out of the slit, and we therefore expect that H19 and H26 have an approximately equal contribution to the flux of what we call ‘R136a6’. We have no way of separating these components and therefore analyse R136a6 as if it were single; however, we exclude it from the analysis of the sample as a whole regarding mass-loss and clumping properties.

The spectra were extracted with multispec, a package tailored to extracting spectra from crowded regions (Maiz-Apellaniz 2005, 2007). Exceptions are the UV spectra of H70 and H141, where the extraction was done with calstis (Bostroem & Proffitt 2011). Since more than two sources are in each pointing the sources are not necessarily centred in the slit, which causes an uncertainty in the absolute wavelength scale that can be as large as ±2 pixels (see also Sect. 2.3). Hα suffers from strong nebular emission, for which we correct by interpolating the Hα emission of off-source spectra and subtracting this from the source spectra. The signal-to-noise ratio (S/N) of the resulting spectra is in the range 7–70, with average values of 23, 19, and 19 for UV, blue-optical and Hα. Figure B.1 shows a distribution of the S/N of the sample stars per wavelength range. A more comprehensive description of the UV data reduction can be found in Crowther et al. (2016), and the optical reduction will be described in more detail in Caballero-Nieves (in prep.).

2.2 GHRS data

For ten sources, we complement the HST/STIS spectra with archival HST/GHRS UV spectra (PIs: Heap & Ebbetts, de Koter et al. 1997, 1998). The wavelength range of these spectra spans from 1150 to 1750 Å, which means that they include the full N iv λ1718 line, contrary to our UV HST/STIS spectra where this line is positioned right on the edge of the grating. We do not use the spectra of R136a1 and a3 because they are possibly contaminated in view of the size of the aperture (0.22′′)³. Furthermore, the sample of de Koter et al. (1998) includes H39 which is not covered by our slits (see Fig. 1) and H42, which is an SB2. We include the N iv λ1718 in our fitting for the ten remaining spectra of their sample, while using HST/STIS data for all other lines. Moreover, the resolving power (λ/Δλ) of these spectra is approximately 5000, which allows us to resolve the interstellar C iv λλ1548–1551 lines and use this for the correction of the HST/STIS data (see Sect. 2.5). For more details of the HST/GHRS observations and data reduction we refer to de Koter et al. (1997, 1998).

2.3 Optical data preparation

Our spectral fitting code needs a set of normalised spectral lines. To this end we have normalised the spectra of the optical and Hα gratings locally around the diagnostic lines, assuming that the continuum can be approximated as a straight line. For each line we obtain the S/N from the continuum selected for the normalisation, and define the uncertainties on the normalised flux as the inverse of the S/N. In some cases we exclude data points, this includes isolated points that deviate a factor 2 or more from the value of the surrounding points, or data points in the centre of Hα where nebular subtraction may have been imperfect⁴.

The radial velocity shift is determined by fitting to the spectral lines in each grating a set of Gaussian functions with centres corresponding to the rest wavelengths of the lines considered. For lines that are affected strongly by the stellar wind (Hα and He ii λ4686) we use the synthetic lines of a small grid of Fastwind models to determine the radial velocity shift.

As described in Sect. 2.1, the absolute wavelength scale of all observations can deviate up to 2 pixels. We correct for this by assuming that the wavelength deviation behaves similar to a Doppler shift. In practice, this means that we correct the spectra for a radial velocity without considering the aforementioned wavelength deviation, that is, we measure the ‘radial velocity’, but this value includes both the true radial velocity, as well as an adjustment for the absolute wavelength deviation. The latter adjustment, typically of the order of 0–100 km s⁻¹ (or 2 pixels), is not physical, but we simply do not have a better model to describe the offset. This approach is thus a pragmatic one, merely to correct the wavelength scale for several effects in order to bring the diagnostic lines in the rest frame of the synthetic spectra. An overview of the derived velocities used for the correction can be found in Appendix B.2.

2.4 UV data preparation

The continuum of hot star UV spectra is hidden by a forest of lines, most notably lines of highly ionised iron group elements such as Fe iv-v. Locating the continuum is thus not trivial, especially since the depth of the so-called pseudo-continuum, formed by the iron lines, depends on stellar properties, in particular on the effective temperature T_eff. To assist the normalisation, we use a grid of cmfgen models in which the iron lines are modelled (Hillier & Miller 1998; Bestenlehner et al. 2014). With the normalised synthetic spectra that these models provide we can recover the shape of the true continuum. Our approach is as follows: Firstly, we select a cmfgen model from the grid that matches best the stellar and wind parameters from Bestenlehner et al. (2020). We may change our choice of model later, see below. Secondly, we mask the wind lines and interstellar lines, so that only the pseudo-continuum is left. Thirdly, we divide the observed UV flux by the normalised flux of the model and fit this quantity with a polynomial, getting a so-called normalisation model. Finally, we divide the observed UV flux by the normalisation model in order to obtain the normalised flux.

Getting a reliable normalisation hinges on the first step. It is especially important that the T_eff of the model matches that of the observed spectrum. In order to assure this, we treat T_eff as a free parameter. We repeat the above steps for all T_eff in the cmfgen grid (ranging from 35–56 kK in nine steps, LMC metallicity; see Bestenlehner et al. 2014) and assess for which temperature the iron pseudo-continuum has the best fit. We also vary the radial velocity v_rad of each model (in steps of 25 km s⁻¹, which is about a tenth of the resolution element) and assess which value fits best. Just as for the optical, this v_rad value includes a possible correction for the wavelength calibration (see also Sect. 2.3). For the micro-turbulent velocity v_micro we assume 10 km s⁻¹ for all models, as this is the only value included in the grid of Bestenlehner et al. (2014) and because the exact value of v_micro has little influence on this specific exercise. An example of the UV normalisation process is shown in Fig. 2.

The final output is a spectrum normalised with the best fitting T_eff, and corrected for v_rad. Before accepting a fit we check it visually; extra care is taken for sources where the fit value of T_eff falls outside the uncertainty margins of the T_eff derived by Bestenlehner et al. (2020). For six sources the model with T_eff as derived from fitting the iron pseudo-continuum did not result in a good fit and in these cases we adopted the value of Bestenlehner et al. (2020) for T_eff (see Table The GHRS data are normalised in a similar way, but instead of fitting T_eff we assume the value found from the HST-STIS iron forest fit.

A by-product of the normalisation process is a measurement of T_eff from the iron lines alone; this measure is independent from the H, He, C, N, O diagnostic lines used for the rest of this work and the analysis of Bestenlehner et al. (2020). These values can be found in Table I.1 and are compared to the H, He, C, N, O temperature measurements in Appendix D.1.

We obtain the S/N of the HST-STIS UV spectra by using the HST-STIS exposure time calculator⁵, assuming the F555W magnitude from De Marchi et al. (2011), and exposure times and A_F555W from Crowther et al. (2016). Using the S/N we get from the calculator, we estimate the uncertainty on the flux points we use for the fit. For the GHRS data we use the provided error spectra.

2.5 Corrections for interstellar absorption lines in the UV

Three interstellar absorption lines blend with important diagnostics: H i at 1215.67 Å (Ly-α), the Si iv doublet at 1393.76–1402.77 Å, and C iv doublet at 1548.20–1550.77 Å. We correct for Ly-α and C iv λλ1548–1551 by recovering the H i and C iv column density in the line of sight (N_H and N_Civ, respectively), by fitting the interstellar profiles with a Voigt-Hjerting function⁶ (Tepper-García 2006, 2007). For the damping factors of the Lorentzian component of the profiles we use the radiative damping constants of the transitions. We fit the interstellar components of multiple spectra and correct the spectra with averaged values rather than the values of the individual fits, since we expect the uncertainties in this fitting process to dominate over the difference in column density from star to star.

In the case of Ly-α, the Voigt-Hjerting profile is fitted to a subset of the points of the normalised UV flux. We do not fit those parts of the spectrum where the Ly-α profile might be blended with the N v 1240 doublet. To estimate where this line ends, we use for each source the edge velocity v_edge from Crowther et al. (2016). Furthermore, for fitting the wings, we select points that trace the stellar continuum: parts of the spectrum that seem free of absorption lines. Figure 3 shows an example of a fit of the Ly-α profile. For finding the average of N_H we fit the Ly-α profiles of the 29 stars brighter than M_V = −5.50 (values from Bestenlehner et al. 2020) and obtain a value of log(N_H [cm⁻²]) = 21.88 ± 0.07, in good agreement with other N_H measurements towards 30 Doradus (e.g. de Boer et al. 1980, who find $\log \left(N_{\text{H}}\left[{\text{cm}}^{-2}\right]\right)=21.85{\pm }_{0.15}^{0.10}$.

For C iv λλ1548–1551, we use the higher resolution GHRS spectra of de Koter et al. (1998). Before fitting the interstellar profile we fit a polynomial through the stellar P Cygni profile of each star and subtract it from the spectrum, so that the interstellar component remains. We resolve two interstellar components with a different velocity in each part of the doublet, so we fit two double Voight-Hjerting profiles, where the ratio of the strength of each of the doublet components and the distance between them is set by the oscillator strengths and the rest wavelength difference, respectively. We fit all spectra of de Koter et al. (1998), except for R136b, where we had problems correcting for the stellar line. From this we find column densities of log(N_Civ[cm⁻²]) = 14.72 ± 0.10 and 14.21 ± 0.07. The individual and mean fits are shown in Fig. 4.

The derived average profiles are used for computing the interstellar line optical depth (as a function of wavelength), which we then subtract from the observed optical depth of each star, obtaining the corrected optical depth, which we convert back to normalised flux. In the case of C iv λλ1548–1551, the average profile is convolved with an instrumental profile corresponding to R = 1250 before it is used for the interstellar corrections of the HST/STIS data. The uncertainty margins on the column densities are used to estimate an uncertainty on the corrected flux, in addition to photon noise.

For all sample stars but one, the Si iv λλ1394–1402 stellar lines are in absorption and can, even in the higher resolution data of de Koter et al. (1998) not be distinguished from the interstellar components. Only in R136b the interstellar component is resolved, however here we are not able to accurately correct for the stellar profile. We therefore cannot correct for the interstellar Si iv λλ394–1402 and do not use this line. The exception is R136b where the line is strongly in emission, and we can clip the interstellar part.

2.6 VLT/SPHERE K_s photometry

In our fitting procedure we calibrate the luminosity with an observed stellar flux (see Sect. 3.3). For this we use the absolute, dereddened K_s-band magnitudes as presented in Bestenlehner et al. (2020). They use VLT/SPHERE K_s magnitudes from Khorrami et al. (2017) and in addition B or U and V magnitudes for the extinction correction (Hunter et al. 1995; De Marchi et al. 2011) and an LMC distance modulus of 18.48 mag (Pietrzyński et al. 2019). The K_s-band is the optimal choice for a luminosity anchor because at these wavelengths (2.2 μm) the extinction is low, while thermal radiation of dust is not yet an issue.

3 Methods

For the analysis we use the model atmosphere code Fastwind (version: V10.3.1) to compute synthetic spectra and the genetic algorithm Kiwi-GA for the fitting. In this section we introduce both tools and also describe our fitting setup and related assumptions.

3.1 Fastwind

Fastwind is a model atmosphere code tailored to hot stars with winds (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2012; Carneiro et al. 2016; Sundqvist & Puls 2018). It solves the NLTE number-density rate equations and takes into account the effects of line-blocking and line-blanketing⁷. The atmosphere consists of a spherically extended photosphere in (pseudo-)hydrostatic equilibrium that is connected to an expanding stellar wind at a velocity transition point v₀ near the sonic point. The stellar wind is parameterised by a mass-loss rate $\dot{M}$, a terminal velocity v_∞, and a wind acceleration parameter β. The wind velocity v_r as a function of radius r is expressed by the classic β-velocity law: $${\upsilon }_{\text{r}}(r)={\upsilon }_{\infty }{(1-b/r)}^{\beta },$$(1)

where b is a radius close to the stellar radius⁸ R_⋆, the exact value of b depending on v₀ (see Santolaya-Rey et al. 1997). Under these assumptions the structure and ionisation and excitation state of the atmosphere and the wind are computed, resulting in a so-called atmosphere model. Using this model, the individual spectral lines are synthesised.

Fastwind stands out in terms of speed, as one model is computed in approximately 15–45 minutes on a single modern CPU. Such a short computation time allows one to compute many models, and thus explore the parameter space thoroughly (see Sect. 3.2). Precision and speed are achieved simultaneously by splitting up the atomic elements in ‘explicit’ and ‘background’ elements. The explicit elements are computed in the co-moving frame using detailed atomic models for the spectral lines, while the background elements are only computed in an approximate way. Individual transitions of the latter are not syn-thesised, but their radiation field is taken into account, which is essential for the treatment of effects of line-blocking and blanketing. To speed up the computation, Fastwind calculates a representative mean radiation background instead of a detailed field (for details, see Puls et al. 2005).

For this work we use Fastwind version V10.3.1 (Sundqvist & Puls 2018), including explicit elements H, He, C, N, O, Si, P. This version is suitable for the analysis of stars with winds that are moderately optically thin in the optical continuum. This condition is met for all sample stars, including the three WNh stars. While the WNh stars do have the densest winds of our sample stars, they are not as dense as those of classical Wolf-Rayet stars. Indeed, when we run a Fastwind model with WNh parameters, we find that the electron scattering optical depth at the sonic point is well below unity at τ_e = 0.28 (for a typical O-star, we find τ_e ≈ 0.02).

3.1.1 Wind clumping, porosity, and vorosity

Clumping is implemented in Fastwind V10.3.1 as detailed in Sundqvist & Puls (2018), employing the two-component formalism introduced in Sundqvist et al. (2014). In this prescription, the clumped wind consists of overdense clumps with a density ρ_cl and an interclump medium with a lower density ρ_ic. The clumps occupy a certain fraction of the total wind volume f_vol, referred to as the volume filling factor. The clumping factor f_cl of the medium can be written as: $$f_{\text{cl}}\equiv \frac{\langle {\rho }^2\rangle }{{\langle \rho \rangle }^2}=\frac{f_{\text{vol}}{\rho }_{\text{cl}}^2+(1-f_{\text{vol}}){\rho }_{\text{ic}}^2}{(f_{\text{vol}}{\rho }_{\text{cl}})+(1-f_{\text{vol}}){\rho }_{\text{ic}}{)}^2}.$$(2)

We note that for the conventional assumption (not adopted here) of a void interclump medium (ρ_ic = 0) this would lead to f_cl = l/f_vol• The parameters describing the clumped medium are illustrated in the left panel of Fig. 5.

The formalism of Sundqvist et al. (2014) allows for the possibility of the clumps becoming optically thick. When clumps become optically thick porosity effects come into play, both in physical and velocity space. Spatial porosity allows photons that would have normally interacted with a slab of gas to escape freely. Some photons will be absorbed by gas that is compressed into dense clumps, but the separation between the clumps allows others to escape without interaction. The velocity field of the wind plays a crucial role in this; the fact that the outflow accelerates results in increasing Doppler shifts throughout the wind and causes the spectral line resonance zones in which the clumps are optically thick for a certain frequency to be very narrow in the radial direction (at least as long as the velocity is not close to υ_∞). If one then assumes that all clumps follow the underlying average velocity field, the amount of leakage for line photons can be directly linked to the spatial volume filling factor f_vol. This effect is illustrated in the middle panel of Fig. 5.

However, the clumps do not necessarily follow the average velocity field⁹. For example, the velocity span of the clumps, δυ, can be larger than that of the underlying smooth field δv_sm (Eq. (1)). In this case the Doppler shifted gas in the clumps spans a wider range of velocities than does the smooth field, which means that effectively the resonance zone in the wind for a given transition becomes larger. In other words, if the clumps are at least somewhat optically thick more gas has the right Doppler shift to absorb a photon of a given frequency, and thus more light is absorbed. This effect of porosity in velocity space, is also called velocity-porosity or ‘vorosity’ (as coined by Owocki 2008). The non-normalised velocity filling factor, f_vor, depends on both f_vol and the relative velocity span of the clumps δv/δv_sm: $$f_{\text{vor}}=f_{\text{vol}}\left|\frac{\delta \upsilon }{\delta {\upsilon }_{\text{sm}}}\right|.$$(3)

In Fastwind a normalised version of this factor is implemented: $$f_{\text{vol}}=\frac{f_{\text{vor}}}{1+f_{\text{vor}}}\frac{f_{\text{vol}}(\delta \upsilon /\delta {\upsilon }_{\text{sm}})}{f_{\text{vol}}(\delta \upsilon /\delta {\upsilon }_{\text{sm}})},$$(4)

taking values from 0 to 1. The parameter f_vel is called the velocity filling factor. It is important to note that this equation reduces to the purely geometrical effect, depending only on f_vol, when δυ follows the underlying smooth field: $$\begin{array}{cc}{f}_{\text{vel}}=\frac{f_{\text{vel}}}{1+f_{\text{vel}}},& [\delta \upsilon /\delta {\upsilon }_{\text{sm}}\to 1].\end{array}$$(5)

The effect of vorosity is illustrated in the middle and right panels of Fig. 5.

These clumping and vorosity effects are implemented in Fastwind by means of an ‘effective opacity formalism’. In this formalism, various properties of the clumps and the interclump medium (such as temperature) are assumed to be similar, and the rate equations, etc., are evaluated for a fiducial clump density P = 〈p〉f_cl. This allows one to approximate the clumpy wind as a one component medium with a certain average ionisation state and a single effective opacity. Essentially, the expensive computation of the NLTE occupation numbers is done only once, obtaining an average opacity 〈χ〉 for a mean clump density, and then re-scaled in order to infer the effective opacity of the two-component clumped wind. The effective opacity χ_eff can be expressed as: $${\chi }_{\text{eff}}=\langle \chi \rangle \frac{1+{\tau }_{\text{cl}}{f}_{\text{ic}}}{1+{\tau }_{\text{cl}}},$$(6)

with τ_cl the clump optical depth (Sundqvist et al. 2014). The interclump density contrast, f_ic, is defined as: $$f_{\text{ic}}\equiv {\rho }_{\text{ic}}/\langle \rho \rangle .$$(7)

The formalism accounts for the vorosity effects by adjusting the clump optical depth. For line opacity the clump optical depth in the rapidly accelerating winds is then computed in the Sobolev approximation: $${\tau }_{\text{cl}}=\frac{{\tau }_{\text{S}}}{f_{\text{vol}}}(1-(1-f_{\text{vol}})f_{\text{ic}}),$$(8)

with τ_S the Sobolev optical depth for the mean wind. In the case of continuum opacity, on the other hand, the clump optical depth will depend on the porosity length h (= l_cl/f_vol, with l_cl the characteristic length scale of clumps). This parameter, describing the spatial porosity, can impact optically thick continua, where it can affect, for example, the ionisation rates. By default in Fastwind a radial variation of this parameter in the form of a so-called velocity-stretch law is assumed: $$h({\upsilon }_{\text{r}})=h_{\infty }{\upsilon }_{\text{r}}/{\upsilon }_{\infty },$$(9)

with h_∞ the porosity length at the terminal wind velocity, given as input by the user. In this work we adopt h_∞ = R_⋆, following Sundqvist & Puls (2018). We refer the reader to Sundqvist et al. (2014) and Sundqvist & Puls (2018) for a more detailed and quantitative explanation of the effective opacity formalism and its implementation in Fastwind.

We conclude our description of the wind structure implementation by noting that we assume a stratified clumping factor, that is, we assume the clumping factor to vary throughout the wind. Several stratifications are implemented in Fastwind. In this work we adopt the implementation used by Sundqvist & Puls (2018) and Hawcroft et al. (2021), where the clumping is described by three parameters: the onset velocity of clumping v_cl,start, the maximum clumping factor f_cl, and the velocity at which this maximum clumping factor is reached, v_cl,max. At the base of the wind the medium is assumed to be unclumped, its structure being only affected by micro-turbulence. Then, from v_r = v_cl,start until v_r = v_cl,max the clumping factor increases linearly with wind velocity from 1 to f_cl, staying constant at f_cl for v_r > v_cl,max. This assumption for the clumping stratification is conform empirical findings in at least the lower and intermediate wind (e.g. Puls et al. 2006; Rubio-Díez et al. 2021). At v_r > v_cl,max the clumping stays constant at the maximum value, f_cl. The values for v_cl,start, v_cl,max and f_cl are specified by the user (see Sect. 3.5). A summary of the wind structure parameters is given in Table 2.

3.1.2 Wind turbulence

The wind structure parameters described in Sect. 3.1.1 are all used in the process of computing the ionisation and excitation structure of the model atmosphere. An additional parameter, the wind turbulence velocity v_windturb, is used only during the synthesis of spectral lines. This parameter introduces a depth-dependence of the micro-turbulent velocity throughout the wind. During the computation of the ionisation and excitation structure the micro-turbulent velocity is assumed to be constant, but when the lines are synthesised, the micro-turbulent velocity increases linearly with wind velocity from v_micro at the base of the wind, to v_windturb, at the point where the wind reaches its terminal velocity (Haser et al. 1995). The wind turbulence velocity is typically of the order of 0.1v_∞ (e.g. Groenewegen et al. 1989), and is used here to mimic the effects of a large wind velocity dispersion upon the spectral line formation. Evidence for such a velocity dispersion is found in both LDI simulations (e.g. Hamann 1980; Puls et al. 1993; Driessen et al. 2019) as well as in observations (e.g. Lucy 1982; Groenewegen et al. 1989; Prinja et al. 1990).

3.1.3 Wind-embedded shocks and X-rays

Instabilities in the winds of massive stars can cause shocks to form in the wind (e.g. Owocki et al. 1988; Feldmeier et al. 1997). These wind-embedded shocks give rise to X-ray emission, which, both by direct and Auger ionisation, can alter the ionisation balance of the wind, as well as the velocity fields of the interclump medium and the clumps. Wind-embedded shocks and associated X-ray emission are implemented into Fastwind, and their characteristics can be tweaked. In our analysis we include X-ray emission by assuming canonical values for each star. Details about the implementation of X-rays in Fastwind and our assumptions regarding the canonical values are detailed in Appendix G.

3.2 A new genetic algorithm: KIWI-GA

In order to find the best fitting Fastwind models we developed a Genetic Algorithm (GA), which we call Kiwi-GA. GAs employ the concept of natural selection or survival of the fittest (Darwin 1859). The goal is to find, for a given dataset (stellar spectrum), the best fitting model. For this, the algorithm starts by computing a group (generation) of random models (individuals). After computing the first generation of models (parents), the fitness of each model is assessed by comparing each model to the data. In our case, we use the χ² value as a fitness measure: $${\chi }^2={\displaystyle \sum _{i=0}^N\left(\frac{{ℱ}_{\text{obs,}i}-{ℱ}_{\mathrm{mod},i}}{{ℰ}_{\text{obs},i}}\right),}$$(10)

where N is the number of data points of the spectrum that is considered in the fit, ℱ_obs,i the observed normalised flux, ℱ_mod,i the normalised flux of the model, and ℰ_obs,i the uncertainty on the observed flux. Generally, models that have parameters that resemble the properties of the observed spectrum will be fitter (have a lower χ²) than models with parameters that are far off. By picking new (offspring) models by combining the parameters of the fittest models of the previous (parent) generation, the offspring models will generally fit the observed spectrum better than the parent generation. For example, a model with a T_eff that is similar to that of the observed star will generally have a better fit than a model with a T_eff that is far off, and this value of T_eff will thus have an increased chance of being selected for models of the new generation. We note that in the process of combining the parameters of two parent models, a fraction of the parameters is altered randomly (we call this mutation), in order to maintain and introduce diversity in the model parameters. By iterating this procedure (the offspring generation becomes the new parent generation), we eliminate parameter values that differ greatly from value that matches the data, while parameter values that match the data well will be kept. This way, the algorithm converges towards models with a better fit to the data. Figure 6 illustrates the workings of the algorithm. Especially for large parameter spaces, this is a very efficient search method. In the past GAs have been successfully used for the analysis of massive star spectra (e.g. Mokiem et al. 2005, 2006; Tramper et al. 2014; Ramirez-Agudelo et al. 2017; Abdul-Masih et al. 2021; Hawcroft et al. 2021).

Kiwi-GA is written in python and uses elements of the algorithm of Mokiem et al. (2005), who in turn use the pikaia algorithm of Charbonneau (1995). The new aspects of the algorithm are introduced after careful assessment of considerations laid out by Pohlheim (2007), who presents an overview of possible structures and operators that can be part of a GA. For Kiwi-GA we selected structures and operators that seemed beneficial for solving our specific optimisation problem; for details, we refer the reader to Appendix C. For the parallel processing we use the schwimmbad package, following Abdul-Masih et al. (2021). Within Kiwi-GA a python command initiates the execution of fastwind, a Fortran executable. Kiwi-GA is publicly available and has a comprehensive documentation in order to be accessible to new users¹⁰.

3.3 Stellar radius and luminosity

In our Kiwi-GA runs the stellar radius is estimated for each model individually using an observed, extinction corrected magnitude (Sect. 2.6), following the procedure described in Mokiem et al. (2005). Based on the temperature of each model, a Planck curve with temperature T = 0.9T_eff is computed and compared to the observed magnitude, after which a radius is chosen such that the Planck curve matches the observed anchor magnitude (Mokiem et al. 2005, who follow Markova et al. 2004). For this, we use a transmission curve of the adopted filter, in our case the VLT/SPHERE K_s filter¹¹. The Planck curve thus serves as an ‘SED estimate’ during the run and the radius is an output of the run, as is luminosity¹². When a run is finalised we compute the real SED of the best fitting model and use this to correct the radius that was estimated during the run. Furthermore we scale the obtained mass-loss rates using the optical depth invariant wind-strength parameter $Q=\dot{M}/{\left(R_{*}{\upsilon }_{\infty }\right)}^{3/2}$ (Puls et al. 1996, Holgado et al. 2018, their Appendix B). It is important to note that there is no need to recompute all models, as the radius has very little impact on the normalised spectrum. The radius corrections for our stars range from 0 to +8% with an average of 3.6% for the O-stars and from –20 to –9% for the WNh stars. We note that for future runs, where the K_s-band is used as an anchor magnitude, a Planck curve with T = 0.83T_eff would be a better guess – the previous estimate of T = 0.9T_eff was tailored to V-band magnitude anchors.

3.4 Best fit parameters and error bars

From the output of each Kiwi-GA run we derive best fit parameters and uncertainties thereof (error bars) with the method of Tramper et al. (2014). For this, we identify the best fitting model (that with the lowest χ²) and use this model to implicitly adjust the uncertainty of each flux point that is fitted, such that the ${\chi }_{\text{red}}^2$ value of the best fitting model equals unity. These adjusted flux uncertainties are then used for recomputing the χ² of each model in the run. This procedure is equivalent to dividing all χ² values of the run by the (original) ${\chi }_{\text{red}}^2$ of the best model. After the flux uncertainty adjustment we find the models that should be considered statistically indistinguishable from the best model, which we call the family of best fitting models. We do this by computing for each model the probability P that the difference between the two models is caused by random fluctuations: $$P\left({\chi }^2,\nu \right)=1-\Gamma \left({\chi }^2/2,\nu /2\right)$$(11)

with Γ(χ²/2, v/2) the incomplete gamma function, representing the cumulative distribution function of the χ² distribution, evaluated at χ²> for v = n_data − n_free the degrees of freedom, where n_data is the number of flux points that is taken into account during the fit and n_free the number of free parameters. The best fitting models are all models where P > 0.32 (i.e. the 68% confidence interval; we call this 1σ) or P > 0.05 (i.e. the 95% confidence interval). From this group we derive error bars by identifying for each parameter what is the lowest and the highest value that is present. In other words, the parameter space spanned by the family of best fitting models determines the size of the error bars. In case the distribution from which we derive the confidence intervals is symmetric and Gaussian, the 68 and 95% confidence intervals translate directly into standard deviations of 1σ and 2σ, respectively. For convenience, we refer to the 68 and 95% confidence intervals as 1 and 2σ uncertainties, even though the confidence intervals we derive are not necessarily symmetric and Gaussian. In all tables we present 1σ uncertainties, unless explicitly stated otherwise. For practical purposes, we generally mark the parameter of the best fitting model, however, we stress that all models in the family of best fitting models should be considered as statistically equivalent.

The normalisation of χ² values that is part of this method relies on the assumption that the best fitting model has a good fit to the data. This condition is satisfied for all stars in our sample except for the three WNh stars, where clear deviations between the best fitting model and data can be seen. In this case the method described above underestimates the error bars and therefore we assume increased error bars for these three stars, such that the error region covers the width of the peak in the χ² distributions. This way, the error bars of the WNh stars are more in line with those of the O-stars of the sample (Sect. 4).

For luminosity, radius and mass-loss rate we increase the error bars given the uncertainty on the magnitudes (as presented in Khorrami et al. 2017). The uncertainty in radius and luminosity is directly related to the uncertainty in the observed flux at the K-band. For the mass-loss rate, we increase the errors propagating the uncertainty on the stellar radius, assuming a scaling of $\dot{M}\propto R_{*}^{3/2}$ (see e.g. Puls et al. 1996).

Lastly, we stress that the uncertainties that we derive from the Kiwi-GA runs, as described in this section, are only statistical uncertainties. Systematic uncertainties, that could arise, for example, due to assumptions regarding extinction, normalisation, or the modelling, are not included in these values.

3.5 Fitting strategy

We fit the full sample two times with Kiwi-GA. The first time we consider only the optical parts of the data, the second time we fit the optical and UV data simultaneously. Ultimately, we are interested in the values of the optical + UV analysis, but the optical fitting still serves a threefold purpose. First, we use the derived values for rotational broadening and helium abundance as fixed values for the optical + UV fitting, reducing the amount of free parameters of those runs. Second, it provides mass-loss rates as derived from recombination lines only, assuming a smooth wind. Third, it allows us to compare our analysis method, fitting with Kiwi-GA, to the spectroscopic analysis with IACOB-GBAT (Simón-Díaz et al. 2011) of the same data by Bestenlehner et al. (2020). The second and third point are addressed in Appendix D. The details of each fitting setup are summarised in Table 3 and explained in detail below.

3.5.1 Optical-only setup

The optical-only runs have 5–7 free parameters n_free, as specified in Table 3. Here, g is the gravitational acceleration, v_eq sin i the projected rotational broadening, Ṁ the mass-loss rate, and x_He the helium surface abundance, where x_He = n_He/n_H, with n_He and n_H the helium and hydrogen number density. If any line is (partially) in emission, we also fit the wind acceleration parameter β and when we see a nitrogen line above the noise we fit the nitrogen abundance n_N (with n_N the number density). The other parameters are held fixed at the values discussed below.

We assume a smooth wind (f_cl = 1) for all stars except for the WNh stars, for which we assume f_cl = 10. Furthermore, we assume v_micro = 10 km s⁻¹, and in case β is fixed we assume β = 0.9. Because the resolution and S/N of the data do not allow us to distinguish between broadening due to rotation versus broadening due to macro-turbulence, we only fit v_eq sin i, assuming v_macro = 0 km s⁻¹. In practice this means that all broadening is captured in a single parameter v_eq sin i and, since for our stars likely v_macro > 0 km s⁻¹ (see e.g. Simón-Díaz et al. 2017), the projected rotational velocities that we find are upper limits of the actual v_eq sin i. The derived v_eq sin i is thus an upper limit. We adopted surface abundances of the CNO-elements using the evolutionary grids of Brott et al. (2011a) and Köhler et al. (2015), based on stellar parameters of Bestenlehner et al. (2020), and other abundances are fixed to Z = 0.5 Z_⊙. We assume the values of Crowther et al. (2016) for the terminal velocities of the winds v_∞. For 12 stars v∞ was not available and in these cases we estimate the velocities by inter- and extrapolating the dependence of v_∞ on luminosity L, that we empirically find using the values of Crowther et al. (2016) and Bestenlehner et al. (2020), for log L/L_⊙ < 5.6.

The optical-only Kiwi-GA runs have a population of 71 to 95 individuals, with the exception of the WNh stars, where we have 191 individuals. The runs of most stars converge in approximately 20 generations. To ensure that all runs are fully converged, we iterate for 30 generations. The runs of sources with strong emission lines converge later and we run them for 40–60 generations.

3.5.2 Optical + UV setup

The optical + UV runs have 6 to 14 free parameters. For 39 stars with relatively high S/N we fit 12 free parameters as listed in Table 3, in which n_C refers to the carbon abundance (by number). For the WNh-stars we fit two additional free parameters (see also below): oxygen abundance n_O (by number) and x_He. The other 17 stars have too low S/N and too weak wind lines to distinguish between 12 free parameters and we therefore only consider 6 free parameters for these stars (Table 3). In this case, the CNO-abundances are fixed to LMC baseline values, for which we assumed log(n_C/n_H) + 12 = 7.75 for carbon, log(n_N/n_H) + 12 = 6.9 for nitrogen and log(n_O/n_H) + 12 = 8.35 for oxygen (Kurt & Dufour 1998, as in Brott et al. 2011a; Köhler et al. 2015). The wind structure parameters are fixed based on typical values we find from the 12-free-parameter runs with lower mass-loss rates: f_ic = 0.05, f_vel = 0.15, v_cl,tart/v_∞ = 0.05, v_windturb/v_∞ = 0.15 (Sect. 5.2). For all optical + UV runs the velocity at which the maximum clumping factor is reached is given by v_cl,max /v_∞ = max(0.3, v_d,start/v_∞).

Oxygen abundance is only a free parameter for the WNh stars. Test runs with free oxygen abundance for the other stars resulted in extremely high values of log(n_O/n_H) + 12 = 9−10. We suspect that this is related to the fact that we have only two oxygen lines in our spectra, both in UV, where there is overlap with various iron lines. We therefore fix it based on the evolutionary grids of Brott et al. (2011a) and Köhler et al. (2015), based on stellar mass, rotation and age as derived by Bestenlehner et al. (2020).

The optical + UV KIWI-GA runs have a population of 95 or 191 individuals (6 or 12 free parameters, respectively). For the WNh stars we have 239 individuals (14 free parameters). The runs of most stars converge in approximately 20 or 40 generations (6 or 12–14 free parameters), so to be on the safe side, we iterate for 30 or 60 generations (6 or 12–14 free parameters). The limits within which each parameter is allowed to vary can be read off from fitness plots shown in the run overview of each star, which can be found on zenodo¹³. We discuss the robustness of this setup in Sect. 4.7.

3.5.3 Diagnostic line selection

We use all strong spectral lines that are present in our spectra that can be synthesised with Fastwind V10.3.1 and where interstellar absorption or emission did not pose a problem. For the optical, these are lines of H, He I, He II, N III, N iv and/or N V. No optical C or O lines were available due to the limited optical wavelength range and moderate S /N and resolution. In some cases data quality was too poor to include a certain line and the line was omitted from fitting; in particular this concerned Hα for 11 sources. We included the following UV lines in the fits of most stars: C iv λ1165, C iii λ1170, N v λ1240, O iv λ1340, O v λ1371, C iv λλ1548–1551, He ii λ1640. It is important to note that:

• For N v λ1240 we only fit the part that could be recovered after the Ly-α correction (see also Sect. 2.5).

• In one case we fitted Si iv λλ1394–1402 (R136b, see also Sect. 2.5).

• Part of O iv λ1340 is clipped because of the presence of the strong interstellar CII 1336 line.

• For stars where O v λ1371 was very weak (cooler stars), we omit it from the analysis completely, as, in those cases, the iron pseudo-continuum dominates the absorption.

• N iv λ1718 was included for about half of the stars. In cases where we had GHRS data we included the full line. In other cases, we included the blue part of the line from the HST data, but only if this was clearly visible in absorption.

An overview of the spectral lines used for the analysis of each individual star is presented in Table I.3.

We note that the UV spectroscopy includes diagnostics for Ṁ, v_∞ as well as for the wind structure parameters f_cl, f_ic, f_vel, so that we can break the degeneracy between these parameters. For example, the strength of Hα depends on the density squared, whereas resonance lines typically depend linearly on density, and so respond to clumping differently (e.g. Puls et al. 2006). Clumps often become optically thick in strong resonance lines, while recombination lines are generally less affected; nonetheless, vorosity effects can sometimes also result in extra light leakage in recombination lines, resulting in weaker profiles (Sundqvist et al. 2010, 2011; Oskinova et al. 2007; Šurlan et al. 2012). Bouret et al. (2005) find that O v λ1371 and N iv λ1718 are also indirectly (due to a modified ionisation structure) sensitive to optically thin clumping, where typically the absorption part of the lines get weaker for higher clumping factors. A non-void interclump density further affects line saturation, for example in the case of N v λ1240 (Zsargó et al. 2008; Šurlan et al. 2012, 2013; Sundqvist & Puls 2018). In particular, both the absorption and emission parts of the line profiles get stronger with a larger interclump density, where Šurlan et al. (2012) find that the effect is most pronounced for the strong lines. Lastly, the onset of clumping affects the line shape close to the line centre (Bouret et al. 2003; Šurlan et al. 2012).

4 Results

For 39 stars, we obtain 14 stellar and wind parameters per star, for the remaining 17 stars, 8 parameters. For the WNh stars we additionally obtain oxygen surface abundance, as their oxygen lines are very strong and dominate the iron pseudo-continuum. A representative example of an output summary is presented in Fig. 7. The top half of the figure shows that the agreement between models and data is good: the best fitting model and the family of best solutions (2σ confidence region) cover the error bars on the data both for the optical as well as the UV data. The bottom half of the figure contains the goodness of fit for all computed models. This is illustrative for the way we derive uncertainties on all parameters: if the fitness distribution of a certain parameter is strongly peaked, the uncertainties on that parameter are small; if it is wide, the uncertainties are large. Output summaries for the other stars can be found on zenodo¹³.

The best fit parameters of the optical + UV runs for all stars are presented in Figs. 8 and 9 and Tables A.1 and A.2. Notes on peculiarities of individual sources can be found in Appendix F. In Appendix I we present additional values derived from the optical + UV runs such as stellar masses, ages and ionising fluxes, as well as best fit values of the optical-only runs. We note that we derived several parameters from both the optical-only as well as the optical + UV analysis. In the remainder of the paper, unless specified otherwise, we show and interpret only one set of values: x_He and v_eq sin i were taken from the optical-only analysis¹⁴ while the remaining parameters were taken from the optical + UV analysis. The WNh stars are an exception: here we do measure x_He in the optical + UV fit, so in this case we use this value instead of the optical-only value. Lastly, our values generally agree well with the stellar properties derived by Bestenlehner et al. (2020) based on optical spectroscopy only. A detailed discussion of the comparison of different methods can be found in Appendix D. In the rest of this section we highlight several results that deserve special attention, and conclude with a robustness analysis.

4.1 Hertzsprung-Russell diagram

Figure 10 shows the derived temperatures and luminosities in a Hertzsprung-Russell diagram (HRD). We review the HRD positions of the stars and inspect the fit of all sources by eye to check for irregularities. From this we conclude that our temperature and luminosity measurements are reliable for all stars except H129 and R136a3. For R136a3 we find a low temperature (42 000 K), but see in the spectrum strong lines of higher ionised ions such as N v and O V. These lines are not matched by the best fit model, where they are weak. A higher temperature for this star thus seems more likely and we test this with an additional run where we fix T_eff to 50 000 K, the value found by Bestenlehner et al. (2020). Although this decreases the fitness of several other lines – and in such a way that we obtain a worse overall fitness – a higher temperature does improve the fit to the N v and O v lines, and places the star closer to the other WNh stars in the HRD. Lastly, from the fit of the iron lines in the UV we found a best fitting temperature of 47 000 K, where the fitness of the 50 000 K model is almost as good, but the fit to the 45 000 K model significantly less, even worse for the 42 000 K model (Fig. F.2). Taking all this into account, we consider the higher T_eff for this star more likely and accept the parameters of the fixed-T_eff run (50 000 K) as the parameters which we use for further analysis (for more details and the spectral fits, see Appendix F.3). The best fitting models of H129 seem to fit the data well, however, the S/N of this source is very low and its position in HRD is far left of the main sequence where we do not expect any O-type stars. Bestenlehner et al. (2020) reported for this star a total-to-selective extinction R_V that was 5σ below the average for R136. This could point to NIR-excess, although this would imply an even lower luminosity for this star, keeping it in the improbable HRD region.

For the subsequent analysis we only use temperatures from our optical + UV analysis. In Appendix D.1 we present a detailed comparison of the different temperature measurements (our three different measurements, plus those of Bestenlehner et al. 2020). Generally, the temperatures agree within their uncertainties.

4.2 Stellar mass and age

In order to derive the evolutionary mass Mevol, the initial mass M_ini and the age τ we use Bonnsai¹⁵ (Schneider et al. 2014) combined with the grids of Brott et al. (2011a) and Köhler et al. (2015). Bonnsai is a Bayesian tool that allows us to compare observed stellar parameters to stellar evolution models in order to infer full posterior distributions of key model parameters such as initial mass and stellar age. Our input parameters are observed luminosity, temperature, helium surface abundance and surface gravity. We use standard settings except for the prior of the initial rotational velocity, for which we assume the distribution of Ramírez-Agudelo et al. (2013) instead of a flat distribution. We find a robust output for all stars except three. For those the observed values match poorly with the posterior distribution of the Bonnsai run. In the case of H129 the value for T_eff cannot be reproduced given the observed log L/L_⊙ and log g. For this star, we deemed our derived luminosity measurement unreliable (Sect. 4.1), and we exclude this source from further analysis. In two other cases, R136b and H30, our observed log g value lies in the P < 0.05 tail of the posterior distribution. Therefore both spectroscopic and evolutionary parameters should be treated with care, although the spectroscopic fits of these stars look good. We do include both sources in further analyses that need M_evol as an input, but check whether the results change drastically upon inclusion or exclusion of R136b and H30, which is not the case. The derived ages and masses can be found in Table I.1. We cross-check our Bonnsai output with that of Bestenlehner et al. (2020) and find generally good agreement, see Appendix E for details. In the remainder of the paper, we use the Bonnsai evolutionary masses when we need stellar masses for our analysis.

4.3 Terminal velocity

For all stars we have set the terminal wind velocity v∞ as a free parameter in the optical + UV fit. For 46 sources we were able to accurately constrain v_∞, albeit with large uncertainties for the stars with lower mass-loss rates (see Figs. 8 and 9). For three of the remaining sources (R136a2, R136b, H36) we do find a tightly constrained value, but see that the fit to the blue wing of C iv λλ1548–1551 is not good: the saturated absorption edge of the best model for these stars extends about 400 km s⁻¹ more to the blue as the absorption edge we see in the data. For the remaining seven sources (H73, H121, H135, H141, H159, H162, and H173) the wind lines crucial for determining v_∞, especially C iv λλ1548–1551, turn out to be too weak to get a constraint: the χ² distribution for v_∞ was flat. In these cases, while we do find some best fit value v_∞, this value is not meaningful and we regard v_∞ as unconstrained.

4.4 Wind acceleration parameterß

The wind acceleration parameter β is fitted for all stars. For stars with log Ṁ > −5.7, we find values up to β = 1.50 with an average of 1.08 ± 0.20, whereas for stars with lower mass-loss rates we find that for all but two sources β is consistent with 0.7 within 1σ errors, with an average of 0.72 ± 0.06 (see Figs. 8 and 9). We note that β = 0.7 was the lowest allowed value during the fitting, this is discussed in Sect. 4.7.

4.5 Onset of clumping

We derive the onset of clumping for 39 sources and find a value of ${\upsilon }_{\text{cl,start}}=0.07{\pm }_{0.07}^{0.10}{\upsilon }_{\infty }$, translating into $R_{\text{cl,start}}=1.02{\pm }_{0.02}^{0.07}{R}_{*}$ on average. There is not much variation across the sample: two thirds of the stars have a value of v_cl,start < 0.1 v_∞ or R_cl,start < 1.08 R_⋆, the higher values that we derive have large error bars – within 2σ all sources are consistent with v_cl,start < 0.05 v_∞. This is visible in Fig. 8.

4.6 Ionising flux

We derive H, He i and He ii ionising fluxes Q₀, Q₁, Q₂ for each star¹⁶ based on the best fitting model (Table We estimate the errors on these values by computing, for one star (H35, spectral type O3 V), the ionising flux for each model in a full Kiwi-GA run; afterwards we apply an error analysis based on the χ² value of each model, as we do for all other parameters¹⁷. From this we find 1σ uncertainties on the derived ionising fluxes to be, approximately, 0.07 dex for log Q₀, 0.1 dex for log Q₁, 0.4 dex for log Q₂. We assume that the uncertainties of the ionising fluxes of the other sources scale with their relative uncertainties (compared to those of H35) on effective temperature and luminosity. Lastly, in Table I.1 we provide also the H-i ionising luminosity, that is, the energy of each star emitted by photons capable of ionising hydrogen, a quantity relevant for large scale simulations involving radiative feedback of massive stars.

4.7 Robustness and systematic errors

In order to check whether our results are robust under small changes in our setup we carry out several test runs. We picked the O3 v star H35, a typical source¹⁸, and fit its spectrum many times, each time changing one aspect of the setup. As a reference point, we compare all runs with the ‘fiducial run’ for H35, that is, the run with the setup as used for all optical + UV runs in this paper. In Fig. 11 this comparison is presented.

First, we show the robustness of the genetic algorithm itself by redoing our fiducial run. We must be sure that the initial pool of individuals contains enough variation. If the variation is large enough to cover the full parameter space, then with exactly the same setup but different random initial parameters, we should get the same or very similar results. Indeed, when we do this test we do see small differences, but generally the agreement between the two runs is very good: for each parameter the 1 and 2σ regions are similar and the best fitting parameters of each runs lie in the 1σ error regions of the other run.

Having done this initial test we then vary the setup, changing one aspect at the time. We see that within uncertainties the different setups show consistent results and our setup is robust to most changes. The choice for micro-turbulence v_micro seems to have the largest effect on the derived parameters, especially T_eff. Changing the value from our fiducial fixed value of 10–5 km s⁻¹ does not have much effect, but changing it to either a fixed value of 15 km s⁻¹ or leaving it as a free variable results in a best fit value of T_eff that is ~2000 K higher than that of our fiducial run, just on the edge of the 2σ error bars. From the data that we have we cannot determine what is the actual value of v_micro and thus of T_eff: the run with v_micro as a free parameter resulted in a velocity exceeding the typical sound speed in the atmospheres of these stars (v_micro = 30 km s⁻¹) and thus seemed not reliable (though for a recent finding, see Schultz et al. 2022). Apart from changes in T_eff we note that also abundances change when we assume a different value for v_micro, the largest change being seen in the carbon abundance when v_micro is lowered to 5 km s⁻¹. This is expected as v_micro impacts the equivalent width of lines. Considering the above, we must thus keep in mind that the lack of atmospheric lines from which we can accurately determine micro-turbulence leads to systematic uncertainties in T_eff and abundances, which we estimate to be about 2000 K and 0.5 dex, respectively.

The mass-loss rates and high clumping factors that we derive are robust within the optically thick clumping framework. We consistently find f_cl > 10, but distinguishing between the higher clumping factors proves difficult. This could be due to the fact that the clumping sensitivity saturates towards higher clumping factors for several of the clumping diagnostics. Since leaving oxygen abundance as a free parameter consistently leads to very high oxygen abundances, we decided to keep the oxygen abundance fixed during our fits. Possible causes for the high obtained oxygen abundances could be blends with iron lines, which are not present in our synthesised model spectrum, or specific shortcomings in the Fastwind oxygen model atom, so that ionisation structure of these ions is not well reproduced by Fastwind. Regardless of the cause we checked the robustness of our results given this uncertainty by doing a run with oxygen abundance left free, and one where we left out both oxygen lines. In both cases, the resulting fit parameters change slightly but are within errors consistent with our fiducial run. This also holds for the wind structure parameters f_ic and f_vel, which show to be unaffected by either leaving the abundance free nor leaving out the lines (see Fig. 11).

For stars with comparatively low mass-loss rates of log Ṁ ≲ –5.8 we find an average wind acceleration parameter of β = 0.72 ± 0.06, with many of the derived values at 0.7. Since this is the lowest value allowed in our fits, we test the effect of extending our parameter space: for four stars we do another run with the same setup except that now we extend the allowed range of β values up to as low as 0.5. We find in all cases that the distribution is nearly flat between 0.5 and 0.7. In these runs β = 0.7 remains in the 2σ error range. Fig. 11 shows that for H35 with the lower best fit value of β the other parameters do not change significantly given the uncertainties.

Apart from the aspects discussed so far, we do also change the prescription for v_cl,max and the X-ray setup. Lastly, we do a run with only UV data. The results seem robust to all these changes. The run with UV data only shows the diagnostic power of these relatively few spectral lines. The optical data adds most to the accuracy of the gravity, though one has to keep in mind that in this ‘UV-only run’ rotation and helium abundance are fixed to values derived from optical data. In conclusion, our fitting setup seems generally robust to the assumptions we made. However, one should be aware of possible systematic errors, especially with regard to uncertainty due to the micro-turbulence that seems to affect the derived T_eff and abundances.

5 Discussion

We discuss our results in the context of theoretical predictions and evolutionary models. In Sect. 5.1 we compare the mass-loss rates that we obtained to the predictions of Vink et al. (2000, 2001); Krticka & Kubát (2018) and Björklund et al. (2021). Here, we also compare the observed and predicted terminal velocities and the modified wind momenta. We conclude this section with a comparison to the CAK-type mass-loss theory of Bestenlehner (2020), and provide an equation for mass-loss rate as a function of the Eddington parameter for electron scattering. The mass-loss rates used in this section rely on the simultaneous fit of the wind structure (clumping) parameters. We observe weak trends in these parameters as a function of mass-loss rate, which is discussed in Sect. 5.2.

The stellar evolution, mass and age of the sources based on the optical data is already discussed in detail by Bestenlehner et al. (2020). After briefly reviewing consistency with their results (Sect. 5.3.1), we add to their discussion based on additional clues we can get from the abundances based on UV spectroscopy (Sect. 5.3.2). Furthermore, Sect. 5.3.1 contains a discussion of the surface gravities and mass estimates that we find for the three WNh stars. We conclude our discussion with a more technical topic, namely the comparison of the terminal velocity measurements from by-eye fitting (Crowther et al. 2016) to our optical + UV spectral analysis (Sect. 5.4).

5.1 Mass-loss rates, wind momentum and terminal velocities

In Fig. 12 we compare our observed mass-loss rates to theoretical predictions of Vink et al. (2001, hereafter Vink01), Krtička & Kubát (2018, hereafter Krtič18) and Björklund et al. (2021, hereafter Björk21). All three predict mass-loss rates of hot luminous stars by computing the line force based on the density and ion-isation structure of a model atmosphere: Vink01 use isa-wind (de Koter et al. 1993), Krtič18 use Metuje (Krtička & Kubát 2010, 2017) and Björk21 use Fastwind (version 11: Puls et al. 2020). Where Vink01 uses a Monte-Carlo method in combination with the Sobolev approximation for the line computation, Krtič18 and Björk21 perform radiative transfer in the co-moving frame, though the former effectively places their critical point upstream from that of Björk21 (see Sundqvist et al. 2019). Furthermore the codes differ in their approach regarding the wind dynamics; Krtič18 and Björk21 solve numerically the equations of motion and Vink01 use a prescribed velocity structure assuming conservation of total radiative and kinetic energy. We note that for this comparison we use the Vink01, Krtič18 and Björk21 prescriptions as presented in their respective papers, even though they assume different values for the solar abundances: Vink01 on the one hand, and Krtič18 and Björk21 on the other hand assume solar metal mass fractions of Z_⊙ = 0.019 and Z_⊙ = 0.013, respectively. For a mass-loss prediction at Z = 0.5 Z_⊙ in the Vink01 prescription one thus implicitly assumes a metal content of a factor 1.46 higher than one would under the same assumption (Z = 0.5 Z_⊙) in the Björk21 or Krtič18 prescriptions. Correcting for this would bring the Vink01 relation approximately 0.13 dex closer to the Björk21 prescription, on average. Other differences in assumptions between the approaches may also result in systematic differences, for example, related to micro-turbulence, line lists and temperature structure.

We assess the goodness of fit of each prescription to our observations for three luminosity regimes: low (log L/L_⊙ < 5.3), intermediate (5.3 < log L/L_⊙ < 6.2) and high (log L/L_⊙ > 6.2). Results are shown in Fig. 12. We note that the highest luminosities that we observed lie outside the grids of Vink01, Krtič18 and Björk21, which extend only to log L/L_⊙ = 6.00, log L/L_⊙ = 6.07 and log L/L_⊙ = 6.25, respectively. Furthermore, it is important to note that for four sources (H73, H116, H121, and H162), for which we derived mass-loss rates in the range log Ṁ = –8.68 to –8.42, the 1σ error bars range (nearly) to the edge of our parameter space, log Ṁ = –9.50. In our fitting and with the derivation of the values, we treated the derived values and uncertainties the same as those of the other points, even though these are technically upper limits. Removing these four points from the fit does not alter the results significantly. The same holds for fitting and χ² values of the modified wind momentum (see below); here we have one additional source with only an upper limit, H135, for which we derive an upper limit only for v_∞.

Comparing observed and theoretical mass-loss rates (Fig. 12), we see that, overall, the predictions of Krtič18 fit best to our observations, outperforming the other two in the regime where mass loss impacts evolution the most (log L/L_⊙ > 5.3), and matching almost perfectly with the linear fit through our observations. In the low-luminosity regime the Björk21 rates match better with our observations: the prediction is in good agreement with the group of low-luminosity observations, reflected in the low ${\chi }_{\text{red}}^2$. The Vink01 rates there are on average an order of magnitude too high in this regime, though several individual points lie well below that average, close to the observed points. Observing weak winds in this luminosity regime is not unusual (e.g. Puls et al. 1996, 2008; Martins et al. 2005; de Almeida et al. 2019). In the intermediate regime, the performance of the three predictions is similar. The Krtič18 rates perform best; their predictions are, on average, only a factor 1.3 lower than the observed rates. The Björk21 rates are a close second, being too low with, on average, a factor 1.5. The Vink01 rates in the intermediate regime are, on average, a factor 2.0 too high, consistent with previous findings (e.g. Bouret et al. 2012; Šurlan et al. 2013; Cohen et al. 2014). However, it is important to note that the overestimation of the Vink01 rates in this regime decreases to only a factor 1.5 if one would apply an average downward shift to correct for the different assumptions of Vink01 versus Björk21 and Krtic18 regarding LMC metallicity. Inspecting the high-luminosity end we see that the Krtič18 rates best match the observations. Both Björk21 and Vink01 overestimate the mass-loss rates of the most luminous stars; of those two, the prediction of Vink01 lies closest to the observed rates.

We conclude the mass-loss rate comparison with noting that recently Vink & Sander (2021) updated the Vink01 Monte Carlo mass-loss recipe with dynamically consistent computations of the terminal wind velocity. The mass-loss rates they predict are similar to those of Vink01, but typically lie a bit higher; for our sample, the Vink & Sander (2021) rates lie on average 0.11 and 0.17 dex above the Vink01 rates, with the average absolute difference being 0.19 and 0.17 dex, using the predicted and observed terminal velocities, respectively. Given that the Vink01 rates generally overpredict the observed mass-loss rates, the updated recipe does so even more and therefore the Vink01 recipe outperforms the Vink & Sander (2021) recipe for our sample.

Another way of comparing the predictions to our observations is by their modified wind momentum, $D_{\text{mom}}=\dot{M}{\upsilon }_{\infty }\sqrt{R_{*}/R_{\odot }}$ (with M and v_∞ in cgs-units), shown in Fig. 13. Björk21 provide an explicit equation for the wind momentum as a function of luminosity and metallicity. Krtič18 do not provide this, instead we compute D_mom for all the models in their LMC grid, and obtain a linear fit through these points as a function of luminosity. Since Vink01 do not predict terminal velocities, their relation for wind momentum (Vink et al. 2000, their Eq. (15)) is semi-empirical; they use observed values for v∞ and R_⋆. The relation plotted in Fig. 13 is corrected for lowered Ṁ and v_∞ as a result of the lower metallicity in the LMC compared to the Milky Way (Vink01; Leitherer et al. 1992). For the modified wind momentum the Björk21 predictions match best our observations in all luminosity regimes. The Vink01 prediction is too high over the full luminosity range, which is to be expected given their overprediction of the mass-loss rates. The Krtič18 predictions lie close to that of Björk21, but their prediction is less steep, translating in underpredicting D_mom for the higher luminosities, and overpredicting it for the less bright stars.

Figure 14 shows observed and predicted terminal velocities as a function of both log L/L_⊙ and T_eff. Neither the predictions of Björk21 nor those of Krtič18 fit the observed values well. Both predict decreasing terminal velocities as a function of T_eff and log L/L_⊙, whereas observations show the opposite. Also in absolute sense the predicted velocities deviate from the observed values. Looking at the left side of the figure, we see that the difference between the observed velocities and those predicted by Björk21 is especially large in the lower luminosity regime, where the Björk21 and Vink01 mass-loss rates diverge the most. The predicted terminal velocities of Krtič18 in this regime match reasonably well our observations, but their predictions are too low in the intermediate regime.

Of all terminal velocity predictions considered here, those of Vink & Sander (2021) are most in line with observations¹⁹; both the absolute values as well as the trend as a function of luminosity are in line with observations (see the green dashed lines in the left panels of Fig. 14), except around 35 000 K, where they predict extremely high velocities. In Fig. 14 we do not show their predictions for T_eff = 35 000 K, and do not include them in the fit.

It is clear that none of the theoretical models considered here can reproduce both the observed mass-loss rates and the terminal velocities that we observe. It has to be noted, however, that while the luminosity range spanned by the hydrodynamic models covers our observations reasonably well, this is not the case for the temperatures of the models: the maximum model T_eff in the grids of Björk21, Krtič18 and Vink & Sander (2021) is 45 000 K, while about 50% of our sample has a T_eff higher than that. The Vink01 grid extends further to 50 000 K, however they do not predict terminal velocities.

We now look more quantitatively at the modified wind momentum, and fit a powerlaw to our observations as a function of luminosity (dark blue dashed line): $$\log D_{\text{mom}}=\log D_0+x\log L/L_{\odot }.$$(12)

For this and other fits where we have to take into account uncertainties in both coordinates we use the orthogonal distance regression (ODR) routine of scipy. Our best fit yields x = 2.00 ± 0.11 and log D₀ = 17.05 ± 0.65. A comparison to other LMC studies shows that this is relatively high: observations of Mokiem et al. (2007); Bestenlehner et al. (2014); Ramírez-Agudelo et al. (2017); Sabín-Sanjulián et al. (2017, grey circles in Fig. 13) give slopes in the range 1.45–1.87. We note, however, that not all these analyses take into account uncertainties in D_mom and/or L in the same way, which might affect the derived slopes (see Markova et al. 2004). Furthermore, beware that the relatively high wind momenta for low-luminosity stars found by other LMC studies (grey circles in the leftmost part of the plot) are likely only upper limits, as these studies lack UV coverage and are therefore rely only on Hα and He ii λ4686 for their mass-loss determinations.

As shown by Puls et al. (1996), the slope x of the modified wind momentum-luminosity relation can be interpreted as a measure for the distribution of line strengths of the spectral lines contributing to the wind driving: it is the inverse of the force multiplier α_eff in (modified) CAK-theory (Castor et al. 1975, with subsequent modifications by Abbott 1982; Pauldrach et al. 1986): x = 1/α_eff. Here, α_eff captures both α, the slope of the line strength distribution, as well as the force multiplier δ, that accounts for the ionisation state of the wind in an approximate way (Abbott 1982), that is, x⁻¹ = α_eff = α − δ in Eq. (12). Our slope of x = 2.00 ± 0.11 then translates into a mean value of α_eff = 0.50 ± 0.03 or α = 0.60 ± 0.03 if we assume the typical value of 0.1 for δ (Abbott 1982; Puls et al. 2008). This is in line with typical values expected for O-stars, α ≈ 0.5−0.6 (Abbott 1982; Pauldrach et al. 1994; Puls et al. 2000, 2008).

We can obtain α_eff from our data in a different manner by inspecting the dependence of mass-loss rate on the Eddington parameter for electron scattering. Bestenlehner (2020) extended the CAK mass-loss prescription (Castor et al. 1975) from the regime of optically thin to that of optically thick winds, by expressing the stellar mass in terms of the electron scattering Eddington parameter Γ_e. The mass-loss rate can then be expressed as a function of Γ_e, the transition mass-loss rate Ṁ_e,trans and the force multiplier α_eff: $$\begin{array}{c}\log \dot{M}=\log {\dot{M}}_{\text{e,trans}}\\ +\underset{\text{I.Dominateswhen\hspace{0.17em}}{\Gamma }_{\text{e}}\ll 1}{\underbrace{\left(\frac{1}{{\alpha }_{\text{eff}}}+0.5\right)\log ({\Gamma }_{\text{e}})}}-\underset{\text{II.Dominates\hspace{0.17em}when\hspace{0.17em}}{\Gamma }_{\text{e}}\to 1}{\underbrace{\left(\frac{1-{\alpha }_{\text{eff}}}{{\alpha }_{\text{eff}}}+2\right)\log (1-{\Gamma }_{\text{e}}).}}\end{array}$$(13)

Our α_eff is equal to what Bestenlehner (2020) calls α, and what Pauldrach et al. (1994) and Puls et al. (1996, Puls et al. 2008) call α′. The transition mass-loss rate corresponds to the transition Eddington parameter Γ_e,trans. At Γ_e = Γ_e,trans the first (I) and the second (II) terms in Eq. (13) are equal. At this point the mass-loss dependency changes from being dominated by the first term to being dominated by the second term.

We compute Γ_e for all our sources²⁰ using the Bonnsai evolutionary masses and fit Eq. (13) to our observations, in order to derive a mean value for α_eff for the sample, as well as the transition mass-loss rate. We find α_eff = 0.46 ± 0.04 (α ≈ 0.56) and log Ṁ_e,trans = −5.19±0.10 (see Fig. 15, red solid line). Our value for Ṁ_e,trans matches well with the mass-loss rate of log Ṁ = −5.2 ± 0.2 that Vink & Gräfener (2012) find for transition objects in the Arches cluster. This agreement was not strictly expected as the definitions for the transition mass-loss rates of Bestenlehner (2020) and Vink & Gräfener (2012) differ, with the transition of Bestenlehner (2020) relating to different terms in Eq. (13), and Vink & Gräfener (2012) deriving their rate based on the inference that stars with a spectral type O4-6If+ correspond to the transition from optically thin to optically thick winds. The α_eff we find from the fit with the Bestenlehner (2020) prescription is lower than we found before using the wind momentum relation, however, taking into account errors on both α_eff and Ṁ_e,trans we see that α_eff = 0.50 found from the wind momentum is just within the 1σ uncertainty range (Fig. 15, shaded area and dashed line). In the high mass-loss regime the relations using the different values for α_eff barely differ, for the low mass-loss regime one sees that the rates match better the low value of α_eff = 0.46.

Filling in the best fit values for Eq. (13), we obtain the following empirical mass-loss rate dependence on the Eddington parameter for electron scattering²¹: $$\log \dot{M}=-5.19+2.69\log ({\Gamma }_{\text{e}})-3.19log(1-{\Gamma }_{\text{e}}).$$(14)

This equation could be used as a mass-loss rate prescription for stellar evolutionary computations for massive stars in the LMC. While our relation is derived based on stars of M > 15 M_⊙, the scatter up to M = 40 M_⊙ (Γ_e,trans = 0.2) is large and the best results will be obtained for stars with M > 40 M_⊙.

5.2 Wind structure

Upon investigating possible trends in wind structure parameters, we plot the obtained values against mass-loss rate. The results are shown in Fig. 16. For computing the averages ℳ and the errors on the averages $\mathcal{U}$ quoted in this section we weigh the best fit value for each individual star, x_i, with the inverse of its 2σ uncertainty (w_i): $$\begin{array}{ccc}ℳ=\frac{{\sum }_{i=0}^{N_x}{w}_i{x}_i}{{\sum }_{i=0}^{N_x}{w}_i}& \text{and}& \mathcal{U}=\frac{{\sum }_{i=0}^{N_x}{w}_i{(x_i-ℳ)}^2}{{\sum }_{i=0}^{N_x}{w}_i}\end{array},$$(15)

where N_x is the number of measurements of the quantity under consideration.

The left panel a of Fig. 16 shows the derived clumping factors. We find that all but five stars have a best fitting clumping factor of f_cl > 10, with an average for all stars of 〈f_cl〉 = 29 ± 15. When dividing the sample in two groups based on their mass-loss rates being lower or higher than log Ṁ = −6, we find typically lower values for f_cl for the stars in the higher mass-loss rate group, although this difference is barely significant; average values and uncertainties are displayed at the top of Fig. 16. Hawcroft et al. (2021), who analyse a sample of eight Galactic O-supergiants with mass-loss rates log Ṁ > −6 in a similar fashion, find 〈f_cl〉 = 25 ± 4 for their sample. This is consistent with our findings.

The middle panels b and c show our values for f_ic and f_vel. We find averages of 〈f_fic〉 = 0.13 ± and 〈f_vel〉 = 0.46 ± 0.39. As before, we also divide the points in two groups based on log Ṁ and compute the averages of each group (displayed at the top of Fig. 16). For the interclump density contrast we find a significant difference between the two groups. The stars with a lower mass-loss rate typically have a lower fic than stars with a higher mass-loss rate (values and errors at the top of Fig. 16). In other words for stars with lower mass-loss rate we find a stronger density contrast between the clumps and the interclump medium, or, equivalently, a relatively lower interclump medium density. We also find a trend in the velocity filling factor, where for the groups with a lower and higher mass-loss rate we find lower and higher values of fvel respectively. We compare our values of fic and fvel to those of Hawcroft et al. (2021) and find that these are generally in agreement with our findings. However, given the small range of mass-loss rates of the stars of the sample that Hawcroft et al. (2021) consider, this comparison cannot give any confirmation of the differences we find between the low- and high mass-loss rate groups.

The last parameter that is related to the wind structure is the wind turbulence v_windturb, shown in the rightmost panel d. Also here, we observe a weak trend where the turbulence seems to be less strong in the stars with higher mass-loss rates. Hawcroft et al. (2021) do not measure this parameter.

For all wind structure parameters, the observations show tentative trends as a function of mass loss. Overall, it appears that the stars with higher mass-loss rates typically have smoother winds than the stars with lower mass-loss rates. All wind structure diagnostics indicate this: the stars with higher mass-loss rates have on average lower clumping factors, a lower contrast between the density in the interclump and clump medium, less wind turbulence, and higher velocity filling factors. The latter may sound like evidence for stronger clumping effects, however a high velocity filling factor too can indicate smooth wind. Namely, as f_vel → 1, this ‘erases’ the density contrast, so that the wind is fully consisting of absorbing material and there are no gaps in velocity space through which the light can escape. In other words, the velocity-porosity effects are no longer present, as it would in either a smooth wind or a wind with only optically thin clumps. This follows explicitly from the equations for the effective opacity (Sundqvist & Puls 2018), summarised in our Sect. 3.1.1: if f_vel → 1, then f_vor → ∞ (Eq. (4)), leading to τ_cl → 0 (Eq. (8)), which then means that χ_eff → (x) (Eq. (6)), such that the limit for either optically thin clumping (if f_cl > 1) or a smooth wind (if f_cl = 1) is recovered. We do stress that, although the wind structure parameters of the higher mass-loss rate stars imply that their winds are on average smoother than those of the low mass-loss rate stars, their clumping factors are still significant (≥4 with an average of 21 ± 15). The high measured velocity filling factors are thus pointing to optically thin clumps, rather than a smooth wind.

We note that while we discuss these trends only as a function of mass-loss rate, we obtain very similar results when we plot the wind structure parameters against other stellar or wind properties, such as luminosity. This is due to strong correlations between mass-loss rate and stellar properties. The interested reader can find the wind structure parameters plotted against luminosity, temperature, Eddington parameter for electron scattering or wind acceleration in Figs. H.1–H.4.

5.3 Evolution

Since the stellar parameters we derive are generally consistent with those of Bestenlehner et al. (2020), we do expect similar results for the age and initial mass distributions of the cluster. Indeed this is what we find when analysing the HRD positions of our sources (see Fig. 10 and Sect. 4.2): we derive a cluster age of 1–2.5 Myr (median: 1.46 Myr) and an initial mass function with a power law slope of γ = 1.99 ± 0.11 in the range 30-200 M_⊙. Detailed comparisons between this work and Bestenlehner et al. (2020) can be found in Appendices D and E. The rest of this section will focus on mass determination of the WNh stars and the evolution of the stars in the context of surface abundances.

5.3.1 Mass and age of the WNh stars

Figure 17 shows the HRD positions of these stars as we find from our analysis, and compares these to the positions derived by Crowther et al. (2010, UV/optical/near-IR spectroscopy) and Bestenlehner et al. (2020, optical spectrocopy). The position of R136a3 as found from our analysis is indicated twice as its effective temperature is hard to constrain; the diamond marks the higher temperature HRD position that we adopt in this discussion, the cross the less likely alternative (see Sect. 4.1 and Appendix F.3 for details).

Looking at the left hand side of Fig. 17 we see that there is a considerable spread in observed temperatures and luminosities for the WNhs stars. However, while the temperature we derive for R136a1 is 7000 K lower than that of Crowther et al. (2010), the initial masses from our and their analyses agree rather well: $273{\pm }_{36}^{25}{M}_{\odot }$ and $320{\pm }_{40}^{100}{M}_{\odot }$, respectively. This is can be explained by differences in the evolutionary tracks used to derive these masses. The right hand side of Fig. 17 shows how different assumptions for the evolutionary computations can lead to divergent theoretical predictions. Nonetheless, it is clear that regardless uncertainty in both observations and theory, the previously accepted initial mass-limit (Figer 2005) challenged by Crowther et al. (2010) is indeed well exceeded by a1 and a2, and likely by a3 too; all different tracks in the left of Fig. 17 point to an initial mass of ≥250 M_⊙ for R136a1, the most massive star in our sample, and, conservatively, ≥150 M_⊙ for R136a2 and R136a3. Of course, these results might not hold if the sources turn out to be in multiples, something that is currently being investigated with radial velocity measurements using HST-observations (Shenar et al. 2019). Furthermore, we note that our masses not only rely on the adequate determination of effective temperatures, and the used evolutionary tracks, but also on the flux calibration and reddening of the anchor magnitude used for our analysis. In this context, we note that Rubio-Díez et al. (2017) focus in particular on the infrared (K-band) flux calibration of R136a1, R136a2 and R136a3, and find considerably lower initial masses for these stars compared to our analysis and that of Crowther et al. (2010), their highest derived initial mass being an upper limit of 194 M_⊙ for R136a1. For their analysis, Rubio-Díez et al. (2017) use VLT/SINFONI K-band spectrophotometry, effective temperatures of Crowther et al. (2010) and evolutionary tracks of Köhler et al. (2015).

We measure the current masses of the WNh stars in a second manner. By measuring the surface gravity, we can obtain the spectroscopic mass $M_{\text{spec}}(=g_e{R}_{*}^2/G,\text{\hspace{0.17em}with\hspace{0.17em}}{g}_e=g+{\left({\upsilon }_{\text{eq}}\sin i\right)}^2/R_{*}$ surface gravity corrected for centrifugal accelerations, and G the gravitational constant). Contrary to the analysis of Bestenlehner et al. (2020), in our fits log g is a free parameter also for the WNh stars. The winds of the WNh stars are so dense that we do not expect to see a very strong signature of the surface gravity in the spectrum. Still, from our Kiwi-GA fits we do constrain log g of both R136a1 and R136a2, for which we find a 2σ range of log g = 3.35–3.9 and log g = 3.55–3.75, respectively. For R136a3 we only find a lower limit, log g > 3.4. Figure 18 shows that the fitness distribution of the gravity of R136a1 clearly favours a value lower than 4.0. When using the measured log g values in combination with the derived radii to derive spectroscopic masses, these compare well with the evolutionary masses (see Table 4), further supporting very high masses for these stars.

5.3.2 CNO abundances

As CNO surface abundances are expected to change over the course of stellar evolution (e.g. Köhler et al. 2015; Eggenberger et al. 2021), they could be used to set apart the more evolved stars from the rest. Before we assess this for our sample, we check whether the derived abundances are consistent with the theory of CNO-processing. This is an especially important check because without exception the C and N abundances are derived from lines that are (mostly) formed in the stellar wind, and thus their strength and shape not only depends on abundance and temperature, but also on a handful of wind properties. Since all those properties can vary independently from each other, consistency with CNO-processing theory is not guaranteed intrinsically and needs to be checked. The diagram in Fig. 19 allows for such a consistency check, by comparing the ratios n_N/n_C to n_N/n_O. Maeder et al. (2014) derive analytically the limits of these ratios given CNO and CN equilibrium: we expect all observations to fall somewhere between the region bordered by those. Furthermore, evolutionary tracks of different masses and initial rotation (Brott et al. 2011a; Köhler et al. 2015) predict quite a narrow range in which we expect points to lie. This method has also been applied by Martins et al. (2015) and Carneiro et al. (2019).

The WNh stars – the only sources for which we measure oxygen abundance – lie in the area predicted by Maeder et al. (2014), and moreover match well the evolutionary predictions depicted in Fig. 19. For the other sources, we do not measure oxygen abundance and we therefore assumed for all stars the LMC baseline value, log n_O/n_H + 12 = 8.35 (Kurt & Dufour 1998, as in Brott et al. 2011a; Köhler et al. 2015), so no enrichment or depletion. With this assumption and considering the uncertainties on the measurements, our observations are generally consistent with CNO-processing, except for R136a7. For this source evolutionary models predict oxygen depletion to n_O = 7.20, assuming an initial mass of 100 M_⊙ and an initial rotation of 300 km s⁻¹. This would move this source into the region consistent with CNO-processing. If we also assume oxygen depletion for other sources, as was done in the optical + UV fitting²², the shift along the log n_N/n_O axis is small (0.1–0.2 dex to the right). We note that there is an uncertainty in the observed LMC baseline abundances of CNO (see, e.g. Dopita et al. 2019) that could affect Fig. 19. However, we find that changes in the diagram that could result from this are smaller than the typical error bars on our observations. Overall, we conclude that no violation of CNO-processing is observed and that abundances and wind parameters can be disentangled from a set of 9–11 (metal) wind lines, albeit with large uncertainties.

Six sources stand out in Fig. 19. These are the WNh stars, R136a5, R136b, and H30. Bestenlehner et al. (2020) show, based on the helium abundance, that the enrichment of these stars is mainly driven by mass loss (see also Bestenlehner et al. 2014). Here, we further investigate the nitrogen enrichment by placing our sources in a Hunter diagram. This diagram, introduced by Hunter et al. 2008, shows the nitrogen surface abundance of stars versus their projected rotational velocity. By comparing these quantities, one can gain insight into mixing processes that occur within the star; rotational mixing being one of the mechanisms that can bring processed elements to the surface. In the Hunter diagrams in Fig. 20 we compare our observed values with evolutionary tracks of different rotation rates (Köhler et al. 2015). In order to account for the different masses of the sources, we compare three subgroups of sources with tracks of three different masses.

All WNh stars (M_ini = 213–273 M_⊙, top panel of Fig. 20) show strong nitrogen enrichment; we find a similar nitrogen surface abundance and age for all three. Given the large uncertainty on the rotation rate, the nitrogen abundance and age of R136a1 match the tracks quite well. For R136a2 and R136a3 the obtained nitrogen abundances are, given the obtained ages, slightly too high for any of the tracks. However, the difference is not large and within uncertainties the single star models and observations match. We emphasise that, while the best fit values for the υ_eq sin i are prominently marked (based on the best fitting model of the optical only run for each star), one should not overlook the large error bars. In the most extreme case of R136a3 these span all displayed υ_eq sin i values, implying that we cannot put any constraints on the current υ_eq sin i of this source. Adding to the uncertainty (not captured in the statistical error bars shown here) is the fact that we derive all rotation rates by convolution of the line profiles, hereby implicitly assuming that the emergent radiation is emitted from one single rotating layer. This assumption likely breaks down for all available spectral lines in the WNh spectra, casting further doubt on the validity of the WNh rotation rates that we derive. A more sophisticated approach would be to include the effects of rotation on the velocity field of the wind into the formal integral (Shenar et al. 2014), however this is not within the current capabilities of Fastwind. Regardless of the observed υ_eq sin i values, a high initial rotation is suggested for all WNh stars if one compares the tracks to the observed ages and nitrogen abundances. This was also found from the BONNSAI runs for these sources based on a comparison of the observed luminosity, temperature, surface gravity and helium surface abundance to the Köhler et al. (2015) tracks.

The middle panel of Fig. 20 shows stars in the mass range M_ini = 92–127 M_⊙. The age and position of the sources generally show good agreement with the evolutionary tracks. The three points that lie close together (R136a4, R136a6 and H36), all seem to have started out with a moderate initial rotation rate. R136b is highly nitrogen enriched, which suggests that the star must have had an initial rotation in the range of 300–450 km s⁻¹ and is older than the other stars. Indeed, the age of the star, derived from the observed luminosity, temperature, surface gravity and helium abundance, is 2 Myr. Looking at the physical position of this star in the cluster (Fig. 1) we see that it is located a bit on the outskirts of the cluster, further away from the centre than most other very massive stars (Fig. 1). Possibly this could be related to its somewhat higher age, however, we do not have any evidence for this; Bestenlehner et al. (2020) considers the ages and positions of all sources and does not find a correlation of age and position (we confirm this finding, see Appendix E). The outlier in this panel is R136a5, that, within its 1σ errors, does not seem to fall on any of those tracks. We note that this star has a slightly higher initial mass than that of the tracks shown here (116 M_⊙). Tracks with higher mass would have more enrichment (see top panel), and thus would bring the tracks and the observations of R136a5 closer together. R136a7 could also be considered an outlier as it is the only source in this mass-range not showing significant enrichment. Yet, it is consistent with the tracks, as, with an age of $0.5{\pm }_{0.46}^{0.37}$ Myr, R136a7 is one of the youngest, if not the youngest, star in the sample²³.

The bottom panel of Fig. 20 shows stars in the mass range M_ini = 52–68 M_⊙. Most of them are not yet nitrogen enriched. For a few we do observe a slight enrichment, though with very large uncertainties. Within errors, the observations are consistent with the 60 M_⊙ tracks of Köhler et al. (2015).

It is interesting that while for the WNh stars we require high initial surface rotation rates (υ_eq,ini > 300 km s⁻¹) in order to match the tracks, this is not the case for most O-stars. The Ostars in Fig. 20, except R136a7 and R136b, suggest initial surface rotation rates υ_eq,ini < 300 km s⁻¹. The inferred initial rotation rates are, for all stars, consistent with the values inferred using BONNSAI based on different observables (luminosity, temperature, surface gravity and helium abundance); with this we find 100 km s⁻¹ for all with the exception of the WNh stars, R136a7 and R136b, for which we find υ_eq,ini > 300 km s⁻¹. The difference between the WNh stars and R136b and R136a7 on the one hand, and the other O-stars on the other hand might indicate that the two groups of stars are formed through a different channel. We note that only for sources for which we observed helium enrichment (x_He ≥ 0.14), a high initial rotation rate was derived using BONNSAI.

In summary, the collection of Hunter diagrams shows the enrichment of young (very) massive stars. At an age of ~1.5 Myr the ~60 M_θ stars are barely nitrogen enriched, while the ~100 M_θ and ~200 M_⊙ stars show enrichment of about one and two orders of magnitude, respectively. This is roughly consistent with the single star models of Köhler et al. (2015). We note that Hunter et al. (2008); Brott et al. (2011b); Grin et al. (2017) report a population of slowly spinning nitrogen enriched stars, which we do not identify in our study. This may potentially point to binary interaction as a source of such stars; given the young age of our population it could be expected that such interactions have not yet occurred frequently.

5.4 Terminal velocity measurements

With spectral fitting we are able to break the degeneracy between υ_∞ and υ_windturb. While based on the edge velocity alone one cannot disentangle υ_∞ and υ_windturb, the shape of the absorption component of C iv λλ1548–1551 is affected differently by the two parameters, and we are able to distinguish between the two. We note that narrow absorption components, unresolved in our spectra, may have contributed to absorption near the terminal velocity of non-saturated line profiles (Prinja et al. 1990). As the equivalent width of these absorption components is small, we do not expect a significant effect on our measurements.

The terminal velocities we find from spectral fitting are systematically larger than those of Crowther et al. (2016), who used the exact same dataset but obtained v∞ by inspecting the blue wing of the C iv λλ1548–1551 P Cygni line (Fig. 21). For 29 sources they identify the maximum blueward extent of (near) saturated absorption profiles, υ_black, and assume υ_∞ = υ_black (following Prinja et al. 1990). For these 29 sources they also identify the velocity at which the violet absorption meets the continuum, υ_edge, and with this derive υ_∞ = 0.8υ_edge. With the latter relation they estimate υ_∞ by identifying υ_edge for the remaining 15 sources, that have less strong wind lines.

Figure 21 shows the υ_∞ values derived from spectral fitting versus those found from υ_black or υ_edge, which we call υ_∞,blue. We find an average difference of 283 ± 30 km s⁻¹ between the two methods, taking into account uncertainties on the measurements, and excluding three sources for which we could not obtain a good fit to the C iv λλ1548–1551 profile (see Sect. 4.3). Moreover, a trend is visible in the difference: on average the difference between υ_∞ derived with the two methods increases for lower terminal velocities. Figure 22 shows C iv λλ1548–1551 profiles for three sources with different mass-loss rates and terminal velocities. The location of v_∞ according to this work and that of Crowther et al. (2016) is indicated. According to the spectral models, the terminal velocity does not coincide with the point of strongest absorption, but lies more towards the blue. Some desaturation thus occurs close to the terminal velocity already in H48, the star with the strongest P Cygni profile of these three. If the approach of Crowther et al. (2016) is followed for obtaining υ_∞, the obtained values υ_∞,blue could be corrected for this effect by using the relation shown as the blue solid line in Fig. 21: $${\upsilon }_{\infty }={\upsilon }_{\infty ,\text{blue}}+0.27\times (3700-{\upsilon }_{\infty ,\text{blue}}),$$(16)

(velocities in km s⁻¹) for 1000 km s⁻¹ > υ_∞,blue > 3500 km s⁻¹. One could use this equation for getting a more accurate first order estimate of υ_∞ from reading off υ_blue from spectra with relatively poor resolution.

Upon comparing the edge velocities of Crowther et al. (2016), obtained by reading off the velocity at the wavelength where the blue edge of C iv λλ1548–1551 meets the continuum, to the edge velocities we obtain by spectral modelling and assuming υ_edge = υ_∞ + υ_windturb, we find that they are generally in good agreement. For our sample we find the following average: $$\begin{array}{cc}{\upsilon }_{\text{windturb}}=0.14\pm 0.06{\upsilon }_{\infty }& \text{or,\hspace{0.17em}equivalently:}\\ {\upsilon }_{\infty }=0.88\pm 0.04{\upsilon }_{\text{edge}}.& \end{array}$$(17)

Had Crowther et al. (2016) assumed this value instead of υ_∞ = 0.8υ_edge, our respective values for υ_∞ would lie closer together for the weaker wind sources (triangles in Fig. 21), although part of the discrepancy would still remain. We stress that Eq. (17) is based on fits of the sources for which we carried out 12-free parameter fits, that is, the sources with stronger winds. Data of higher S/N and resolution is required to disentangle υ_∞ and υ_windturb for the sources with weaker winds.

6 Summary and outlook

We have simultaneously analysed optical and UV spectroscopy of a population of 56 stars in the core of the R136 star cluster, nine members of which have masses M ≳ 100 M_⊙. For the first time we investigate the wind structure parameters of a large range of spectral types while fitting the interclump density, the wind turbulence, and the effects of optically thick clumps such as velocity-porosity. By taking these effects into account, we have improved the accuracy of mass-loss determinations of the most massive stars. Moreover, the derived mass-loss rates are no longer affected by the well-known mass-loss and clumping dichotomy, but they are actual values. Our main findings are the following:

• The HRD positions of the sources suggest a cluster age of 1–2.5 Myr, which is in line with the findings of Bestenlehner et al. (2020). The ages of the highly nitrogen enriched WNh stars are in line with the age of the rest of the population.

• Our conservative estimate for the initial mass of R136a1, which is the most massive star in our sample, is 250 M_⊙. The initial masses of R136a2 and R136a3 well exceed 150 M_⊙. The spectroscopic masses of these sources, which we measure here for the first time, further support this conclusion.

• We compared the theoretical predictions of Vink et al. (2001), Krtička & Kubát (2018), Björklund et al. (2021), and Vink & Sander (2021) to the observed mass-loss rates and terminal velocities and found that none of the predictions satisfactorily reproduce both quantities. The largest discrepancies for the terminal velocities were found for stars with log L/L_⊙ ≲ 5.3 for the Björklund et al. (2021) predictions, and stars with log L/L_⊙ ≲ 5.3 for the Krtička & Kubát (2018) predictions.

• Overall, the mass-loss recipe of Krtička & Kubát (2018) best matches the observed mass-loss rates of the stars in our sample. The predictions of Björklund et al. (2021) match almost as nicely, performing better in the low-luminosity regime (log L/L_⊙ ≲ 5.3), but worse for the higher luminosities. The prescriptions of Vink et al. (2001) and Vink & Sander (2021) overpredict the mass-loss rates for all luminosity regimes.

• The stellar winds of the stars in our sample are highly clumped, with an average clumping factor of f_cl = 29 ± 15.

• We find tentative trends in the wind structure parameters as a function of mass loss, where the stars with the highest mass-loss rates seem to have smoother, albeit still clumpy, winds (but see below).

• We provide a prescription for the mass-loss rates of the most massive stars as a function of the electron scattering Eddington parameter, following the work of Bestenlehner (2020). For this, we used our best fit values of the CAK force multiplier parameter α_eff = α − δ = 0.46 ± 0.04 and the transition mass-loss rate log ${\dot{M}}_{\text{e,trans}}=-5.19\pm 0.10$.

• The point with the strongest absorption in a P-Cygni profile does not typically correspond with the terminal velocity, which lies more bluewards. We provide an equation that quantifies this effect.

This is the first investigation of the trends in the wind structure parameters of massive stars. While the measurements of Hawcroft et al. (2021) are not in contradiction with our results, they cannot confirm the trends we observed either, given the limited span of mass-loss rates in their sample. Further investigation of these trends is thus necessary. Hawcroft (in prep.) are undertaking such a study, analysing optical and UV spectroscopy of a sample of about 30 LMC and SMC stars covering most O-type subclasses (O3–O9). Furthermore, our study can be considered a pilot for a larger investigation making use of the high-quality HST UV spectra of the ULLYSES project²⁴ (Roman-Duval et al. 2020). The ULLYSES sample, when complete, will consist of ~250 massive stars (mostly in the Magellanic Clouds), including ~150 O-stars, covering all O-star subtypes and luminosity classes. Complemented by the optical XshootU program²⁵, ULLYSES will provide an excellent opportunity for further study of the structure of massive star winds.

Acknowledgements

We like to thank the anonymous referee, whose constructive criticism has helped us to improve the presentation of this work. This publication is part of the project ‘Massive stars in low-metallicity environments: the progenitors of massive black holes’ with project number OND1362707 of the research TOP-programme, which is (partly) financed by the Dutch Research Council (NWO). Observations were taken with the NASA/ESA HST, obtained from the data archive at the Space Telescope Institute. Computations were carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. S.A.B. would like to thank the late Lykle Voort for his help and patience. F.A.D. and J.O.S. acknowledge support from the Odysseus program of the Belgian Research Foundation Flanders (FWO) under grant G0H9218N. T.S. acknowledges support from the European Union’s Horizon 2020 under the Marie Sklodowska-Curie grant agreement No 101024605. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 945806). This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster). G.G. acknowledges support from Deutsche Luft- und Raumfahrt (DLR) Grant No. 50 OR2009.

Appendix A Best fit parameters

The best fit parameters and 1σ error bars for the optical + UV fits can be found in Table A.1 to A.3.

Appendix B Data quality

B.1 S/N

Distributions of S/N for all sources in our sample are shown in Fig. B.1. We see that generally the UV spectra have the highest S/N, but we note that typically C iv λ1165-CIII λ1170, which is very close to the blue grating edge and is also slightly affected by Ly-α absorption, has a lower S/N than the values shown in the figure. Also, in the most red part of the UV grating, He ii λ11640 and N iv λ1718 (if used at all) are typically a bit more noisy than the other lines.

B.2 Wavelength correction

The radial velocity shifts inferred for all stars, for each grating, are shown in Fig. B.2. If our radial velocity measurements were accurate, we would expect that they have a mean velocity of 267.7 ± 25 km s⁻¹, as measured by Hénault-Brunet et al. (2012), who analysed high-quality spectra of stars around the R136 core. Due to the calibration related wavelength offset however, the dispersion we find is expected to be higher. Indeed, the standard deviation on the radial velocities is of the order of ±2 pixels for most gratings. The exception is Hα, where the typical offset is 3 to 4 pixels.

Appendix C Technical details of Kiwi-GA

We now discuss how KIWI-GA differs from the algorithm of Mokiem et al. (2005) and Abdul-Masih et al. (2021), who follow the approach of the former. First of all, KIWI-GA applies operations of recombination (mixing of parameters of two models) and mutation (addition of random variations to a subset of the parameters) directly on the model parameters. Mokiem et al. (2005) store the parameter values in the form of one string per model, where each character of the string can have a value from 0… 9 and different parts of the string indirectly represent the values of all parameters. On these strings the mutation and recombination operations are applied. The fact that this concept is abandoned in KIWI-GA has consequences for the way that mutation and recombination are implemented. With this in mind, we list here the most significant changes of our algorithm compared to that of Mokiem et al. (2005):

• Recombination is carried out per parameter. Given two sets of parameters (parents), the new parameters (offspring) are chosen by, for each parameter, picking the value of one of the parents. The genome can thus be split in more locations than only one or two, but not in the middle of a parameter. This implies that mutation due to recombination can never occur.

• No parents are cloned, that is, recombination always occurs²⁶. All new models are thus the product of both recombination and mutation, and none of mutation alone. This increases the diversity of the population.

• Two types of mutation are introduced, each having a different rate of occurrence. In each case the size of the mutation, or the amount with which the parameter value changes, follows a Gaussian distribution. Small mutations, where the value of the mutated parameter lies close to the original value, occur frequently: after recombination, each parameter has a large chance to undergo a small mutation. Larger mutations occur with a lower frequency. The exact frequencies can be set by the user, in this work we assumed a rate of 0.5 for the small mutations, and a maximum rate of 1/n_free for the larger mutations. The latter means that on average one parameter per model undergoes a large mutation. This manner of mutating parameters is very different from the scheme used in Mokiem et al. (2005, see Appendix C.1 for details).

• Reinsertion of new individuals into the parent population is done according to an elitist-fitness scheme, where the parent population can be larger than the amount of individuals in one generation. This means that after each generation, the fitness of both parent and offspring models is assessed and the top models of each sub-population will form the new parent population. An elitist-fitness scheme leads to a fitter parent population, while diversity is ensured by always including a certain fraction of offspring models²⁷.

Together, these adaptions result in a faster parameter space exploration, that is, less generations need to be computed.

Especially the exploration around the best fit, done in order to assess the uncertainty regions, benefits from the algorithm adaptions.

Apart from these major points there are other additions, such as the fact that H-I, He-I and He-II ionising fluxes and uncertainties thereof can be given as an output, the possibility to make the first generation larger compared to the rest of the generations²⁸, the option to treat X-ray parameters as free parameters, and the option to estimate an appropriate volume filling fraction ƒ_X for each model such that L_X/L_bol ≈ 10⁻⁷ (see Appendix G.2). Furthermore, KIWI-GA can be used in combination with Fastwind version 11, which can treat radiative transfer of the full spectrum in the co-moving frame (Puls et al. 2020).

C.1 Mutation in Kiwi-GA

Figure B.3 shows several mutation schemes implemented in the algorithms of Mokiem et al. (2005) and Abdul-Masih et al. (2021), next to the mutation scheme used in KIWI-GA. The scheme used in KIWI-GA results in a mutation distribution that covers the parameter space more regularly. In practice, this gives the user more control over the ratio of small to large mutations.

On the one hand, parameters close to fit models (created by small mutations) are expected to be more successful (because the original model was already fit), but on the other hand parameters very different from the fit models can cause larger improvement in fitness. This is the reason that we chose for small mutations with a high probability, and large mutations with a low probability.

Optimising hyper parameters of a genetic algorithm, such as mutation rate, is not trivial. An extensive study costs a lot of resources, so we limited ourselves to educated guesses and a few test runs. Nonetheless, the choices we have made seem to work well in practice, when the parameter space is well behaved. In the end, the goal is to map for each parameter the envelope of the fitness distribution (i.e. as a function of each free parameter, find the lowest χ²), but without computing every model in the parameter space. In the grid-fitting approach each model in the parameter space is computed, and one can be certain that the envelope of the fitness distribution is completely mapped. With a genetic algorithm, the mapping is only approximate; one might miss certain models in the parameter space. In Fig. B.4 we compare fitness distributions of KIWI-GA to those of pyga, using runs with identical setup. The KIWI-GA distributions are smoother (have less gaps), meaning that that set of models provide a closer representation of the true fitness envelope. We see that both algorithms find very similar best fit solutions, that are in both cases close to the true values, that is, the ones that the original model have (mind that deviations from this could be caused by the simulated data: a S/N of 30 was imposed on the synthetic spectra). In practice, this thus means that the algorithms work the same and give similar outputs, but that with KIWI-GA one needs to compute less models, that is, needs a lower amount of computation time to trace the envelope of the fitness distribution. We suspect that the reason for this lies mostly in the mutation scheme, but also other changes in the algorithm could have contributed to the improved performance. We note that KIWI-GA has the option to use the mutation operators used in the algorithms of Charbonneau (1995); Mokiem et al. (2005, demonstrated in Fig. B.3). Furthermore, we note that, regardless of the genetic algorithm or mutation operators used, the amount of generations must be increased with higher quality data.

Appendix D Comparison of methods

In this section we compare the results of different methods. We stress that the error bars have to be taken into account when comparing the different methods; with a few exceptions, the agreement is good. To illustrate typical differences in the best fit models of the different methods, we show in Figure D.1 a selection of optical lines for two stars: one star where there is some discrepancy between parameters, another star for which the agreement between parameters is good; the latter being representative for the sample.

We derive the parameters T_eff, log g, $\dot{M}$, log L/L_⊙, and in some cases also β and nitrogen abundance, in two different manners: the optical-only and optical + UV runs. Additionally, we obtain a third measure for T_eff by fitting cmfgen models to the iron pseudo-continuum in the UV. Lastly, Bestenlehner et al. (2020) have already carried out an optical-only analysis of the same spectra, but using a different method. For the O-stars, Bestenlehner et al. (2020) use the IACOB-GBAT tool Simón-Díaz et al. (2011) instead of Fastwind and KIWI-GA, while for the WNh stars, they use cmfgen. For a few stars, values where obtained both with Fastwind and cmfgen, in which case we adopt, in this section, the result of their Fastwind analysis. Furthermore, in our analysis nitrogen abundance was a free parameter that could range from 6.9 to 10.0, while Bestenlehner et al. (2020) use only three grid values, namely 7.1, 8.2 and 8.5. In this section we compare the outcomes of the different analyses.

D.1 Temperatures

Figure D.2 shows a comparison of the three T_eff measurements from this paper to each other and to those of the optical-only analysis of Bestenlehner et al. (2020). For the hotter stars (T_eff ≳ 45000 K) we find systematically higher temperatures from our optical-only analysis than does Bestenlehner et al. (2020, top left panel), though within errors also there the agreement between the methods is good. Inspecting the top right and lower left panels, we see that these higher temperatures we find from our optical-only analysis are not found from the optical + UV fit. Assuming that the optical + UV analysis is the most reliable method because it takes into account the most spectral information, the temperatures of Bestenlehner et al. (2020) seem to be more reliable than our optical-only values, although the difference is small. It may be the case that the nitrogen lines do not affect strongly the goodness of fit that KIWI-GA uses for fitting, while this should be the case. Namely, for measuring the temperature for stars hotter than ~ 45000 K one relies mostly on (weak) nitrogen lines, as He-I/He-II ratio can no longer be used since all He-I is gone at those temperatures. Upon inspecting our fits and especially the fits to the nitrogen lines, we do not see cases where the nitrogen lines have an especially bad fit. The reason for the slight, albeit significant overestimation of our optical-only temperatures is thus not understood.

The bottom right panel of Fig. D.2 contains a comparison of our optical + UV temperatures and the temperatures we get from fitting the iron lines with the cmfgen grid. Clearly, the iron (Fe) lines show a strong temperature dependence that is consistent with the temperatures of the H-He-C-N-O lines. At the lower end the Fe lines indicate a systematically higher temperature. This can be explained by the fact that the lowest temperature of the cmfgen grid is 35481 K. At the higher temperature end the points are scattered. This could be due to the fact that in the wavelength range from 1600 – 1700 Å, where most of the Fe iv lines lie, the S/N of our spectra is poor. Lastly, we note that both Fe lines as well as H-He-C-N-O lines are sensitive to micro-turbulence, but that a fixed value for micro-turbulence was assumed in both the cmfgen grid and in the optical + UV runs (10 km s⁻¹ in both cases).

D.2 Mass-loss rates & β

We see a systematic offset offset between the mass-loss rates obtained by our optical-only analysis, and that of Bestenlehner et al. (2020). This can be explained by the different assumptions and or fit-values for the wind acceleration parameter β. For most sources, in the optical-only analysis we assume β = 0.9, while Bestenlehner et al. (2020) leaves this parameter free, and finds for most sources a value of 1.0 (16 sources) or 1.2 (28 sources)²⁹. When basing the mass loss only on optical lines without strong emission features, the parameter β is degenerate with mass-loss rate, where a lower mass-loss rate or a lower β are having the same effect on the wind line. This explains the observed discrepancy in measured mass-loss rates: we find or assume lower values for β, resulting in higher values for mass-loss rate, compared to Bestenlehner et al. (2020). For deducing additional uncertainties on $\dot{M}$ related to lack of knowledge of the value of β, one should vary β and assess for each value the corresponding $\dot{M}$.

Markova et al. (2004) and Repolust et al. (2004) carry out such an analysis and find typical uncertainties of 0.1 – 0.3 dex.

The discrepancy between the optical-only and optical + UV measurements of the mass-loss rate is related predominantly to clumping. We investigate whether the observed differences is in line with expectations by using the individual f_cl values obtained from our optical + UV fits to do a clumping correction for the mass-loss rates from our optical-only fits, where we assumed a smooth wind, that is, f_cl = 1. In other words, we check whether $${\dot{M}}_{\text{opt+UV}}={\dot{M}}_{\text{opt}}/\sqrt{f_{\text{cl,opt+UV}}}$$(D.1)

is satisfied³⁰. In theory, the density of the ρ²-sensitive optical recombination lines from which the optical-only mass-loss rates are derived should scale with $\sqrt{f_{\text{cl}}}$ and the above should hold, assuming the adopted clumping factor is spatially constant. Upon applying the clumping correction to ${\dot{M}}_{\text{opt}}$ we find that on average the shifted values are too low compared to ${\dot{M}}_{\text{opt+UV}}$. Closer inspection reveals that this discrepancy is related to different values of β that are assumed or derived for the optical-only and the optical + UV runs. Again, the degeneracy of β and $\dot{M}$ plays a role: For sources where β_opt < β_opt+UV the correction with $\sqrt{f_{\text{cl,opt+UV}}}$ is too large, while for sources where β_opt ≥ β_opt+UV Eq. D.1 is satisfied (see Fig. D.3). The latter group includes the sources where Hα is in emission and β could be measured from the optical-only run. Lastly, we note that after taking into account the different values of β there remains a group where the correction with $\sqrt{f_{\text{cl,opt+UV}}}$ is not large enough. An explanation for this could be that the mass-loss rates found from the optical-only analysis are higher than the true mass-loss rates. For these sources we probably derive an upper limit rather than the actual mass-loss rate from the optical-only data.

D.3 Other parameters

Figure D.4 shows a comparison of the optical-only analysis of this paper to the optical-UV analysis, and the analysis of Besten-lehner et al. (2020). The most striking differences can be seen in the nitrogen abundance. The match of the Bestenlehner et al. (2020) analysis with our optical-only analysis is not particularly good, but it is important to note that Bestenlehner et al. (2020) stress that the nitrogen abundance value that they provide is an indication rather than a precise measurement³¹ and that the error bars on our optical-only analysis are very large. Furthermore, upon closer inspection of individual sources we see that agreement is good for all six sources that have significant overabundance (R136a1, R136a2, R136a3, R136a4, R136a5, R136b, H36, H48). This is also the case for these sources if one compares the optical-only versus optical + UV analysis, though there is a slight systematic towards higher abundances from the optical-only analysis.

For log g, log L/L_⊙, υ_eq sin i, x_He there is generally good agreement between the different runs. The exception are a few υ_eq sin i values of our analysis that seem to be very high, however this concerns a few very low S/N stars for which υ_eq sin i is essentially unconstrained. Another point that shows a clear mismatch is the one with the highest He abundance in the bottom plot: this is a3, where we from the optical only analysis found a temperature of 42750 K where with the UV data we assessed that 50000 K is likely closer to the true value, given the NV and OV lines in the UV. With a higher temperature, we find a lower He abundance.

The above shows that, even if the same code is used for the analysis, the outcome (at least the best fit parameters) can strongly depend on the analysis method, at least when the data quality is not superb.

Appendix E Stellar masses and ages

The power law slope of the initial mass function that we derive, γ = 1.99 ± 0.11 in the range 30 – 200 M_⊙, is consistent with the value of $1.90{\pm }_{0.26}^{0.37}$ derived for massive stars in the wider 30 Doradus region (Schneider et al. 2018b). The sample of Schneider et al. (2018b) excludes the core of R136 and the most massive star in their sample is $203{\pm }_{44}^{40}$ M_⊙. The highest initial mass that we derive, $273{\pm }_{36}^{25}$ M_⊙ for R136a1, is consistent with the stochastic sampling analysis of Schneider et al. (2018b), that excludes maximum initial masses of over 500 M_⊙ in 30 Doradus with 90% certainty. We note that for our BONNSAI runs we used the Salpeter (1955) initial mass-function, that, with a slope of 2.35, prefers lower masses. This can impact the slope that we derive, albeit in a conservative way, that is, without this prior we would have found an even shallower slope. Regarding the ages; while we derive a cluster age of 1 – 2.5 Myr, there is a tail of older stars, suggesting that the stars in R136 are not coeval; this could be either the result of a projected view onto this dense central region in 30 Doradus, or is related to the formation of R136 (see also Crowther et al. 2016; Schneider et al. 2018a; Bestenlehner 2020).

A comparison of the evolutionary and spectroscopic masses (see also table Table I.1) can be found in Fig. E.1. The masses agree well within errors, but for M ≳ 70M_⊙ the spectroscopic masses are, with a few exceptions, systematically higher than the evolutionary masses. Bestenlehner et al. (2020) find the same (their Figure 3), although in some cases (e.g. R136a5, R136b) the current analysis gives larger values. Furthermore, we note that in addition we derive spectroscopic masses for the WNh stars R136a1-a3. For an in depth discussion of the mass discrepancy, see Bestenlehner et al. (2020).

The age distribution we derive for R136 based on the BONNSAI values is shown in Fig. E.2. We find a strong peak at 1.1–1.3 Myr, consistent with the results of Bestenlehner et al. (2020). A second peak peak may be suggestive of a second burst of star formation around 2.0 Myr, though given the uncertainties the evidence for this is not strong.

Figure E.3 shows the initial mass distribution we derive from the BONNSAI values. We fit a power law in the form ξ(M) ∝ M^−γ to the initial mass distribution in the mass range 30–200M_⊙ and find a slope of γ = 1.99 ± 0.11. This is consistent with the results of Bestenlehner et al. (2020) of the R136 core (γ = 2.0 ± 0.3), and the findings of Schneider et al. (2018b, $\gamma =1.90{\pm }_{0.26}^{0.37}$) for the rest of the Tarantula Nebula, but steeper than the standard Salpeter initial mass function.

Lastly, we compare the age of our sources with the spatial distribution (Fig. E.4). While many young sources are found close to the centre, several of them are also found in the outskirts. Moreover, the old sources are scattered throughout the cluster. The fact that we do not observe them in the very centre might be an observational bias: the very bright WNh stars could dominate the older, lower mass stars in that region. We thus conclude that we do not see a strong correlation between the spatial position of the stars and their age. This is in line with what Bestenlehner et al. (2020) conclude (their Figure 11).

Appendix F Notes on individual sources

F.1 R136a1

Possibly contaminated by H17, which lies north of R136a1 with an angular separation of 75 mas and a flux ratio of 0.112 in the K-band (Khorrami et al. 2017). Due to the low flux ratio and the fact that H17 lies mostly outside the slit we assume the spectral analysis of R136a1 is not affected by H17.

For this star O v λ1371 is poorly reproduced – something that is seen also for the other WNh stars of the sample (R136a2, R136a3). In particular, the line is not broad enough compared to the data, implying that there is not enough O v in the outer wind. We explored whether this could be related to our (canonical) X-ray assumptions by probing the X-ray parameter space, and ruled this out as a cause. Another possibility is that something might go wrong with the Fastwind treatment particularly for very dense winds, either due to indirect line overlap effects in the outer wind, or something else. Lastly, it could be that the wind structure that is assumed (i.e. the clumping stratification) is not representative for what actually happens in the wind, and it could be that this discrepancy especially shows, or, is largest for, the WNh stars. The exact cause of the issue is currently unidentified. We note that O v λ1371 is sensitive to clumping, however we do not believe the poor fit of O v λ1371 results in spurious clumping factors for the WNh stars. Namely, for this particular line a higher clumping factor leads to a weaker profile which is the opposite behaviour of some other of the clumping factor sensitive lines (e.g. C iv λλ1548–1551). In the WNh stars we find rather high clumping factors, where a better fit of O v λ1371 (deeper profile) would be obtained with a lower clumping factor.

It thus seems that other clumping sensitive lines dominate the value that we obtain.

F.2 R136a2

O v λ1371 was poorly reproduced; see explanation in Appendix F.1.

F.3 R136a3

The best fit of the optical + UV run of R136a3 results in a T_eff of 42 kK, as opposed to 50 kK found by Bestenlehner et al. (2020) and 53 kK found by Crowther et al. (2010). Upon inspecting the fit we see that while in general the fit is matching the data well, this does not hold for the higher ionisation lines C iv λ1165, O v λ1371, and N v λλ4604–1620, suggesting the derived T_eff of 42 kK is too low. We carried out another fit where we fixed T_eff to 50 kK and it is important to note that the fits to O v λ1371 and N v λλ4604–4620 improved considerably, but the fit to many other lines got worse (Fig. F.1). The mass-loss rate and helium abundance we find from the two runs differ significantly, we find log $\dot{M}=-4.2{\pm }_{0.04}^{0.07}$ and log $\dot{M}=-4.5{\pm }_{0.07}^{0.04}$, and $x_{\text{He}}=0.49{\pm }_{0.11}^{0.01}$ and $x_{\text{He}}=0.37{\pm }_{0.10}^{0.10}$ for the low and high temperature, respectively. Because there is no way to reproduce the strong O v λ1371 and N v λλ4604–4620 lines with a lower temperature, and the fact that when assuming a low effective temperature R136a3 lies far away from the other two WNh stars in the HRD, we deem the higher temperature more likely. Further support for that comes from the fit to the iron lines in the UV presented in Fig. F.2, resulting in a temperature of around 47000 – 50000 K. For these reasons, we adopt a temperature of T_eff = 50kK for R136a3. The other stellar parameters we obtain from the optical + UV run with T_eff fixed to that value (blue model in Fig. F.1). For the υ_eq sin i we stick to the value of the optical-only run (best fit: T_eff = 42750 K) as we do for the other stars.

We note that Khorrami et al. (2017) resolve R136a3 into two sources, however the flux ratio is 0.124 in the J-band (0.044 in the K-band), so in practice one source will dominate our optical and UV spectroscopy. We do not see evidence for unresolved binarity of R136a3 where the source would consist of a hot and a cooler component, at least not directly: no He i lines are visible.

Furthermore, we note that, while we adapted the higher T_eff for R136a3 for the rest of our analysis for reasons explained in this section, we have to keep in mind that the low values we obtained from our original fit may point to the limits of the applicability of this Fastwind version, which was not designed for situations where, for example, Hγ is in strong emission. Lastly, we note that, also in the higher T_eff fit, O v λ1371 was poorly reproduced; see explanation in Appendix F.1.

F.4 R136a6

R136a6 consists of two stars with an approximately equal flux contribution (see also Sect. 2.1). We carry out the spectral analysis but exclude it from the analyses of the sample as a whole regarding mass-loss and clumping properties.

F.5 R136a8

No optical spectra are available and thus they are not included in the analysis of Bestenlehner et al. (2020). The optical + UV run of this source was thus in fact carried out on UV data only. The setup was the same as for the other sources, but instead of fixing the υ_eq sin i and helium abundance on values derived of the optical, we fixed them at 150 km s⁻¹ and x_He = 0.10, respectively; these are typical values given the optical-only fits of the other sample stars.

F.6 H36

From the optical-only run we find T_eff = 42000 K for this star. This value is low compared to what we find from the optical + UV run, 49500 K, and inconsistent with the spectral type of the star (O2 If*). The optical-only fit is good (i.e. profiles of models covering the uncertainty ranges coincide with the scatter of the data points), with the exception of N v λλ4604–4620, which is too weak in our models. Upon comparison of the optical lines of the optical-only and optical + UV fits we see that the higher temperature of the optical + UV run does not improve the fit of most optical lines, not even for N v λλ4604–4620, which has a similar best fit profile in both our runs. Most notable is the worsening of the Balmer line fits, which generally are too deep for the optical + UV fit, but match the data well for the optical-only fit. In this temperature regime and with this data quality, the He i and He ii are of not much help for constraining the temperature, between our optical-only and optical + UV runs these lines are, within the noise level, unchanged. Bestenlehner et al. (2020) find, based on optical data only, a high temperature for this star (52000 K). Their N v λλ4604–4620 seems to have a better fit than we get, this seems to be at the cost of the fitness of He ii λ4686, and, as for our optical + UV run, of the Balmer lines, which for their best model are too deep compared to what is observed. However, the higher ionisation ions, both the N iv lines in the optical and UV as well as the strong O v λ1371, do strongly support a high temperature, which we therefore deem more likely for this source.

F.7 H129

The position of H129 in the HRD is on the left of main sequence. We do not use the parameters derived for this star for the further analysis, but exclude it from the analyses of the sample as a whole regarding mass-loss and clumping properties (see also Sect. 4.1). For the optical-only run we estimated υ∞ based on L/L_⊙ from Bestenlehner et al. (2020), but the extrapolated velocity was only 100 km s⁻¹ we assumed a velocity of 500 km s⁻¹ instead.

Appendix G X-ray

G.1 X-rays implementation in Fastwind

Wind-embedded shocks and associated X-ray emission are implemented into Fastwind by assuming a very small fraction of the stellar wind to be very hot and emit X-rays³² (Carneiro et al. 2016). The shocks are then described by five parameters: the volume filling fraction of the X-ray emission f_X, the maximum jump velocity of the shocks u∞ and the exponent γ_X that relates the outflow velocity to the jump velocity (u_∞ and γ_X together describing the shock temperature, typically being of the order 10⁶ K), and two parameters related to the onset radius of the X-ray emission³³: $R_{\min }^{\text{input}}$ and a factor m_X. The total X-ray luminosity is thus not a free parameter, but follows as output from the model. For further details on the implementation of wind-embedded shocks in Fastwind see Carneiro et al. (2016).

G.2 Assumptions regarding X-rays

For both optical-only and optical + UV fits, we include wind-embedded shocks and resulting X-rays. The corresponding parameters are not fitted, but instead fixed at certain values for each star. For all stars, we assume γ_X = 0.75, m_X = 30, and $R_{\min }^{\text{input}}$. The value of γ_X sets the gradient of the shock strength relative to the wind velocity (Pauldrach et al. 1994); a value of γ_X < 1 means that the relative increase in shock jump velocity as a function of radius is higher than that of the wind velocity. Our assumed value γ_X = 0.75 lies in between the higher value assumed by Krticka & Kubát (2009); Carneiro et al. (2016) and the lower values adopted in models of Pauldrach et al. (1994). For m_X we use the best fit value from Pauldrach et al. (2001), and in our choice of $R_{\min }^{\text{input}}$ we follow, for example, Pauldrach et al. (1994).

The X-ray parameters having the most profound influence on the ionisation fractions are the maximum jump velocity of the shocks u_∞ and the volume filling fraction of the X-ray emission f_X (Carneiro et al. 2016), because they directly relate to the temperature distribution and total X-ray flux. We tailor these parameters per star; for the first, we assume u_∞ = 0.3υ_∞ (after Carneiro et al. 2016, who follow Krtička & Kubát 2009). Here we take υ_∞ values from Crowther et al. (2016), or use the estimated values as discussed in Sect. 3.5.1. For two stars with υ_∞ ≈ 3500 (H36 and H46) we assume u_∞ = 0.2υ_∞, based on LDI simulations of winds of luminous stars with high terminal velocities (F. Driessen, priv. comm.). The typical maximum jump velocities we assume translate to maximum shock temperatures ranging from 0.6 – 28.0 ··10⁶ K.

Then, given these assumptions, we choose f_X such that the total X-ray luminosity in the ROSAT band (0.1–2.5 keV) equals L_X/L = 10⁻⁷, the canonical ratio for O-stars (e.g. Chlebowski et al. 1989; Berghoefer et al. 1997; Sana et al. 2006; Nazé et al. 2011; Rauw et al. 2015. for Galactic O-stars; Crowther in prep. confirm that this is a reasonable assumption for our sample stars, finding an average of log L_X/L = −6.6 ± 0.3 for nine X-ray sources associated with R136a.). However, L_X is not a direct input parameter of Fastwind. Instead, we estimate f_X as to give L_X close to L_X/L = 10⁻⁷. This is done for each model individually during the KIWI-GA run. Based on $\dot{M}$ and υ_∞ of the model, f_X is computed to satisfy: $$\log (f_{\text{X}})=-5.45-1.05log({\dot{M}}_6/{\upsilon }_{\infty }),$$(G.1)

where ${\dot{M}}_6=\dot{M}/{10}^{-6}$, with $\dot{M}$ in M_⊙ yr⁻¹. This equation is derived from the observational values of Kudritzki et al. (1996, their Figure 6). We note that we extrapolate the relation that they find to weaker winds (lower ${\dot{M}}_6/v_{\infty }$). The input and output X-ray parameters for each source can be found in Table G.1.

G.3 Validity of the X-ray assumptions

We check the validity of the assumptions for our X-ray setup by comparing the output log L_X/L of each optical + UV run to the canonical value, log L_X/L = −7. Figure G.1 shows that indeed, generally the output values are close to canonical. In the six cases where it is not, we underestimated the assumed X-ray flux. Exact values of X-ray input and output can be found in Table G.1.

Appendix H Trends in wind structure

Figure H.1 to H.4 show the trends in the wind structure parameters as a function of effective temperature, luminosity, β, and the Eddington parameter for electron scattering. As these parameters are strongly correlated with each other, in all these plots similar trends are visible (see Sect. 5.2).

Appendix I Additional tables

Table I.2 contains the values of the optical-only runs. Table I.1 contains stellar masses, ages, and ionising fluxes based on the parameters derived from optical + UV runs. Furthermore, it contains the temperatures derived based on the iron continuum. We note that the uncertainties we quote for M_evol, M_ini, and the age are only statistical uncertainties that result from the BONNSAI tool. Systematic uncertainties, that is, those resulting from the chosen input-physics of the evolutionary model, for example, the assumed mixing scheme, or mass-loss rate prescription, are not included. We discuss the systematic uncertainties on evolutionary masses for the WNh stars in Sect. 5.3.1. Table I.3 lists for each star the diagnostic lines that were used in the analysis.

References

---

All Tables

All Figures