11institutetext:  Anton Pannekoek Institute for Astronomy, University of
Amsterdam, 1090 GE Amsterdam, The Netherlands 11email<EMAIL_ADDRESS>22institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001
Leuven, Belgium 33institutetext: Department of Physics and Astronomy
University of Sheffield, Sheffield, S3 7RH, UK 44institutetext: LMU München,
Universitätssternwarte, Scheinerstr. 1, 81679 München, Germany
55institutetext: Department of Aerospace, Physics and Space Sciences, Florida
Institute of Technology, 150 W. University Boulevard, Melbourne, FL 32901, USA
66institutetext: European Southern Observatory, Alonso de Córdova 3107,
Vitacura, Santiago, Chile 77institutetext: Centro de Astrobiología, CSIC-INTA.
Crtra. de Torrejón a Ajalvir km 4. 28850 Torrejón de Ardoz (Madrid), Spain
88institutetext: Argelander-Institut für Astronomie, Universität Bonn, Auf dem
Hügel 71, 53121 Bonn, Germany 99institutetext: Heidelberger Institut für
Theoretische Studien, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany
1010institutetext: Astronomisches Rechen-Institut, Zentrum für Astronomie der
Universität Heidelberg, Mönchhofstr. 12-14, 69120 Heidelberg, Germany
1111institutetext: Armagh Observatory, College Hill, Armagh BT61 9DG, UK

# The R136 star cluster dissected with Hubble Space Telescope/STIS. III. The
most massive stars and their clumped winds

Sarah A. Brands 11 Alex de Koter 1122 Joachim M. Bestenlehner 33 Paul A.
Crowther 33 Jon O. Sundqvist 22 Joachim Puls 44 Saida M. Caballero-Nieves 55
Michael Abdul-Masih 2266 Florian A. Driessen 22 Miriam García 77 Sam Geen 11
Götz Gräfener 88 Calum Hawcroft 22 Lex Kaper 11 Zsolt Keszthelyi 11 Norbert
Langer 88 Hugues Sana 22 Fabian R.N. Schneider 991010 Tomer Shenar 11 Jorick
S. Vink 1111

###### Abstract

Context. The star cluster R136 inside the Large Magellanic Cloud hosts a rich
population of massive stars, including the most massive stars known. The
strong stellar winds of these very luminous stars impact their evolution and
the surrounding environment. We currently lack detailed knowledge of the wind
structure that is needed to quantify this impact.

Aims. To observationally constrain the stellar and wind properties of the
massive stars in R136, in particular the wind-structure parameters related to
wind clumping.

Methods. We simultaneously analyse optical and ultraviolet spectroscopy of 53
O-type and 3 WNh-stars using the Fastwind model atmosphere code and a genetic
algorithm. The models account for optically thick clumps and effects related
to porosity and velocity-porosity, as well as a non-void interclump medium.

Results. We obtain stellar parameters, surface abundances, mass-loss rates,
terminal velocities and clumping characteristics and compare these to
theoretical predictions and evolutionary models. The clumping properties
include the density of the interclump medium and the velocity-porosity of the
wind. For the first time, these characteristics are systematically measured
for a wide range of effective temperatures and luminosities.

Conclusions. We confirm a cluster age of 1.0-2.5 Myr and derive an initial
stellar mass of $\geq 250$ M⊙ for the most massive star in our sample, R136a1.
The winds of our sample stars are highly clumped, with an average clumping
factor of $f_{\rm cl}=29\pm 15$. We find tentative trends in the wind-
structure parameters as a function of mass-loss rate, suggesting that the
winds of stars with higher mass-loss rates are less clumped. We compare
several theoretical predictions to the observed mass-loss rates and terminal
velocities and find that none satisfactorily reproduces both quantities. The
prescription of Krtička & Kubát (2018) matches best the observed mass-loss
rates.

###### Key Words.:

– Stars: massive – Stars: fundamental parameters – Stars: winds, outflows –
Stars: mass-loss – Galaxies: star clusters: individual: R136 – Magellanic
Clouds

##  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\ {\rm M}_{\odot}$ and
a few even surpass $\sim 150\ {\rm M}_{\odot}$ (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 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 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, that are estimated to have a
characteristic mass of tens to hundreds 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 is residing 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\gtrsim 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, 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
spectroscopy of this dataset, Crowther et al. (2016) derive a cluster age of
$1.5\pm^{0.3}_{0.7}$ Myr, and find that the He ii $\lambda$1640 emission is
completely dominated by stars with initial masses $\gtrsim 100$ M⊙.
Bestenlehner et al. (2020) focus on the optical spectra and derive 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 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
inhomogeneities 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 the 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 out-flowing 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$\alpha$ line, but can also affect lines indirectly due to changes in the
ionisation/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, 2000; Owocki et
al. 2004). These velocity-porosity or ‘vorosity’ effects impact mostly
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 was 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 in detail the mass-loss rates and wind structure. 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 will yield the most accurate mass-loss rate
measurements possible with current model atmosphere codes, and furthermore
will allow us for the first time to investigate 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 Section 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
parameterisation of the wind structure parameters (Section 3.1.1). The results
of our analysis are presented in Sect. 4. This section is concluded with
several tests of robustness (Section 4.7). We discuss our results in the
context of theoretical predictions and evolutionary models in Sect. 5; Section
5.1 is dedicated to mass-loss rates and wind momentum; Section 5.2 to the
potential trends that we observe in the wind structure parameters, and in
Section 5.3 we consider our findings in the context of stellar evolution. Two
methods for measuring terminal velocities are compared in Section 5.4. We
conclude with a summary and outlook (Sect. 6).

Figure 1: HST/WFC3 V-band (F555W) photometry of the core of R136 (O’Connell
2010). Positions of the stars in our sample are indicated (yellow to red
closed circles) with respect to the slits of the HST/STIS observations (light
blue lines). The colour of the circles indicates the current (evolutionary)
mass of each source as derived in this paper. Identifications starting with
‘H’ from Hunter et al. (1995); those starting with ‘a’ and ‘b’ from Weigelt &
Baier (1985). H68 and H129 (located $4.61"$ and $\sim 4.10"$ from a1,
respectively) are part of our sample but fall outside the region shown here.
Three stars identified here (blue open circles) are not in our sample: H42 and
H77 (SB2s), and H39 (analysed in de Koter et al. 1998, outside our slit
coverage). Twelve stars of our sample overlap with that of de Koter et al.
(1997, 1998, black open circles). Not indicated in the image are H17, north of
a1, and the two components of a6, H19 (north) and H26 (south); see Section 2.1
for more details.

##  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 5 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 dataset111For one star, R136a8,
there are no optical spectra available except the H$\alpha$ line. This star
was not included in the sample of Bestenlehner et al. (2020), but is included
in our optical + UV analysis. 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

Table 1: Observational setup of the HST/STIS long-slit spectroscopy.  Grating | Grating positions | Wavelength | $\lambda/\Delta\lambda$
---|---|---|---
G140L | 1425 | 1150 – 1717 Å | 1250
G430M | 3936, 4194, 4451, 4706 | 3795 – 4743 Å | 7700
G750M | 6581 | 6297 – 6866 Å | 5850

For our analysis we use blue-optical, H$\alpha$, and far-UV HST-STIS
spectroscopy (PI: Crowther, Crowther et al. 2016). For this, 17 HST-STIS long-
slit (52”x0.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°/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
H.2). 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 mas222These sources are clearly visible in
the extreme adaptive optics VLT/SPHERE K-band images of Khorrami et al. (2021,
their Fig. 2), but were already identified by Hunter et al. (1995) with
HST/WFC2. (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 $\pm$2
pixels (see also Section 2.3). H$\alpha$ suffers from strong nebular emission,
for which we correct by interpolating the H$\alpha$ 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$\alpha$. Figure A.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

Figure 2: UV normalisation and $T_{\rm eff}$ fitting procedure for H35. The
left panels show _a)_ the observed flux in units of
$10^{-14}\leavevmode\nobreak\ {\rm ergs\ cm}^{-2}{\rm s}^{-1}$ Å-1, _b)_ the
observed flux divided by the normalised cmfgen model (dark blue points, masked
points are shown semi-transparent) and a fit through the non-masked points
(red solid line), and _c)_ the normalised spectrum (dark blue solid line), the
normalised, Ly-$\alpha$ corrected spectrum (blue dashed line) and the cmfgen
model used for obtaining the normalisation model (solid orange line). The
yellow shaded areas indicate the regions used for determining the goodness of
fit of the iron pseudo-continuum of the model; roman numerals indicate the
dominant iron ion in each region. The right panel shows for each model $T_{\rm
eff}$ and $\varv_{\rm rad}$ the dimensionless goodness of fit to the
normalised iron pseudo-continuum of the data (defined as the inverse of the
squared sum of the difference between the normalised fluxes of the model and
the data in the iron continuum). At the peak lie the $T_{\rm eff}$ and
$\varv_{\rm rad}$ that give the best fit: these are the values that are used
for the final normalisation shown on the left.

For 10 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 $\lambda$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”)333For the same reason, de Koter et al. (1997,
1998) do not analyse the spectrum of R136a2.. 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 $\lambda$1718 in our fitting
for the 10 remaining spectra of their sample, while using HST/STIS data for
all other lines. Moreover, the resolving power ($\lambda/\Delta\lambda$) of
these spectra is approximately $5000$, which allows us to resolve the
interstellar C iv $\lambda\lambda$1548-1551 lines and use this for the
correction of the HST/STIS data (see Section 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$\alpha$ 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$\alpha$
where nebular subtraction may have been imperfect 444The removal of points in
the core of H$\alpha$ may increase uncertainties on the mass-loss rate
determination. Given the UV wind lines that we take into account in our fits,
we expect this effect to be of minor importance..

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$\alpha$ and He ii $\lambda$4686) we use the synthetic
lines of a small grid of Fastwind models to determine the radial velocity
shift.

As described in Section 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 on the order of $0-100$ km s-1 (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 Section A.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_{\rm 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:

  * •

Select a cmfgen model from the grid that matches best the stellar and wind
parameters from Bestenlehner et al. (2020). _But see also below_.

  * •

Mask the wind lines and interstellar lines, so that only the pseudo-continuum
is left.

  * •

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.

  * •

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_{\rm eff}$ of the model matches that of the observed
spectrum. In order to assure this, we treat $T_{\rm eff}$ as a free parameter.
We repeat the above steps for all $T_{\rm eff}$ in the cmfgen grid (ranging
from 35-56 kK in 9 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 $\varv_{\rm rad}$ of each model (in steps of 25
km s-1, which is about a tenth of the resolution element) and assess which
value fits best. Just as for the optical, this $\varv_{\rm rad}$ value
includes a possible correction for the wavelength calibration (see also
Section 2.3). For the micro-turbulent velocity $\varv_{\rm micro}$ we assume
10 km s-1 for all models, as this is the only value included in the grid of
Bestenlehner et al. (2014) and because the exact value of $\varv_{\rm 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_{\rm eff}$,
and corrected for $\varv_{\rm rad}$. Before accepting a fit we check it
visually; extra care is taken for sources where the fit value of $T_{\rm eff}$
falls outside the uncertainty margins of the $T_{\rm eff}$ derived by
Bestenlehner et al. (2020). For six sources the model with $T_{\rm 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_{\rm eff}$ (see Table H.2). The GHRS data are normalised in a similar way,
but instead of fitting $T_{\rm 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_{\rm 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 H.2 and are
compared to the H, He, C, N, O temperature measurements in Section C.1.

We obtain the S/N of the HST-STIS UV spectra by using the HST-STIS exposure
time calculator555https://etc.stsci.edu/etc/input/stis/spectroscopic/,
assuming the F555W magnitude from De Marchi et al. (2011), and exposure times
and $A_{\rm 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-$\alpha$), the Si iv doublet at 1393.76-1402.77 Å, and C iv
doublet at 1548.20-1550.77 Å. We correct for Ly-$\alpha$ and C iv
$\lambda\lambda$1548-1551 by recovering the H i and C iv column density in the
line of sight ($N_{\rm H}$ and $N_{\rm C\sc{iv}}$, respectively), by fitting
the interstellar profiles with a Voigt-Hjerting function666Often simply called
Voigt function (e.g., Mihalas 1978, Hubeny & Mihalas 2014) (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-$\alpha$, 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-$\alpha$ 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_{\rm 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-$\alpha$ profile. For finding the average of $N_{\rm H}$ we fit
the Ly-$\alpha$ profiles of the 29 stars brighter than $M_{V}=-5.50$ (values
from Bestenlehner et al. 2020) and obtain a value of $\log(N_{\rm H}\ [{\rm
cm}^{-2}])=21.88\pm 0.07$, in good agreement with other $N_{\rm H}$
measurements towards 30 Doradus (e.g., de Boer et al. 1980, who find
$\log(N_{\rm H}\ [{\rm cm}^{-2}])=21.85\pm^{0.10}_{0.15}$).

For C iv $\lambda\lambda$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_{\rm C\sc{iv}}[{\rm cm}^{-2}])=14.72\pm 0.10$ and $14.21\pm 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
$\lambda\lambda$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 $\lambda\lambda$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 $\lambda\lambda$1394-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.

Figure 3: The Ly-$\alpha$ fitting procedure for H35. The upper panel shows the
original data (blue: flux points indicated with a large dot are considered in
the fit, small ones are not) and the best fit Ly-$\alpha$ profile (orange
dashed line). The lower panel shows the Ly-$\alpha$ corrected flux and the
values of $N_{\textrm{H}}$ and $\lambda_{Ly-\alpha}$ used for the correction.
Indicated with vertical lines are the rest wavelengths of transitions of
diagnostic lines (yellow dotted) and the position of the edge of the N v
$\lambda$1240 line (yellow dashed).

###  2.6 VLT/SPHERE $K_{s}$ photometry

In our fitting procedure we calibrate the luminosity with an observed stellar
flux (see Section 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\leavevmode\nobreak\ \mu\textrm{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-blanketing777See, e.g., Pauldrach et al. 2001, for an explanation of
these concepts.. 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 $\varv_{\rm 0}$ near the sonic point. The
stellar wind is parameterised by a mass loss rate $\dot{M}$, a terminal
velocity $\varv_{\infty}$, and a wind acceleration parameter $\beta$. The wind
velocity $\varv_{\rm r}$ as a function of radius $r$ is expressed by the
classic $\beta$-velocity law:

$\varv_{\rm r}(r)=\varv_{\infty}(1-b/r)^{\beta},$ (1)

where $b$ is a radius close to the stellar radius888As defined in Santolaya-
Rey et al. (1997), their Eq. (10). $R_{\star}$, the exact value of $b$
depending on $\varv_{\rm 0}$ (see Santolaya-Rey et al. 1997). Under these
assumptions the structure and ionisation/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.

Figure 4: Best fits of the interstellar C iv $\lambda\lambda$1548-1551 lines
for all sources of de Koter et al. (1998). For each source we show the
normalised flux of the interstellar lines (black dots) and the best fit
(orange lines). The last panel contains the best fit of all sources (blue
lines).  Figure 5: Illustration of the effects of clumping and porosity in
velocity space. Each sketch shows the star (yellow) and its clumpy wind
(different shades of grey). In the left figure we indicate different
parameters that describe the wind structure. In the middle and right sketches
we illustrate how the effective transparency of spectral lines in the wind
depends on volume filling fraction as well as on the velocity span of the
clumps. Depicted are photons of different frequencies ($\nu_{1}$ in red and
$\nu_{2}$ in blue, with $\nu_{1}<\nu_{2}$) that are intercepted by a strong
absorption line if they are Doppler shifted with the right velocity. The
corresponding resonance zones lie at 0.5$v_{\infty}$ and 0.6$v_{\infty}$
(shaded red and blue, respectively).

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 Section 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
synthesised, but their radiation field is taken into account, which is
essential for the treatment of effects of line-blocking/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 $\tau_{\rm
e}=0.28$ (for a typical O-star, we find $\tau_{\rm e}\approx 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 $\rho_{\rm cl}$ and an interclump medium with a lower density
$\rho_{\rm ic}$. The clumps occupy a certain fraction of the total wind volume
$f_{\rm vol}$, referred to as the volume filling factor. The clumping factor
$f_{\rm cl}$ of the medium can be written as:

$f_{\rm
cl}\equiv\frac{\langle\rho^{2}\rangle}{\langle\rho\rangle^{2}}=\frac{f_{\rm
vol}\rho_{\rm cl}^{2}+(1-f_{\rm vol})\rho_{\rm ic}^{2}}{(f_{\rm vol}\rho_{\rm
cl}+(1-f_{\rm vol})\rho_{\rm ic})^{2}}.$ (2)

We note that for the conventional assumption (not adopted here) of a void
interclump medium ($\rho_{\rm ic}=0$) this would lead to $f_{\rm cl}=1/f_{\rm
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
$\varv_{\infty}$). 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_{\rm vol}$. This effect is
illustrated in the middle panel of Fig. 5.

However, the clumps do not necessarily follow the average velocity field999The
velocity field can, for example, be severely altered by shocks (see Section
3.1.3).. For example, the velocity span of the clumps, $\delta\varv$, can be
larger than that of the underlying smooth field $\delta\varv_{\rm 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_{\rm vor}$, depends on both $f_{\rm vol}$ and the
relative velocity span of the clumps $\delta\varv/\delta\varv_{\rm sm}$:

$f_{\rm vor}=f_{\rm vol}\left|\frac{\delta\varv}{\delta\varv_{\rm
sm}}\right|.$ (3)

In Fastwind a normalised version of this factor is implemented:

$f_{\rm vel}=\frac{f_{\rm vor}}{1+f_{\rm vor}}=\frac{f_{\rm
vol}(\delta\varv/\delta\varv_{\rm sm})}{1+f_{\rm
vol}(\delta\varv/\delta\varv_{\rm sm})},$ (4)

taking values from 0 to 1. The parameter $f_{\rm vel}$ is called the velocity
filling factor. Note that this equation reduces to the purely geometrical
effect, depending only on $f_{\rm vol}$, when $\delta\varv$ follows the
underlying smooth field:

$f_{\rm vel}=\frac{f_{\rm vol}}{1+f_{\rm vol}},\hskip
28.45274pt[\delta\varv/\delta\varv_{\rm sm}\rightarrow 1].$ (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 $\rho=\langle\rho\rangle f_{\rm 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
$\langle\chi\rangle$ 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 $\chi_{\rm eff}$ can be expressed as:

$\chi_{\rm eff}=\langle\chi\rangle\frac{1+\tau_{\rm cl}f_{\rm ic}}{1+\tau_{\rm
cl}},$ (6)

with $\tau_{\rm cl}$ the clump optical depth (Sundqvist et al. 2014). The
interclump density contrast, $f_{\rm ic}$, is defined as:

$f_{\rm ic}\equiv\rho_{\rm 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_{\rm cl}=\frac{\tau_{\rm S}}{f_{\rm vor}}(1-(1-f_{\rm vol})f_{\rm ic}),$
(8)

with $\tau_{\rm 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$ ($\equiv l_{\rm cl}/f_{\rm vol}$, with $l_{\rm 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(\varv_{\rm r})=h_{\infty}\varv_{\rm r}/\varv_{\infty},$ (9)

with $h_{\infty}$ the porosity length at the terminal wind velocity, given as
input by the user. In this work we adopt $h_{\infty}=R_{\star}$, 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 $\varv_{\rm cl,start}$, the maximum
clumping factor $f_{\rm cl}$, and the velocity at which this maximum clumping
factor is reached, $\varv_{\rm 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 $\varv_{\rm r}=\varv_{\rm cl,start}$ until $\varv_{\rm
r}=\varv_{\rm cl,max}$ the clumping factor increases linearly with wind
velocity from 1 to $f_{\rm cl}$, staying constant at $f_{\rm cl}$ for
$\varv_{\rm r}>\varv_{\rm 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
$\varv_{\rm r}>\varv_{\rm cl,max}$ the clumping stays constant at the maximum
value, $f_{\rm cl}$. The values for $\varv_{\rm cl,start}$, $\varv_{\rm
cl,max}$ and $f_{\rm cl}$ are specified by the user (see Section 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 Section 3.1.1 are all used in the
process of computing the ionisation/excitation structure of the model
atmosphere. An additional parameter, the wind turbulence velocity $\varv_{\rm
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/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 $\varv_{\rm micro}$ at the base of the wind, to $\varv_{\rm
windturb}$, at the point where the wind reaches its terminal velocity (Haser
et al. 1995). The wind turbulence velocity is typically on the order of
$0.1\varv_{\infty}$ (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).

Table 2: Parameters that describe the wind structure in Fastwind. All these
parameters can directly affect the shape of a spectral line, except for $h$,
which is only impacting continuum opacity. Note that $f_{\rm cl}$ is in fact
the maximum value of clumping that is reached at $\varv_{\rm cl,max}$, and
$\varv_{\rm windturb}$ the maximum value of the wind turbulence reached at the
terminal velocity. Furthermore, note that not all parameters listed here are
treated as free parameters in our analysis (see Section 3.5).  Wind structure
parameters in Fastwind
---
$f_{\rm cl}$ |  Clumping factor
$f_{\rm ic}$ |  Interclump density contrast
$f_{\rm vel}$ |  Velocity filling factor
$\varv_{\rm cl,start}$ |  Onset of clumping
$\varv_{\rm cl,max}$ |  Clumping reaches maximum ($f_{\rm cl}$)
$h$ |  Porosity length
$\varv_{\rm windturb}$ |  Wind turbulence

####  3.1.3 Wind-embedded shocks & 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 F.

###  3.2 A new genetic algorithm: Kiwi-GA

Figure 6: Flowchart portraying the workings of Kiwi-GA. In each generation,
depicted by the green circle with arrows, a large amount of models is computed
(typically $50-250$). The reproduction step is a means of natural selection
towards fitter models, which eventually leads to convergence of the algorithm
towards the best fitting solution (typically after $30-100$ generations).

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 $\chi^{2}$ value as a fitness measure:

$\chi^{2}=\sum_{i=0}^{N}\left(\frac{\mathcal{F}_{{\rm
obs},i}-\mathcal{F}_{{\rm mod},i}}{\mathcal{E}_{{\rm obs},i}}\right)^{2},$
(10)

where $N$ is the number of data points of the spectrum that is considered in
the fit, $\mathcal{F}_{{\rm obs},i}$ the observed normalised flux,
$\mathcal{F}_{{\rm mod},i}$ the normalised flux of the model, and
$\mathcal{E}_{{\rm 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 $\chi^{2}$) 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_{\rm eff}$ that is similar
to that of the observed star will generally have a better fit than a model
with a $T_{\rm eff}$ that is far off, and this value of $T_{\rm 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; Ramírez-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 B.
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
users101010https://github.com/sarahbrands/Kiwi-GA.

###  3.3 Stellar radius & luminosity

In our Kiwi-GA runs the stellar radius is estimated for each model
individually using an observed, extinction corrected magnitude (Section 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_{\rm 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}$
filter111111The filter is specified by the user and any filter of which a
transmission curve is available can be chosen. Currently the following filters
are implemented: Johnson V, HST F555W, and VLT/SPHERE $K_{s}$. The
transmission curves are taken from the SVO Filter Profile Service:
http://svo2.cab.inta-csic.es/theory/fps/. The Planck curve thus serves as an
‘SED estimate’ during the run and the radius is an output of the run, as is
luminosity121212Kiwi-GA also has the option to set the radius to a fixed value
for all models.. 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}/(R_{\star}\varv_{\infty})^{3/2}$
(Puls et al. 1996, Holgado et al. 2018, their Appendix B). 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_{\rm eff}$ would be a better
guess – the previous estimate of $T=0.9T_{\rm 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 $\chi^{2}$)
and use this model to implicitly adjust the uncertainty of each flux point
that is fitted, such that the $\chi^{2}_{\rm red}$ value of the best fitting
model equals unity. These adjusted flux uncertainties are then used for
recomputing the $\chi^{2}$ of each model in the run. This procedure is
equivalent to dividing all $\chi^{2}$ values of the run by the (original)
$\chi^{2}_{\rm red}$ 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 $\Gamma(\chi^{2}/2,\nu/2)$ the incomplete gamma function, representing
the cumulative distribution function of the $\chi^{2}$ distribution, evaluated
at $\chi^{2}$, for $\nu=n_{\rm data}-n_{\rm free}$ the degrees of freedom,
where $n_{\rm data}$ is the number of flux points that is taken into account
during the fit and $n_{\rm free}$ the number of free parameters. The best
fitting models are all models where $P>0.32$ (i.e., the 68% confidence
interval, we will call this $1\sigma$) 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\sigma$
and $2\sigma$, respectively. For convenience, we will refer to the 68% and 95%
confidence intervals as $1\sigma$ and $2\sigma$ uncertainties, even though the
confidence intervals we derive are not necessarily symmetric and Gaussian. In
all tables we present 1$\sigma$ 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 $\chi^{2}$ 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 $\chi^{2}$ 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_{\star}^{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 C. The details of each fitting setup are summarised in Table 3 and
explained in detail below.

Table 3: Free parameters in the optical-only and optical + UV fits. Parameters in brackets are free only in a subset of the runs, see text. The last column contains for each fit a reference to the table where the best fit parameters are presented. Parameter names are defined in the text.  Fit |  Free parameters |  Results
---|---|---
Optical-only |  $T_{\rm eff}$, $g$, $\varv_{\rm eq}\sin i$, $\dot{M}$, $x_{\rm He}$, ($n_{\rm N}$, $\beta$) |  Table H.3
Optical + UV (high S/N) |  $T_{\rm eff}$, $g$, $\dot{M}$, $\beta$, $\varv_{\infty}$, $f_{\rm cl}$, $\varv_{\rm cl,start}$, $\varv_{\rm windturb}$, $f_{\rm ic}$, $f_{\rm vel}$, $n_{\rm N}$, $n_{\rm C}$, ($n_{\rm O}$, $x_{\rm He}$) |  Table 4, 6
Optical + UV (low S/N) |  $T_{\rm eff}$, $g$, $\dot{M}$, $\beta$, $\varv_{\infty}$, $f_{\rm cl}$ |  Table 5

####  3.5.1 Optical-only setup

The optical-only runs have 5 to 7 free parameters $n_{\rm free}$, as specified
in Table 3. Here, $g$ is the gravitational acceleration, $\varv_{\rm eq}\sin
i$ the projected rotational broadening, $\dot{M}$ the mass-loss rate, and
$x_{\rm He}$ the helium surface abundance, where $x_{\rm He}=n_{\rm He}/n_{\rm
H}$, with $n_{\rm He}$ and $n_{\rm H}$ the helium and hydrogen number density.
If any line is (partially) in emission, we also fit the wind acceleration
parameter $\beta$, and when we see a nitrogen line above the noise we fit the
nitrogen abundance $n_{\rm N}$ (with $n_{\rm N}$ the number density). The
other parameters are held fixed at the values discussed below.

We assume a smooth wind ($f_{\rm cl}=1$) for all stars except for the WNh
stars, for which we assume $f_{\rm cl}=10$. Furthermore, we assume $\varv_{\rm
micro}$ $=10$ km s-1, and in case $\beta$ is fixed we assume $\beta=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 $\varv_{\rm eq}\sin i$, assuming $\varv_{\rm macro}$ $=0$ km s-1.
In practice this means that all broadening is captured in a single parameter
$\varv_{\rm eq}\sin i$ and, since for our stars likely $\varv_{\rm macro}$
$>0$ km s-1 (see e.g., Simón-Díaz et al. 2017), the projected rotational
velocities that we find are upper limits of the actual $\varv_{\rm eq}\sin i$.
The derived $\varv_{\rm 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\leavevmode\nobreak\
{\rm Z}_{\odot}$. We assume the values of Crowther et al. (2016) for the
terminal velocities of the winds $\varv_{\infty}$. For 12 stars
$\varv_{\infty}$ was not available and in these cases we estimate the
velocities by inter- and extrapolating the dependence of $\varv_{\infty}$ on
luminosity $L$, that we empirically find using the values of Crowther et al.
(2016) and Bestenlehner et al. (2020), for $\log L/L_{\odot}<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_{\rm C}$ refers to the carbon abundance (by number). For the WNh-stars we
fit two additional free parameters (see also below): oxygen abundance $n_{\rm
O}$ (by number) and $x_{\rm 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}_{\rm C}/{n}_{\rm H})+12=7.75$ for carbon, $\log({n}_{\rm N}/{n}_{\rm
H})+12=6.9$ for nitrogen and $\log({n}_{\rm O}/{n}_{\rm 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_{\rm ic}$ $=0.05$,
$f_{\rm vel}$ $=0.15$, $\varv_{\rm cl,start}$/$\varv_{\infty}=0.05$,
$\varv_{\rm windturb}$/$\varv_{\infty}=0.15$ (Section 5.2). For all optical +
UV runs the velocity at which the maximum clumping factor is reached is given
by $\varv_{\rm cl,max}/\varv_{\infty}={\rm max}(0.3,\varv_{\rm
cl,start}/\varv_{\infty})$.

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}_{\rm O}/{n}_{\rm 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 in Appendix
I. We discuss the robustness of this setup in Section 4.7.

Figure 7: Kiwi-GA output summary for the optical + UV run of H35 (60
generations). The top of the figure shows the line profiles that were
considered in the fit. For each spectral line we show the observed spectrum
(black bars), the best fit model (green solid line), and the family of best
fitting models, that is, the region spanned by all models in the $2\sigma$
confidence interval (light green shaded area). The bottom of the figure shows
for each free parameter the goodness-of-fit (expressed as $1/$$\chi^{2}_{\rm
red}$) of each model of the run represented by a dark blue dot). The position
of the best model, as well as the 1$\sigma$ and 2$\sigma$ error regions (dark
and light shaded yellow, respectively) are indicated. Output summaries for the
other runs can be found in Appendix I.

Figure 8: Best fit parameter ranges for the 39 stars where a 12-free-parameter
optical + UV fit was done, ordered by decreasing $\dot{M}$. First four columns
(blue colours) show basic stellar parameters $T_{\rm eff}$ and $\log g$ as
well as $\varv_{\rm eq}\sin i$ and $x_{\rm He}$ derived from the optical-only
fit. The next eight columns (green to red colours) show the parameters that
describe the wind and wind structure. The last two columns (pink and purple)
concern the C and N abundances. The darker and lighter shaded regions
correspond to $1\sigma$ and $2\sigma$ uncertainties, respectively. Note that
R136a8 has no optical data and thus no optical-only fit: in this case
$\varv_{\rm eq}\sin i$ and $x_{\rm He}$ are indicated with a $\bullet$.

Figure 9: As Fig. 8, but for the best fit parameter ranges for the 17 stars
where a 6-free-parameter optical + UV fit was done. The columns correspond to
the first eight columns of Fig. 8.

####  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$\alpha$ for 11 sources.
We included the following UV lines in the fits of most stars: C iv
$\lambda$1165, C iii $\lambda$1170, N v $\lambda$1240, O iv $\lambda$1340, O v
$\lambda$1371, C iv $\lambda\lambda$1548-1551, He ii $\lambda$1640\. Note
that:

  * •

For N v $\lambda$1240 we only fit the part that could be recovered after the
Ly-$\alpha$ correction (see also Section 2.5).

  * •

In one case we fitted Si iv $\lambda\lambda$1394-1402 (R136b, see also Section
2.5).

  * •

Part of O iv $\lambda$1340 is clipped because of the presence of the strong
interstellar Cii 1336 line.

  * •

For stars where O v $\lambda$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 $\lambda$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 H.4.

We note that the UV spectroscopy includes diagnostics for $\dot{M}$,
$\varv_{\infty}$ as well as for the wind structure parameters $f_{\rm cl}$,
$f_{\rm ic}$, $f_{\rm vel}$, so that we can break the degeneracy between these
parameters. For example, the strength of H$\alpha$ 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 $\lambda$1371 and N iv $\lambda$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 $\lambda$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\sigma$ 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 in Appendix I.

The best fit parameters of the optical + UV runs for all stars are presented
in Figs. 8 and 9 and Tables 4 and 5. Notes on peculiarities of individual
sources can be found in Appendix E. In Appendix H 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. Note that we
derive 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_{\rm He}$ and $\varv_{\rm
eq}\sin i$ were taken from the optical-only analysis131313We do not fit this
in the optical + UV analysis. while the remaining parameters were taken from
the optical + UV analysis. The WNh stars are an exception: here we do measure
$x_{\rm 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 C. 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: Positions of the R136 sources (optical + UV analysis, dark blue
points) in the Hertzsprung-Russell diagram. Yellow dashed and dotted yellow
lines are LMC evolutionary tracks of Brott et al. (2011a) and Köhler et al.
(2015), respectively. Red solid lines are isochrones. All tracks have an
initial rotation of $\approx$ 160 km s-1, representative for the O-stars in
the sample. Note that the tracks shown in Fig. 17 have a higher initial
rotation.

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 (42000 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_{\rm eff}$ to 50000 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 47000 K, where the
fitness of the 50000 K model is almost as good, but the fit to the 45000 K
model significantly less, even worse for the 42000 K model (Fig. E.14). Taking
all this into account, we consider the higher $T_{\rm eff}$ for this star more
likely and accept the parameters of the fixed-$T_{\rm eff}$ run (50000 K) as
the parameters which we use for further analysis (for more details and the
spectral fits, see Section E.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\sigma$ 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 Section C.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.

Table 4: Best fit parameters and $1\sigma$ error barsa) for the optical + UV fits of 39 stars where we fitted 12 free parameters.  Source | $\log L/{\rm L}_{\odot}$ |  $T_{\rm eff}$(K) | $\log g$ | $R_{\star}/{\rm R}_{\odot}$ | $\log\dot{M}$ |  $v_{\infty}$ (km/s) | $v_{\rm windturb}/v_{\infty}$ | $f_{\rm cl}$ | $f_{\rm ic}$ | $f_{\rm vel}$ | $\beta$ | $v_{\rm cl,start}/v_{\infty}$ |  $\log(\frac{n_{\rm C}}{n_{\rm H}})$+12 |  $\log(\frac{n_{\rm N}}{n_{\rm H}})$+12
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
R136a1b) |  6.86$\pm$ | ${}^{0.04}_{0.07}$ |  46000$\pm$ | ${}^{1250}_{2375}$ |  3.65$\pm$ | ${}^{0.20}_{0.20}$ |  42.7$\pm$ | ${}^{1.6}_{0.9}$ |  -4.57$\pm$ | ${}^{0.13}_{0.03}$ |  3150$\pm$ | ${}^{300}_{250}$ |  0.04$\pm$ | ${}^{0.05}_{0.03}$ |  43$\pm$ | ${}^{7}_{20}$ |  0.48$\pm$ | ${}^{0.03}_{0.12}$ |  0.97$\pm$ | ${}^{0.05}_{0.25}$ |  1.18$\pm$ | ${}^{0.23}_{0.17}$ |  0.03$\pm$ | ${}^{0.05}_{0.03}$ |  7.60$\pm$ | ${}^{0.28}_{0.60}$ |  8.75$\pm$ | ${}^{0.20}_{0.25}$
R136a2b) |  6.71$\pm$ | ${}^{0.03}_{0.03}$ |  47000$\pm$ | ${}^{1000}_{625}$ |  3.55$\pm$ | ${}^{0.25}_{0.05}$ |  34.7$\pm$ | ${}^{0.7}_{0.8}$ |  -4.48$\pm$ | ${}^{0.12}_{0.06}$ |  2900$\pm$ | ${}^{300}_{200}$ |  0.07$\pm$ | ${}^{0.03}_{0.06}$ |  29$\pm$ | ${}^{4}_{13}$ |  0.50$\pm$ | ${}^{0.01}_{0.42}$ |  1.00$\pm$ | ${}^{0.03}_{0.10}$ |  1.18$\pm$ | ${}^{0.17}_{0.23}$ |  0.03$\pm$ | ${}^{0.07}_{0.03}$ |  7.20$\pm$ | ${}^{0.45}_{0.23}$ |  9.05$\pm$ | ${}^{0.25}_{0.15}$
R136a3b,c) |  6.70$\pm$ | ${}^{0.02}_{0.02}$ |  50000$\pm$ | ${}^{2500}_{8000}$ |  4.05$\pm$ | ${}^{1.00}_{0.40}$ |  30.2$\pm$ | ${}^{0.6}_{0.6}$ |  -4.64$\pm$ | ${}^{0.08}_{0.05}$ |  2700$\pm$ | ${}^{300}_{400}$ |  0.07$\pm$ | ${}^{0.10}_{0.07}$ |  46$\pm$ | ${}^{5}_{16}$ |  0.41$\pm$ | ${}^{0.10}_{0.25}$ |  1.00$\pm$ | ${}^{0.03}_{0.03}$ |  1.30$\pm$ | ${}^{0.17}_{0.25}$ |  0.02$\pm$ | ${}^{0.09}_{0.02}$ |  7.75$\pm$ | ${}^{0.42}_{0.78}$ |  9.15$\pm$ | ${}^{0.25}_{0.25}$
R136a4 |  6.28$\pm$ | ${}^{0.02}_{0.05}$ |  50000$\pm$ | ${}^{500}_{2000}$ |  4.15$\pm$ | ${}^{0.20}_{0.15}$ |  18.5$\pm$ | ${}^{0.5}_{0.3}$ |  -5.84$\pm$ | ${}^{0.08}_{0.10}$ |  3150$\pm$ | ${}^{75}_{150}$ |  0.14$\pm$ | ${}^{0.03}_{0.07}$ |  49$\pm$ | ${}^{1}_{14}$ |  0.00$\pm$ | ${}^{0.03}_{0.01}$ |  0.45$\pm$ | ${}^{0.07}_{0.07}$ |  0.72$\pm$ | ${}^{0.06}_{0.05}$ |  0.03$\pm$ | ${}^{0.03}_{0.03}$ |  7.38$\pm$ | ${}^{0.23}_{0.12}$ |  7.90$\pm$ | ${}^{0.25}_{0.45}$
R136a5 |  6.32$\pm$ | ${}^{0.03}_{0.03}$ |  48000$\pm$ | ${}^{750}_{750}$ |  4.35$\pm$ | ${}^{0.17}_{0.35}$ |  21.1$\pm$ | ${}^{0.4}_{0.4}$ |  -5.09$\pm$ | ${}^{0.10}_{0.10}$ |  3250$\pm$ | ${}^{200}_{200}$ |  0.06$\pm$ | ${}^{0.04}_{0.03}$ |  16$\pm$ | ${}^{5}_{7}$ |  0.44$\pm$ | ${}^{0.04}_{0.09}$ |  0.30$\pm$ | ${}^{0.29}_{0.10}$ |  1.10$\pm$ | ${}^{0.07}_{0.10}$ |  0.11$\pm$ | ${}^{0.03}_{0.05}$ |  7.03$\pm$ | ${}^{0.17}_{0.05}$ |  8.30$\pm$ | ${}^{0.30}_{0.40}$
R136a6 |  6.24$\pm$ | ${}^{0.03}_{0.07}$ |  52000$\pm$ | ${}^{1000}_{3000}$ |  4.00$\pm$ | ${}^{0.10}_{0.20}$ |  16.4$\pm$ | ${}^{0.5}_{0.3}$ |  -5.60$\pm$ | ${}^{0.11}_{0.08}$ |  3200$\pm$ | ${}^{150}_{150}$ |  0.18$\pm$ | ${}^{0.01}_{0.05}$ |  34$\pm$ | ${}^{8}_{12}$ |  0.19$\pm$ | ${}^{0.13}_{0.14}$ |  0.70$\pm$ | ${}^{0.17}_{0.11}$ |  0.78$\pm$ | ${}^{0.12}_{0.10}$ |  0.05$\pm$ | ${}^{0.04}_{0.05}$ |  7.53$\pm$ | ${}^{0.09}_{0.50}$ |  7.80$\pm$ | ${}^{0.33}_{0.45}$
R136a7 |  6.36$\pm$ | ${}^{0.05}_{0.07}$ |  54000$\pm$ | ${}^{2000}_{3000}$ |  4.30$\pm$ | ${}^{0.23}_{0.40}$ |  17.5$\pm$ | ${}^{0.5}_{0.4}$ |  -5.52$\pm$ | ${}^{0.18}_{0.10}$ |  2900$\pm$ | ${}^{150}_{150}$ |  0.15$\pm$ | ${}^{0.04}_{0.07}$ |  12$\pm$ | ${}^{6}_{6}$ |  0.39$\pm$ | ${}^{0.10}_{0.18}$ |  0.25$\pm$ | ${}^{0.42}_{0.12}$ |  0.93$\pm$ | ${}^{0.15}_{0.15}$ |  0.02$\pm$ | ${}^{0.06}_{0.02}$ |  7.70$\pm$ | ${}^{0.25}_{0.30}$ |  7.10$\pm$ | ${}^{0.85}_{0.25}$
R136a8 |  6.17$\pm$ | ${}^{0.03}_{0.03}$ |  49500$\pm$ | ${}^{1250}_{1250}$ |  4.25$\pm$ | ${}^{0.30}_{0.25}$ |  16.6$\pm$ | ${}^{0.4}_{0.3}$ |  -5.82$\pm$ | ${}^{0.10}_{0.10}$ |  3000$\pm$ | ${}^{75}_{150}$ |  0.11$\pm$ | ${}^{0.04}_{0.02}$ |  37$\pm$ | ${}^{8}_{12}$ |  0.13$\pm$ | ${}^{0.05}_{0.06}$ |  0.42$\pm$ | ${}^{0.06}_{0.25}$ |  0.70$\pm$ | ${}^{0.06}_{0.03}$ |  0.02$\pm$ | ${}^{0.06}_{0.02}$ |  7.72$\pm$ | ${}^{0.35}_{0.07}$ |  6.95$\pm$ | ${}^{0.55}_{0.10}$
R136b |  6.35$\pm$ | ${}^{0.03}_{0.03}$ |  35500$\pm$ | ${}^{750}_{750}$ |  3.55$\pm$ | ${}^{0.25}_{0.12}$ |  40.0$\pm$ | ${}^{0.8}_{0.8}$ |  -5.15$\pm$ | ${}^{0.05}_{0.10}$ |  1850$\pm$ | ${}^{75}_{50}$ |  0.06$\pm$ | ${}^{0.01}_{0.03}$ |  20$\pm$ | ${}^{6}_{2}$ |  0.45$\pm$ | ${}^{0.04}_{0.07}$ |  0.95$\pm$ | ${}^{0.05}_{0.15}$ |  1.50$\pm$ | ${}^{0.17}_{0.03}$ |  0.03$\pm$ | ${}^{0.06}_{0.02}$ |  7.28$\pm$ | ${}^{0.10}_{0.15}$ |  9.20$\pm$ | ${}^{0.10}_{0.15}$
H30 |  5.76$\pm$ | ${}^{0.05}_{0.03}$ |  40000$\pm$ | ${}^{1500}_{1000}$ |  4.20$\pm$ | ${}^{0.35}_{0.15}$ |  15.9$\pm$ | ${}^{0.3}_{0.4}$ |  -6.09$\pm$ | ${}^{0.11}_{0.29}$ |  2650$\pm$ | ${}^{250}_{150}$ |  0.13$\pm$ | ${}^{0.07}_{0.09}$ |  31$\pm$ | ${}^{14}_{14}$ |  0.01$\pm$ | ${}^{0.03}_{0.01}$ |  0.10$\pm$ | ${}^{0.10}_{0.05}$ |  0.70$\pm$ | ${}^{0.12}_{0.03}$ |  0.28$\pm$ | ${}^{0.04}_{0.27}$ |  7.15$\pm$ | ${}^{0.16}_{0.17}$ |  8.50$\pm$ | ${}^{0.25}_{1.00}$
H31 |  5.98$\pm$ | ${}^{0.03}_{0.03}$ |  47500$\pm$ | ${}^{1000}_{750}$ |  4.00$\pm$ | ${}^{0.10}_{0.15}$ |  14.6$\pm$ | ${}^{0.3}_{0.3}$ |  -6.15$\pm$ | ${}^{0.11}_{0.08}$ |  3050$\pm$ | ${}^{100}_{75}$ |  0.15$\pm$ | ${}^{0.04}_{0.02}$ |  42$\pm$ | ${}^{9}_{12}$ |  0.14$\pm$ | ${}^{0.06}_{0.04}$ |  0.47$\pm$ | ${}^{0.09}_{0.09}$ |  0.70$\pm$ | ${}^{0.07}_{0.03}$ |  0.02$\pm$ | ${}^{0.08}_{0.02}$ |  7.67$\pm$ | ${}^{0.20}_{0.12}$ |  6.90$\pm$ | ${}^{0.38}_{0.05}$
H35 |  5.82$\pm$ | ${}^{0.04}_{0.03}$ |  47500$\pm$ | ${}^{1500}_{1000}$ |  4.08$\pm$ | ${}^{0.28}_{0.07}$ |  12.1$\pm$ | ${}^{0.2}_{0.3}$ |  -5.97$\pm$ | ${}^{0.05}_{0.25}$ |  3050$\pm$ | ${}^{150}_{150}$ |  0.13$\pm$ | ${}^{0.06}_{0.04}$ |  21$\pm$ | ${}^{15}_{7}$ |  0.08$\pm$ | ${}^{0.04}_{0.05}$ |  0.35$\pm$ | ${}^{0.09}_{0.10}$ |  0.70$\pm$ | ${}^{0.07}_{0.03}$ |  0.01$\pm$ | ${}^{0.07}_{0.01}$ |  7.58$\pm$ | ${}^{0.23}_{0.15}$ |  7.00$\pm$ | ${}^{0.25}_{0.15}$
H36 |  6.27$\pm$ | ${}^{0.03}_{0.03}$ |  49500$\pm$ | ${}^{750}_{1000}$ |  4.10$\pm$ | ${}^{0.35}_{0.20}$ |  18.6$\pm$ | ${}^{0.4}_{0.4}$ |  -5.29$\pm$ | ${}^{0.05}_{0.10}$ |  3900$\pm$ | ${}^{150}_{75}$ |  0.08$\pm$ | ${}^{0.03}_{0.02}$ |  31$\pm$ | ${}^{12}_{3}$ |  0.49$\pm$ | ${}^{0.02}_{0.10}$ |  0.25$\pm$ | ${}^{0.38}_{0.15}$ |  1.05$\pm$ | ${}^{0.14}_{0.10}$ |  0.03$\pm$ | ${}^{0.04}_{0.03}$ |  8.07$\pm$ | ${}^{0.12}_{0.36}$ |  8.05$\pm$ | ${}^{0.20}_{0.35}$
H40 |  5.93$\pm$ | ${}^{0.05}_{0.09}$ |  47500$\pm$ | ${}^{2000}_{3250}$ |  3.90$\pm$ | ${}^{0.20}_{0.17}$ |  13.8$\pm$ | ${}^{0.6}_{0.3}$ |  -6.12$\pm$ | ${}^{0.14}_{0.18}$ |  3250$\pm$ | ${}^{100}_{150}$ |  0.04$\pm$ | ${}^{0.06}_{0.03}$ |  44$\pm$ | ${}^{7}_{27}$ |  0.02$\pm$ | ${}^{0.14}_{0.01}$ |  0.60$\pm$ | ${}^{0.17}_{0.16}$ |  0.70$\pm$ | ${}^{0.10}_{0.03}$ |  0.27$\pm$ | ${}^{0.07}_{0.27}$ |  8.18$\pm$ | ${}^{0.03}_{0.47}$ |  6.95$\pm$ | ${}^{0.70}_{0.10}$
H45 |  5.80$\pm$ | ${}^{0.07}_{0.03}$ |  41500$\pm$ | ${}^{2500}_{1000}$ |  4.15$\pm$ | ${}^{0.17}_{0.30}$ |  15.5$\pm$ | ${}^{0.3}_{0.5}$ |  -6.27$\pm$ | ${}^{0.08}_{0.35}$ |  3100$\pm$ | ${}^{200}_{250}$ |  0.17$\pm$ | ${}^{0.03}_{0.09}$ |  28$\pm$ | ${}^{18}_{13}$ |  0.03$\pm$ | ${}^{0.10}_{0.02}$ |  0.05$\pm$ | ${}^{0.14}_{0.05}$ |  0.70$\pm$ | ${}^{0.24}_{0.03}$ |  0.39$\pm$ | ${}^{0.01}_{0.19}$ |  7.47$\pm$ | ${}^{0.38}_{0.23}$ |  7.20$\pm$ | ${}^{0.85}_{0.30}$
H46 |  6.10$\pm$ | ${}^{0.02}_{0.06}$ |  47500$\pm$ | ${}^{500}_{2500}$ |  3.90$\pm$ | ${}^{0.55}_{0.15}$ |  16.7$\pm$ | ${}^{0.5}_{0.3}$ |  -5.15$\pm$ | ${}^{0.18}_{0.06}$ |  3650$\pm$ | ${}^{200}_{200}$ |  0.09$\pm$ | ${}^{0.05}_{0.06}$ |  4$\pm$ | ${}^{3}_{4}$ |  0.39$\pm$ | ${}^{0.11}_{0.20}$ |  0.82$\pm$ | ${}^{0.17}_{0.45}$ |  1.05$\pm$ | ${}^{0.17}_{0.14}$ |  0.36$\pm$ | ${}^{0.05}_{0.15}$ |  7.40$\pm$ | ${}^{0.39}_{0.17}$ |  7.80$\pm$ | ${}^{0.23}_{0.95}$
H47 |  5.98$\pm$ | ${}^{0.05}_{0.03}$ |  43500$\pm$ | ${}^{1750}_{1000}$ |  4.45$\pm$ | ${}^{0.30}_{0.50}$ |  17.3$\pm$ | ${}^{0.4}_{0.4}$ |  -5.25$\pm$ | ${}^{0.08}_{0.11}$ |  3450$\pm$ | ${}^{150}_{150}$ |  0.03$\pm$ | ${}^{0.04}_{0.01}$ |  6$\pm$ | ${}^{7}_{4}$ |  0.30$\pm$ | ${}^{0.11}_{0.18}$ |  0.05$\pm$ | ${}^{0.49}_{0.05}$ |  0.95$\pm$ | ${}^{0.15}_{0.12}$ |  0.40$\pm$ | ${}^{0.01}_{0.18}$ |  7.33$\pm$ | ${}^{0.45}_{0.33}$ |  7.05$\pm$ | ${}^{0.35}_{0.15}$
H48 |  5.96$\pm$ | ${}^{0.03}_{0.04}$ |  46500$\pm$ | ${}^{1000}_{1500}$ |  3.90$\pm$ | ${}^{0.23}_{0.07}$ |  14.9$\pm$ | ${}^{0.4}_{0.3}$ |  -5.60$\pm$ | ${}^{0.09}_{0.10}$ |  3200$\pm$ | ${}^{150}_{100}$ |  0.17$\pm$ | ${}^{0.01}_{0.05}$ |  9$\pm$ | ${}^{6}_{3}$ |  0.18$\pm$ | ${}^{0.18}_{0.05}$ |  0.60$\pm$ | ${}^{0.20}_{0.28}$ |  0.93$\pm$ | ${}^{0.24}_{0.15}$ |  0.02$\pm$ | ${}^{0.20}_{0.02}$ |  7.25$\pm$ | ${}^{0.38}_{0.07}$ |  7.60$\pm$ | ${}^{0.35}_{0.33}$
H49 |  5.76$\pm$ | ${}^{0.09}_{0.04}$ |  44000$\pm$ | ${}^{3500}_{1500}$ |  3.85$\pm$ | ${}^{0.55}_{0.30}$ |  13.1$\pm$ | ${}^{0.3}_{0.5}$ |  -6.22$\pm$ | ${}^{0.10}_{0.33}$ |  3300$\pm$ | ${}^{150}_{200}$ |  0.12$\pm$ | ${}^{0.07}_{0.06}$ |  38$\pm$ | ${}^{13}_{15}$ |  0.18$\pm$ | ${}^{0.23}_{0.10}$ |  0.20$\pm$ | ${}^{0.40}_{0.15}$ |  0.72$\pm$ | ${}^{0.20}_{0.05}$ |  0.02$\pm$ | ${}^{0.17}_{0.02}$ |  7.12$\pm$ | ${}^{0.36}_{0.15}$ |  7.15$\pm$ | ${}^{0.88}_{0.30}$
H50 |  5.85$\pm$ | ${}^{0.05}_{0.04}$ |  47000$\pm$ | ${}^{2000}_{1250}$ |  4.15$\pm$ | ${}^{0.15}_{0.15}$ |  12.8$\pm$ | ${}^{0.3}_{0.3}$ |  -6.12$\pm$ | ${}^{0.03}_{0.23}$ |  2850$\pm$ | ${}^{200}_{75}$ |  0.20$\pm$ | ${}^{0.01}_{0.04}$ |  15$\pm$ | ${}^{32}_{2}$ |  0.07$\pm$ | ${}^{0.07}_{0.03}$ |  0.25$\pm$ | ${}^{0.17}_{0.07}$ |  0.70$\pm$ | ${}^{0.09}_{0.03}$ |  0.06$\pm$ | ${}^{0.11}_{0.04}$ |  7.83$\pm$ | ${}^{0.20}_{0.33}$ |  7.05$\pm$ | ${}^{0.40}_{0.15}$
H52 |  5.70$\pm$ | ${}^{0.04}_{0.04}$ |  45500$\pm$ | ${}^{1500}_{1500}$ |  4.05$\pm$ | ${}^{0.24}_{0.12}$ |  11.5$\pm$ | ${}^{0.3}_{0.2}$ |  -6.22$\pm$ | ${}^{0.10}_{0.20}$ |  2900$\pm$ | ${}^{150}_{100}$ |  0.20$\pm$ | ${}^{0.01}_{0.07}$ |  36$\pm$ | ${}^{15}_{19}$ |  0.08$\pm$ | ${}^{0.06}_{0.05}$ |  0.17$\pm$ | ${}^{0.15}_{0.10}$ |  0.75$\pm$ | ${}^{0.10}_{0.07}$ |  0.02$\pm$ | ${}^{0.07}_{0.02}$ |  7.45$\pm$ | ${}^{0.33}_{0.20}$ |  7.75$\pm$ | ${}^{0.70}_{0.45}$
H55 |  5.77$\pm$ | ${}^{0.08}_{0.04}$ |  47500$\pm$ | ${}^{3000}_{1500}$ |  3.95$\pm$ | ${}^{0.25}_{0.12}$ |  11.5$\pm$ | ${}^{0.3}_{0.4}$ |  -6.27$\pm$ | ${}^{0.14}_{0.20}$ |  3150$\pm$ | ${}^{75}_{200}$ |  0.10$\pm$ | ${}^{0.05}_{0.04}$ |  28$\pm$ | ${}^{16}_{10}$ |  0.22$\pm$ | ${}^{0.13}_{0.13}$ |  0.30$\pm$ | ${}^{0.21}_{0.15}$ |  0.72$\pm$ | ${}^{0.07}_{0.05}$ |  0.02$\pm$ | ${}^{0.15}_{0.02}$ |  7.75$\pm$ | ${}^{0.38}_{0.39}$ |  7.05$\pm$ | ${}^{0.47}_{0.20}$
H58 |  5.87$\pm$ | ${}^{0.08}_{0.06}$ |  47500$\pm$ | ${}^{3000}_{2250}$ |  4.40$\pm$ | ${}^{0.23}_{0.35}$ |  12.8$\pm$ | ${}^{0.4}_{0.4}$ |  -6.52$\pm$ | ${}^{0.18}_{0.25}$ |  3000$\pm$ | ${}^{100}_{200}$ |  0.15$\pm$ | ${}^{0.06}_{0.05}$ |  32$\pm$ | ${}^{17}_{14}$ |  0.11$\pm$ | ${}^{0.18}_{0.04}$ |  0.30$\pm$ | ${}^{0.15}_{0.23}$ |  0.70$\pm$ | ${}^{0.11}_{0.03}$ |  0.22$\pm$ | ${}^{0.10}_{0.22}$ |  8.18$\pm$ | ${}^{0.05}_{0.28}$ |  7.00$\pm$ | ${}^{0.70}_{0.15}$
H62 |  5.63$\pm$ | ${}^{0.08}_{0.07}$ |  45000$\pm$ | ${}^{3000}_{2250}$ |  4.05$\pm$ | ${}^{0.28}_{0.20}$ |  10.8$\pm$ | ${}^{0.3}_{0.4}$ |  -6.02$\pm$ | ${}^{0.10}_{0.15}$ |  3050$\pm$ | ${}^{125}_{150}$ |  0.20$\pm$ | ${}^{0.01}_{0.05}$ |  15$\pm$ | ${}^{28}_{4}$ |  0.05$\pm$ | ${}^{0.07}_{0.02}$ |  0.38$\pm$ | ${}^{0.10}_{0.19}$ |  0.72$\pm$ | ${}^{0.28}_{0.05}$ |  0.09$\pm$ | ${}^{0.14}_{0.06}$ |  7.15$\pm$ | ${}^{0.41}_{0.17}$ |  7.15$\pm$ | ${}^{0.93}_{0.30}$
H64 |  5.85$\pm$ | ${}^{0.05}_{0.07}$ |  46000$\pm$ | ${}^{2000}_{2500}$ |  4.30$\pm$ | ${}^{0.25}_{0.40}$ |  13.3$\pm$ | ${}^{0.4}_{0.3}$ |  -6.67$\pm$ | ${}^{0.25}_{0.25}$ |  2250$\pm$ | ${}^{250}_{150}$ |  0.20$\pm$ | ${}^{0.01}_{0.10}$ |  37$\pm$ | ${}^{13}_{14}$ |  0.02$\pm$ | ${}^{0.05}_{0.02}$ |  0.30$\pm$ | ${}^{0.10}_{0.15}$ |  0.70$\pm$ | ${}^{0.25}_{0.03}$ |  0.01$\pm$ | ${}^{0.14}_{0.01}$ |  8.05$\pm$ | ${}^{0.17}_{0.39}$ |  7.25$\pm$ | ${}^{0.88}_{0.40}$
H65 |  5.73$\pm$ | ${}^{0.06}_{0.04}$ |  42000$\pm$ | ${}^{2000}_{1500}$ |  3.85$\pm$ | ${}^{0.38}_{0.15}$ |  13.9$\pm$ | ${}^{0.3}_{0.4}$ |  -6.07$\pm$ | ${}^{0.15}_{0.28}$ |  2700$\pm$ | ${}^{150}_{150}$ |  0.14$\pm$ | ${}^{0.06}_{0.05}$ |  25$\pm$ | ${}^{15}_{12}$ |  0.10$\pm$ | ${}^{0.14}_{0.06}$ |  0.03$\pm$ | ${}^{0.28}_{0.03}$ |  0.72$\pm$ | ${}^{0.20}_{0.05}$ |  0.14$\pm$ | ${}^{0.14}_{0.07}$ |  7.30$\pm$ | ${}^{0.36}_{0.20}$ |  7.45$\pm$ | ${}^{0.68}_{0.55}$
H66 |  5.66$\pm$ | ${}^{0.04}_{0.04}$ |  47500$\pm$ | ${}^{1500}_{1500}$ |  4.10$\pm$ | ${}^{0.20}_{0.10}$ |  10.1$\pm$ | ${}^{0.2}_{0.2}$ |  -6.22$\pm$ | ${}^{0.10}_{0.09}$ |  2700$\pm$ | ${}^{125}_{100}$ |  0.20$\pm$ | ${}^{0.01}_{0.05}$ |  42$\pm$ | ${}^{7}_{5}$ |  0.07$\pm$ | ${}^{0.04}_{0.03}$ |  0.10$\pm$ | ${}^{0.14}_{0.05}$ |  0.70$\pm$ | ${}^{0.09}_{0.03}$ |  0.01$\pm$ | ${}^{0.05}_{0.01}$ |  7.70$\pm$ | ${}^{0.04}_{0.35}$ |  7.00$\pm$ | ${}^{0.95}_{0.07}$
H68 |  5.68$\pm$ | ${}^{0.07}_{0.07}$ |  42000$\pm$ | ${}^{2250}_{2500}$ |  3.95$\pm$ | ${}^{0.30}_{0.17}$ |  13.2$\pm$ | ${}^{0.4}_{0.4}$ |  -6.62$\pm$ | ${}^{0.18}_{0.30}$ |  2650$\pm$ | ${}^{350}_{250}$ |  0.18$\pm$ | ${}^{0.02}_{0.07}$ |  30$\pm$ | ${}^{18}_{11}$ |  0.01$\pm$ | ${}^{0.06}_{0.01}$ |  0.05$\pm$ | ${}^{0.11}_{0.05}$ |  0.75$\pm$ | ${}^{0.11}_{0.07}$ |  0.31$\pm$ | ${}^{0.05}_{0.14}$ |  8.15$\pm$ | ${}^{0.05}_{0.56}$ |  8.15$\pm$ | ${}^{0.40}_{1.25}$
H69 |  5.47$\pm$ | ${}^{0.06}_{0.04}$ |  41000$\pm$ | ${}^{2000}_{1500}$ |  4.15$\pm$ | ${}^{0.23}_{0.30}$ |  10.9$\pm$ | ${}^{0.3}_{0.3}$ |  -6.83$\pm$ | ${}^{0.40}_{0.20}$ |  3050$\pm$ | ${}^{200}_{250}$ |  0.20$\pm$ | ${}^{0.01}_{0.10}$ |  50$\pm$ | ${}^{1}_{35}$ |  0.04$\pm$ | ${}^{0.04}_{0.04}$ |  0.05$\pm$ | ${}^{0.23}_{0.05}$ |  0.78$\pm$ | ${}^{0.19}_{0.10}$ |  0.35$\pm$ | ${}^{0.03}_{0.21}$ |  7.58$\pm$ | ${}^{0.31}_{0.40}$ |  7.00$\pm$ | ${}^{1.35}_{0.15}$
H70 |  5.71$\pm$ | ${}^{0.06}_{0.05}$ |  45500$\pm$ | ${}^{2250}_{2000}$ |  4.15$\pm$ | ${}^{0.50}_{0.17}$ |  11.7$\pm$ | ${}^{0.3}_{0.3}$ |  -6.32$\pm$ | ${}^{0.13}_{0.15}$ |  2600$\pm$ | ${}^{300}_{125}$ |  0.16$\pm$ | ${}^{0.04}_{0.13}$ |  49$\pm$ | ${}^{2}_{18}$ |  0.09$\pm$ | ${}^{0.11}_{0.06}$ |  0.03$\pm$ | ${}^{0.17}_{0.03}$ |  0.70$\pm$ | ${}^{0.23}_{0.03}$ |  0.02$\pm$ | ${}^{0.15}_{0.02}$ |  7.15$\pm$ | ${}^{0.36}_{0.17}$ |  7.05$\pm$ | ${}^{0.85}_{0.20}$
H71 |  5.47$\pm$ | ${}^{0.07}_{0.07}$ |  45000$\pm$ | ${}^{2500}_{2500}$ |  3.90$\pm$ | ${}^{0.20}_{0.20}$ |  9.1$\pm$ | ${}^{0.3}_{0.3}$ |  -6.57$\pm$ | ${}^{0.10}_{0.50}$ |  2650$\pm$ | ${}^{150}_{150}$ |  0.15$\pm$ | ${}^{0.04}_{0.05}$ |  21$\pm$ | ${}^{29}_{8}$ |  0.05$\pm$ | ${}^{0.08}_{0.03}$ |  0.23$\pm$ | ${}^{0.09}_{0.17}$ |  0.70$\pm$ | ${}^{0.21}_{0.03}$ |  0.01$\pm$ | ${}^{0.24}_{0.01}$ |  7.83$\pm$ | ${}^{0.40}_{0.39}$ |  7.75$\pm$ | ${}^{0.65}_{0.60}$
H75 |  5.47$\pm$ | ${}^{0.05}_{0.05}$ |  46000$\pm$ | ${}^{2000}_{2000}$ |  4.45$\pm$ | ${}^{0.28}_{0.25}$ |  8.6$\pm$ | ${}^{0.2}_{0.2}$ |  -6.47$\pm$ | ${}^{0.10}_{0.35}$ |  2700$\pm$ | ${}^{200}_{150}$ |  0.15$\pm$ | ${}^{0.05}_{0.07}$ |  26$\pm$ | ${}^{25}_{10}$ |  0.11$\pm$ | ${}^{0.24}_{0.06}$ |  0.05$\pm$ | ${}^{0.23}_{0.05}$ |  0.72$\pm$ | ${}^{0.12}_{0.05}$ |  0.04$\pm$ | ${}^{0.13}_{0.04}$ |  7.53$\pm$ | ${}^{0.35}_{0.40}$ |  7.35$\pm$ | ${}^{0.82}_{0.45}$
H78 |  5.52$\pm$ | ${}^{0.03}_{0.11}$ |  46000$\pm$ | ${}^{1000}_{4000}$ |  4.00$\pm$ | ${}^{0.20}_{0.20}$ |  9.1$\pm$ | ${}^{0.4}_{0.2}$ |  -6.67$\pm$ | ${}^{0.30}_{0.18}$ |  2700$\pm$ | ${}^{125}_{150}$ |  0.15$\pm$ | ${}^{0.05}_{0.07}$ |  44$\pm$ | ${}^{4}_{29}$ |  0.16$\pm$ | ${}^{0.14}_{0.06}$ |  0.07$\pm$ | ${}^{0.12}_{0.07}$ |  0.70$\pm$ | ${}^{0.11}_{0.03}$ |  0.06$\pm$ | ${}^{0.21}_{0.03}$ |  7.58$\pm$ | ${}^{0.23}_{0.42}$ |  7.25$\pm$ | ${}^{0.33}_{0.40}$
H80 |  5.12$\pm$ | ${}^{0.03}_{0.06}$ |  35500$\pm$ | ${}^{1000}_{1750}$ |  3.90$\pm$ | ${}^{0.25}_{0.25}$ |  9.6$\pm$ | ${}^{0.3}_{0.2}$ |  -8.29$\pm$ | ${}^{0.20}_{0.48}$ |  2400$\pm$ | ${}^{300}_{250}$ |  0.07$\pm$ | ${}^{0.09}_{0.08}$ |  32$\pm$ | ${}^{19}_{23}$ |  0.08$\pm$ | ${}^{0.21}_{0.08}$ |  0.03$\pm$ | ${}^{0.42}_{0.03}$ |  0.70$\pm$ | ${}^{0.31}_{0.03}$ |  0.34$\pm$ | ${}^{0.07}_{0.17}$ |  8.18$\pm$ | ${}^{0.05}_{0.39}$ |  7.60$\pm$ | ${}^{0.70}_{0.45}$
H86 |  5.38$\pm$ | ${}^{0.06}_{0.03}$ |  45500$\pm$ | ${}^{2500}_{1000}$ |  3.85$\pm$ | ${}^{0.23}_{0.20}$ |  8.0$\pm$ | ${}^{0.2}_{0.3}$ |  -6.35$\pm$ | ${}^{0.05}_{0.19}$ |  2750$\pm$ | ${}^{200}_{100}$ |  0.15$\pm$ | ${}^{0.05}_{0.07}$ |  26$\pm$ | ${}^{22}_{7}$ |  0.15$\pm$ | ${}^{0.16}_{0.04}$ |  0.15$\pm$ | ${}^{0.16}_{0.10}$ |  0.72$\pm$ | ${}^{0.28}_{0.03}$ |  0.03$\pm$ | ${}^{0.06}_{0.03}$ |  7.00$\pm$ | ${}^{0.35}_{0.03}$ |  7.15$\pm$ | ${}^{0.57}_{0.30}$
H90 |  5.35$\pm$ | ${}^{0.11}_{0.06}$ |  42000$\pm$ | ${}^{4000}_{2000}$ |  3.95$\pm$ | ${}^{0.40}_{0.12}$ |  9.0$\pm$ | ${}^{0.3}_{0.4}$ |  -6.62$\pm$ | ${}^{0.20}_{0.40}$ |  2800$\pm$ | ${}^{300}_{250}$ |  0.20$\pm$ | ${}^{0.01}_{0.10}$ |  22$\pm$ | ${}^{28}_{16}$ |  0.04$\pm$ | ${}^{0.17}_{0.03}$ |  0.07$\pm$ | ${}^{0.19}_{0.07}$ |  0.70$\pm$ | ${}^{0.21}_{0.03}$ |  0.21$\pm$ | ${}^{0.15}_{0.17}$ |  7.62$\pm$ | ${}^{0.47}_{0.33}$ |  7.40$\pm$ | ${}^{0.60}_{0.55}$
H92 |  5.24$\pm$ | ${}^{0.04}_{0.03}$ |  40000$\pm$ | ${}^{1250}_{750}$ |  4.15$\pm$ | ${}^{0.17}_{0.30}$ |  8.8$\pm$ | ${}^{0.2}_{0.2}$ |  -7.69$\pm$ | ${}^{0.35}_{0.30}$ |  2700$\pm$ | ${}^{150}_{300}$ |  0.20$\pm$ | ${}^{0.01}_{0.10}$ |  50$\pm$ | ${}^{1}_{26}$ |  0.09$\pm$ | ${}^{0.20}_{0.04}$ |  0.03$\pm$ | ${}^{0.30}_{0.03}$ |  0.78$\pm$ | ${}^{0.17}_{0.10}$ |  0.19$\pm$ | ${}^{0.16}_{0.19}$ |  8.18$\pm$ | ${}^{0.05}_{0.44}$ |  7.75$\pm$ | ${}^{0.65}_{0.85}$
H94 |  5.38$\pm$ | ${}^{0.08}_{0.05}$ |  44500$\pm$ | ${}^{3000}_{1750}$ |  3.95$\pm$ | ${}^{0.28}_{0.20}$ |  8.3$\pm$ | ${}^{0.2}_{0.3}$ |  -7.19$\pm$ | ${}^{0.25}_{0.25}$ |  2550$\pm$ | ${}^{150}_{150}$ |  0.18$\pm$ | ${}^{0.03}_{0.09}$ |  43$\pm$ | ${}^{8}_{16}$ |  0.12$\pm$ | ${}^{0.22}_{0.05}$ |  0.05$\pm$ | ${}^{0.41}_{0.05}$ |  0.75$\pm$ | ${}^{0.17}_{0.07}$ |  0.03$\pm$ | ${}^{0.18}_{0.03}$ |  8.20$\pm$ | ${}^{0.03}_{0.28}$ |  6.95$\pm$ | ${}^{0.70}_{0.10}$
H143 |  5.16$\pm$ | ${}^{0.06}_{0.10}$ |  40000$\pm$ | ${}^{2000}_{3000}$ |  3.75$\pm$ | ${}^{0.40}_{0.20}$ |  8.0$\pm$ | ${}^{0.4}_{0.2}$ |  -7.79$\pm$ | ${}^{0.40}_{0.25}$ |  1950$\pm$ | ${}^{250}_{250}$ |  0.20$\pm$ | ${}^{0.01}_{0.20}$ |  18$\pm$ | ${}^{31}_{9}$ |  0.10$\pm$ | ${}^{0.20}_{0.09}$ |  0.07$\pm$ | ${}^{0.21}_{0.07}$ |  0.72$\pm$ | ${}^{0.31}_{0.05}$ |  0.01$\pm$ | ${}^{0.34}_{0.01}$ |  8.20$\pm$ | ${}^{0.03}_{0.53}$ |  6.95$\pm$ | ${}^{1.25}_{0.10}$
a) Note that when the error bar reaches the edge of the parameter space, the
best fit value is in fact an upper or lower limit. Figure 8 shows when this is
the case.
b) Run with 14 free parameters. Values for the helium and oxygen surface
abundance can be found in Table 6.
c) Formal uncertainties from the Kiwi-GA run with fixed, estimated $T_{\rm
eff}$ of 50 kK (see text). In reality uncertainties on all parameters for this
source are larger due to the uncertain $T_{\rm eff}$.

Table 5: Best fit parameters and $1\sigma$ error barsa) for the optical + UV fits of the 17 stars where we fitted 6 free parameters.  Source | $\log L/{\rm L}_{\odot}$ |  $T_{\rm eff}$(K) | $\log g$ | $R_{\star}/{\rm R}_{\odot}$ | $\log\dot{M}$ |  $v_{\infty}$ (km/s) | $f_{\rm cl}$ | $\beta$
---|---|---|---|---|---|---|---|---
H73 |  5.18$\pm$ | ${}^{0.09}_{0.12}$ |  32500$\pm$ | ${}^{2500}_{3000}$ |  4.10$\pm$ | ${}^{0.35}_{0.23}$ |  12.4$\pm$ | ${}^{0.7}_{0.5}$ |  -8.48$\pm$ | ${}^{0.28}_{0.95}$ |  4000$\pm$ | ${}^{600}_{2800}$ |  1$\pm$ | ${}^{25}_{1}$ |  0.90$\pm$ | ${}^{0.85}_{0.20}$
H108 |  4.89$\pm$ | ${}^{0.10}_{0.18}$ |  40000$\pm$ | ${}^{3500}_{5000}$ |  4.05$\pm$ | ${}^{0.42}_{0.38}$ |  5.9$\pm$ | ${}^{0.5}_{0.2}$ |  -8.09$\pm$ | ${}^{0.35}_{0.50}$ |  1350$\pm$ | ${}^{525}_{350}$ |  31$\pm$ | ${}^{20}_{25}$ |  0.70$\pm$ | ${}^{0.61}_{0.03}$
H112 |  5.11$\pm$ | ${}^{0.11}_{0.12}$ |  35500$\pm$ | ${}^{3250}_{3250}$ |  4.10$\pm$ | ${}^{0.47}_{0.38}$ |  9.5$\pm$ | ${}^{0.5}_{0.4}$ |  -8.40$\pm$ | ${}^{0.40}_{0.40}$ |  1300$\pm$ | ${}^{875}_{425}$ |  37$\pm$ | ${}^{14}_{31}$ |  0.70$\pm$ | ${}^{0.61}_{0.03}$
H114 |  5.22$\pm$ | ${}^{0.07}_{0.11}$ |  44500$\pm$ | ${}^{2500}_{4000}$ |  4.20$\pm$ | ${}^{0.38}_{0.35}$ |  6.9$\pm$ | ${}^{0.3}_{0.2}$ |  -7.48$\pm$ | ${}^{0.30}_{0.40}$ |  2100$\pm$ | ${}^{725}_{400}$ |  46$\pm$ | ${}^{5}_{32}$ |  0.70$\pm$ | ${}^{0.59}_{0.03}$
H116 |  4.87$\pm$ | ${}^{0.13}_{0.12}$ |  36000$\pm$ | ${}^{4000}_{3250}$ |  4.05$\pm$ | ${}^{0.40}_{0.30}$ |  7.0$\pm$ | ${}^{0.4}_{0.4}$ |  -8.68$\pm$ | ${}^{0.50}_{0.75}$ |  1900$\pm$ | ${}^{1350}_{600}$ |  45$\pm$ | ${}^{6}_{41}$ |  0.70$\pm$ | ${}^{0.66}_{0.03}$
H120 |  4.85$\pm$ | ${}^{0.05}_{0.27}$ |  40000$\pm$ | ${}^{1500}_{8000}$ |  4.45$\pm$ | ${}^{0.23}_{0.65}$ |  5.6$\pm$ | ${}^{0.7}_{0.1}$ |  -7.99$\pm$ | ${}^{0.30}_{0.50}$ |  1600$\pm$ | ${}^{975}_{300}$ |  37$\pm$ | ${}^{13}_{35}$ |  0.72$\pm$ | ${}^{0.42}_{0.05}$
H121 |  4.79$\pm$ | ${}^{0.10}_{0.14}$ |  33500$\pm$ | ${}^{3000}_{3500}$ |  4.05$\pm$ | ${}^{0.30}_{0.35}$ |  7.4$\pm$ | ${}^{0.5}_{0.3}$ |  -8.42$\pm$ | ${}^{0.23}_{1.10}$ |  4300$\pm$ | ${}^{300}_{3400}$ |  1$\pm$ | ${}^{49}_{1}$ |  0.75$\pm$ | ${}^{0.64}_{0.07}$
H123 |  4.93$\pm$ | ${}^{0.06}_{0.04}$ |  40000$\pm$ | ${}^{2000}_{1250}$ |  4.00$\pm$ | ${}^{0.30}_{0.25}$ |  6.1$\pm$ | ${}^{0.1}_{0.2}$ |  -7.69$\pm$ | ${}^{0.20}_{0.30}$ |  2000$\pm$ | ${}^{350}_{250}$ |  40$\pm$ | ${}^{11}_{23}$ |  0.70$\pm$ | ${}^{0.26}_{0.03}$
H129 |  4.48$\pm$ | ${}^{0.12}_{0.10}$ |  42000$\pm$ | ${}^{4500}_{3250}$ |  4.25$\pm$ | ${}^{0.55}_{0.40}$ |  3.3$\pm$ | ${}^{0.2}_{0.2}$ |  -8.09$\pm$ | ${}^{0.35}_{0.75}$ |  1300$\pm$ | ${}^{1350}_{300}$ |  9$\pm$ | ${}^{41}_{7}$ |  0.70$\pm$ | ${}^{0.66}_{0.03}$
H132 |  5.04$\pm$ | ${}^{0.10}_{0.07}$ |  40000$\pm$ | ${}^{3500}_{2500}$ |  4.20$\pm$ | ${}^{0.40}_{0.30}$ |  6.9$\pm$ | ${}^{0.2}_{0.3}$ |  -8.34$\pm$ | ${}^{0.55}_{0.65}$ |  1300$\pm$ | ${}^{1650}_{500}$ |  43$\pm$ | ${}^{8}_{37}$ |  0.70$\pm$ | ${}^{1.25}_{0.03}$
H134 |  4.75$\pm$ | ${}^{0.16}_{0.10}$ |  36000$\pm$ | ${}^{5000}_{2750}$ |  4.05$\pm$ | ${}^{0.42}_{0.35}$ |  6.1$\pm$ | ${}^{0.3}_{0.4}$ |  -8.09$\pm$ | ${}^{0.48}_{0.38}$ |  2000$\pm$ | ${}^{750}_{500}$ |  36$\pm$ | ${}^{15}_{25}$ |  0.70$\pm$ | ${}^{0.54}_{0.03}$
H135 |  4.71$\pm$ | ${}^{0.15}_{0.17}$ |  29000$\pm$ | ${}^{4000}_{4000}$ |  3.90$\pm$ | ${}^{0.30}_{0.35}$ |  9.0$\pm$ | ${}^{0.7}_{0.5}$ |  -8.25$\pm$ | ${}^{0.60}_{0.63}$ |  400$\pm$ | ${}^{1950}_{300}$ |  41$\pm$ | ${}^{10}_{35}$ |  1.30$\pm$ | ${}^{0.72}_{0.60}$
H139 |  4.91$\pm$ | ${}^{0.06}_{0.10}$ |  40000$\pm$ | ${}^{2000}_{3000}$ |  4.15$\pm$ | ${}^{0.40}_{0.30}$ |  6.0$\pm$ | ${}^{0.3}_{0.2}$ |  -8.20$\pm$ | ${}^{0.23}_{0.63}$ |  1300$\pm$ | ${}^{500}_{300}$ |  42$\pm$ | ${}^{9}_{32}$ |  0.70$\pm$ | ${}^{0.66}_{0.03}$
H141 |  4.69$\pm$ | ${}^{0.19}_{0.12}$ |  30000$\pm$ | ${}^{5500}_{3000}$ |  3.80$\pm$ | ${}^{0.70}_{0.40}$ |  8.3$\pm$ | ${}^{0.5}_{0.6}$ |  -7.56$\pm$ | ${}^{0.60}_{0.65}$ |  4300$\pm$ | ${}^{300}_{2500}$ |  47$\pm$ | ${}^{4}_{42}$ |  0.70$\pm$ | ${}^{0.66}_{0.03}$
H159 |  4.71$\pm$ | ${}^{0.07}_{0.10}$ |  31000$\pm$ | ${}^{2000}_{2500}$ |  3.95$\pm$ | ${}^{0.35}_{0.30}$ |  7.9$\pm$ | ${}^{0.4}_{0.3}$ |  -7.71$\pm$ | ${}^{0.35}_{0.35}$ |  4100$\pm$ | ${}^{500}_{1500}$ |  42$\pm$ | ${}^{9}_{25}$ |  1.30$\pm$ | ${}^{0.72}_{0.57}$
H162 |  4.81$\pm$ | ${}^{0.15}_{0.25}$ |  36000$\pm$ | ${}^{4750}_{6500}$ |  4.30$\pm$ | ${}^{0.75}_{0.65}$ |  6.6$\pm$ | ${}^{0.8}_{0.4}$ |  -8.47$\pm$ | ${}^{1.05}_{1.05}$ |  100$\pm$ | ${}^{4400}_{100}$ |  46$\pm$ | ${}^{5}_{46}$ |  1.02$\pm$ | ${}^{1.00}_{0.35}$
H173 |  4.46$\pm$ | ${}^{0.18}_{0.23}$ |  26500$\pm$ | ${}^{4500}_{4250}$ |  3.85$\pm$ | ${}^{0.50}_{0.42}$ |  8.1$\pm$ | ${}^{0.9}_{0.6}$ |  -8.15$\pm$ | ${}^{0.43}_{0.90}$ |  3700$\pm$ | ${}^{900}_{3400}$ |  49$\pm$ | ${}^{2}_{31}$ |  0.70$\pm$ | ${}^{1.32}_{0.03}$
a) Note that when the error bar reaches the edge of the parameter space, the
best fit value is in fact an upper or lower limit. Figure 9 shows when this is
the case.
Table 6: Best fit parameters and $1\sigma$ error bars for the oxygen and helium abundances of the optical + UV runs of the WNh stars. Only for these stars oxygen and helium abundance were fitted in the optical + UV runs. Source | $x_{\rm He}$ |  $\log n_{\rm O}/n_{\rm H}$+12
---|---|---
R136a1 | 0.22$\pm$ | $0.05$ | 8.30$\pm$ | ${}^{0.05}_{0.65}$
R136a2 | 0.39$\pm$ | ${}^{0.11}_{0.07}$ | 7.80$\pm$ | ${}^{0.50}_{0.20}$
R136a3 | 0.37$\pm$ | $0.10$ | 7.45$\pm$ | ${}^{0.70}_{0.75}$

###  4.2 Stellar mass & age

In order to derive the evolutionary mass $M_{\rm evol}$, the initial mass
$M_{\rm ini}$ and the age $\tau$ we use Bonnsai141414The BONNSAI web-service
is available at https://www.astro.uni-bonn.de/stars/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_{\rm eff}$ cannot be reproduced given the observed $\log L/{\rm
L}_{\odot}$ and $\log g$. For this star, we deemed our derived luminosity
measurement unreliable (Section 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_{\rm evol}$ as an input, but check
whether the results change drastically upon inclusion/exclusion of R136b and
H30, which is not the case. The derived ages and masses can be found in Table
H.2. We cross-check our Bonnsai output with that of Bestenlehner et al. (2020)
and find generally good agreement, see Appendix D for details. In the
remainder of the paper we will 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 $\varv_{\infty}$ as a
free parameter in the optical + UV fit. For 46 sources we were able to
accurately constrain $\varv_{\infty}$, albeit with large uncertainties for the
stars with lower mass-loss rates (see Figs. 8 and 9). For 3 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 $\lambda\lambda$1548-1551 is not good:
the saturated absorption edge of the best model for these stars extends about
400 km s-1 more to the blue as the absorption edge we see in the data. For the
remaining 7 sources (H73, H121, H135, H141, H159, H162, H173) the wind lines
crucial for determining $\varv_{\infty}$, especially C iv
$\lambda\lambda$1548-1551, turn out to be too weak to get a constraint: the
$\chi^{2}$ distribution for $\varv_{\infty}$ was flat. In these cases, while
we do find some best fit value $\varv_{\infty}$, this value is not meaningful
and we regard $\varv_{\infty}$ as unconstrained.

###  4.4 Wind acceleration parameter $\beta$

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

Figure 11: Comparison of parameters from the different test runs, fitting in
each instance the spectrum of H35. For each run and each parameter we indicate
the best fit value (blue dots) and 1$\sigma$ and 2$\sigma$ error bars (dark
and light yellow). In the first column we show the parameters of our ‘fiducial
run’: this is the setup as used throughout the paper. The other columns show
parameters of runs were we changed the setup, one aspect at the time:
‘Fiducial (redo)’ – different initial random population of models,
‘$\varv_{\rm micro}5$ km s-1’ and ‘$\varv_{\rm micro}15$ km s-1’ – assumed
value for $\varv_{\rm micro}$ to 5 and 15 km s-1, respectively, instead of 10
km s-1, ‘$\varv_{\rm micro}$ free’ – $\varv_{\rm micro}$ a free parameter,
‘Oxygen free’ – oxygen abundance a free parameter, ‘No O lines’ – exclude both
O iv $\lambda$1340 and O v $\lambda$1371 from the fitting, ‘No X-rays ’ – do
not include any X-rays, ‘X-rays free’ – $f_{\rm X}$ a free parameter,
‘$\varv_{\rm cl,max}$’ – assume $\varv_{\rm cl,max}=2\varv_{\rm cl,start}$
instead of $\varv_{\rm cl,max}/\varv_{\infty}={\rm max}(0.3,\varv_{\rm
cl,start}/\varv_{\infty})$, ‘Min $\beta=0.5$’ - set lower limit of $\beta$ to
0.5 instead of 0.7 and ‘Only UV’ – only fit the UV spectra. For reference, the
best fit value and $2\sigma$ error region of the fiducial run are shown in
blue throughout all columns.

###  4.5 Onset of clumping

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

###  4.6 Ionising flux

We derive H, He i and He ii ionising fluxes $Q_{0}$, $Q_{1}$, $Q_{2}$ for each
star151515Here, by convention, $Q_{x}=q_{x}4\pi R_{\star}^{2}$, with $q_{x}$
the ionising radiation (number of photons) per unit surface area per second
and $x\in\\{1,2,3\\}$. based on the best fitting model (Table H.2). 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 $\chi^{2}$ value of each model, as we
do for all other parameters161616Uncertainties on ionising fluxes are now a
standard output of Kiwi-GA, but this functionality was implemented only after
we had done all fits. Since carrying out a fit is computationally expensive,
we decided to do a rerun with the new functionality only for one star.. From
this we find $1\sigma$ uncertainties on the derived ionising fluxes to be,
approximately, $0.07$ dex for $\log\ Q_{0}$, $0.1$ dex for $\log\ Q_{1}$,
$0.4$ dex for $\log\ Q_{2}$. 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 H.2
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
source171717We consider H35 a ‘typical source’ because O3 V is the most common
spectral type in our sample and the data quality of H35 is representative:
there are sources with higher S/N, but the H35 data is good enough so that we
can constrain each of our 12 free parameters. , 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\sigma$ and $2\sigma$ regions
are similar and the best fitting parameters of each runs lie in the $1\sigma$
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 $\varv_{\rm micro}$ seems to have the largest effect on the
derived parameters, especially $T_{\rm eff}$. Changing the value from our
fiducial fixed value of 10 km s-1 to 5 km s-1 does not have much effect, but
changing it to either a fixed value of 15 km s-1 or leaving it as a free
variable results in a best fit value of $T_{\rm eff}$ that is $\sim 2000$ K
higher than that of our fiducial run, just on the edge of the $2\sigma$ error
bars. From the data that we have we cannot determine what is the actual value
of $\varv_{\rm micro}$ and thus of $T_{\rm eff}$: the run with $\varv_{\rm
micro}$ as a free parameter resulted in a velocity exceeding the typical sound
speed in the atmospheres of these stars ($\varv_{\rm micro}$ = 30 km s-1) and
thus seemed not reliable (though for a recent finding, see Schultz et al.
2022). Apart from changes in $T_{\rm eff}$ we note that also abundances change
when we assume a different value for $\varv_{\rm micro}$, the largest change
being seen in the carbon abundance when $\varv_{\rm micro}$ is lowered to 5 km
s-1. This is expected as $\varv_{\rm 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_{\rm 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_{\rm 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_{\rm ic}$ and $f_{\rm 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\
\dot{M}\lesssim-5.8$ we find an average wind acceleration parameter of
$\beta=0.72\pm 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 $\beta$ 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 $\beta=0.7$ remains in the 2$\sigma$ error range. Fig. 11 shows
that for H35 with the lower best fit value of $\beta$ the other parameters do
not change significantly given the uncertainties.

Apart from the aspects discussed so far, we do also change the prescription
for $\varv_{\rm 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_{\rm eff}$ and abundances.

##  5 Discussion

We discuss our results in the context of theoretical predictions and
evolutionary models. In Section 5.1 we compare the mass-loss rates that we
obtained to the predictions of Vink et al. (2000, 2001); Krtička & 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 Section 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 (Section 5.3.1), we add to their
discussion based on additional clues we can get from the abundances based on
UV spectroscopy (Section 5.3.2). Furthermore, Section 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 (Section
5.4).

###  5.1 Mass-loss rates, wind momentum and terminal velocities

Figure 12: Mass-loss rates from the optical + UV fits (dark blue solid
diamonds) compared to the mass-loss predictions of Vink et al. (2001, light
blue triangles), Krtička & Kubát (2018, red solid line) and Björklund et al.
(2021, orange dot-dashed line). The prescription of Vink et al. (2001) depends
on more parameters than luminosity; a linear fit through the individual points
(light blue dashed line) is plotted to guide the eye. Thin dotted lines in
corresponding colours show the extrapolation of each prescription beyond the
coverage of their respective model grids. For reference, a linear fit through
the data points is shown (thin dashed darkblue line). For all mass-loss
prescriptions we assess the goodness of fit ($\chi^{2}$-values, top) for three
ranges of luminosity. Figure 13: As Fig. 12 but now for the modified wind
momentum. In grey circles we show the values or upper limits of Mokiem et al.
(2007); Bestenlehner et al. (2014); Ramírez-Agudelo et al. (2017); Sabín-
Sanjulián et al. (2017), which we all lowered by 0.7 dex, assuming $f_{\rm
cl}=25$.

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 ionisation 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_{\odot}=0.019$ and $Z_{\odot}=0.013$, respectively. For a mass-loss
prediction at $Z=0.5\leavevmode\nobreak\ Z_{\odot}$ 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\leavevmode\nobreak\ Z_{\odot}$) 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/{\rm L}_{\odot}$ $<5.3$), intermediate
($5.3\leq$ $\log L/{\rm L}_{\odot}$ $\leq 6.2$) and high ($\log L/{\rm
L}_{\odot}$ $>6.2$). Results are shown in Fig. 12. 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/{\rm L}_{\odot}$ $=6.00$, $\log L/{\rm
L}_{\odot}$ $=6.07$ and $\log L/{\rm L}_{\odot}$ $=6.25$, respectively. Note
furthermore, that for four sources (H73, H116, H121, H162), for which we
derive mass-loss rates in the range log $\dot{M}$ $=-8.68$ to $-8.42$, the
$1\sigma$ error bars range (nearly) to the edge of our parameter space, log
$\dot{M}$ $=-9.50$. In our fitting and with the derivation of the $\chi^{2}$
values, we treat the derived values and uncertainties the same as those of the
other points, even though these are technically upper limits. Removing these 4
points from the fit does not alter the results significantly. The same holds
for fitting and $\chi^{2}$ 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 $\varv_{\infty}$.

Figure 14: Terminal velocity $\varv_{\infty}$ (top) and the ratio
$\varv_{\infty}/\varv_{\rm esc,eff}$ (bottom) against $\log L/{\rm L}_{\odot}$
(left) and $T_{\rm eff}$ (right) for the R136 sample (solid dark blue
diamonds) compared to predicted values (yellow squares, red stars, and green
pentagons Björklund et al. 2021, Krtička & Kubát 2018, and Vink & Sander 2021,
respectively). Light purple triangles denote R136 sources for which we could
not derive $\varv_{\infty}$ values (see Section 4.3), and in the bottom plot
also includes the WNh stars, where we cannot be too confident about
$\varv_{\rm esc,eff}$. Grey circles around the points of the WNh stars allows
to distinguish them from the O-type stars. Dark and light coloured error bars
denote $1\sigma$ and $2\sigma$ uncertainties, respectively. The lines, plotted
in different styles (see legend) show linear fits to both the observed (dark
blue) and predicted values (red, orange and green). In the bottom panel, the
black dashed line shows the empirically derived $\varv_{\infty}/\varv_{\rm
esc,eff}$ = 2.6 (Lamers et al. 1995). The grey shaded regions in the
luminosity plots correspond to those in Fig. 12. Note that we do not show any
points of models with ${T}_{\rm eff}=35000$ K for the Vink & Sander 2021
predictions; in this regime the predicted terminal velocities largely exceed
the velocity scale of this plot. The linear fit through the Vink & Sander 2021
predictions also excludes these points.

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/{\rm L}_{\odot}$ $>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^{2}_{\rm red}$. 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, 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 Krtič18 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_{\rm mom}=\dot{M}\varv_{\infty}\sqrt{R_{\star}/{\rm
R}_{\odot}}$ (with $\dot{M}$ and $\varv_{\infty}$ 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_{\rm 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 $\varv_{\infty}$ and
$R_{\star}$. The relation plotted in Fig. 13 is corrected for lowered
$\dot{M}$ and $\varv_{\infty}$ 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_{\rm 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/{\rm L}_{\odot}$ and $T_{\rm 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_{\rm eff}$ and $\log L/{\rm
L}_{\odot}$, 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 observations181818In Fig. 14 we plot the
$\varv_{\infty}$ of the $Z=1/3\ {\rm Z}_{\odot}$ models, multiplied by a
factor $(0.5/(1/3))^{0.19}=1.08$ so that they correspond to our assumed LMC
metallicity of $Z=0.5\ {\rm Z}_{\odot}$.; 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 35000 K,
where they predict extremely high velocities. In Figure 14 we do not show
their predictions for $T_{\rm eff}$ $=35000$ 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_{\rm eff}$ in
the grids of Björk21, Krtič18 and Vink & Sander 2021 is 45000 K, while about
50% of our sample has a $T_{\rm eff}$ higher than that. The Vink01 grid
extends further to 50000 K, however they do not predict terminal velocities.

We will 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_{\rm mom}=\log\ D_{\rm 0}+x\log\ L/{\rm 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\pm 0.11$ and $\log\ D_{\rm 0}=17.05\pm
0.65$. A comparison to other LMC studies shows that this is relatively high:
observations of Mokiem et al. (2007, grey circles in Fig. 13); Bestenlehner et
al. (2014, grey circles in Fig. 13); Ramírez-Agudelo et al. (2017, grey
circles in Fig. 13); Sabín-Sanjulián et al. (2017, grey circles in Fig. 13)
give slopes in the range 1.45-1.87. Note however that not all these analyses
take into account uncertainties in $\ D_{\rm 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$\alpha$ and He ii $\lambda$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 $\alpha_{\rm eff}$ in (modified) CAK-
theory (Castor et al. 1975, with subsequent modifications by Abbott 1982;
Pauldrach et al. 1986): $x=1/\alpha_{\rm eff}$. Here, $\alpha_{\rm eff}$
captures both $\alpha$, the slope of the line strength distribution, as well
as the force multiplier $\delta$, that accounts for the ionisation state of
the wind in an approximate way (Abbott 1982), that is, $x^{-1}=\alpha_{\rm
eff}=\alpha-\delta$ in Eq. 12. Our slope of $x=2.00\pm 0.11$ then translates
into a mean value of $\alpha_{\rm eff}=0.50\pm 0.03$ or $\alpha=0.60\pm 0.03$
if we assume the typical value of 0.1 for $\delta$ (Abbott 1982; Puls et al.
2008). This is in line with typical values expected for O-stars,
$\alpha\approx 0.5-0.6$ (Abbott 1982; Pauldrach et al. 1994; Puls et al. 2000,
2008).