© The Authors 2024

1. Introduction

An important process during the life of massive stars are stellar winds, which are driven by radiation that is intercepted by a multitude of atmospheric spectral lines (Castor et al. 1975). Over time these radiation driven winds may cause stars to lose a sizeable fraction of their initial mass and angular momentum. This significantly influences the evolution of massive stars and ultimately the nature and characteristics of its end of life products (e.g. Heger et al. 2003; Langer 2012).

While studying luminous B stars that experience such line-driven winds, Pauldrach & Puls (1990) and Lamers et al. (1995) noticed a peculiar behaviour in the ratio of terminal wind speed v_∞ over surface escape speed v_esc. It seemed to decrease sharply (by a factor of about two) for stars with an effective temperature (T_eff) cooler than ∼21 kK. Although these early results already hinted at the possibility of a corresponding increase in Ṁ, such an increase was quantified by Vink et al. (1999, 2001). These authors (possibly) identified a recombination of iron (from Fe IV to Fe III) near the sonic point. The increase in the amount of line transitions for Fe III resulted in more efficient line driving as the underlying physical reason for the mass-loss jump. Based on these models, Vink et al. (2001) provided fit-formulae and corresponding ‘mass-loss recipes’ that are divided into a ‘hot’ and ‘cool’ prescription, where the exact location of the ‘bistability’ divide depends on the metallicity Z and the classical Eddington parameter (Γ_e, see Equation (8))¹. Depending on these parameters, the cool solution in this recipe may have up to an order of magnitude more mass loss than the hot solution, resulting in a big mass-loss jump, typically located somewhere in the range T_eff ∼ 22.5 − 27.5 kK. Since it was first introduced, this prescription has become a standard choice to describe mass loss from line-driven winds in stellar evolution codes, such as in Modules for Experiments in Stellar Astrophysics (MESA) (Paxton et al. 2010, 2013, 2015) and the Geneva stellar evolution code GENEC (e.g., Groh et al. 2019). Studies of the effects of such a bistability mass-loss jump in evolution calculations suggest that if this occurs during the final parts of the main sequence (or during the post-main sequence evolution of stars that do not reach the red supergiant stage) substantial mass and angular momentum loss may occur. This then more easily allows for the formation of Wolf-Rayet stars (Björklund et al. 2023) and also causes bistability braking, which is strong reduction of surface rotation speed at T_eff cooler than ∼21 kK (Vink et al. 2010; Keszthelyi et al. 2017; Britavskiy et al. 2024).

Considerable effort has been invested to understand the physical nature and effects of terminal velocity and mass-loss bi-stability. The Vink et al. (2000) models parametrise the wind velocity law v(r), apply the Sobolev (1960) approximation for line transfer, constrain Ṁ from a global energy balance argument, and employ a modified nebular approximation to compute the state of the gas. In follow-up work, Müller & Vink (2008) and Muijres et al. (2012) used their computed Monte-Carlo line force to fit a parametrised radial function. By means of this function then, the equation of motion can be solved in an iterative manner until Ṁ and v_∞ are converged. Vink (2018) applied this technique to luminous stars of different T_eff and found sharply decreasing v_∞ accompanied by mass-loss rates that increased by more than an order of magnitude when decreasing T_eff from 25 to 18 kK (while keeping metallicity, luminosity, and mass fixed, see their Fig. 1). More recently, steady-state line-driven wind models that are based on co-moving frame (CMF) radiative transfer have been calculated (Krtika & Kubát, 2017; Sander et al. 2017; Sundqvist et al. 2019; Björklund et al. 2021, 2023; Krtika et al., 2024). These simulations solve the equation of motion without parametrisation, and keep momentum conserved locally throughout the atmosphere and wind. Specifically, because of the difficulty of driving wind material through the critical sonic point the models computed by Björklund et al. (2023) for the B-star regime do not show any upward increase in mass-loss rate with decreasing T_eff (again while keeping the metallicity and stellar luminosity to mass ratio fixed). The models by Krtika et al. (2024) do find a localised and small mass-loss increase (followed by a sharp decrease) but at significantly lower T_eff than 21 kK (see also Petrov et al. 2016). These calculations are also based on CMF radiative transfer and solve the full steady-state equation of motion. However, Krtika et al. (2024) scale their radiative accelerations to corresponding calculations based on the Sobolev approximation causing their critical point to shift from the sonic point to the point in the supersonic wind where the radiative-acoustic wave (known as an ‘Abbott wave’, Abbott 1980) speed equals the flow velocity. Additionally, Krtika et al. (2024) compute bound-bound rates in the statistical equilibrium equations using this Sobolev approximation, neglecting velocity curvature effects in near-sonic regions (Owocki & Puls 1999). These choices may be the potential causes for the differences in Ṁ behaviour seen between the (at first glance very similar) models by Krtika et al. (2024) and Björklund et al. (2023). Mass-loss rates from the three approaches outlined above are shown in Fig. 1, illustrating the significant differences in predicted rates for the region under study in this work. We point out that what is shown are fit-formulae to the underlying model-sets that not always capture the full Ṁ behaviour, introducing additional differences. For instance, in the case of the Vink et al. (2000) rates these formulae often enhance the size of the jump (which in the underlying model-set typically is about a factor of 5, Vink et al. 1999) through enforcing a discontinuous mass-loss jump.

Fig. 1. Comparison of mass-loss rates from prescriptions by Vink et al. (2001), Björklund et al. (2023), Krtika et al. (2024) over the temperature range from 60 to 12.5 kK. Key stellar parameters in the left panel are: log10(L⋆/L⊙) = 5.3, M⋆ = 25 M⊙ and, Z = 0.5 Z⊙. This results in Γe = 0.21. In the right panel adopted properties are log10(L⋆/L⊙) = 5.6, M⋆ = 25 M⊙ and, Z = 0.5 Z⊙, yielding Γe = 0.42. The full lines signify that the result of the fit-formulate given by these authors are shown rather than the actual underlying model-sets. For the Vink et al. (2000) rates, for which a mass-loss jump is identified to cover the regime ∼22.5 − 27.5 kK, enforcing a strict discontinuity (i.e. at a single temperature) implies that the size of the Ṁ-jump is somewhat enhanced as compared to the underlying models.

So far, observational evidence of a bistability mass-loss jump has proven difficult to obtain. Trundle et al. (2004) comment that in their sample of small Magellanic cloud (SMC) B-(super)giants, for which they obtained optical spectra, the mass-loss rates of stars at the cool side of the jump are noticeably smaller than the Vink et al. (2001) prescription. Additionally, Markova & Puls (2008) did not find a jump from an optical H α analysis of B-supergiants (BSGs) and neither have Rubio-Díez et al. (2022), when studying continuum infrared and radio data of bright OB stars. Most recently a large scale study on 116 B-supergiants in the Galaxy, based on optical spectra only, found no evidence for a bistability jump (de Burgos et al. 2024). Benaglia et al. (2007), however, claimed to have tentative signs of a mass-loss increase for stars cooler than T_eff ∼ 21 kK in their radio analysis. A common denominator of these studies is that they rely solely on diagnostic features for which the associated opacities depend on the square of the density. This makes these mass-loss constraints prone to degeneracies introduced by wind inhomogeneities (or ‘clumping’, see, e.g., Puls et al. 2008). By combining spectroscopy in the optical and UV (see also Bernini-Peron et al. 2023), we aim to break these degeneracies and derive absolute constraints on mass-loss rates as well as wind clumping properties.

With the completion of the UV data taken in the Hubble UV Legacy Library of Young Stars as Essential Standards (ULLYSES) (Roman-Duval et al. 2020) and XShooting ULLYSES (XShootU) (Vink et al. 2023; Sana et al. 2024) programs, we now a have a unique UV and optical dataset of late O-type and early to mid B type stars in the Large Magellanic Cloud (LMC). This allows us, for the first time, to constrain stellar, mass-loss, and wind-clumping characteristics across the anticipated bistability region for a statistically relevant sample. Hence, it also permits the most direct confrontation to date of empirical results with the predictions of mass loss and clumping (Driessen et al. 2019) in the 14 − 33 kK blue supergiant regime. To fit these new data, we utilised a combination of FASTWIND (v10.5, Puls et al. 2005; Sundqvist & Puls 2018) and Kiwi-GA (Brands et al. 2022). Previous empirical studies utilising this method for spectroscopic analysis include Hawcroft et al. (2021), Brands et al. (2022), Hawcroft et al. (2024a), Backs et al. (2024) Brands et al. (in prep.). From here on, we refer to Hawcroft et al. (2021) as H21 and Brands et al. (2022) as B22.

In this work, we first present the dataset utilised to conduct our investigation in Sect. 2.1 and introduce FASTWIND (Sect. 2.2) and Kiwi-GA (Sect. 2.3). In Sect. 3, we describe the stellar, wind, and clumping parameters obtained from fitting the data, with special attention to mass-loss rates and clumping properties. In Sect. 4, we first discuss the influence of some of the choices made during the fitting process, after which we examine the results. Our main findings are that mass-loss does not seem to ‘jump’ across the bi-stability jump, nor do the clumping properties (which have sizeable uncertainties) show a discontinuous behaviour. We do find a similar downward trend in terminal wind speeds with effective temperature as has been reported previously, though the scatter in especially the ratio v_∞/v_esc is large for our sample. We also find that, on average, about 40% of the total wind mass seems not to reside in dense clumps, bur rather in the more diluted medium around them. Implications of these results to the wider field are discussed in Section 4.5.

2. Observations and methodology

In this section we describe our chosen sample of LMC O- and B-supergiants. We introduce the stellar atmosphere and spectral synthesis code FASTWIND used to model the stars in our sample, and the genetic algorithm (GA) code, Kiwi-GA, used to determine the optimum parameter values and their uncertainties.

2.1. Dataset

We selected 15 BSG stars and two late O I-III stars in the LMC from the Hubble Space Telescope legacy program ULLYSES (Roman-Duval et al. 2020) in which UV spectra for 250 sources in the Magellanic Clouds are taken with the Space Telescope Imaging Spectrograph (STIS) and Cosmic Origins Spectrograph (COS) spectrographs. This program is complemented by the Very Large Telescope (VLT)/X-shooter large program XshootU (Vink et al. 2023; Sana et al. 2024) that collects optical spectra of the same sources. Our dataset covers the spectral range B7 through O8.5, which roughly corresponds to temperatures from 14 kK to 33 kK, which is the range in which the bistability jump in mass loss is predicted and the terminal wind velocity divided by the escape speed has been seen to change by about a factor two (Lamers et al. 1995). The reduced STIS spectra used are from ULLYSES data release DR6; the X-shooter spectra are from XShootU DR2. The photometric K-band data used for absolute calibration of the spectra is from Vink et al. (2023). Star identifications, spectral type and K-band photometry are given in Table 1. The Hubble space telescope (HST)/COS has multiple gratings, the G130M/1096 grating has a spectral resolving power of 6000 ranging from 940 − 1240 Å. The G130M/1291 and the G160M grating have a resolving power ranging from 11 000 to 19 000 and cover the wavelength range from 1141 − 1783 Å. The HST/STIS spectra obtained using the E140M and the E230M grating cover a range of 1141 − 1708 Å and 1608 − 2366 Å respectively. The resolving power of the E140M grating is 48 500, and the E230M grating has a resolving power of 30 000. If multiple datasets would cover a modelled feature, the used data was chosen based on signal-to-noise ratio and resolving power. The optical data have two sections, the VIS and the UVB arm from VLT/X-shooter which have been stitched together by Sana et al. (2024). The UVB arm covers a spectral range from 3100 to 5500 Å with a spectral resolving power of 6700 for the chosen slit width of 0.8″. The VIS arm covers the longer wavelengths from 5500 Å to 8000 Å and has a resolving power of 11 400 for the slit width of 0.7″. The visual spectra have a high S/N that is usually higher than 100. The UV data have a lower S/N at around 20.

2.2. FASTWIND

The unified model atmosphere code FASTWIND solves the NLTE rate equations assuming statistical equilibrium in a spherically symmetric and stationary extended stellar envelope comprising both the photosphere and the outflowing wind. We employ here version 10.5 (see Puls et al. 2005; Rivero González et al. 2012; Sundqvist & Puls 2018; Carneiro et al. 2018), which accounts for the accumulative feedback-effects from the multitudes of metallic spectral lines upon the radiation field and atmospheric structure by means of a computationally efficient statistical method. Some of those chemical elements, that are subsequently used for detailed spectroscopy (here H, He, Si, C, N, O), are treated separately using more precise radiative transfer calculations (for details, Puls et al. 2005).

The atmospheric structure is computed assuming (quasi-)hydrostatic equilibrium in the deeper atmospheric layers, connecting smoothly to an analytic wind outflow with radial velocity field v(r) = v_∞(1 − bR_*/r)^β and average density structure $⟨\rho (r)⟩=\dot{M}/(4\pi v(r)r^2)$. The transition between the (quasi-)hydrostatic atmosphere and the wind is set to occur at a velocity 0.1 times the gas sound speed at effective temperature. Here the input parameters describing the wind are the mass-loss rate Ṁ, wind acceleration parameter β, and terminal wind speed v_∞. bR_* is the radius at the photospheric boundary. These parameters supplement the effective temperature T_eff, surface gravity g, stellar radius R_*, and atmospheric chemical abundances in the list of input stellar and wind parameters describing the model atmosphere. Regarding elemental abundances we here adopt the solar composition by Asplund et al. (2009), and scale the overall metallic abundances (Z) accordingly to 0.5 times solar composition to match the LMC. The helium to hydrogen number abundance ratio, Y_He, as well as some specific metal abundances (carbon, nitrogen, oxygen), are free-parameters in the fitting procedure. The silicon abundance is set to 7.2 dex (defined as log₁₀ϵ_Si = log(N_Si/N_H + 12) as in Brott et al. (2011) (differing slightly from the value given in Vink et al. 2023) and is not used as a fit parameter. This constant Si-abundance has been chosen as we expect the Si-abundance to be nearly constant over the LMC. Additionally, the GA methodology used here is not optimised to constrain abundances as these will only have a small effect on the overall goodness of fit. We re-fitted several of the targets spread over the temperature range to determine if including the Si-abundance (and micro-turbulence, see below) as a free parameter changes the results of the other parameters. All fits of the test stars had a large uncertainty margin for the Si-abundance that enveloped our baseline of 7.2 dex. This suggests that our choice of not including the Si-abundance as a fit parameter has not influenced our other results. Additionally, the new best fit Si-abundance was within 0.2 dex of the assumed value.

Inhomogeneities in the wind are described with the effective opacity formalism from Sundqvist & Puls (2018), incorporating modifications in opacity associated with clumping in physical and velocity-space. This formalism assumes a stochastic two-component wind consisting of over-dense clumps and a rarefied interclump medium. Mean densities are related to the root-mean-square (rms) average via a clumping factor f_cl ≡ ⟨ρ²⟩/⟨ρ⟩² ≥ 1 and the interclump densities ρ_ic set by f_ic ≡ ρ_ic/⟨ρ⟩≥0. Light-leakage effects in spectral lines (‘velocity-porosity’) are accounted for by a velocity filling factor 0 ≤ f_vel ≤ 1, defined as the fraction of the velocity field covered by the dense clumps. Wind clumping is assumed to start at some wind velocity v_cl, start > a, where a is the gas sound speed, and then increase linearly with velocity until the input parameter-values for f_cl, f_ic, and f_vel are reached at v_cl, max. Clump optical depths for spectral lines are calculated from the input parameters using the Sobolev (1960) approximation; that is, we do not assume that clumps are optically thin or thick but instead compute them for all lines² according to the structure parameters in order to evaluate corresponding effects upon the ionisation balance and spectrum formation. The ionisation balance is then calculated for an average effective medium taking into account the clumped and interclumped components. Moreover, the micro-turbulent velocity is assumed to increase linearly with wind velocity from a fixed photospheric value v_mic = 10 km/s to a maximum value set by an input parameter scaled to the terminal wind speed. To quantify the effects of fixing the micro-turbulence to 10 km/s, we refitted a select number of stars over the investigated temperature range while instead fixing v_mic to 5 km/s, 20 km/s, and finally allowing it to be a free parameter. Results of these runs are summarised in Fig. A.2. In short, we found that the best fit effective temperature and mass-loss rate changed only within the error margins. The CNO abundance overall also did not show any strong changes; typically the oxygen and nitrogen abundance stayed the same within their sizeable uncertainties. The carbon abundance also generally showed only variations within the error margin. However, for one source the abundance increased by 0.5 dex when lowering v_mic to 5 km/s, and for another source, when allowing the micro-turbulence to be free, the carbon abundance decreased by 1 dex. When letting the micro-turbulence free, on average we obtained slightly higher values than the assumed 10 km/s. However, as mentioned above, this does not significantly impact the overall results we focus on here.

The full method we employ to account for effects of wind clumping is presented and discussed in detail in Sundqvist & Puls (2018), and recent quantitative spectroscopic applications involve H21 and B22. Finally, we note that in the limit of a void interclump medium (f_ic → 0), neglecting velocity-porosity, and assuming all clumps are optically thin and follow the mean velocity field, our method recovers the alternative clumping descriptions included as defaults in alternative atmospheric codes such as Potsdam Wolf-Rayet Stellar Atmospheres (PO WR) (Gräfener et al. 2002; Hamann & Gräfener 2003), and CMFGEN (Hillier & Miller 1998) (see also discussion in Sect. 4.5). We note, however, that even in this limit, the onset and radial stratification of the clumping parameters we assume in this paper may be different from what is typically assumed in these alternative codes.

We also include X-rays produced by wind embedded shocks, as this turns out to be necessary to reproduce some features in the UV resonance lines, particularly in the surprisingly strong C IV profiles of cooler B stars. The details of the X-ray implementation can be found in Carneiro et al. (2016). In this description the energy emitted by the hot gas is given by:

$$\begin{array}{c}\hfill {ϵ}_{\nu }=f_{\mathrm{X}}(r)n_{\mathrm{p}}(r){\mathrm{\Lambda }}_{\nu }(n_{\mathrm{e}}(r),T_{\mathrm{s}}(r)).\end{array}$$(1)

Here n_p and n_e are the proton and electron density of the stationary pre-shock wind, T_s the shock temperature and f_X is defined as f_X = 16e_s², where e_s is the X-ray volume filling factor and the factor 16 is included to account for the density jump in a strong adiabatic shock (making the allowed range for f_X zero to 16). The term Λ_ν is the frequency-dependent volume emission coefficient per proton and per electron. The shock temperature is approximated in the strong shock limit by

$$\begin{array}{c}\hfill T_{\mathrm{s}}(r)=\frac{3}{16}\frac{\mu m_{\mathrm{H}}}{k_{\mathrm{b}}}{u}^2\end{array}$$(2)

with u being the shock velocity and μ the mean atomic weight. The shock velocity itself is set by

$$\begin{array}{c}\hfill u(r)=u_{\infty }{\left[\frac{v(r)}{v_{\infty }}\right]}^{{\gamma }_{\mathrm{s}}}\end{array}$$(3)

where both the maximum jump speed u_∞ and the hardness-parameter γ_s are input parameters. Finally an onset radius of the X-ray emission is chosen by the minimum of $R_{\min }^{\mathrm{input}}$ and the radius at which v_min = m_Xa is reached. Here, m_X is a parameter that one can tune freely if one would want to change the onset of X-rays. In this paper we keep $R_{\min }^{\mathrm{input}}=1.45R_{\star }$, γ_s = 0.75, m_X = 30.0 constant. $R_{\min }^{\mathrm{input}}$ has been taken from Pauldrach et al. (1994), γ_s is chosen to be in between the values of Krtika & Kubát (2009), Carneiro et al. (2016) and Pauldrach et al. (1994), finally m_x is taken from the best-fit value of Pauldrach et al. (2001). We do fit the filling factor (f_X) and the maximum jump velocity (u_∞); effectively this means we are only fitting the overall strength of the X-rays while keeping the onset and hardness constant for each star.

2.3. Kiwi-GA

We are trying to constrain many parameters simultaneously in complex systems. Quite generally this means that regular fitting routines may have issues getting around local minima to find robust global solutions at a reasonable computation cost. For the problem at hand, making a grid large enough to cover the full parameter space is currently not a practical solution. A genetic algorithm (GA) is a way to try and minimise these issues, as it is capable of testing new models that are not decided by gradient descent; the specific genetic algorithm used in this work is Kiwi-GA³ (Brands et al. 2022).

The base version of Kiwi-GA starts out by computing a first generation of FASTWIND models within the chosen parameter space for all parameters one wants to fit. For all models in the initial generation a goodness of fit is calculated as

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=0}^N{\left(\frac{{\mathcal{F}}_{\mathrm{obs},\mathrm{i}}-{\mathcal{F}}_{\mathrm{mod},\mathrm{i}}}{{\mathcal{E}}_{\mathrm{obs},\mathrm{i}}}\right)}^2.}\end{array}$$(4)

Here N is the number of points in the spectrum, ℱ_obs, i is the observed normalised flux, ℱ_mod, i is the normalised model flux and ℰ_obs, i is the uncertainty in the observed flux. The points in the spectra are modified object to object to minimise the continuum as this might lower χ² artificially. However, one should keep in mind that the absence of specific lines also helps characterise stars. The best fit will remain the same even when including a lot of continuum, but the error range might change slightly as it scales with χ² (see Equations (5), (6)).

The next generation will be formed by combining parameters of the fittest previous generations. For example, the effective temperature of one of the best fitting models may be combined with the surface gravity of another well fitting model to create a model in the new generation. In an attempt to get a good coverage of the parameter space, two possible mutations (change to the individual parameters) can occur on the parameters. There is a 70% chance that a small mutation occurs, and the magnitude of this mutation is determined by a Gaussian centred around the current value with a small width. The goal of these small mutations is to give an in-depth exploration of the region around the best fit value. There is a 30% chance that a large mutation occurs; these mutations are also normally distributed but the width of the Gaussian for this mutation is 30% of the full parameter space. These two mutation changes are not additive, 21% of the time the parameter is kept constant. This method has been compared to the traditional by-eye fitting and a grid-search (Sander et al. 2024) using different spectral synthesis codes (PO WR, Gräfener et al. 2002; Hamann & Gräfener 2003; Oskinova et al. 2007; CMFGEN, Hillier & Miller 1998). They studied the differences in resulting fit parameters when using the different fitting routines on three O stars. Overall the different methods gave comparable results, but differences in T_eff up to 3000 K were noted and mass-loss also saw variation of around 0.3 dex. We note that since, in addition to the different fitting techniques, also clumping descriptions vary between these codes, consequently that part of these differences may arise because of that (see also discussion in Sander et al. 2024).

2.4. Two-step process

Described above is the base version of Kiwi-GA. However, in an attempt to make the fitting of the 18 parameters more consistent, we split the above routine into 2 steps. The motivation being that, even though a GA attempts to not get stuck in local minima, in the normal approach we could clearly see that T_eff was not consistently correctly estimated. This was apparent from He II lines not being present in the data but being quite strong in the best fit. To solve this issue, we aimed to start the GA in a good estimate of the T_eff. To obtain this estimate, we decided to first perform a simple run on a selection of optical lines only (all optical H, He, and Si lines). This first optical-only run only includes 6 fit parameters (T_eff, log g, Y_He, max v_rotsin i, β, Ṁ) and a smooth wind outflow is assumed, which makes this part of the fitting routine very efficient while giving consistent results after only 20 generations with 107 models for each generation.

The results of the first run are then used as a simulated first generation for the full run using all 18 parameters (T_eff, g_eff, Y_He, CNO-abundance, upper limit on v_rotsin i, Ṁ, v_∞, β, f_cl, f_ic, f_vel, v_cl, start, v_cl, max, wind turbulence, f_X, and u_∞) and all lines (Table B.1). Here g_eff is the measured surface gravity that is reduced by the rotation (g_eff = g_⋆ − (v_rotsin i)²/R_⋆). This means the first generation of the full run all have the same 6 parameters of the optical-only run, the others are distributed randomly across the parameter space. We limit the variation of Y_He and the maximum of v_rotsin i to the uncertainty region of the first fit to keep the influence of the UV lines on these parameters limited. The full run with 18 parameters is ran for an additional 30–50 generations resulting in a total of 50–70 generations with 107 models for each generation. The values and uncertainty margins of Y_He and max v_rotsin i are also derived from the optical-only run. This approach ensures that the full run starts with predisposed values at approximately the correct temperature, surface gravity, helium abundance, mass-loss rate, and maximum rotational velocity. Additionally, the first run also gives a good upper limit for Ṁ before adding all complexities of wind clumping and means we know something might be wrong if we find a higher Ṁ in the full fits.

The consequence of using this method is that, if for some reason the UV lines are formed with very different parameters than the optical lines we could be biasing the fits more heavily. However, this should in general not be the case and the opposite is a problem of the one-step routine where the UV lines, as they have a lot more data points in the spectral window, carry a bigger weight in the fit while they have a bigger uncertainty. Another concern could be that we are biasing ourselves to a local minimum, which is harder to escape, and is predisposed to a lower clumping factor as our best fit from the optical starts at a smooth wind. However, after some tests we found consistently lower or equal χ² values for the two-step approach compared to the normal approach after the same total generations. We note that, although both were run for approximately the same number of generation (50–70), the computation time of the two-step approach is considerably lower as the simple optical-only run takes less time as we do not include such elements as X-rays and clumping.

2.5. Uncertainty estimate

We used the new method developed by Brands et al. (in prep.) to estimate the error margins of the fits, which uses the root mean square error of approximation (RMSEA) (Steiger 2016). Starting from the χ² value (Eq. (4)) we compute the RMSEA as:

$$\begin{array}{c}\hfill \mathrm{RMSEA}=\sqrt{max(\frac{{\chi }^2-n_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}}{n_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}(N-1)},0)}.\end{array}$$(5)

Here n_d.o.f. are the degrees of freedom and N the total number of data points. As the lines are cut at slightly different wavelengths depending on the object the n_d.o.f. varies around a value of 5000. When the RMSEA is close to 0, the fits are good. There is no formal cut-off point to decide when a model is no longer acceptable; in this work the 1-σ error cut-off for models is set to:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{RMSEA}}=1.04·min\left(\mathrm{RMSEA}\right).\end{array}$$(6)

The RMSEA values of the best fitting model are typically on the order of 0.1–0.01. The 2-σ error is set to 1.09. These values are calibrated to give similar results to the error margins in Brands et al. (2022). As the error cut-off scales with the best fit, bad fitting stars will have a larger error margin than good fits.

When evaluating the results of this paper we do not place too much weight on the individual best fit value, but instead focus more on the uncertainty region. There are many parameters in our fits that can be disentangled only up to a certain point, which the uncertainty takes into account naturally. Additionally, FASTWIND (like all current atmospheric models), is most likely not able to account for all multidimensional phenomena (see, e.g., Schultz et al. 2023; Debnath et al. 2024) in its parameterised 1D descriptions of the various structure parameters. For instance, wind clumping is handled as a statistical two-component medium which is unlikely to be a perfect description of reality (see also discussion in Sect. 4.5). The genetic algorithm approach is by its nature a statistical approach and although it does try to converge to a ‘best-fitting’ model the real value of this approach is that we can show the range of parameters for which the deviations in the spectra are small enough that the goodness of fit does not change significantly. This gives us global errors that automatically take into account all degeneracies.

2.6. Diagnostic line selection

The initial line selection in the optical and UV was based on the lines used in Hawcroft et al. (2021) and Brands et al. (2022). These lines include all prominent H, He, C, N, O, and Si lines in the observed spectra and can be synthesised using FASTWIND V10.5. As these papers focused primarily on hotter O stars; however, a somewhat modified line list has been used here to focus more on the ionisation stages of C, N, O, and Si that are prominent in the cooler BSGs. For example, in this temperature range there are fewer clear lines from these atoms in the UV domain. For the optical line list we also considered the line list used in Trundle et al. (2004). In the appendix we list all lines used in the Kiwi-GA fitting routine (Table B.1).

2.7. Stellar parameter determinations

Not all parameters we are interested in are direct fit parameters. There are a couple of extra parameters we can derive that are of particular interest, including: the mass, the radius, and Γ_e. The mass can be obtained using

$$\begin{array}{c}\hfill M_{\star }=\frac{g_{\star }{R}_{\star }^2}{G},\end{array}$$(7)

where G is Newton’s gravitational constant and g_⋆ is the surface gravity corrected for centrifugal effects (g_⋆ = g_eff + (v_rotsin i)²/R_⋆). The stellar radius R_⋆ is derived by using the fitted temperature and photometrically obtained stellar luminosity in the K-band (Vink et al. 2023), which we de-reddened. Shifting this apparent magnitude to absolute magnitude is then readily done using the distance modulus to the LMC set to 18.48 (Pietrzyński et al. 2019). The de-reddening will influence this luminosity anchor very little as the K-band extinction is low and thermal emission by dust only becomes relevant at longer wavelengths. When investigating how the K-band reddening influences the mass-loss rate, we computed that a change of 0.1 mag in the K-band luminosity results in only an 8% change in the mass-loss rate. This is thus lower than the expected error due to the depth of the LMC, which is on the order of 16% (Subramanian & Subramaniam 2009), and the typical errors of our fits are on the order of a factor 2. With these values we only need to find the radius for which the luminosity in the chosen band matches the observations.

Obtaining good estimates of stellar masses is also important for the wind analysis, since wind strength scales strongly with the classical Eddington parameter:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{e}}=\frac{{\kappa }_{\mathrm{e}}{L}_{\star }}{4\pi cG{M}_{\star }},\end{array}$$(8)

with κ_e the electron scattering opacity and c as the speed of light. The term Γ_e gives the ratio of radiative to gravitational acceleration in the case that the only opacity source is κ_e. When including these parameters on the various plots illustrating our results, we adopted

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{e}}=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{H}}}\left(\frac{1+I_{\mathrm{He}}·Y_{\mathrm{He}}}{1+4·Y_{\mathrm{He}}}\right).\end{array}$$(9)

Here σ_T is the Thomson cross section, m_H is the mass of hydrogen, I_He is the ionisation of He giving either one or two free electrons for each He-atom (set to one for all B stars in our sample and two for the O stars); hydrogen is assumed to be fully ionised, and Y_He is the helium to hydrogen number abundance ratio. Note that Γ_e does not take into account the line-driving effect of the radiation acceleration, but instead gives a generic radius independent value as a measure of the closeness of the star to the classical Eddington limit.

3. Results

In this section we present the stellar and wind parameters determined by the fitting method described above, applied to the sample of LMC stars. All the derived parameters and corresponding error estimates are given in Tables 2, 3, and 4.

Figure 2 illustrates fits for a selection of important lines for four typical stars in our sample. The stars are organised from hot to cold, as clearly visible in the He II lines that are present for the hottest star, but disappear when moving to the cooler objects. For the cooler stars the Si II lines instead become stronger, while the Si IV line that is present in the red wing of the Hδ line vanishes. The H α line for all stars shows clear signatures of wind emission filling in the photospheric absorption component, indicating this line is a good mass-loss indicator throughout our sample (though see discussion below on some complicating factors). All four stars show P-Cygni like features in the UV Si IV and C IV lines, again indicating the presence of significant wind outflows. In this respect, we note that C IV is clearly present in the wind of objects at around 17kK (although the doublet line there is narrower and weaker than in the other two stars); we discuss this further in Sect. 3.3, as it has consequences for our modelling of the wind ionisation balance and derivation of terminal wind speeds in this regime.

Fig. 2. Comparison of diagnostic lines for four stars. The stars (Sk−67 ° 107, Sk−68 ° 41,Sk−69 ° 52,RMC-109) are organised from high (left) to low (right) effective temperature (in spectral type: O8.5 II, B0.7 Ia, B2 Ia, B5 Ia). The lines are selected due to their roles as good diagnostics for stellar or wind parameters. From top to bottom the lines are sorted by wavelength. We note that the line names here correspond to the names in Table B.1. The green line shows the best fit and the green shaded region shows the 1-σ uncertainty interval.

3.1. Stellar parameters

Although our analysis and discussion will be focused mostly on wind parameters, we start by presenting the stellar parameters in Table 1. We fit seven stellar parameters in our models, namely the effective temperature (T_eff), the effective gravity at the stellar surface (log g), a maximum (see further below) projected equatorial rotational velocity (v_rotsin i), the surface helium abundance (Y_He) and the surface carbon, oxygen, and nitrogen (CNO) abundances. Figure 3 shows stellar parameter distributions of the inverse of RMSEA for models of a characteristic star in our sample. The best fit values are marked with a red line and the 1- and 2-σ confidence intervals are the dark and light grey shaded regions. The colour scheme marks which generation in the GA the model stems from. We observed that the distributions for effective temperature, surface gravity, and maximum v_rotsin i all have well-defined and fairly sharp peaks in their distributions, suggesting good statistical constraints on the parameters. The abundances are generally more difficult to constrain, but also these sometimes display relatively peaked distributions, as exemplified for carbon in the figure.

Fig. 3. Distribution of the quality of the fits for stellar parameters of star Sk−68 ° 41. In dark grey the 1-σ interval is shown while the light grey shows the 2-σ interval. The red line shows the best fit value. We only display the C-abundance as the other elements have similar fitting curves. As mentioned in Section 2.3 and the following sections we derive the maximum projected rotational velocity from the optical-only fit which is why there are only 21 generations.

The sensitivity to effective temperature of the GA is best constrained from the silicon lines, as well as helium lines for stars on the hotter end of our sample. We note that other lines are sensitive to T_eff and that due to the nature of a GA all of these are taken into account when deriving the best fit and the uncertainty margins. For all stars in the sample we see a strongly peaked distribution for the effective temperature meaning that the combination of stellar lines is adequate for determining T_eff. The derived effective temperatures for the sample range from about 15 kK–30 kK with luminosities ranging from 4.7 to 5.8 log₁₀(L/L_⊙). Figure 4 shows the stars in the Hertzsprung–Russell (HR) diagram with their spectroscopic mass indicated by their colour. We plotted the evolutionary tracks from the MIST data base (Dotter 2016; Choi et al. 2016; Paxton et al. 2010, 2013, 2015) for a 15, 26, and 40 M_⊙ star in the same window. Comparing these models to our spectroscopic results, the majority of our sample appears to be in the post-main sequence phase of evolution. We note, however, that the evolutionary tracks shown in the figure assume a convective core overshoot parameter f_ov, core = 0.2 for a step overshoot model (Choi et al. 2016); if stronger overshooting would be assumed this might extend the main sequence so that some of the stars in the sample could still be on it (Castro et al. 2014).

Fig. 4. Hertzsprung–Russell diagram showing all stars in the sample and with the colour showing the derived spectral masses. The diamonds mark all the stars with spectral types Ia, while the ‘v’ symbols mark the lower luminosity classes. The stellar evolution tracks are from the MIST data base without any rotation (Dotter 2016; Choi et al. 2016; Paxton et al. 2010, 2013, 2015).

The derived surface gravities are presented in Table 2. For which, in particular Hγ, Hδ are sensitive. The fit values for log₁₀(g cm/s²) range from 2.2 to 3.6. Only four stars in our sample have a surface gravity 3.0 < log₁₀(g cm/s²).

Surface gravity is an important stellar parameter because of its direct influence on the derived spectroscopic mass of the star as described in Section 2.7. The resulting derived spectral masses from these surface gravities lie between between 15 and 75 M_⊙, with corresponding stellar radii between 19 and 70 R_⊙. The latter are shown by the colour scheme in the HR-diagram of Figure 4. Γ_e in this sample ranges from 0.10 to 0.46. This parameter is important for the analysis below as it represents another key axis over which wind behaviour changes, in addition to the temperature axis.

The maximum projected stellar rotation (max v_rotsin i) is the only macroscopic broadening mechanism included in our line fits, meaning we do not attempt here to differentiate between different macroscopic broadening mechanisms (e.g. rotation and ‘macro-turbulence’). As such, the quoted v_rotsin i values are upper limits, and not the actual projected equatorial rotation velocity (as its measurement is influenced by the broadening mechanisms we do not handle explicitly). Figure 5 shows the derived upper limits of v_rotsin i as function of the effective temperature. We note that, with the exception of three of the hottest objects, our sample has max (v_rotsin i) < 100 km/s. Indeed, it is for these relatively low values that we expect it to be difficult to disentangle rotation broadening from other macroscopic broadening mechanisms. Specifically, Sundqvist et al. (2013) showed that when applying standard techniques to disentangle effects from macro-turbulence v_mac and v_rotsin i for massive stars with known very long rotation periods, corresponding to v_rotsin i ≲ 1 km/s, one still ends up deriving v_rotsin i ∼ v_mac whenever significant macro-turbulence is present in the observational data. This demonstrates the difficulty of present methods to disentangle these two effects, motivating our methodological choice here (indeed, visual inspection of our observed strategic line-profiles often show more Gaussian-like shapes than the rounded shape predicted by broadening due to rotation). That is, our quoted values for max v_rotsin i may as well reflect the presence of large Gaussian-like turbulent motions and not stellar rotation.

Fig. 5. Derived maximum projected rotational velocity as a function of effective temperature. We note that this maximum value in practice reflects the combined effect of all sources of macroscopic line broadening (see text).

The average helium abundance is Y_He = 0.1 (or alternatively ϵ_He = 11.0, see below) with the highest and lowest values of the helium abundance 1 dex apart (ϵ_He ≈ 11.4 − 10.4) with relative small error margins on the order of 0.2 dex. As one may have noted Table 2 shows values that reach helium abundances below Y_He = 0.08. This might seem problematic as this is below the expected baseline for the LMC at $\mathit{y}=0.{091}_{-0.012}^{+0.014}$ (Russell & Dopita 1990). However, the fitting range in all cases reaches values above this 0.08 baseline and as mentioned above, this range is the important property to assess. We refitted all problematic stars as a test, while not allowing the helium abundance to be lower than Y_He = 0.08 and no values changed appreciatively from Table 2.

Finally, to make the GA sensitive to CNO-abundances, we added atmospheric CNO lines of various ionisation stages. However, because this paper makes an effort to focus on mass-loss rates and not on abundances there are comparatively few data points which are sensitive to the abundance. As such, small changes to a specific abundance only modify the quality of the fit of these specific lines slightly (for instance we include only one NII quintuplet) and typically do not affect significantly the overall fit quality used to define our ‘best-fit’ model. Consequently, we see high uncertainty regions in CNO abundance as compared to the other stellar parameters. The abundances of CNO in this paper are given on the standard scale ϵ_x = 12 + log₁₀(n_x/n_H) where n_H is the number abundance of hydrogen abundance and n_x is the number abundance of the element. The average CNO abundances we find for our sample are ϵ_C = 8.0, ϵ_N = 7.2, ϵ_O = 8.3. Although these are fitted using different lines they do display similar behaviour for their error margins, which are typically on the order of 0.5–1 dex for all elements. Abundances and 1-σ error margins for each star in our sample are given in Table 2. Vink et al. (2023) gathered abundances of common elements in the LMC and averaged the abundances found in different stellar populations to find ϵ_C = 8.01, ϵ_N = 7.03, ϵ_O = 8.4; in comparison to these results, we appear to have relative agreement. However, the uncertainty of our derived abundances is usually on the order of 0.5–1 dex. 5 of the stars in our sample have been studied in the optical regime by Menon et al. (2024) focusing on the CNO abundances to determine if they are possible binary merger products. All of the stellar parameters (e.g., T_eff and log g) have overlap in the error-margins, however for some of the CNO abundances we find drastically different best-fit solutions. Due to the very big uncertainty intervals on our CNO-abundances we do still overlap in most cases, but as mentioned above getting close constraints on CNO-abundances has not been our focus. For a small expansion on our derived CNO abundance see Figure A.1.

3.2. Mass-loss rates

One of the preferred lines for determining mass-loss rates is the H α line, as it is a sensitive mass-loss diagnostic. However, it is important to realise here that since H α is a recombination line in this regime, the line is also sensitive to the wind-clumping factor f_cl. Additionally, it has been suggested by Petrov et al. (2014) that the H α line behaviour and morphology may change with temperature around the supposed bistability jump. If f_cl would be kept constant throughout the complete line-forming region, and clumps would be optically thin, the scaling invariant would simply be $\sqrt{f_{\mathrm{cl}}}\dot{M}$. However, since H α often is formed in the region where wind clumping starts to become effective in our models (i.e. in regions around v_cl, start and between v_cl, start and v_cl, max), deviations from this basic scaling may occur. Even so, the degeneracy between f_cl and Ṁ causes problems when trying to determine absolute mass-loss rates using only the strength of H α. Figure 6 illustrates H α’s strong dependency on mass loss and the corresponding degeneracy with f_cl. Both stars in Figure 6 have the same input parameters (specifically T_eff = 20 kK, log₁₀g = 2.7, v_∞ = 750 km/s, f_ic = 0.1, f_vel = 0.5), except for the mass-loss rate and clumping factor. The orange line is for a smooth outflow and Ṁ = 5 · 10⁻⁷ M_⊙/yr, and the blue line has f_cl = 40 and Ṁ = 8 · 10⁻⁸ M_⊙/ yr. The mass-loss rate is thus different for both stars, but the product $\sqrt{f_{\mathrm{cl}}}\dot{M}$ is the same. One can see that the height of the peak of the Hα line is similar in the two models. The match is not perfect, however, and the overall shape is different because the H α line is partly formed in regions where f_cl has not reached its maximum value (here we assumed v_cl, start = 0.05 and v_cl, max = 0.1 for both model stars). In the optically thin limit and assuming the clumping is constant over the entire velocity range the 2 H α profiles would match perfectly, but the models we use here relax these assumptions. In an attempt to break this mass loss-clumping degeneracy, we fit not only the sensitive H α line but also wind resonance lines in the UV. These often depend differently on wind clumping as their extinction coefficients typically scale as ∼⟨ρ⟩, rather than ∼⟨ρ²⟩ as in the case of recombination lines. We note that additional dependencies from the ionisation balance and velocity-porosity effects (as well as saturation and interclump density effects) can modify these natural scalings of recombination and resonance lines, making it complicated to sort out all inter-dependencies. We show this in the lower panel of Figure 6, displaying the Si IV resonance doublet for the same stars and parameter combinations as the H α line discussed above; it is clear that the dependence on f_cl is very different for these lines than for H α. Specifically for this case, we then observed how the absorption troughs in the doublet lines become weaker for an increasing f_cl, which is here because the ionisation balance of silicon shifts when strong clumping is introduced. In our multi-diagnostic automated fits, we utilised the differences of how lines react to derive simultaneous constraints on f_cl and Ṁ (and the other stellar and wind parameters).

Fig. 6. Clumping influence on the spectral lines. The figures show the influence a change in fcl has on the line profile of on one hand the H α recombination line and on the other hand the Si IV 1400 resonance doublet. Both model spectra are using almost the same input, they only differ in mass-loss rate and clumping. The blue line both has Ṁ = 8 · 10−8 M⊙/yr, and is highly clumped (fcl = 40). The orange line is a smooth outflow but has a significantly higher mass-loss rate (Ṁ = 5 · 10−7 M⊙/yr); the product $\sqrt{f_{\mathrm{cl}}}\dot{M}$ is thus the same for both lines. If we assume an optically thin clumping over the complete line forming region, the 2 H α lines would show perfect agreement. However, this is not the case as we use models that relax the optically thin assumption and are only clumped above a fitted onset-velocity.

Overall, our multi-diagnostic GA approach is able to isolate the mass-loss dependence rather well and so derive absolute values of Ṁ. This is demonstrated in Figure 7, showing the distribution of the inverse of RMSEA for models of the same star but with different mass-loss rates. The more peaked this distribution is, the better the GA can distinguish the different parameter values from each other. The GA method takes into account the interactions of all degeneracies automatically when deriving uncertainty intervals. In the cases where the clumping factor stays somewhat degenerate with the mass-loss rate we see a broadness in both the mass loss and derived clumping factor. This is more common in the low temperature stars. Best-fitting mass-loss rates for all stars in our sample are shown in Figure 12 and vary from ∼10⁻⁶ M_⊙/yr for the hotter stars to ∼10⁻⁸ M_⊙/yr for the cooler stars. Figure 8 further displays the correlation between Ṁ and f_cl for the fit of star Sk−68 ° 41, showing 1/RMSEA for each combination of clumping factor and mass-loss rate (note that one such combination may have several models due to variation in other parameters). Due to the sparse sampling there is some stochasticity in the RMSEA values. The figure illustrates that although the uncertainty in f_cl is rather high, no severe degeneracy issues between f_cl and Ṁ seem to be present in our fits as there is a clear optimal region.

Fig. 7. Quality of all models in a GA run ranked based on the inverse of the RMSEA value as a function of the mass-loss rate. The models for these plot are based on a run trying to fit the star Sk−68 ° 41 (SpT B0.7 Ia). Each point here is a separate model where the colour defines the generation in which this model was run. The dark grey region highlights the 1-σ uncertainty interval and the light grey region is the 2-σ uncertainty.

Fig. 8. Correlation between the clumping factor and the mass-loss rate for star sk−68 ° 41. The colour shows the 1/RMSEA value.

3.3. Wind speeds, accelerations, and turbulence

The lines that provide a good diagnostic of terminal wind speed to the GA fits are Si IV 1400 and C IV1550. These lines often appear as distinct ‘P-Cygni’ profiles, allowing for fairly straightforward visual verification of the fitted v_∞ from the broadness of their blue-shifted absorption. The dependence of the strength of these lines on temperature can be seen in Figure 2, where a few characteristic lines are shown for four stars ranging from high (left) to low (right) effective temperature. The observed C IV and Si IV lines for higher T_eff are well-developed, strong or even saturated towards the blue absorption edge. The features become less broad when reducing the temperature, which generally points to lower terminal velocity. In the standard modelling of the cooler stars in our sample, C IV is actually no longer present in the wind. Nonetheless, as seen in Figure 2, some of these cool stars still display significant observed P-Cygni profiles for C IV1550 (albeit weaker than for the hotter objects). To reproduce this behaviour in the modelling, we have included X-rays in the FASTWIND simulations also in this temperature region (previous versions did not consider X-ray ionisation below T_eff ∼ 25 kK), which then can raise the wind ionisation balance and produce enough C IV to explain the observed lines. The influence of these X-rays upon the determination of v_∞ in this parameter range is further discussed below. The terminal wind speeds resulting from our fits of the full sample are displayed in panel 1 of Figure 9 as a function of effective temperature. For the coolest star in our sample (Sk−67 ° 195) all UV lines that are sensitive to the terminal wind speed are no longer present. As a result, we do not get good constraints on the terminal wind speed this is also visible in the fit distribution shown in figure 17 of the appendix found on Zenodo.

Fig. 9. Derived terminal velocities, wind acceleration parameter β, and maximum wind micro turbulence, as a function of effective temperature. The colour for all panels shows the Γe value. The first panel also shows the linear fit to the terminal velocities found in this work in red compared to the fit (in blue) from Hawcroft et al. (2024b).

The wind acceleration is parametrised by the β value in the analytic velocity law (see previous section). As can be seen in panel 2 of Figure 9, we find relatively high values ≳1.5, with a sample average ⟨β⟩ = 2.1. This is significantly higher than values derived for hotter O stars by B22, who found an average β ≲ 1.0 in that regime, but comparable to the optical only BSG study in the Galaxy by Haucke et al. (2018) and the UV+optical SMC BSG study by Bernini-Peron et al. (2024). Overall, the slower wind acceleration we find for BSGs compared to the faster wind acceleration found for early O stars agree well with a range of previous studies (overview in Puls et al. 2008). The error margins on these values are appreciable but still keep the 1-σ error interval clear above 1. A larger scatter is seen for the cooler stars in the sample.

The inferred values of maximum wind micro-turbulence range from 0.05 to 0.4 in units of v_∞, with a sample average 0.2. The uncertainty interval of this fit parameter is generally high with typical errors around 0.1, but with outliers reaching error margins as high as 0.4 on this scale. We note that this value then also has a significant effect of the derived terminal wind speed, as turbulence adds to the blueward extension of the modelled line-profile.

3.4. Clumping

To recapitulate, in addition to the clumping factor (f_cl), we fit the interclump-density parameter f_ic, the velocity filling factor f_vel, the onset of clumping v_cl, start, and the wind speed at which maximum values for the clumping parameters are reached v_cl, max. All parameter-values displayed in figures and quoted in text are the maximum values, applied in our models from v_cl, max to the outer boundary. As also found in the previous studies by H21 and B22 it is generally challenging to obtain simultaneous constraints on these parameters.

Figure 10 shows the determined clumping parameters as function of effective temperature. The plots are generally characterised by significant scatter and large error bars, illustrating the general challenge of obtaining these parameter values or finding a trend with temperature. Within 1-σ errors, clumping factors range between the lower limit 1 up to ∼40 for our sample and velocity filling factors between very low values ≲0.1 all the way up to the upper limit unity. While the scatter in interclump density is high, this parameter does seem to show some preference for values ≳0.1. v_cl, start is generally preferred to have wind clumping start rather close to our lower allowed bound at v_cl, start ≈ a, although the scatter in this parameter is high with typical values within 1-σ ranging from v_cl, start/v_∞ ≲ 0.01 up to ∼0.1. v_cl, max typically is about ∼0.2v_∞, with similarly broad ranges including a few outliers at the high end reaching values of 0.35; the 1-σ interval for this parameter ranges from 0.05 to 0.15 v_∞. Additionally, panel four in Figure 10 uses colour to show the difference between v_cl, max and v_cl, start, where we find that most stars have values v_cl, max − v_cl, start ≈ 0.2 in units of v_∞. A typical distribution of fit quality for the clumping factors can be in seen in Fig. 11.

Fig. 10. Clumping parameters with respect to the effective temperature. The panels from top left to bottom right show the clumping factor, interclump density (in units of mean density), velocity filling factor and the onset clumping velocity (vcl, start) and the velocity of maximum clumping (vcl, max). The three first panels colour according to Γe for the stars, while in the last plot the colour indicates the difference between vcl, start and vcl, max. For each star in this plot the two points are connected by a dashed line (not visible in some objects as the error margins are covering the dashed line).

Fig. 11. Distribution of the quality of the fits for clumping parameters. We show a typical distribution of the models at end of a GA fitting run. The colour indicates the generation of the model and the grey areas indicate the 1σ and 2σ uncertainty intervals.

4. Discussion

4.1. The absence of a generic mass-loss jump in the T_eff range 15–35 kK

A key focus of the present analysis is to investigate whether there are empirical signs of a general increase in Ṁ when moving from hotter to cooler stars in the region T_eff ∼ 15.0 − 27.5 kK. Previous empirical studies attempting to tackle this question have often been hampered by the mass loss-clumping degeneracy of optical recombination lines like H α (e.g., Markova & Puls 2008). Although Markova & Puls (2008) do not find a jump comparable to the predicted jump (Vink et al. 2001) as there is some indication that the clumping factor might reversely decrease with temperature as well (Driessen et al. 2019) this was considered to not be fully conclusive. As outlined in the previous section, here we break this degeneracy by a detailed account of wind clumping and by considering a multitude of diagnostic lines in the optical and UV. This allows us to derive absolute empirical constraints on Ṁ.

Figure 12 shows the empirically derived mass-loss rates of our sample as function of effective temperature. We stress that this figure is showing stars which not only differ in T_eff but also in their L_⋆/M_⋆ ratios, where the latter have a big influence on the mass-loss rate. This ratio is shown efficiently by colour coding the stars according to their Γ_e values. Two stars with the same T_eff may have different Ṁ depending on their individual values of Γ_e, where typically the star with the higher value also has a higher mass-loss rate. Therefore, we indicate Γ_e with a colourbar in Figure 12. The figure demonstrates that there is a clear general downward trend in Ṁ when moving from the hotter to cooler objects in our sample. If we inspect only stars with similar Γ_e values such a trend persists. For objects with higher Eddington ratios (green and yellow in figure), there is a fairly small but noticeable downward trend in the region T_eff ∼ 25 − 15 kK. For objects with somewhat lower Eddington ratios (purple in figure) we also observe that the collection of stars above T_eff ∼ 30 kK have higher mass-loss rates than corresponding objects around and below T_eff ∼ 20 kK. That is, in our current sample there are no empirical signs of an upward jump (or upward trend) in Ṁ with decreasing T_eff in this region.

Fig. 12. Derived mass-loss rates for our sample as function of effective temperature. The top plot shows the derived mass-loss rates and he colour shows the derived classical Eddington parameter for each star. The bottom plot shows the ratio between the theoretical prescriptions by Vink et al. (2001), Björklund et al. (2023), Krtika et al. (2024) and the empirically found values. The green line indicates where these are equal.

We next compare our empirically derived rates (Figure 12) to the predictions by Vink et al. (2001), Krtika et al. (2024) and Björklund et al. (2023). As discussed in the introduction, the Ṁ(T_eff) behaviour for a fixed mass and luminosity differs between these prescriptions. In the context of the present study, the main difference in the predictions of these schemes is that Vink et al. (2001) predict significantly higher mass-loss rates for stars below T_eff ∼ 22 − 25 kK than for hotter stars, whereas Björklund et al. (2023) do not predict such a jump but instead a steady decrease of mass-loss rates with decreasing temperature. The Krtika et al. (2024) rates do not exhibit such a monotonic decrease of mass loss with decreasing temperature, but rather has a characteristic bump at T_eff ∼ 13 − 15 kK. This bump is very localised and predicts a rather modest mass loss increase (followed by a steep decrease), which also lies at the edge of our observational coverage range. As such, it is not possible to scrutinise this particular low-temperature feature in the mass-loss behaviour for the present set of stars (spectral types O8.5-B7).

The comparisons between Vink et al. (2001), Björklund et al. (2023), and Krtika et al. (2024) are shown in the bottom panel of Figure 12, where the ordinate displays predicted Ṁ divided by the empirical Ṁ derived in this paper, again plotted as function of effective temperature. All predicted rates have been computed by means of the ‘mass-loss recipes’ provided in the corresponding papers, using the stellar parameters derived for each star (L_⋆, spectroscopically derived M_⋆, T_eff, v_∞), and a metallicity reflecting that of the LMC. Here we corrected for the higher solar metallicity used in Vink et al. (2001) compared to the Björklund et al. (2023) rates as described in Sundqvist et al. (2019). As such, additional dependencies on for instance Γ_e are naturally accounted for in these comparisons. At temperatures above ∼22.5 kK all three prescriptions are, on average, fairly aligned with the empirical rates. However, at lower temperature the predictions by Vink et al. (2001) jump up as a result of the bistability jump causing their mass-loss recipe to over-predict the empirical rates by more than an order of magnitude. The recipe by Björklund et al. (2023) shows no such jump, and instead displays a rather constant trend with T_eff, aligned with the trend inferred for the empirical rates. For the coolest star in our sample the Björklund et al. (2023) rates are much lower than our best-fit values. This might be an outlier or due to the extrapolation of the Björklund et al. (2023) model grid. For stars hotter than 26 kK, two stars agree very closely with the observation while two underestimate by about a factor 8. When comparing the observations to the Krtika et al. (2024) prescription we see that they also agree fairly well. However, in contrast to the Björklund et al. (2023) prescription the Krtika et al. (2024) prescriptions oscillate between over- and under-predictions.

Although the prescription by Björklund et al. (2023) reproduces the empirically found trend rather well, we observe from the bottom panel of Figure 12 an almost constant underestimation with a relative offset of $⟨|{\dot{M}}_{\mathrm{P}}-{\dot{M}}_{\mathrm{E}}|/{\dot{M}}_{\mathrm{E}}⟩\sim 0.5$ from our best-fit values, where Ṁ_P denotes the prescription results and Ṁ_E is the empirical mass-loss rate. This average is weighted with the error margin. For the prescriptions by Vink et al. (2001) we split such average comparisons into stars with T_eff above and below the jump. For stars above this threshold we find $⟨|{\dot{M}}_{\mathrm{P}}-{\dot{M}}_{\mathrm{E}}|/{\dot{M}}_{\mathrm{E}}⟩\sim 0.8$ and for stars below $⟨|{\dot{M}}_{\mathrm{P}}-{\dot{M}}_{\mathrm{E}}|/{\dot{M}}_{\mathrm{E}}⟩\sim 20$. The Krtika et al. (2024) prescription fluctuates between over and under estimating but there is no clear cut-off temperature to split the prescription. If we compute the average relative offset for the full sample, we find $⟨|{\dot{M}}_{\mathrm{P}}-{\dot{M}}_{\mathrm{E}}|/{\dot{M}}_{\mathrm{E}}⟩\sim 0.5$. Which is exactly the same as the Björklund et al. (2023) offset and close to the Vink et al. (2001) offset for the hot stars.

Again these comparisons and simple averages show quite clearly that the sudden increase in mass-loss rate predicted at the supposed bistability jump by the Vink et al. (2001) prescription is not present in these observational results.

We may further compare our results to other empirical studies. Bernini-Peron et al. (2023) show a study of 4 BSGs (B2-B5) in the Galaxy with access to UV spectra and therefore also try to constrain the clumping factor. The mass-loss rates acquired in that study are also higher for similar effective temperatures as expected due to the higher metallicity. They do find some localised increase when compared to 4 literature stars at hotter temperature, but this is hard to verify due to not covering the supposed bistability jump with the stars from their sample. In contrast Markova & Puls (2008) assume a smooth wind outflow and derive mass-loss rates for Galactic BSGs using H α. Within our current spectroscopic formalism for wind clumping, these mass-loss values should be interpreted as upper limits. Assuming effects of clumping do not change the inferred rates of mass loss significantly across the region, we may still examine the mass-loss trend with effective temperature. Over their investigated temperature range the derived mass-loss rates decrease monotonically, showing no signs of an increase around 22.5 kK. The derived Ṁ in these Galactic stars are generally higher than those we infer here for LMC stars. As the luminosity ranges are comparable, this is due to the higher metallicity in the Milky Way. More recently, Rubio-Díez et al. (2022) examined winds of OB-stars, including 13 Galactic BSGs, by means of radio and far-infrared continuum observations. They derive values of Ṁ for the BSGs in their sample that are significantly higher than here. This is to be expected as, first, their measured mass-loss rate is an upper limit; second, they targeted very luminous objects, and, third, they consider a Galactic sample. Compared to this study, their general trend with effective temperature is similar. A comparison between their rates and the predicted Vink et al. (2001) rates therefore too shows a large discrepancy at temperatures below 22.5 kK, with the predictions higher than the empirical rates by between 1 and 2 orders of magnitude (for their estimates of maximum mass loss).

Very recently, de Burgos et al. (2024) used optical spectra of 116 Galactic BSGs in the temperature range 15–35 kK to derive upper limits to mass-loss rates (upper limits because they used smooth wind models); they did not find any signs of an upward trend in mass-loss rate with decreasing temperature in the predicted bistability region.

4.2. Terminal wind speeds

As can be seen in the top panel of Figure 9, the resulting terminal wind speeds suggest a linear dependence on effective temperature. Fitting a linear function v_∞ = aT_eff − b, we find best fit values a = (9.7 ± 0.6)10⁻² km (sK)⁻¹ and b = 1360 ± 130 km s⁻¹. As mentioned in the previous section, the linear behaviour as well as the numerical fit values are in line with the recent findings of Hawcroft et al. (2024b), who fitted a = (8.5 ± 0.5)10⁻² and b = 1150 ± 170 km s⁻¹ for their larger sample. Compared to the prescription from the modelling efforts by Vink & Sander (2021) we notice that the discontinuity at the supposed bistability range does not agree with the observed values. Instead we see a continues increase of v_∞ with T_eff. The classical study of UV P-Cygni lines by Lamers et al. (1995) found a sudden decrease in the ratio v_∞/v_esc around T_eff ∼ 21 kK, a behaviour that has often been associated with the presence of a bistability jump, whereas other studies have found a more gradual decrease (e.g., Markova & Puls 2008). Here v_esc² = 2GM_⋆(1 − Γ_e)/R_⋆ is the effective (i.e. reduced by Thomson scattering) escape speed from the stellar surface. Figure 13 shows our measured v_∞/v_esc against effective temperature. The figure is dominated by large scatter in the inferred values, reflecting also uncertainties in the stellar parameters used to compute v_esc. We do indeed observe some indication of a downward trend in the range T_eff ∼ 25 − 17 kK, although also in our relatively small sample at least some outliers to this trend clearly exists on the cold end. Actually, one may detect a similar behaviour also in the above-mentioned larger study by Hawcroft et al. (2024b), if one zooms in on the region around T_eff ∼ 20 kK for LMC stars (see their Figure 5). However, in that study this quite subtle ‘by-eye’ detected trend is statistically washed away by the large number of stars with hotter effective temperatures and the scarcity of stars with T_eff < 20 kK.

Fig. 13. Terminal wind speed divided by escape speed plotted over temperature. The mass used to determine the escape speed is the effective spectral mass.

Nonetheless, although uncertainties are still large, there may be signs of the general trends observed by Crowther et al. (2006) present in this sample as well. This trend is similar to the trend by Lamers et al. (1995) but instead of a sudden downward jump upwards in v_∞/v_esc a more gradual decrease with decreasing T_eff is seen instead.

4.3. X-rays

Without X-ray ionisation the later-type stars in our sample show very little triple ionised carbon in their winds and outer atmospheres, resulting in absent or very weak C IV1550 lines, without any P-Cygni like wind signatures. Observations, in contract, show that some cooler stars in the sample have strong C IV1550 lines with clear wind signatures. To reproduce this, we added X-rays produced due to shocks in the wind also in this temperature region, as outlined in Sect. 2.2. To ensure all results are comparable, even the stars that fit well without X-rays were fit with X-rays, and the results are noted here, with the exception of the coolest star in the sample for which the X-rays caused convergence issues. This also allowed us to see if there is a noticeable improvement in fitting by including X-rays. The hotter stars that do not need X-rays to excite carbon or the cool stars with no noticeable C IV give X-ray values with much higher error margins or very low values and the fit quality stays comparable.

Figure 14 shows the difference the X-rays make to the ability to fit the C IV and S IV wind resonance lines of sk−70 ° 16. The best-fit v_∞ of the model including X-ray ionisation is 1000 km/s while the one without gives v_∞ = 500 km/s. Notice how this big change in v_∞ and the addition of X-rays as a whole did not influence the Si IV line noticeably while transforming the C IV line substantially, showing that the Si IV line is a poor indicator of v_∞. Although introducing X-rays forces us to consider even more free input-parameters in our modelling, the upside of the addition is that it does not significantly affect any of the other parameters within our sample. For example, T_eff is typically slightly lowered when X-rays are included, but the reduction is always within the 1-σ uncertainty of the originally found effective temperature. The issue is often that while weak, unsaturated UV resonance lines may still show clear signatures of line formation in the lower wind, they can be insensitive to v_∞ due to very low optical depths in the outermost wind. To properly capture this one thus needs to model also the detailed wind ionisation balance, which is typically not done in the simplified large-scale approaches that work well on strong and (almost) saturated UV wind lines (e.g. Hawcroft et al. 2024b; Lamers et al. 1995).

Fig. 14. Three versions of the same doublet line per panel. In black we show the data from Sk−70 ° 16 a B2 II star, in blue the best fit when not allowing for X-rays, and in red is the best fit when allowing for X-rays. The top panel shows the C IV 1550 line, and the bottom panel shows the Si IV1400 line.

We found an average log₁₀(L_X/L_bol) =  − 7.3 for our sample, but with very large error margins of ±1 to ±2.5. For individual stars the large errors cause the best fit range to reach down to log₁₀(L_X/L_bol) <  − 8. Our average value is on the same order as the often-quoted empirical relation log₁₀(L_X/L_bol) =  − 7.2 ± 0.2 found from X-ray observations of O stars (Rauw et al. 2015). Crowther et al. (2022) have studied the X-ray properties of specifically O-stars in the LMC and found a similar relation log₁₀(L_X/L_bol) =  − 6.90 ± 0.65. It also agrees well with the X-ray luminosity range found by Bernini-Peron et al. (2023) (log₁₀(L_X/L_bol)≈[−7, −8.3]). When computing the range of shock temperatures using Equation (2), we note typical shock temperatures ranging from 10⁵ K, at the X-ray onset radius, to 5 × 10⁶ − 10⁷ K at the terminal wind speed. Note again that the reasoning behind adding X-rays here was solely to reproduce the C IV1550 doublet in a few cooler stars but implemented in all for standardised comparison (see above). For the rest of the stars in our sample including X-rays does not improve the best fit quality. As such, the corresponding fit parameters become redundant and the uncertainties of X-ray luminosities extremely high.

4.4. Clumping

Figure 10 shows that clumping is ubiquitous. However, the scatter and error margins of the clumping parameters are large, and our sample is rather small, making robust quantitative interpretations challenging. The Pearson correlation coefficient to find correlations between all clumping factors, the β parameter, T_eff, and mass-loss rate, showed no significant correlations for the clumping parameters with any of the other parameters. However, the aforementioned trends, including the mass-loss rate, T_eff correlation, and the v_∞ − T_eff correlation, are reflected by the Pearson coefficient.

As by eye we perceived possible trends for the subsample of stars less than 30 kK, we performed a similar test for those stars. However, this did not reveal significant trends in any of the combinations save for the (f_cl, T_eff) and (f_cl, Ṁ) relations. First off we observed a tenuous increase in f_cl for a decreasing T_eff⁴. From line de-shadowing instability (LDI) simulations (Driessen et al. 2019) we expect a downward jump of f_cl as we go from O stars to B stars, but assuming our observed trend is real it could explain why this jump is difficult to observe. Nonetheless, we want to highlight again that this observed trend is rather tenuous. Because of the clear correlation between Ṁ and T_eff, we also find a slightly stronger correlation between f_cl and Ṁ⁵. This is inline with the trends observed in O stars by Brands et al. (2022), although they also highlight the delicate nature of this result.

Instead of correlations, we can study the general trends of the parameters not in respect to another. Figure 15 shows an approximate kernel distribution where each separate fit result is interpreted as a Gaussian with the uncertainty margin as the width of the Gaussian and each result is scaled by its 1/χ² value. Adding all these Gaussians together yields a distribution that is normalised to unity. This Gaussian allows us to talk about the behaviour of this parameter over the full sample when no clear trend is present. The approximate kernel distribution of the clumping factor peaks at a value 20. For low clumping factors, the likelihood stays strong, due to a single good fit, the tail towards the higher clumping factors decreases rapidly. The interclump density is less clearly peaked, but still has a visible peak around 0.4. The apparent long tail towards lower values also has to do with the displayed logarithmic axis, which visually prioritises small values and compresses the higher end of the axis. The velocity filling factor shows one of the possible pitfalls of this method, where one data point is dominating the distribution due to fitting well and having a low error-margin. The rest of the distribution is still useful, however, where we see a broad distribution with a relatively steep drop-off towards the lowest values.

Fig. 15. Approximate kernel distribution of the clumping parameters. In full lines we show the normalised approximate kernel distribution of the clumping parameter (fcl), the interclump density (fic), and the velocity filling factor (fvel). In dashed lines the contributions of the separate fits contributing to the distribution.

When taking a weighted average of the three clumping parameters we find: ⟨f_cl⟩ = 18 ± 10, ⟨f_ic⟩ = 0.38 ± 0.23 and, ⟨f_vel⟩ = 0.66 ± 0.18. These values are representing the part of the wind where it is most structured; at low velocities (below the gas sound speed) the atmosphere is always assumed smooth. How the three clumping parameters on average change with velocity is shown in Figure 16. The high mean value inferred for the interclump density sticks out, suggesting that we can expect that the interclump component comprises close to 40% of the average wind density. This is in stark contrast to alternative descriptions of wind clumping that are based on an effectively void interclump component. Non-zero f_ic have also been found in the previous studies (of O stars) by H21 and B22. These best fits, although quite a bit lower to the value quoted here (H21 f_ic = 0.13 ± 0.08, B22 $f_{\mathrm{ic}}=0.{13}_{-0.13}^{+0.15}$), seem at least not-pointing towards zero. As further discussed in the next section, this changes significantly the interpretation of the clumping parameters in terms of the relations between clumping factors, clump over-densities, and volume filling factors, and also invalidates the commonly applied assumption that the complete wind mass is compressed into dense and small clumps.

Fig. 16. Change of the 3 clumping parameters as a function of v/v∞. The dark blue line shows the averages, with the blue shaded region showing the 1-σ region of these averages.

The inferred mean value of the (maximum) clumping factor is higher than in models by Driessen et al. (2019) and also higher than the Galactic observational study by Bernini-Peron et al. (2023) where they find clumping factors of at most 2. We do note that both the fitting method and clumping description in their empirical study are different to ours, but the large difference is still noteworthy. Our derived clumping factors are, in contrast, lower than the sample averages for Galactic O-supergiants by H21 (f_cl ≈ 25) and for LMC O and WNh stars by B22 (f_cl ≈ 29 ± 15) who use the same clumping model. A downward trend with T_eff would be consistent with the linear analysis and non-linear instability simulations by Driessen et al. (2019), but establishing such a potential trend empirically would require a (significantly) larger sample than investigated here. In view of these empirical uncertainties and the fact that the theoretical clumping-study by Driessen et al. (2019) only used 1D simulations neglecting the influence of a turbulent stellar surface (Jiang et al. 2015; Debnath et al. 2024) on the clumpy outflow, it is at this point not meaningful to discuss more quantitative comparisons of empirical and theoretical clumping factors in this regime.

The empirical data suggest that clumping on average starts at low wind velocities within our sample, with several of the 1σ error ranges extending down towards our lower allowed bound at v ≈ a. While this is qualitatively consistent with previous empirical studies for hotter O stars (Puls et al. 2006; Cohen et al. 2011; Hawcroft et al. 2021; Brands et al. 2022), corresponding 1D instability simulations (Sundqvist & Owocki 2013), and 3D radiation transfer studies (Šurlan et al. 2013) also here we caution against more quantitative interpretations due to the fact that current wind clumping implementations may not be ideally suited for turbulent layers near the photosphere (see discussion in Debnath et al. 2024).

Finally, the fits for the velocity filling factor suggest that velocity-porosity effects indeed play a role in the spectrum formation of these winds, with a mean value that seems reasonable in view of general theoretical expectations (see discussion in Sundqvist & Puls 2018). The scatter is high, however, and no significant trends can be identified from the present data sample. We note in this context that there exists an independent observational indication for velocity-porosity (parametrised in the present study by a velocity filling factor) in BSG winds, namely via direct comparison of the depth of the blue and red absorption dips of unsaturated UV resonance doublets (Prinja & Massa 2013). Parsons et al. (2024) investigated this for stars in the ULLYSES sample and found similarly that there are clear indications for velocity-porosity across the sample but no significant trend with temperature. They also find very large scatter in their inferred optical depth ratios, as well as a strong temporal variation. Such temporal variations are not taken into account here as we combine two observations (UV+OPT) not taken simultaneously.

Comparing the obtained optical depths from Parsons et al. (2024) to the velocity filling factors derived here quantitatively is challenging, but equation 23 from Sundqvist et al. (2014) does allow us to approximate the vorosity factor (f_vor) from the optical depth of the blue and red absorption dips. This f_vor is then related to f_vel by equation 14 in Sundqvist & Puls (2018). The empirically derived blue and red optical depths from Parsons et al. (2024) then provides an independent estimate of f_vel, which can be directly compared to the values derived in this paper. The key qualitative point from this analysis is that velocity-porosity clearly is required in order to explain the blue to red optical depth ratios measured by Parsons et al. (2024). However, due to high uncertainties in these measured ratios, we still find significant quantitative error margins in f_vel (see Fig. A.3). Nonetheless, within these (admittedly very large) error margins, the two methods seem to provide results that overall are in reasonable agreement. This lends some support to the ability of the fast method utilised by Parsons et al. (2024) for identifying stars where velocity-porosity is important in their UV resonance wind line formation.

4.5. On the interpretation and validity of current wind clumping methods for spectroscopic studies

The generally high values found for the interclump density challenge some common assumptions and interpretations regarding wind clumping in massive stars. Namely, in most previous studies it has been assumed that the interclump medium is effectively void, so that all wind mass is contained within clumps and the relation between the clumping factor and clump volume filling factor f_vol is simply $f_{\mathrm{cl}}=f_{\mathrm{vol}}^{-1}$, accompanied by a characteristic clump over-density D = ρ_cl/⟨ρ⟩=f_cl. But for the typical values of f_ic ≡ ρ_ic/⟨ρ⟩ found here, these simple relations no longer hold. Instead we had to use the full relations for the stochastic two-component medium (Sundqvist & Puls 2018):

$$\begin{array}{c}\hfill f_{\mathrm{vol}}=\frac{{(f_{\mathrm{ic}}-1)}^2}{f_{\mathrm{cl}}+f_{\mathrm{ic}}^2-2f_{\mathrm{ic}}}\end{array}$$(10)

$$\begin{array}{c}\hfill D=\frac{{\rho }_{\mathrm{cl}}}{⟨\rho ⟩}=\frac{1-f_{\mathrm{ic}}(1-f_{\mathrm{vol}})}{f_{\mathrm{vol}}}\end{array}$$(11)

for a mean wind density

$$\begin{array}{c}\hfill ⟨\rho ⟩=f_{\mathrm{vol}}{\rho }_{\mathrm{cl}}+(1-f_{\mathrm{vol}}){\rho }_{\mathrm{ic}}\end{array}$$(12)

and clumping factor

$$\begin{array}{c}\hfill f_{\mathrm{cl}}\equiv \frac{⟨{\rho }^2⟩}{{⟨\rho ⟩}^2}=\frac{f_{\mathrm{vol}}{\rho }_{\mathrm{cl}}^2+(1-f_{\mathrm{vol}}){\rho }_{\mathrm{ic}}^2}{{(f_{\mathrm{vol}}{\rho }_{\mathrm{cl}}+(1-f_{\mathrm{vol}}){\rho }_{\mathrm{ic}})}^2}.\end{array}$$(13)

We point out first that while the above involves the quantities f_vol, f_cl, f_ic, and D, only two of these are independent. This is readily seen from the equation for mean density, which we may write in normalised form 1 = f_volD + (1 − f_vol)f_ic. Thus, if one chooses to set, say, f_ic and D, f_vol follows accordingly, and thereby also f_cl from the relations above. That is, while in our modelling we chose the input parameters f_cl and f_ic we may as well have chosen any pair of the above four quantities and then calculated out the others. Similarly, in the previously standard method of simply assuming f_ic = 0, only one of the above parameters is necessary to set; e.g. for an input f_vol we see directly that indeed $f_{\mathrm{vol}}^{-1}=D=f_{\mathrm{cl}}$.

To illustrate the situation when f_ic ≠ 0, we take our derived best-fit values for star Sk−67 ° 14; f_ic = 0.86 and f_cl = 8.0. This yields a clump volume filling factor f_vol = 0.003 that is now much lower than $f_{\mathrm{cl}}^{-1}=0.13$. Similarly we find for the characteristic clump over-density D = 51, whereas interpreting this the standard way would yield D = f_cl = 8. Moreover, identifying mass-contributions from the clumps and the interclump medium from the equation for ⟨ρ⟩ one gets for the former Df_vol = 0.15 and the latter (1 − f_vol)f_ic = 0.85; that is, we find here that most of the wind mass is actually not contained in the dense clumps, but rather in the interclump medium. Although this star lies on the extreme end of our sample, similar re-interpretations are necessary also when taking the weighted averages of our sample f_cl = 18 ± 10 and f_ic = 0.38 ± 0.23, whereby we obtain f_vol = 0.022 ± 0.020 (vs. $f_{\mathrm{cl}}^{-1}=0.06\pm 0.03$), D = 28 ± 27 (vs. f_cl = 18 ± 10) and that the interclump medium contributes 37 ± 22% of the wind mass. These general tendencies arise because we typically find a significant interclump density together with a rather high rms over-density $\sqrt{⟨{\rho }^2⟩}/⟨\rho ⟩$ from the observational data. This combination forces clumps to be confined into very small volumes, so that relative contributions from the clump and interclump media to the total wind mass become approximately 1 − f_ic and f_ic, respectively. Additionally, since typically f_vel is significantly higher than f_vol, it means that these spatially very confined clumps within our formalism are quite spread out in velocity space.

H21 find values f_ic ∼ 0.15 − 0.3 from their sample of Galactic O supergiants. B22 split their sample of O and WNh stars into two categories according to stellar luminosity, and find sample averages for their low and high luminosity stars f_ic = 0.07 ± 0.06 and f_ic = 0.28 ± 0.21, respectively. Our sample average here is thus noticeably higher, but still in overall qualitative agreement with these previous results for high luminosity O-stars regarding the possibility of a significant non-clump wind mass (we note in this respect that, originally, the discussion in Sundqvist & Puls 2018, based largely on 1D simulations, argued for low “standard” values on order ∼0.01–0.2 for this parameter). The increase of f_ic in our study compared to H21 and B22 could be due to the different temperature region or due to slight differences in the fitting routine such as different fitting assumptions. On the other hand, in an independent study Šurlan et al. (2013) used 3D Monte-Carlo radiation transfer, in combination with the 1D code PoWR, for a Galactic O star analysis and showed that the inclusion of interclump density was needed to reproduce the P V 1118 doublet. They got best correspondence for relatively high interclump density (0.1 − 0.4⟨ρ⟩), however including the interclump density reduced the effect of the clumping factor on the line profiles. These counteracting effects probably also increase the error-margin of our clumping parameters, but due to the global fits of the GA method such inter-dependencies are all taken into account.

Overall, the general picture that emerges is different from that which has been assumed in spectroscopic studies of massive-star wind clumping previously; a significant fraction of the wind mass seems actually not to be contained in small-scale clumps, but rather in the medium in between these. While this is in stark contrast with the assumptions underlying the vast majority of previous spectroscopic studies, it actually correlates rather well with recent theoretical results from multi-dimensional radiation-hydrodynamic simulations. The 3D models of Wolf-Rayet stars by Moens et al. (2022) indeed find that approximately half of the wind mass is contained in parcels having densities lower than the mean density of the wind, illustrating again the general issue with assuming a wind completely dominated by dense clumps.

However, in these recent simulations it is also found that the density distributions do not resemble a two-component medium, but rather have Gaussian-like distributions (likely log-normal ones, see Owocki & Sundqvist 2018 and Schultz et al. 2020) where the most probable density is quite close to the mean and the dispersion is large. Since similar results are suggested also by the 2D O-star simulations by Debnath et al. (2024) and Driessen (2022), the latter including effects of the (LDI, Owocki & Rybicki 1984) (see their Fig. 8.3), this may indicate issues also with the generic assumption applied here of an effective two-component medium consisting of clumps and an interclump medium. Additionally, the models by Jiang et al. (2018) and Debnath et al. (2024) that consider deep sub-surface layers clearly indicate that also the photosphere is very variable and structured, whereas in the clumping description applied here we assume an inner boundary for possible structure formation at the wind sonic point (see Fig. 16). In this respect, we note further that (as mentioned in Sect. 2.1) in our current approach the ionisation balance is only derived for an effective medium, and thus not for the interclump and clumped components separately. Overall, in view of these results a general re-calibration of methods used to account for wind clumping effects in spectroscopic studies may thus be needed (see Owocki & Sundqvist 2018 for a first attempt that focuses on transfer effects arising from the turbulent density structures typically seen in multi-D simulations of hot star winds).

5. Summary and outlook

We have analysed the optical X-shooter and UV ULLYSES spectra of 15 OB supergiants in the LMC using the model atmosphere code FASTWIND and the GA fitting approach Kiwi-GA, resulting in 18 consistently determined stellar (seven) stellar and (11) wind parameters. Derived spectroscopic masses in our sample range from 15–75 M_⊙ and effective temperatures lie in the range 35–14 kK with uncertainty margins of around ±1500 K per star. Our sample of stars have been selected with the goal of determining if the so-called bistability jump, an upward jump in mass loss towards cooler photospheres within the observed temperature regime, is observable when the parameter degeneracy of wind clumping is broken by the availability of UV data. As such, we focused especially on determining accurate mass-loss rates. Our derived rates range from −7.7 to −6 log₁₀(M_⊙/yr), taking into account inhomogeneities of the wind in the most detailed way that current 1D atmosphere codes allow.

By determining the mass-loss rates of all stars within ±0.3 dex for most of them, we could see that no sudden mass-loss increase in this effective temperature regime is present in the empirical analyses. From comparisons to different theoretical predictions we observed that the jump described by Vink et al. (2001) (which is the current standard recipe to include in applications such as stellar evolution modelling) drastically overestimates the mass-loss rates on the cool side of their predicted jump. For our sample, we derived, similarly to how we derive the uncertainty margins on the cumulative clumping parameters in Figure 15, an overestimation by a factor of ∼24. In contrast, on the hotter side of the predicted jump these same predictions are typically rather well aligned with the empirically derived rates, on average the offset is ∼0.8 times the empirically derived rates. When comparing to the alternative theoretical rates by Björklund et al. (2023), which decreases monotonically with temperature, and we did not find any mass-loss jump (see discussions in previous sections). We noticed that the trend with T_eff of these predictions is rather well aligned with the empirical data, showing a relatively constant underestimation and an offset of ∼0.5 times the empirically derived rates. If a jump were present in the observed sample it would not be possible to find such a constant offset with a prescription that is only decreasing with decreasing T_eff. We note that for the computation of these ratios, we did not take into account the uncertainty margins of the prescriptions, which are always present as they are, on purpose, simple fits of a grid of models, resulting in some scatter between the fits and the actual model results. As such, the quoted error margins are likely to be slightly underestimated. Further, the Krtika et al. (2024) prescriptions also align well overall with the observed mass-loss rates, but our current dataset is not optimised for the testing of the weaker mass-loss bump present at cooler temperatures, as it is also very localised( decreasing aggressively when moving towards even cooler temperatures; see Fig. 1). As such, we can neither confirm nor exclude a presence of the modest bump seen in the models by Krtika et al. (2024). The localised behaviour in temperature and the comparatively weaker strength of the bump would make this feature relatively hard to verify in general, and such verification is impossible with our sample due to its temperature position.

Our derived terminal wind speeds are very much in line with the findings of Hawcroft et al. (2024b). The terminal wind speeds show a very linear dependence on the effective temperature even when masses and luminosities vary heavily. When trying to derive v_∞ from the C IV1550 doublet, which is typically the main diagnostic for terminal wind speeds, we noticed that stars with optical spectra that point clearly towards cool B stars sometimes still had strong C IV lines, although they are too cool to ionise a sufficient amount of carbon atoms to C IV in order to have strong C IV lines. To solve this we included ionisation due to additional X-rays in the wind also for stars down to T_eff ≈ 15 kK. This produced a sufficient amount of C IV atoms to model the C IV1550 doublet, whereas, interestingly enough, it neither changed the Si IV wind lines nor had a noticeable change on other fit parameters such as mass loss or clumping parameters. Such effects from X-ray ionisation for BSGs mimic the well known so-called wind ‘super-ionisation’ for O-stars (often discussed in the context of OVI, Cassinelli & Olson 1979; Macfarlane et al. 1994; Krtika & Kubát 2009). However, even in B-stars this effect has been noted although less commonly mentioned (Macfarlane et al. 1994).

We derived constraints on a collection of clumping parameters that describe inhomogeneities in the wind (f_cl, f_ic, f_vcl, start, f_vcl, max, f_vel). For f_cl we obtained a sample average 18 ± 10, where the large error margin is primarily due to the large spread of best-fit clumping factors. Perhaps most strikingly in our empirical study of wind clumping are the high values of f_ic. Here we obtained a sample average of 0.38 ± 0.23, where again the large uncertainty is primarily due to the large spread of the best-fit values; indeed, for individual stars we found values up to $0.{92}_{-0.22}^{+0.02}$. These high interclump densities may be viewed as surprising in the sense that a widely used assumption in 1D atmospheric and spectroscopic modelling has been that all wind mass is contained in the clumps. In contrast, we find that typically about half of the wind mass is actually within the interclump regions. On the other hand, in view of recent multi-dimensional radiation-hydrodynamical simulations of hot star atmospheres with winds (Moens et al. 2022; Debnath et al. 2024), this may not be so surprising after all, since these simulations tend to display Gaussian-like density distributions centred quite close to the mean rather than bimodal distributions with dominant over-dense clumps (which had been the prevailing thought based largely on 1D LDI simulations; e.g. Sundqvist & Owocki 2013; Driessen et al. 2019). In the simulations by Moens et al. (2022), the accumulative density distributions indeed show that almost half of the wind mass is contained in parcels that have densities lower than the mean density of the wind. Since similar results seem to be indicated also by multi-D simulations including the LDI (Driessen 2022) but are yet to be further quantified, our results indicate a need to rethink how the radiation-driven wind is described in 1D codes used for spectroscopic modelling. This very basic assumption greatly influences not only the density structure of the wind but also its ionisation structure.

An interesting follow-up to these results will be the study of similar objects in the SMC, for which there are equivalent data thanks to the ULLYSES program (Roman-Duval et al. 2020) and the X-shooting ULLYSES (Vink et al. 2023; Sana et al. 2024). The goal will be to find if these empirical mass-loss and clumping trends continue to a lower metallicity. Additionally, with a very comparable methodology for SMC and LMC samples, it should be possible to also study the mass-loss behaviour scaling with metallicity. To this end, one would be able to compare this study with SMC studies on the same data in other temperature regimes (Parsons et al. 2024; Backs et al. 2024; Bernini-Peron et al. 2024) and other LMC samples studying the O-star regime (Brands et al. 2022; Hawcroft et al. 2024a) (Brands et al., in prep.) and comparing to studies in the Milky Way (Hawcroft et al. 2021; Bernini-Peron et al. 2023; de Burgos et al. 2024). A first look at how the sample range can easily be expanded by studies using similar fitting routines and methodology is shown in Figure 17. Brands et al. (in prep.) studied O stars in the LMC and therefore complements our sample to the hotter side. Figure 17 shows our sample again but the thinner crosses show the mass-loss rates and effective temperatures of the study by Brands et al. (in prep.). We can see that the trend of increasing mass-loss rates with increasing temperatures and Γ_e is still present in this larger sample. As expected the error margins in this high temperature regime are smaller, as we already saw in our sample, where the hotter stars have lower error in mass-loss rate. Object sk−68 ° 155 is present in both samples and have very comparable fit-parameters in both studies (with the exception of the carbon abundance).

Fig. 17. Expanding the present sample with study by Brands et al. (in prep.).

Finally, another interesting step would be to further study the effects of not including a strong bistability jump in evolution calculations. The lack of a prominent jump in observations means that stellar evolution codes currently being used as standards to describe single massive star evolution are also likely to significantly overestimate the loss of mass and angular momentum in this temperature and luminosity regime, perhaps then even influencing the viability of single stripped star formation (see, e.g., discussion in Björklund et al. 2023) or changing surface abundances (Josiek et al. 2024), and rotation Keszthelyi et al. (2017), Britavskiy et al. (2024). However, in this respect we note that for even higher luminosity to mass ratios than those examined in this paper, multi-dimensional models in this temperature regime suggest that strong winds can already be ignited from hot and optically thick sub-surface layers (Jiang et al. 2018), which may enhance mass-loss rates beyond the standard luminosity scalings of the stars with line-driven winds studied in this paper. And indeed, empirically this then approaches the region where we find the so-called luminous blue variables, which are well known to undergo vigorous and variable mass loss (see overview in Vink 2022). We argue that it should be a key focus for future (both empirical and theoretical) studies to try to better constrain the mass-loss properties of this (very) high-luminosity region in the HR diagram.

Data availability

The fits of all stars discussed here and the distribution of the GA models are available at https://zenodo.org/records/13948998

Acknowledgments

The resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government. We would also like to thank professor Leen Decin for her contribution to this work. O.V. and J.S. acknowledge the support of the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N and FWO grant G077822N and, KU Leuven C1 grant MAESTRO C16/17/007. O.V. also acknowledges the FWO travel grant under grant number K210124N. J.S., F.B., and P.S. further acknowledge the support of the European Research Council (ERC) Horizon Europe under grant agreement number 101044048. B.K. acknowledges the support from the Grant Agency of the Czech Republic (GAČR 22-34467S) and RVO:67985815. AACS and MBP are supported by the Deutsche Forschungsgemeinschaft (DFG – German Research Foundation) in the form of an Emmy Noether Research Group – Project-ID 445674056 (SA4064/1-1, PI Sander). AACS and MBP further acknowledge funding from the Federal Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments. RK acknowledges financial support via the Heisenberg Research Grant funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant no. KU 2849/9, project no. 445783058. F.N., acknowledges support by grant PID2022-137779OB-C41 funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe”.

References

Appendix A: Remarks on the results

The goal of this paper was not to focus on the CNO abundances in psarticular. However, as they are needed fit parameters to derive the other parameters we shall briefly highlight what we can learn from them. The CNO abundance as a whole is expected to stay constant during the evolution of the star as one element gets converted into the other. To this effect to check the validity of our CNO abundances we have made figure A.1 showing the sum of CNO massfractions and propagating the uncertainties for each of them. As expected, the sum of CNO abundances with the uncertainty stays close to the base sum.

Fig. A.1. Cumulative CNO mass fraction; Here the CNO abundance is plotted for each star, (C in green, N in yellow, and O in red) in the sample with in black-lines the uncertainty on the sum of the CNO abundance. The blue dots indicates the expected cumulative CNO massfraction from Vink et al. (2023)

Fig. A.2. Results of the micro-turbulence variation tests. We show for three stars, four different fits. These 4 fits are very similar except for the photospheric micro turbulence, which is changed from a fixed 10 km/s we used in the paper (blue), to 5 km/s (red), to 20 km/s (purple), and left as a fit parameter (green). We plot the 1-σ uncertainty of 5 parameters besides the vmic. These are: effective temperature, mass loss rate, and the CNO abundance. The three different stars each get a column in the 6 different panels and the 4 different fits are slightly overlapping in this column.

Fig. A.3. Velocity filling factor in this study compared to Parsons et al. (2024). This figure shows a comparison between the velocity filling factor we derived here, and the velocity filling factor from Parsons et al. (2024). We compute the velocity filling factor from optical depths of the red and blue absorption dips obtained by Parsons et al. (2024) using equation 23 from Sundqvist et al. (2014). When the ratio of the red and blue absorption dips reached values above 2 we excluded it here as this results in non-physical fvel.

Appendix B: Line list

Note: The first column shows the atom and its ionisation stage which is responsible for the transition. The second column shows the corresponding wavelength with possible multiplets. The third column shows where to find this line in the fit summary which is available on Zenodo.

All Tables

All Figures

Fig. 1. Comparison of mass-loss rates from prescriptions by Vink et al. (2001), Björklund et al. (2023), Krtika et al. (2024) over the temperature range from 60 to 12.5 kK. Key stellar parameters in the left panel are: log10(L⋆/L⊙) = 5.3, M⋆ = 25 M⊙ and, Z = 0.5 Z⊙. This results in Γe = 0.21. In the right panel adopted properties are log10(L⋆/L⊙) = 5.6, M⋆ = 25 M⊙ and, Z = 0.5 Z⊙, yielding Γe = 0.42. The full lines signify that the result of the fit-formulate given by these authors are shown rather than the actual underlying model-sets. For the Vink et al. (2000) rates, for which a mass-loss jump is identified to cover the regime ∼22.5 − 27.5 kK, enforcing a strict discontinuity (i.e. at a single temperature) implies that the size of the Ṁ-jump is somewhat enhanced as compared to the underlying models.

In the text

Fig. 2. Comparison of diagnostic lines for four stars. The stars (Sk−67 ° 107, Sk−68 ° 41,Sk−69 ° 52,RMC-109) are organised from high (left) to low (right) effective temperature (in spectral type: O8.5 II, B0.7 Ia, B2 Ia, B5 Ia). The lines are selected due to their roles as good diagnostics for stellar or wind parameters. From top to bottom the lines are sorted by wavelength. We note that the line names here correspond to the names in Table B.1. The green line shows the best fit and the green shaded region shows the 1-σ uncertainty interval.

In the text

Fig. 3. Distribution of the quality of the fits for stellar parameters of star Sk−68 ° 41. In dark grey the 1-σ interval is shown while the light grey shows the 2-σ interval. The red line shows the best fit value. We only display the C-abundance as the other elements have similar fitting curves. As mentioned in Section 2.3 and the following sections we derive the maximum projected rotational velocity from the optical-only fit which is why there are only 21 generations.

In the text

	Fig. 4. Hertzsprung–Russell diagram showing all stars in the sample and with the colour showing the derived spectral masses. The diamonds mark all the stars with spectral types Ia, while the ‘v’ symbols mark the lower luminosity classes. The stellar evolution tracks are from the MIST data base without any rotation (Dotter 2016; Choi et al. 2016; Paxton et al. 2010, 2013, 2015).
In the text

	Fig. 5. Derived maximum projected rotational velocity as a function of effective temperature. We note that this maximum value in practice reflects the combined effect of all sources of macroscopic line broadening (see text).
In the text

Fig. 6. Clumping influence on the spectral lines. The figures show the influence a change in fcl has on the line profile of on one hand the H α recombination line and on the other hand the Si IV 1400 resonance doublet. Both model spectra are using almost the same input, they only differ in mass-loss rate and clumping. The blue line both has Ṁ = 8 · 10−8 M⊙/yr, and is highly clumped (fcl = 40). The orange line is a smooth outflow but has a significantly higher mass-loss rate (Ṁ = 5 · 10−7 M⊙/yr); the product $\sqrt{f_{\mathrm{cl}}}\dot{M}$ is thus the same for both lines. If we assume an optically thin clumping over the complete line forming region, the 2 H α lines would show perfect agreement. However, this is not the case as we use models that relax the optically thin assumption and are only clumped above a fitted onset-velocity.

In the text

Fig. 7. Quality of all models in a GA run ranked based on the inverse of the RMSEA value as a function of the mass-loss rate. The models for these plot are based on a run trying to fit the star Sk−68 ° 41 (SpT B0.7 Ia). Each point here is a separate model where the colour defines the generation in which this model was run. The dark grey region highlights the 1-σ uncertainty interval and the light grey region is the 2-σ uncertainty.

In the text

	Fig. 8. Correlation between the clumping factor and the mass-loss rate for star sk−68 ° 41. The colour shows the 1/RMSEA value.
In the text

	Fig. 9. Derived terminal velocities, wind acceleration parameter β, and maximum wind micro turbulence, as a function of effective temperature. The colour for all panels shows the Γe value. The first panel also shows the linear fit to the terminal velocities found in this work in red compared to the fit (in blue) from Hawcroft et al. (2024b).
In the text

Fig. 10. Clumping parameters with respect to the effective temperature. The panels from top left to bottom right show the clumping factor, interclump density (in units of mean density), velocity filling factor and the onset clumping velocity (vcl, start) and the velocity of maximum clumping (vcl, max). The three first panels colour according to Γe for the stars, while in the last plot the colour indicates the difference between vcl, start and vcl, max. For each star in this plot the two points are connected by a dashed line (not visible in some objects as the error margins are covering the dashed line).

In the text

	Fig. 11. Distribution of the quality of the fits for clumping parameters. We show a typical distribution of the models at end of a GA fitting run. The colour indicates the generation of the model and the grey areas indicate the 1σ and 2σ uncertainty intervals.
In the text

Fig. 12. Derived mass-loss rates for our sample as function of effective temperature. The top plot shows the derived mass-loss rates and he colour shows the derived classical Eddington parameter for each star. The bottom plot shows the ratio between the theoretical prescriptions by Vink et al. (2001), Björklund et al. (2023), Krtika et al. (2024) and the empirically found values. The green line indicates where these are equal.

In the text

	Fig. 13. Terminal wind speed divided by escape speed plotted over temperature. The mass used to determine the escape speed is the effective spectral mass.
In the text

	Fig. 14. Three versions of the same doublet line per panel. In black we show the data from Sk−70 ° 16 a B2 II star, in blue the best fit when not allowing for X-rays, and in red is the best fit when allowing for X-rays. The top panel shows the C IV 1550 line, and the bottom panel shows the Si IV1400 line.
In the text

	Fig. 15. Approximate kernel distribution of the clumping parameters. In full lines we show the normalised approximate kernel distribution of the clumping parameter (fcl), the interclump density (fic), and the velocity filling factor (fvel). In dashed lines the contributions of the separate fits contributing to the distribution.
In the text

	Fig. 16. Change of the 3 clumping parameters as a function of v/v∞. The dark blue line shows the averages, with the blue shaded region showing the 1-σ region of these averages.
In the text

	Fig. 17. Expanding the present sample with study by Brands et al. (in prep.).
In the text

	Fig. A.1. Cumulative CNO mass fraction; Here the CNO abundance is plotted for each star, (C in green, N in yellow, and O in red) in the sample with in black-lines the uncertainty on the sum of the CNO abundance. The blue dots indicates the expected cumulative CNO massfraction from Vink et al. (2023)
In the text

Fig. A.2. Results of the micro-turbulence variation tests. We show for three stars, four different fits. These 4 fits are very similar except for the photospheric micro turbulence, which is changed from a fixed 10 km/s we used in the paper (blue), to 5 km/s (red), to 20 km/s (purple), and left as a fit parameter (green). We plot the 1-σ uncertainty of 5 parameters besides the vmic. These are: effective temperature, mass loss rate, and the CNO abundance. The three different stars each get a column in the 6 different panels and the 4 different fits are slightly overlapping in this column.

In the text

Fig. A.3. Velocity filling factor in this study compared to Parsons et al. (2024). This figure shows a comparison between the velocity filling factor we derived here, and the velocity filling factor from Parsons et al. (2024). We compute the velocity filling factor from optical depths of the red and blue absorption dips obtained by Parsons et al. (2024) using equation 23 from Sundqvist et al. (2014). When the ratio of the red and blue absorption dips reached values above 2 we excluded it here as this results in non-physical fvel.

In the text