© The Authors 2024

1 Introduction

Feedback from massive stars plays an important role in the thermal and dynamical evolution of the interstellar medium and galaxy halos (e.g., Hopkins et al. 2011, 2018; Andersson et al. 2020). It is provided through jets, (ionizing) radiation, powerful stellar outflows, and, when their cores ultimately collapse, through supernovae (e.g., Kuiper & Hosokawa 2018; Dale & Bonnell 2011; Luisi et al. 2016; Geen et al. 2020; Efstathiou 2000; McLeod et al. 2019).

The mass and angular momentum that are lost through their stellar outflows (processes that persist throughout their lives) strongly influence the evolution of massive stars (e.g., Puls et al. 2008; Langer 2012; Renzo et al. 2017; Vink 2022). At Galactic metallicities, mass loss causes the mass at the end of carbon burning to be 25–40% of the initial mass for stars that started their lives with 20 M_⊙ or more (Ekström et al. 2012). The cumulative mass loss experienced by a massive star throughout its evolution is found to be a function of the initial stellar metal content (see below).

Metal-poor hot luminous stars have weaker winds, as has been shown empirically for the Large Magellanic Cloud (LMC, e.g., Brands et al. 2022), the Small Magellanic Clould (SMC, e.g., Ramachandran et al. 2019), and for sub-SMC (e.g., Garcia et al. 2014; Bouret et al. 2015). One consequence of this is that properties of core-collapse supernovae are expected and are found to differ for galaxies with low and high metal content. Long-duration gamma-ray bursts (e.g., Modjaz et al. 2008; Heger et al. 2003; Savaglio et al. 2009), superluminous supernovae (e.g., Chen et al. 2017), and pair-instability supernovae (e.g., Young et al. 2010), for instance, favor lower-metallicity environments. It also implies that the first stars to form in the Universe (e.g., Hirano et al. 2015; Hosokawa et al. 2016) must have lost relatively little mass (e.g., Kudritzki 2002), with important consequences for their end-of-life products (e.g., Marigo et al. 2003).

An accurate understanding of the mass loss from massive stars is therefore very important, including its dependence on the metal content. On the theory side, the hot-phase driving mechanism has been identified to be radiation pressure on the metal lines of mostly CNO and iron-group elements that are clustered in the ultraviolet, where the stars emit the bulk of their radiation (e.g., Lucy & Solomon 1970; Castor et al. 1975; Abbott & Lucy 1985; Pauldrach et al. 1986; Puls et al. 1996). Combined stellar atmosphere and hydrodynamical methods that take into account several 10⁵ of these lines have been developed to make quantitative predictions (de Koter et al. 1997; Vink et al. 2001; Sander et al. 2017; Krtička & Kubát 2018; Björklund et al. 2021; Gormaz-Matamala et al. 2022a).

The behavior of mass loss versus metallicity was probed by studies of the solar metallicity and of the LMC (Z = 0.5 Z_⊙), and SMC (Z = 0.2 Z_⊙) metallicities. When they were based on relatively high-luminosity sources (log L/L_⊙ ≳ 5.4), where wind signatures in Hα and He II 4686 are clearly visible, the studies inferred a luminosity-independent behavior $\dot{M}$ ∝ Z^0.5–0.8 (Mokiem et al. 2007; Marcolino et al. 2022). To constrain the mass-loss properties for lower-luminosity stars, the spectral analysis relies strongly on the ultraviolet spectral range, where UV resonance lines such as C IV 1548, 1550 Å provide more sensitive probes of the wind strength. Marcolino et al. (2022) used a compilation of UV and optical studies to first address this parameter space and tentatively found a weaker or even vanishing dependence on Z: $\dot{M}$ ∝ Z^0.1 at log L/L_⊙ = 4.5. They stressed that their findings should be tested with larger samples that include more accurate determinations of the terminal wind velocity.

The combination of the Hubble Space Telescope (HST) UV Legacy Library of Young Stars as Essential Standards (ULLYSES) (Roman-Duval et al. 2020) and Very Large Telescope (VLT) X-Shooting ULLYSES (XShootU) (Vink et al. 2023) programs provides high-quality spectra of OB stars in the Magellanic Clouds with an extensive wavelength coverage from the UV to the optical. We analyze a sample of O4 to O9.5 dwarfs, giants, bright giants, and supergiants in the SMC from these programs. The spectra allow us to constrain the photospheric conditions and a range of wind properties, including inhomogeneities in the outflow and the terminal velocity. To gain insight into the metallicity dependence of our findings, we compare them with similar studies of Galactic (Hawcroft et al. 2021) and LMC stars (Hawcroft et al. 2024a, Brands et al., in prep.).

A major challenge in the empirical determination of mass-loss rates is the inhomogeneity in the wind, which is also referred to as clumping. These clumps are overdense regions that are surrounded by lower density material. Clumps affect the strength of spectral lines both directly and indirectly. They affect the lines directly in the case of recombination lines, such as Hα and He II 4686, the strength of which depends on the mean of the square of the density of the medium; and they affect them indirectly in the case of lines that trace ionic species whose ionization balance shifted significantly because the different density. Scattering lines of abundant ionic species, such as C IV 1550 Å, may be unaffected by clumping depending on the stellar temperature, however.

Because of these different effects, clumping has been suggested and used as solution for the discrepancies between inferred mass-loss rate from different diagnostic lines (e.g. Bouret et al. 2003; Fullerton et al. 2006). Whether wind clumping itself is a function of metallicity and/or luminosity has been little explored so far. If this is the case, it may affect the derived mass-loss properties and should thus be accounted for when the empirical $\dot{M}$(Z) or $\dot{M}$(L, Z) relations are derived. This challenging problem is a goal of this study as well. So far, Driessen et al. (2022) reported a weak relation between the metallicity and the clumping factor in 2D line-deshadowing instability simulations. Parsons et al. (2024) reported a similar trend for structures that appeared in single-epoch observations of wind lines in the UV spectra of B supergiants.

This paper is organized as follows. We describe the data and the sample in Section 2, along with the data preparation and normalization process. The model and fitting routine are presented in Section 3. Section 4 shows and describes the best-fit line profiles and the corresponding stellar and wind parameters. In Section 5, we briefly discuss the evolutionary stage of the sample stars and present the implications for $\dot{M}$(L, Z), and the wind structure parameters. Finally, Section 6 lists the main conclusions.

2 Data and sample

Our sample consists of 13 SMC O-type giants to supergiants that are available in ULLYSES DR5 (Roman-Duval et al. 2020) and XShootU eDR1 (Vink et al. 2023; Sana et al. 2024, XShootU I and II). The ULLYSES data consist of a mix of FUSE¹, HST/COS², and HST/STIS³ spectra, with both archival and new observations. The ULLYSES sample was selected to include presumably single stars (Vink et al. 2023). However, the presence of binaries in our sample cannot be excluded. We used the high-level science products as reduced by Roman-Duval et al. (2020). The XShootU data consist of VLT/X-shooter spectra covering the UVB and VIS arms (~3100–10 000 Å). The sample was selected based on data availability in the optical and UV wavelengths, where we required that the UV data cover at least the wavelength range 1150–1700 Å. The sample of stars with their spectral type and available UV observations is listed in Table 1. The literature spectral types we used for target selection are included, as are updated spectral types using the XShootU data from Bestenlehner et al., in prep.

All FUSE data were taken with the LWRS aperture, which has a spectral resolution of ~17 500 and covers 905–1180 Å. For AV 307 alone, the data were taken with HST/COS using the G130M/1096 grating. This grating covers the wavelength range from 940–1240 Å and has a resolution of R ~ 6000 around λ ~ 1120 Å. The HST/COS observations using the G130M/1291 and G160M/1611 gratings cover the wavelength range 1141–1783 Å and have a resolution from 11000 to 19000. The HST/STIS data that were obtained using the E140M and E230M gratings have a resolution of 45800 and 30000, respectively. The E140M grating covers the wavelength range from 1141 to 1708 Å and the E230M grating covers the wavelengths 1608–2366 Å. When multiple data sets covered a modeled feature, the data were chosen based on the signal-to-noise ratio, resolving power, and possible systematic effects. The optical data consist of the UVB and VIS arms of VLT/X-shooter. The UVB arm covers 3100–5500 Å and has a resolution of 6700 for the chosen slit width of 0.″8. The VIS arm covers 5500–8000 Å with a resolution of 11400 with a slit width of 0.″7.

The signal-to-noise ratio (S/N) of the data varies depending on the observation and wavelength range. The optical data typically have a high signal-to-noise ratio of >100. The UV data have an S/N ~ 20 on average. In some cases, a higher S/N was obtained. The lowest signal-to-noise ratio is found in the FUSE data, where it can be as low as ~5.

2.1 Photometry

Photometric information of the stars was used to determine the line-of-sight extinction. The extinction and the distance were used to determine the absolute K_s-band magnitude (listed in Table 1), which was used as luminosity anchor in the analysis. All photometric data were taken from the Bonanos et al. (2010) catalog. From this catalog, we used the U, V, B, and I photometry of the Magellanic Cloud Photometric Survey (Zaritsky et al. 2002), and the J, H, and K_s bands from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) and the InfraRed Survey Facility (Kato et al. 2007). For the luminosity anchor (see Section 2.2), we prioritized the 2MASS photometry. Finally, we used the extinction curves of Fitzpatrick (1999) to parameterize the line-of-sight extinction. We adopted R_V = 3.1 and determined the A_V by fitting the photometry to Castelli-Kurucz model spectral energy distributions (SED) (Castelli & Kurucz 2004). We opted to fix R_V as only marginal differences are found when it is made a free parameter (Sander et al. 2024). For this part of the analysis, the temperature of the model was selected based on the literature spectral type, as listed in Table 1. These spectral types were then used to obtain an approximate temperature following Martins et al. (2005).

2.2 Data preparation

The spectrum-fitting routine employed in this work (see Section 3) requires normalized line profiles. Therefore, the observed spectra need to be normalized. Incorrectly identifying the continuum may affect the strength and shape of spectral lines, so it is of great importance to determine it as accurately as possible. To do this, we attempted to improve on the normalization done by Sana et al. (2024, XShootU II) in eDR1 by pursuing the following normalization routine.

We divided the observed flux-calibrated spectrum by a normalized CMFGEN model (Hillier & Miller 1998). Ideally, that is, if the model fits perfectly, this results in a featureless pseudocontinuum with noise. This pseudo-continuum was fit with a polynomial. The normalized spectrum was then obtained by dividing the observed spectrum by the polynomial. The CMFGEN model parameters were selected from a limited grid of models, in which we used the model with the lowest χ² when we compared the normalized data to the normalized model spectrum. In this process, we masked diverging features that are not properly covered in the model, such as wind-sensitive lines and interstellar features, so that they did not affect the polynomial fit or χ² determination. The normalization procedure was done locally, that is, around spectral lines of interest, in the optical, and per grating in the UV spectral range (see below).

By using normalized spectra to fit our models and by using a polynomial to fit the pseudo-continuum, we minimized the possible effect of uncertainties in the extinction on the analysis. This was aided by choosing the K_s-band magnitude as the luminosity anchor, as the effect of extinction toward our targets is limited at longer wavelengths.

For the UV data, the normalization was performed simultaneously with the determination of the radial velocity, using the same χ² analysis as in the selection of the model. The optical X-shooter data were corrected for their radial velocity using a cross-correlation with a selection of hydrogen and helium lines. This should give corrections accurate to ~10 km s⁻¹, depending on the signal-to-noise ratio and projected rotational velocities. The UV spectrum of hot stars in the SMC is rich in Fe lines. These lines cover the whole spectral range, making it hard to identify the continuum between the lines. Therefore, the normalization was applied to the full spectrum separately for each grating.

The differences between the normalization performed in XShootU eDR1 (Sana et al. 2024, XShootU II) and this work are typically small (<1%). We ascribe them to some broad local features that may have affected the global normalization of Sana et al..

Strong interstellar Lyα absorption is present in all spectra. The wings of this absorption overlap with the N V λ1240 feature and are often wide enough to affect the C III λ1176 and C IV λ1169 features. Since interstellar Lyα is not part of our modeling, we corrected for it. We fit a Voigt-Hjerting function to the affected data (Tepper-García 2006). This resulted in good fits, which allowed us to set the continuum around these features at unity. Features in the spectral regimes of interest that were not modeled, such as missing lines or elements, diffuse interstellar bands, and interstellar absorption lines, were clipped from the observed spectrum. In this way, they did not affect the fitting efforts.

3 Methods

We aim to determine detailed stellar and wind parameters of our targets by calculating model spectra and comparing those to the observations. The model spectra were calculated using the model atmosphere code FASTWIND. As the number of models required to fully explore the parameter space is substantial, we used a genetic algorithm to efficiently explore the parameter space and converge to the optimal parameters and their associated uncertainties. Sect. 3.1 describes FASTWIND in more detail. The genetic algorithm is covered in Sect. 3.2. Finally, in Sect. 3.3, we describe the fitting approach and free parameters.

3.1 Fastwind

FASTWIND⁴ is a stellar atmosphere and radiative transfer code that is optimized for hot stars and their winds (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2012; Carneiro et al. 2016; Sundqvist & Puls 2018). The model assumes nonlocal thermodynamic equilibrium (NLTE) and includes line-blanketing. FASTWIND aims to minimize computational cost. To this end, the elements are split into explicit and background species. The explicit elements are treated in detail in the comoving frame, while the background elements are only used to account for line-blocking and line-blanketing. For this reason, only explicit elements can be used to produce diagnostic line profiles. The explicit elements used here are H, He, C, N, O, Si, and P. FASTWIND uses a pseudo-hydrostatic photosphere that smoothly transitions into a trans-sonic wind. The wind is then described in a parameterized way, with a specified mass-loss rate and a radially increasing smooth-wind β-velocity law that asymptotically approaches the terminal velocity v_∞. Additionally, we used the clumping prescription that we refer to as “optically-thick”, as described in Sundqvist & Puls (2018), to account for small-scale inhomogeneities in the outflow. In this prescription, the wind and clumps are not assumed to be optically thin, which allows for both optically thin and thick parts. The clumped wind is described by a set of six wind-structure parameters, f_cl, f_ic, f_vel, v_cl,start, v_cl,max, and v_windturb, which we introduce briefly (for more detailed descriptions, see Sundqvist & Puls 2018; Brands et al. 2022).

The medium was assumed to consist of two components: regions in which the density ρ_cl is relatively high (clumps), and regions in which the density ρ_ic is relatively low (inter-clump medium). The clumps fill a fraction f_vol of the total volume, such that the mean density $$⟨\rho ⟩=f_{\mathrm{v}\mathrm{o}\mathrm{l}}{\rho }_{\mathrm{c}\mathrm{l}}+(1-f_{\mathrm{v}\mathrm{o}\mathrm{l}}){\rho }_{\mathrm{i}\mathrm{c}}.$$(1)

The density of the rarefied medium in between the clumps is set by the inter-clump density contrast, $$f_{\text{ic }}\equiv \frac{{\rho }_{\text{ic }}}{⟨\rho ⟩}.$$(2)

The clumping factor f_cl relates the mean density to the meansquare density as $$f_{\text{cl}}\equiv \frac{⟨{\rho }^2⟩}{⟨\rho {⟩}^2}=\frac{f_{\mathrm{v}\mathrm{o}\mathrm{l}}{\rho }_{\mathrm{c}\mathrm{l}}^2+(1-f_{\mathrm{v}\mathrm{o}\mathrm{l}}){\rho }_{\mathrm{i}\mathrm{c}}^2}{{[f_{\mathrm{v}\mathrm{o}\mathrm{l}}{\rho }_{\mathrm{c}\mathrm{l}}+(1-f_{\mathrm{v}\mathrm{o}\mathrm{l}}){\rho }_{\mathrm{i}\mathrm{c}}]}^2},$$(3)

such that for a void inter-clump medium (ρ_ic = 0), f_cl = 1 / f_vol or f_vol = 1 / f_cl. For the more general case of a nonvoid inter-clump medium, $$f_{\mathrm{v}\mathrm{o}\mathrm{l}}=\frac{{(1-f_{\mathrm{i}\mathrm{c}})}^2}{f_{\mathrm{c}\mathrm{l}}-2f_{\mathrm{i}\mathrm{c}}+f_{\mathrm{i}\mathrm{c}}^2}.$$(4)

All the above parameters are a function of radial distance, which we left out for simplicity of notation. The model inherently assumes the clumping to result from the line-deshadowing instability (Owocki et al. 1988). Therefore, the clumping only starts when the wind starts to accelerate significantly. Here, the onset of clumping is a free parameter (see Section 3.3), with the lowest allowed value v_cl,start = 0.01v_∞. From this starting velocity, the clumping increases linearly with velocity until it reaches its maximum value of f_cl at v_cl,max.

Each clump is assumed to have an internal velocity dispersion δv. To put this in perspective, as clumps have a physical size, it may be related to the velocity span of the clumps as a result of the underlying radially increasing smooth outflow velocity, δv_sm. If the clumps have small internal δv and are located relatively far from one another, the gas in optically thick clumps can only absorb a modest amount of light as Doppler shifts are small and do not spread out the line opacity over a wide velocity range. However, if δv is larger than δv_sm, light will be much more effectively blocked. This effect of porosity in velocity space (termed velocity-porosity or vorosity in Owocki 2008) is quantified using a normalized velocity-filling factor f_vel that takes values between 0 and 1. Finally, v_windturb describes the turbulence in the wind.

The wind-structure parameters f_ic, f_vel, v_cl,start, and v_windturb are defined as a function of the clumping factor f_cl. Therefore, if no inhomogeneities are present in the flow, that is, when the wind is truly smooth, the clumping factor f_cl equals unity. In this case, the other wind-structure parameters no longer have any effect or meaning. This implies that when the fitting procedure (see below) selects f_cl = 1 for a model calculation, none of these five detailed structure parameters affects the line profiles.

3.2 Genetic algorithm

Given the large number of parameters, it was not feasible to fully explore the parameter space with a standard grid-based approach. However, it was essential that all parameters were explored simultaneously, as many parameters are connected and correlated. Determining parameter values sequentially would then likely result in an underestimation of the uncertainties and possibly in a suboptimal final fit. Therefore, an efficient fitting algorithm was required to fit all parameters simultaneously. A genetic algorithm (GA) was found to be efficient, while remaining robust enough to not become stuck in local minima (e.g., Mokiem et al. 2005; Tramper et al. 2014). We used the genetic algorithm Kiwi-GA⁵ (Brands et al. 2022) to find the optimal fit parameters, and we refer to Brands et al. (2022) for an in-depth discussion of the method. Kiwi-GA functions both as fitting routine and as FASTWIND wrapper. The algorithm starts with a set of models, the parameters of which result from a random uniform sampling of the parameter space.

Subsequent generations of models were then generated by combining two models of the parent population. The parameters of these models were mixed, giving the new model the parameter value from either parent. Additionally, the parameters were allowed to mutate, which changed the value. The parent models were selected semi-randomly based on their χ² value, with a lower χ² giving a higher probability to be selected. The χ² value was determined by comparing the normalized model spectra to the observed spectra and their uncertainty. Only a selected set of spectral features was compared (see Section 3.3).

We determined the confidence intervals of the fit parameters using the root mean square error of approximation (RMSEA, Steiger 1998). This statistics was chosen over a regular χ² statistics because the latter requires the residuals of the data to follow a standard normal distribution. Because of the large number of high-quality (i.e., with a high signal-to-noise ratio) data points and inherently imperfect models, too many residuals diverged from a normal distribution. This resulted in high χ² values and too large Δχ² when it deviates from the best-fit model. The RMSEA corrects for this deviation from a normal distribution by rescaling the χ² distribution as follows: $$\text{RMSEA}=\sqrt{max(\frac{{\chi }^2-n_{\mathrm{d}\mathrm{o}\mathrm{f}}}{n_{\mathrm{d}\mathrm{o}\mathrm{f}}(N-1)},0)},$$(5)

with N the number of data points, and n_dof the number of degrees of freedom. The best-fitting model is given by the lowest RMSEA value, the 1σ confidence interval is given by the models that have RMSEA < 1.04 × min(RMSEA), and the 2σ intervals are given by RMSEA < 1.09 × min (RMSEA). These values of 1.04 and 1.09 were calibrated such to obtain similar confidence intervals as the χ² analysis performed in Brands et al. (2022). The RMSEA statistics does not affect the sampling performed in the genetic algorithm.

Using Kiwi-GA, we were able to fit 15 free parameters and converged to a robust solution in 80 generations. Each generation consisted of 128 models, resulting in ~10 000 models per star in total. For comparison, this is already less effort than calculating all parameter value permutations given only two values per parameter (2¹⁵ = 32 768). A grid consisting of five values per parameter, which is unrealistically coarse for most parameters, would constitute ~3 × 10¹⁰ models, and would only be computational-cost effective if more than 3 × 10⁶ stars were scrutinized to the level pursued here.

3.3 Fitting approach

The parameters we aimed to constrain are listed in Table 2. We fit essential atmosphere parameters such as the effective temperature T_eff, the surface gravity log g, and the projected equatorial rotational velocity v sin i, along with the surface abundance of helium, carbon, nitrogen, oxygen, silicon, and, when available, phosphorus. The wind (structure) parameters we fit are the mass-loss rate $\dot{M}$, the velocity law index β, the terminal velocity v_∞, the (maximum) clumping factor f_cl, the inter-clump density contrast f_ic, the velocity porosity f_vel, the onset velocity of clumping v_cl,start, and the turbulent wind velocity v_windturb. We assumed that maximum clumping was reached at v_cl,max = 2v_cl,start after increasing linearly from v_cl,start.

The analysis of each star was split into an optical-only fit to constrain the projected rotational velocity v sin i and the helium abundance y_He of the star. We did not distinguish between v sin i and additional broadening from macroturbulence. This parameterization of the broadening profile is not expected to critically affect other parameters such as the temperature and surface gravity (Holgado et al. 2018). Then, another fit addressing both optical and UV diagnostics was preformed using the v sin i and y_He value from the optical fit. This two-step approach was used to ensure that we found the correct value of the rotation. From test calculations, we found that the v sin i value was often found to be degenerate with, notably, v_windturb in wind lines, which allows slightly better fits in the UV at the cost of worse fits to many optical lines. The helium abundance was fixed in the optical and UV fits as there was little to no sensitivity to this parameter in the UV and because it would lead to more uncertain and less accurate values. Additionally, the two-step approach made it possible to check for systematic changes or biases caused by including the UV diagnostics.

In the optical-only fits, the only free wind parameter was the mass-loss rate. All other wind parameters were fixed to fiducial values of β = 1, f_cl = 10, f_ic = 0.1, f_vel = 0.5, v_cl,start = 0.05v_∞, and v_windturb = 0.1v_∞, and we used the v_∞ from ULLYSES spectroscopy (Hawcroft et al. 2024b, XShootU III).

The spectral features scrutinized in this work are listed in Table 3. The top and bottom parts of the table list UV and optical diagnostics, respectively. All lines that are part of the same complex were fit together, and the formal solution of their radiative transfer was solved together.

Wind-embedded shocks, caused by instabilities in the wind (Owocki et al. 1988, e.g.), can cause the emission of X-rays. The diagnostic features we selected to scrutinize have limited sensitivity to these X-rays. However, for completeness, we included the emission in our modeling with FASTWIND (Carneiro et al. 2016). Our implementation follows the prescription of Brands et al., in prep.

Some lines such as the OV1371 lines and N IV 1718 lines were excluded from the analysis. These did not show clear wind signatures, and the photospheric absorption lines are blended with metal lines that are not included in the model. The N V 1240 line was excluded because it strongly depends on the assumptions regarding X-rays. This results in systematic changes compared to the optical-only fits (see also Appendix D).

The mix of spectral features of different elements with varying ionization stages allowed us to accurately determine the temperature of these stars. The Balmer lines allowed us to infer the surface gravity. The large selection of wind-sensitive lines, including resonance lines and recombination lines, allowed us to constrain the mass-loss rate and terminal velocity, and if possible, also the wind-clumping parameters. We found that the latter required a relatively strong outflow (see Section 4).

Additionally, we fit the CNO-cycle element abundances and silicon and phosphorus abundances, the latter of which was only fit in the optical and UV fits. However, not many diagnostics were available for studying the Si and P abundances, with the added complexity that these diagnostics are also sensitive to wind.

The radius of the stars, and therefore, their luminosity, was determined using the absolute K_s-band magnitude as an anchor. This means that the stellar parameters determine the shape of the SED, which was then scaled with the appropriate radius to match the observed absolute magnitude. The K_s-band magnitude was chosen as it is widely available and is not strongly affected by the interstellar extinction or the thermal radiation of dust. The uncertainty on the absolute magnitude is a result of the uncertainty on the photometry, the extinction toward the star, and the distance to the star. This uncertainty affects the radius, luminosity, spectroscopic mass, ionizing flux, and mass-loss rate. The latter is affected because we assumed $\dot{M}/R_{\star }^{3/2}$ to be constant (Puls et al. 1996). However, intrinsic uncertainties on the mass-loss rate are typically significantly higher than the uncertainty on the radius. We assumed a distance of 62.44±2 kpc to the sources in the SMC (Graczyk et al. 2020; Subramanian & Subramaniam 2009). Here, the uncertainty is an estimate of the depth of the SMC, as it is unclear where the source is in the direction of the line of sight through the galaxy.

3.4 Comparison with evolutionary models

We determined the initial mass, current evolutionary mass, and age using BONNSAI⁶ (Schneider et al. 2014). BONNSAI is a Bayesian framework that can compare observed stellar properties with those of evolutionary models to obtain the posterior distribution of additional stellar properties. As input, we used the luminosity, temperature, and an upper limit of v sin i in combination with the SMC evolutionary tracks of Brott et al. (2011). We used an upper limit for the rotation as we did not include macroturbulence in our modeling (e.g., Simón-Díaz et al. 2017). Therefore, the rotational velocity could be lower than what we found.

The inferred temperature and luminosity of AV 15, AV 70, and AV 80 are not covered by the SMC model grid of Brott et al. (2011), and we therefore used the LMC grids of Brott et al. (2011) and Köhler et al. (2015) instead. We used the default BONNSAI settings with the exception of the prior for the initial rotation velocity, for which we used the distribution of Ramírez-Agudelo et al. (2013).

4 Results

Tables 4 and 5 show the best-fit stellar atmosphere parameters and derived parameters along with their uncertainty. For completeness, we include the spectroscopic mass M_spec, the number of H I and He I ionizing photons (Q₀ and Q₁), and the Eddington parameter for electron scattering Γ_Edd,e in Table 5. We do not discuss these parameters further. The helium surface abundance y_He and the projected rotational velocity v sin i are the result of the optical-only fit. All other parameters are from the optical and UV fit. The remainder of the best-fit parameters of the opticalonly fits, which are not the focus of this paper, can be found in Appendix E. Table 6 shows the best-fit mass-loss and windstructure parameters and their uncertainties. Figs. H.1 and H.2 show an overview of the diagnostic lines and line complexes and their fits. More detailed figures of the line profiles and fit results can be found online⁷.

Below, we describe the general fits for the UV and wind lines in Sect. 4.1, and we describe the fit for the optical lines in Sect. 4.2. In Sect. 4.3 we comment on some poorer fits.

4.1 UV and wind lines

Overall, the UV lines are well reproduced by the models, with the exception of the phosphorus lines. The optical wind-sensitive lines are generally reproduced reasonably well. The fits are described per line group below.

Phosphorus. The P V lines at 1118 and 1128 Å are in absorption in all of the sources. The two components of the doublet are very sensitive to the mass-loss properties and temperature of the star. For AV 15, AV 95, AV 207, and AV 307, the best-fit line profiles match the observations very well. For the other stars, the best-fitting results show a stronger in-filling of the line core or even P-Cygni profiles, while photospheric absorption profiles were observed.

Carbon. The C IV 1550 doublet resonance lines are well reproduced by the models, with clear P-Cyngi wind profiles for most stars. The two components of the C IV 1169 and six components of the C III 1176 lines are generally photospheric, with the C IV lines typically being relatively weak. The C III line complex only shows clear wind signatures for AV 83 and AV 372. The photospheric profiles are typically well reproduced, with some exceptions where the model shows some emission from the wind. In these cases, notably AV 479, AV 70, and to a limited extent, AV 469, the produced absorption feature is not strong enough.

Oxygen. Only O IV 1340 was modeled in the UV. The O V 1371 line is either not present or too contaminated by blended metal lines to be included in our analysis. The O IV 1340 line, consisting of three transitions, is dominated by the photospheric component for all stars in the sample, with only some blueshifted absorption for some of the stars. Typically, the lines are well reproduced, but in some cases, the best-fit model results in a too weak absorption profile, particularly for the blue component.

Silicon. The Si III 1113 line is photospheric for all stars, and in most cases, it is well reproduced by the model. The Si IV 1400 doublets typically show both photospheric and wind features. The line profile resulting from the wind is well reproduced. However, the photospheric component is often underestimated by the model. This could be due to contamination by the superimposed interstellar absorption feature, which is difficult to distinguish from the stellar feature and was therefore not clipped from the data.

Helium. The He II 1640 line is blended with iron lines on the blue side of the profile, and therefore, only the red wing can be properly fit with FASTWIND. Only AV 80 and AV 83 show P-Cygni profiles for this line, and these are reasonably well reproduced. The stars showing a photospheric profile are also well reproduced.

Nitrogen. When available, the N III 1750 lines are photospheric in nature. They are typically well reproduced. For some of the cooler stars, we find slightly weaker lines in the models than in the observations.

Optical emission lines. The main optical wind diagnostics are He II 4686 and Hα. Whereas the main UV wind lines are resonances lines, these are recombination lines. The combination of the two types of lines is essential for constraining the clumping parameters of the wind. The He II 4686 line has proven hard to reproduce in this work. The exact strength of the line is often not recovered, with a slightly too strong central absorption component. The most extreme case of this is AV 479, where we find a significantly stronger absorption component than is observed. 2dFS 163 shows only a very weak and broad absorption, while the best-fit model features a significant emission profile (see also Section 4.3). The Hα lines are typically reproduced quite well within the uncertainties. We note that AV 372 and AV 70 show complex line shapes for the hydrogen and helium lines that are not fully recovered by the model, but the strength of the features are recovered within the uncertainty. The shape of these complex profiles can be reproduced by a 1D FASTWIND model, but they remain very sensitive to many parameters, as can also be seen in the large uncertainty region. The complex of C III and N III lines from 4630 to 4660 Å is well reproduced. These lines are very sensitive to various other stellar parameters, such as temperature and gravity, and therefore, the uncertainties on the model profiles are substantial, as indicated by the shaded region.

4.2 Optical absorption lines

Below, we describe the overall fit of several groups of absorption lines in the optical spectra. Their widths are generally well reproduced, suggesting that the v sin i was properly determined using the optical-only fit.

Hydrogen. The hydrogen lines Hβ through Hε are accurately reproduced by the models, including the wings, which are sensitive to the surface gravity of the stars. Exceptions apply for stars with poor fits (see Section 4.3). For some stars, the wind emission was strong enough to (partially) fill in the absorption profile. This was accurately reproduced by the models.

Helium. The helium lines, which are good indicators of the accuracy of the temperature determination, fit the observations well, with the exception of 2dFS 163 and AV 327. The observations of both stars show stronger He II lines than the best-fit models, suggesting that the temperature is too low. However, the remaining diagnostics prevent a higher temperature from resulting in a better overall fit.

Nitrogen, carbon, oxygen, and silicon. The optical nitrogen and carbon lines are typically well reproduced. The same holds for oxygen, where the O III 5592 line is consistent with the data, and the weak O III and O II lines in blends with other lines are reproduced within the uncertainties. These uncertainties are relatively large, however, which we tentatively ascribe to the strong temperature sensitivity of these lines. These considerations, combined with the fact that only a few lines could be used in the analysis, result in rather large uncertainties on the oxygen abundance.

The silicon lines, located in the wings of Hδ, are reproduced with varying success. The lines are very sensitive to temperature, gravity, and mass loss. As a result the lines may appear in as emission or absorption profiles depending on these parameter values. The observed strength of the line therefore varies significantly. This resulted in rather poor fits for AV 83, AV 207, and AV 479.

4.3 Anomalous fits

FASTWIND successfully reproduced the observed spectra for most stars, but for a few stars, some spectral features could not be fully matched with the model. For three stars (AV 70, AV 83, and 2 dFS 163), the mismatch between model and observation was too significant for us to consider the inferred stellar parameters to be representative of the source. These stars were therefore left out in the further analysis, are colored gray in the tables, and have a red outline in the figures.

2dFS 163. The observed N IV 3480 lines of 2dFS 163 are significantly stronger than those calculated by FASTWIND. Additionally, the strength of the He II lines is systematically underpredicted by the model, except for He II 4686. This wind-sensitive line is almost absent in the observed spectrum, but the model predicts a strong emission. The emission is the result of the high clumping factor and mass loss required to reproduce the other wind-sensitive lines. This likely also prevented the temperature from increasing further to allow the strength of the He II lines to match the observations, as a higher temperature would result in even stronger emission in He II 4686. Furthermore, the hydrogen and He I lines appear to be slightly blueshifted, while this does not seem to be the case for the other lines. The He I lines also appear to be slightly narrower than the other lines, suggesting a different rotation rate. The Hα line (see Fig. H.1 and online⁸) is double peaked, as are the Paschen lines, which we do not display here. Taken together, we suspect that 2dFS 163 is a composite source. Ramachandran et al. (2024, XShootU VIII) identified this star as a post-interaction binary system. Therefore, this source was left out in the further analysis, is marked with a red outline in the figures, and is listed in gray in the tables.

AV 83. The spectrum of AV 83 is likely composite. It shows a feature in the red wing of the He II 5411 line that is also visible in the spectrum used by Hillier et al. (2003), but shifted in wavelength. This results in a broad line. The optical He II, C IV, and N IV lines all appear to be slightly redshifted compared to the other lines. Some He I lines and the optical Si IV lines show emission that is not reproduced by the model. As a result, the best-fit parameters of this star are likely not representative of its properties. This source was left out in the further analysis, is marked with a red outline in the figures, and is listed in gray in the tables.

AV 70. The He I lines in the spectrum of AV 70 show a broad component that is not visible in other lines, nor is it reproduced by the model. The broad component could suggest the presence of a rapidly rotating cooler companion. This is backed by the apparent flat-bottomed C IV 1550 absorption trough of the P-Cygni profile, for which the flux remains at 30% of the continuum at the bottom. For this reason, the star was left out in the further analysis, is marked red in the figures, and is listed in gray in the tables.

AV 307. Additional absorption lines can be seen in the C III-N  III 4630-4655 complex of AV 307. These additional lines coincide with the transitions of O II, but they are not observed in any of the other stars in the sample, with the exception of AV 327. In this star, the O II lines are considerably weaker, however. A closer inspection shows the spectrum of AV 307 to be rich in O II lines. AV 307 and AV 327 are the coolest stars in the sample, both with T_eff =29 500 K as the best-fit value, and they feature the lowest v sin i, with 55 and 80 km s⁻¹. The abundance of oxygen is found to be higher in AV 327, mainly based on the strength of O IV 1340.

AV 479. The He I and He II lines of AV 479 fit well within the uncertainties, with the exception of He II 4686, for which the best-fit model displays a significantly stronger absorption than the observations. The Hα line shows a slightly too strong absorption feature in the model. A consistent spectral fit for both the UV and optical wind signatures therefore does not appear to be feasible: The Si IV 1400 and C IV 1550 lines are fit well, at the expense of Hα and He II 4686. Therefore, a higher mass-loss rate would be at the expense of the fit quality of the two ultraviolet lines. Similarly, the fit to the optical wind lines could benefit from a higher clumping factor. However, this would affect the ionization structure of Si and C, resulting in a worse fit for the Si IV 1400 and C III 1176 lines. Alternatively, an earlier onset of clumping, that is, a lower value for v_cl,start, could also improve the line fits of Hα and He II 4686. Again, this is disfavored for the ultraviolet wind lines.

5 Discussion

The main part of this discussion focuses on the mass-loss properties of the SMC sample in contrast to comparable samples in the Milky Way and LMC, such as to probe the dependence of mass loss on metallicity. However, we start with a brief discussion of other derived properties to characterize the target stars and their evolutionary state better.

Many of the stars in our sample have been the subject of other studies. We compare our findings to those studies in Appendix A. Most parameters are consistent with previous findings, although we find slightly higher temperatures and luminosities.

5.1 Evolutionary state

Figure 1 presents our sample of giants, bright giants, and supergiants in the Hertzsprung-Russell diagram (HRD). The figure shows that the sample covers a substantial part of the (T_eff, L) plane. All objects have moved away from the zero-age main sequence, but are still on the main sequence when considering single-star evolutionary tracks (Brott et al. 2011), and their ages range from approximately 2–6 Myr. For three stars (AV 83, 2dFS 163, and AV 70; see also Section 4.3), we were not able to determine reliable parameter values. While we present their bestfit values, these sources were not considered in the mass-loss rate and clumping analysis that is the topic of the next subsections. In Fig. 1 and also in Figs. 3–5 below, these sources are present and marked with a red border, but they are not included in any fits or analyses.

Figure 2 shows the derived surface mass fractions for He, C, N, O, Si, and P. We also plot the SMC baseline mass fraction following Dopita et al. (2019) for all elements save phosphorus, for which we scaled the Asplund et al. (2009) solar abundance by a factor of 1/5. We also show indications of the expected surface mass fractions, using an initially 40 M_⊙ SMC model that rotates at the end of formation at 389 km s⁻¹ and has a typical age of 3 Myr (Brott et al. 2011). For the rotational velocity we chose, the range of helium mass fraction is well reproduced, but we remark that for a lower initial spin, rotation-induced mixing is less efficient. X_He is therefore only marginally enriched (X_He < 0.3). The range of nitrogen mass fractions is well reproduced. The C and O abundances are slightly above the baseline, albeit with large uncertainties, while depletion is expected. In Appendix B we compare the total combined C, N, and O surface abundances to the baseline abundance and to the study of Martins et al. (2024; XShootU V). Their sample overlaps with ours. They analyzed the same data with a different method that is more focused on surface abundances.

Mixing is not expected to affect silicon and phosphorus. The surface abundances of these elements are poorly constrained; for silicon, they are typically above the baseline, and for phosphorous, they cluster around the baseline.

The latter may be an effect of the sampling. Because the effect of the surface abundance on the strength of the P V doublet is limited, the sampling is close to uniform. For example, AV 307 has an excellent fit to the P V lines, but the phosphorus abundance is unconstrained due to the dependence on temperature, among other parameters. This correlation might be too strong to determine the abundance reliably with only one available ionization stage.

We repeat that 2dFS 163, AV 83, and AV 70 might be binaries (see Sect. 4.3) and that their best-fit parameter values are likely not representative of their true values. They were left out in the further analysis.

Fig. 1 Hertzsprung-Russell diagram of our SMC sample of stars with temperatures and luminosities resulting from the fit. Points with a red border indicate stars with unreliable parameters. We overplot the evolutionary tracks (thin solid gray lines) and isochrones (thin dashed gray lines) of Brott et al. (2011). The solid black line indicates the zero-age main sequence. The points are labeled with the name of the object in a small font that is visible when zoomed in.

5.2 Comparison of empirical and predicted mass-loss rates

We briefly show how the mass-loss rates in our sample compare to theoretical predictions for SMC stars. Figure 3 shows the mass-loss rate as a function of luminosity, compared to the mass-loss rate predictions of Vink et al. (2001); Krtička & Kubát (2018); Björklund et al. (2021); for the SMC metallicity, we assumed 20% of solar for all predictions, following Mokiem et al. (2007). The dotted line indicates a linear fit to the empirical mass-loss rates using orthogonal distance regression (ODR Boggs et al. 1987), which takes the uncertainties of all variables into account, that is here, mass-loss rate and luminosity, given by $$\log \dot{M}=(-20.81\pm 1.48)+(2.48\pm 0.26)\times \log \left(L/L_{\odot }\right).$$(6)

The luminosity dependence that we find is stronger than that of the theoretical predictions that are shown in the figure. Overall, the Björklund et al. (2021) mass-loss rate-luminosity relation agrees best with the empirical mass-loss rate from our sample. Vink et al. (2001), Gormaz-Matamala et al. (2022b), and Krtička & Kubát (2018) predicted generally higher mass-loss rates at low luminosity, with the discrepancy decreasing toward higher luminosity. For the comparison to the prediction of Vink et al. (2001), we used the best-fit stellar properties and the spectroscopic mass. Slightly different results can be expected when the evolutionary mass as a discrepancy between the two masses remains (see Appendix C).

Appendix A.1 compares the mass-loss rates found here to those found in previous optical and UV stellar wind analyses in the SMC. We mostly find mass-loss rates consistent with those of previous studies.

For two stars, AV 372 and AV 469, a significantly higher mass-loss rate is found than for the rest of the sample given their luminosity. It is unclear why these stars require such a mass-loss rate for an optimal fit. Both stars agree within the uncertainties with the Vink et al. (2001) and Gormaz-Matamala et al. (2022b) predictions, while the rest of the sample agrees with those of Björklund et al. (2021). Both stars are classified as supergiants, but other supergiants in the sample follow the same trend as the rest of the sample. The He abundance for both stars is significantly enriched from the baseline, showing the highest abundances in the sample. Additionally, we find a strong nitrogen enrichment in AV 469. The diverging values are still within the typical scatter found when determining mass-loss rates using optical and UV spectroscopy (see e.g., Brands et al. 2022; Bouret et al. 2021). It remains possible, however, that these stars are physically different from the rest of the sample in ways in that are currently unknown, for instance, they may be post-interaction binaries. Alternatively, these stars could be in a transient episode of higher mass loss.

Fig. 2 Distribution of the surface mass fractions of He, C, N, O, Si, and P from the FASTWIND/GA fitting of SMC stars. The solid green line indicates an approximate kernel density estimate of the distribution taking the uncertainties on the best-fit values into account. The thin dashed lines show the contribution of the individual sources. The vertical black line shows the baseline mass fraction from Dopita et al. (2019) for all elements except phosphorus, for which we show the scaled (by 0.2) solar abundance from Asplund et al. (2009). The area marked in red shows the surface depletion or enrichment that can be expected due to mixing based on the evolution of a 40 M⊙ star with an initial rotation of 389 km s−1 (Brott et al. 2011). For phosphorus, only stars with a reasonable fit are included. The helium mass fractions are taken from the optical-only fit. All other fractions are from the optical and UV fit.

Fig. 3 Mass-loss rate as a function of luminosity. The fit to the empirical mass-loss rates and its 1σ uncertainty are indicated with the dotted green line and the shaded region, respectively. The mass-loss rate – luminosity relations of Krtička & Kubát (2018) and Björklund et al. (2021) are shown with the dashed blue and yellow lines, respectively. For the mass-loss predictions of Vink et al. (2001), we used their equation (24) to obtain the predicted rates for each star based on their properties. The magenta V indicate these rates, and the dashed magenta line corresponds to the fit to these points. The light blue stars indicate the mass-loss rates of equation 7 from Gormaz-Matamala et al. (2022b). Points with a red border are unreliable fits that are likely not representative of their physical properties; these are not included in the fit.

5.3 Dependence of the modified wind momentum on luminosity and metallicity

To empirically test the relation between mass loss, luminosity, and metallicity, we compared the mass-loss rates of samples of Galactic (Hawcroft et al. 2021), LMC⁹ (Hawcroft et al. 2024a, Brands et al. in prep.), and SMC (this work) O stars. In these studies, the winds were scrutinized using a FASTWIND GA-analysis of both their optical and UV spectrum in which the inhomogeneous wind structure accounts for optically thick clumping. The uniform nature of these analyses – they all applied the same formalism for the wind inhomogeneities – allows for the best possible comparison of wind mass loss in different metallicity environments. However, it is important to note that the despite the similar approaches, the properties of the samples are not identical. Hawcroft et al. (2021) focused on a set of Galactic supergiants with only one star with a spectral type later than O7 and all stars more luminous than 10^5.6 L_⊙. The LMC stars are also on average of an earlier type than the SMC sample, with 14 stars having a spectral type earlier that O7, while only 2 such stars appear in the SMC sample.

To isolate and quantify the impact of metallicity on mass loss from the limited size samples available, we resorted to the mechanical momentum of the stellar wind, modified by a term depending on stellar radius. This property, the so-called modified wind momentum, $$D=\dot{M}{v}_{\mathrm{\infty }}\sqrt{\frac{R_{\star }}{R_{\odot }}}\left[\mathrm{g}\mathrm{r}\mathrm{c}\mathrm{m}{\mathrm{s}}^{-2}\right],$$(7)

is primarily a function of stellar luminosity and often predicted to be almost independent of stellar mass (Kudritzki et al. 1995; Kudritzki & Puls 2000). We note that for this to be valid, however, the CAK-α parameter (see below) needs to be constant at a value of 2/3. This may not be the case throughout our parameter space or everywhere in the wind (e.g., Muijres et al. 2012). To compare our sample with other observations and theoretical predictions, we fit the following relation through the observed values: $$\log D=\log D_0+x\log L/L_{\odot }.$$(8)

From now on, we refer to x as the slope, and to log D₀ as the offset of the modified wind momentum-luminosity relation. Figure 4 shows the relation for the three samples, using ODR to constrain the fit variables x and log D₀. The empirical relations are distinguished by color: purple for Galactic stars, orange for LMC stars, and green for SMC stars. Although the uncertainties for the Galactic sample are sizeable, the overall trend is indeed that the wind momentum increases with metallicity. Furthermore, we find a tentative trend in the slope of the modified wind momentum-luminosity relation, with steeper slopes toward lower metallicities. Our fit values are presented in Table 7.

To account for a metallicity dependence of the modified wind momentum, we fit all three samples simultaneously with the function $$\log D(L,Z)=(a+b\log \frac{Z}{Z_{\odot }})\log \frac{L}{{10}^6{L}_{\odot }}+c\log \frac{Z}{Z_{\odot }}+d,$$(9)

following Krtička & Kubát (2018) and Björklund et al. (2021), where a through d are the fitting parameters. We assumed Z_MW = Z_⊙, Z_LMC = 0.5 Z_⊙, and Z_SMC = 0.2 Z_⊙, consistent with Mokiem et al. (2007), Krtička & Kubát (2018), Björklund et al. (2021), and Marcolino et al. (2022). We took the uncertainty in the wind momentum, luminosity, and metallicity into account. For the latter, we assumed the uncertainty to be Δ log Z/Z_⊙ = 0.1, which is conservative considering the typical uncertainties on the measured abundances (e.g. Dopita et al. 2019). Figure 5 shows the best fit of Equation (9) to the data evaluated at Z = 1, 0.5 and 0.2 Z_⊙. The fit shows that over the wide range in luminosity probed by the samples, the modified wind momentum is only poorly represented by a luminosity-independent power-law dependence on metallicity. For relatively dim stars, the metallicity dependence is stronger than for relatively bright stars. Although the fit agrees nicely with the LMC and SMC data points, it slightly overestimates the wind momentum of the MW stars. The best-fit values of the parameters a, b, c, and d in Eq. (9) are listed in Table 8. These parameters are strongly correlated, and therefore, their individual variances are a poor representation of their true uncertainty. The full covariance matrix of the fit can be found in Appendix G.

Comparing the two fitting procedures, we note that the gradient in slope between the direct linear fits to the individual samples is smaller than in the combined fit, as can be seen by comparing Figures 4 and 5. This is a result of the difference in wind momentum between the Milky Way and the LMC, which is smaller than the same difference between the LMC and SMC stars, while in both cases, the difference in metallicity is similar. As a result, the metallicity-offset term, c in Equation (9), cannot accommodate both differences at once. To mitigate this, the metallicity slope, b in Equation (9), was adjusted. This works because of the lack of relatively low-luminosity points in the Milky Way sample. Therefore, caution is advised: our fit result for Equation (9) may overestimate the slope in the metallicity dependence. Additional studies using similar methods to model the clumped winds of, in particular, Galactic stars with luminosities log L/L_⊙ ≤ 5.6 are needed to scrutinize this further. Alternatively, the discrepancy between the two fitting results may signal that the parameterization in Equation (9) is not sufficiently suitable. This might be mitigated by adding a nonlinear (in log-space) metallicity term.

Fig. 4 Modified wind momentum as function of luminosity for samples of SMC (this work), LMC (Hawcroft et al. 2024a, Brands et al., in prep.), and Galactic stars (Hawcroft et al. 2021). The dotted lines indicate the linear fits to each subsample, and the shaded regions show their 1σ uncertainties. The points with a red border represent unreliable fits that are likely not representative of their physical properties; these are not included in the fit.

Fig. 5 Modified wind momentum as function of luminosity and metallicity, as in Figure 4, but now fit with Equation (9) to the three samples, instead of separate linear fits for each sample. The dashed magenta lines and pink shaded regions indicate the best fit and the 1σ confidence interval, respectively. Equation (9) has been evaluated at Z = 1.0, 0.5, and 0.2 Z⊙. The points with a red border are unreliable fits that are likely not representative of their physical properties; these are not included in the fit.

5.3.1 Comparison to other empirical studies

Figure 6 repeats the fitting results from Fig. 4 along with the best relations from Mokiem et al. (2007) and Marcolino et al. (2022), who also empirically investigated the metallicity dependence of mass loss. Mokiem et al. (2007) used data from optical studies of Galactic and Magellanic Clouds stars. Mass-loss rates follow from FASTWIND analyses without considering clumping effects. Because ultraviolet diagnostics are not considered, quite a few of their relatively low-luminosity sources yield only upper limits for $\dot{M}$. Marcolino et al. (2022) collected results from CMFGEN (Hillier & Miller 1998) studies that relied on both ultraviolet and optical spectra. In these MW and SMC studies, an optically thin clumping prescription was used that sometimes only yielded rough estimates of the clumping factor.

Mokiem et al. (2007) estimated the metallicity dependence at only a single luminosity, that is, log L/L_⊙ = 5.75, a typical mean for their sample. Near this luminosity, their slope is not too different from ours. Their wind momentum is higher than ours, however, which likely is a result of Mokiem et al. not considering clumping. Marcolino et al. (2022) reported a significantly steeper slope in the D(L) diagram for both the MW and SMC. Interestingly, their MW and SMC curves diverge for higher luminosities, while our curves converge. The luminosity range considered in Marcolino et al. (2022) extends to ~10^4.6 L_⊙, which is significantly lower than our sample. This might lead to a different slope.

We quantitatively compared our results to the pioneering work of Mokiem et al. (2007) by evaluating Equation (8) at log L/L_⊙ = 5.75 with the SMC, LMC, and MW coefficients as listed in Table 7. Under the assumptions that the mass-loss rate and terminal wind velocity scale as $$\dot{M}\propto Z^m$$(10)

and $$v_{\mathrm{\infty }}\propto Z^n,$$(11)

such that $$m+n=\frac{\mathrm{d}\log D_{\mathrm{m}\mathrm{o}\mathrm{m}}}{\mathrm{d}\log Z},$$(12)

with n = 0.13 (Leitherer et al. 1992), we find from an ODR fitting of the three points on the D(Z, L = 10^5.75 L_⊙)-plane $$m=\frac{\mathrm{d}\log \dot{M}}{\mathrm{d}\log Z}=1.02\pm \mathrm{0.30.}$$(13)

This is slightly higher than but consistent with the value of Mokiem et al. (m = 0.83 ± 0.16).

In the top panel of Figure 7, we repeat this procedure for all luminosities (dashed blue line and light blue uncertainty region). At log L/L_⊙ = 5.75, m = 1.02, as discussed above. However, we do not find a constant value for m, but higher m for lower luminosities. The luminosity-independent result m = 0.83 of Mokiem et al. is plotted as the red point. We plot with a dotted green line the results of Marcolino et al. (2022). As their metallicity dependence only relies on MW and SMC results, they did not require a three-point fitting in order to obtain m. At the reference luminosity log L/L_⊙ = 5.75, they found m ~ 0.6, which is lower than Mokiem et al. (2007) and just outside our 1σ confidence interval. However, Marcolino et al. (2022) reported a positive trend of m with luminosity, which is opposite to our findings and not in line with theoretical predictions and expectations (see Section 5.3.2). It is important to note, however, that the sample of Marcolino et al. (2022) extends to luminosities of ~10^4.6 L_⊙, where the winds are significantly weaker for both Galactic and SMC stars. Ramachandran et al. (2019) and Rickard et al. (2022) also studied the atmosphere properties and mass-loss rates of massive stars in the SMC. Their samples consisted predominantly of dwarf stars and mostly covered lower luminosities. Ramachandran et al. (2019) reported a value of ~2, as shown in Fig. 7 with the dash-dotted dark green line, which is consistent with our findings at lower luminosities. Rickard et al. (2022) found a similar trend with luminosity, that is, higher values of m at lower luminosity, with the metallicity dependence scaling as m(L) = −1.2 log(L/L_⊙) + 10.3. Their value of m would range from 4.3 to 3 in our luminosity span. This significantly higher metallicity dependence might be caused by the so-called weak-wind problem in SMC dwarf stars.

Fig. 6 Fits to the empirically determined modified wind momentum as a function of luminosity. The solid colored lines show the fits to the MW, LMC, and SMC samples presented here. The dashed and dash-dotted black and gray lines indicate the modified wind momentum relations from Mokiem et al. (2007) and Marcolino et al. (2022), respectively. The values representing the fits are listed in Table 7.

Fig. 7 Slope of the mass-loss rate as a function of metallicity, m, for different values of the luminosity. The top panel shows the empirical determinations of Mokiem et al. (2007), Ramachandran et al. (2019), and Marcolino et al. (2022), along with our determination based on the linear fits as shown in Figs. 4 and 6. The empirical values of m from Rickard et al. (2022) lie above the vertical range of this figure. The bottom panel shows the predicted values of m of Krtička & Kubát (2018) and Björklund et al. (2021), along with our empirical determination. The value of m = 0.83 from Mokiem et al. (2007) is also shown for reference. In all cases, d log v∞/d log Z = n = 0.13 is assumed (Leitherer et al. 1992).

5.3.2 Comparison to predictions

Table 8 lists our fit parameters a, b, c, and d of Equation 9, along with the predicted values from Krtička & Kubát (2018) and Björklund et al. (2021). Additionally, these relations are shown graphically in Appendix G. In order to obtain the values of Krtička & Kubát, we calculated the modified wind momentum rates using the values from Tables 1 from both Krtička & Kubát (2017) and Krtička & Kubát (2018) up to 42 500 K. We find a significantly stronger metallicity dependence on the modified wind momentum-luminosity slope (parameter b) than both of the two predictions. Consequently, this implies a stronger metallicity dependence of the mass loss at low luminosity.

The bottom panel of Figure 7 again shows our empirical value of m as a function of luminosity and the result from Mokiem et al. (2005) of m = 0.83, but now compared to the two sets of theoretical results. Both predictions show a decreasing trend of m with luminosity as in our empirical findings, albeit much weaker, and opposite to that of the work of Marcolino et al. (2022). We note that especially the Björklund et al. (2021) result is close to the single m value of Mokiem et al..

In passing, we note that we performed the same analysis for the derived mass-loss rates directly, instead of using the modified wind momentum. This can be found in Appendix F. The results from that analysis are consistent with the wind momentum results presented here. This suggests that the value of d log v_∞/d log Z = n = 0.13 from Leitherer et al. (1992) that we adopted is reasonable, and a more recent determination of n = 0.19 by Vink & Sander (2021) would also be consistent with our findings. A weak dependence of the terminal wind speed on metallicity was also reported by Hawcroft et al. (2024b, Paper III).

5.3.3 Line-driven wind theory considerations

The question arises whether the stronger dependence of the mass loss on metallicity for lower-luminosity stars can be explained by line-driven wind theory. The original theory by Castor et al. (1975, henceforth CAK), and subsequent CAK-like theories such as Puls et al. (2000); Gayley (1995), and re-formulations of CAK-like theories, such as Lamers & Cassinelli (1999), essentially have a principal scaling of $$\dot{M}\propto {\overline{Q}}^{(1-\alpha )/\alpha }\cdot L^{1/\alpha }\cdot [M(1-\mathrm{\Gamma }){]}^{-(1-\alpha )/\alpha },$$(14)

where M is the stellar mass, Γ is the Eddington parameter for Thomson scattering, $\overline{Q}$ represents the maximum value of the multiplication factor by which the total continuum opacity in the outflow is to be multiplied to also account for line opacities (Gayley 1995), and α may be interpreted as the ratio of the line force from optically thick lines to the total force (Puls et al. 2000); hence 0 < α < 1. For stars that lie not too close to the Eddington limit, it thus holds that $\dot{M}$ ∝ L^1/α. Assuming no direct dependence of terminal velocity on luminosity, we arrive at $$D\propto L^{1/\alpha }=L^x,$$(15)

where the last equality is what we used in Eq. (8).

The results in Table 7 would thus be matched by a decreasing value of α with metallicity: α = 0.72, 0.53, 0.41 for Z = 1, 0.5, 0.2 Z_⊙, respectively. Puls et al. (2000) explored the metallicity (and effective temperature and radial depth) dependence of α in OBA-star outflows in the regime 0.1–3 Z_⊙ (their figure 27, where $\hat{\alpha }\sim \alpha$). A CAK-optical depth parameter t ~ 10⁻⁵ indeed shows a roughly similar behavior to our findings. However, for higher values of optical depth (i.e., closer to the stellar surface), this dependence disappears. The overall effect may therefore be modest (such as found by Krtička & Kubát 2018; Björklund et al. 2021), at least in the metallicity regime probed here. Still, further theoretical studies, probing a wider metallicity range, would be extremely helpful in assessing the reasons for the luminosity dependence of the modified wind momentum versus metallicity dependence reported here and by Marcolino et al. (2022). These studies may also address the weak empirical dependence of D(Z) at high luminosities found here. This behavior might indicate that for very strong winds, the line force is dominated by optically thick lines, such that a Z-dependence is lost. Alternatively, at high luminosity, the Eddington limit may be exceeded already quite deep down in the photosphere, initiating a mass loss that the star is simply not able to cope with, leading to fall-backs and a maximum sustainable mass loss.

In addition to theoretical studies, empirical studies should be pursued to improve the coverage of the parameter space (especially in the low-luminosity regime of Galactic stars). Alternatives to the parameterization Eq. (9) should be explored.

5.4 Wind structure and metallicity

The dependence of the mass-loss properties on metallicity as described above is based on analyses that accounted for the presence of (optically thick and thin) clumps in the outflow as well as for the associated vorosity. The derived wind structure parameters for our target SMC stars are listed in Table 6. They show a wide range of clumping factors f_cl, of 2–30, and a density of the inter-clump medium that is relatively low, with f_ic ~ 10⁻² − 10⁻¹. The velocity-porosity effect, characterized by f_vel, shows a wide range of values, from 0.04 to 0.74. The wind begins to form inhomogeneities close to the onset of the wind at v_cl,start ~ 0.01 v_∞ to 0.23 v_∞. For the turbulent velocity, our analysis only probes the range of values from 0–0.3 v_∞, following earlier empirical results by Groenewegen et al. (1989). We find that the results fully span this allowed range, with v_windurb ~ 0.01–0.29 v_∞.

The obtained clumping constraints allowed us to probe for the first time for O-type stars whether the wind-clumping properties themselves are metallicity dependent and thus contribute to the empirical behavior of D(L, Z). If this dependence were strong, an in-depth treatment of clumping would be a prerogative for any empirical study of the dependence of mass loss on metal content. We continue here with caution as these first results should be considered as only exploratory for three important reasons. First, a comparison of clumping properties at different metallicities requires a similar treatment of the wind structure to within great detail. For this reason, we opted to only discuss the analyses of the sample of SMC stars described here and of LMC stars described in Brands et al. in prep., where very similar methods are applied to very similar data sets. We thus excluded the study by Hawcroft et al. (2021), who also accounted for optically thick clumping, but where slightly different assumptions regarding its treatment were adopted. Second, the clumping prescription adopted here still lacks a quantitative comparison to 2D or 3D hydro-simulations of radiation-driven outflow, in which wind structure naturally develops (e.g., Sundqvist et al. 2018). This comparison would help us to assess its appropriateness. Third, the O-star SMC and LMC samples scrutinized here are still relatively small and show some differences in intrinsic properties, with generally hotter and brighter stars in the LMC sample. The latter therefore implies overall somewhat more massive stars, with stronger stellar winds.

Figure 8 shows the distribution of the five wind-structure parameters (f_cl, f_ic, f_vel, v_cl,start, v_windturb) and β, describing the acceleration rate of the smooth wind flow at LMC and SMC metallicities. The histograms also show the approximate kernel density estimation (KDE) based on the (often sizable) uncertainty on the parameter values. For the clumping factor, the square root of its value is displayed, as for optically thin clumping and a void inter-clump medium the optical wind-diagnostic lines are essentially invariant for the property $\dot{M}\cdot \sqrt{f_{\mathrm{c}\mathrm{l}}}$.

To assess the significance of the differences between the two distributions, we performed a two-sample Kuiper test. We found that the differences are consistent with resulting from a random sampling for f_vel, v_cl,start, v_windturb and β. For f_cl and f_ic, we found that the probabilities of the difference result purely from random sampling of 0.08 and 0.11, suggesting that the underlying distributions may differ. We note, however, that the Kuiper test does not take the uncertainty on the parameter values into account, which may make the difference less significant.

f_cl and f_ic appear to show higher values at higher metallicity. At face value, this would suggest that the inter-clump medium in the LMC is relatively dense (i.e., high f_ic) and that the clumps have a high density (i.e., high f_cl), but only fill a relatively small fraction of the wind volume. In the wind medium of the SMC stars, the inter-clump medium is relatively rarefied, the clumps have a modest density, but fill a relatively large fraction of the wind volume. Driessen et al. (2022) reported a similar trend in the clumping factor with metallicity. They reported that the clumping factor scales with metallicity as f_cl ∝ Z^0.15. This would imply an increase in f_cl of ~15%, however, we find a difference of more than a factor of 5. Whether our analysis indeed implies these differences between the properties of the wind structure in the LMC and SMC sample requires additional studies (of larger samples). The differences might be associated with correlations between the clumping properties and other stellar parameters (e.g., mass loss). However, a search for these dependences did not yield any convincing correlation. As the clumping properties impact the strength and shape of spectral lines in complex ways, which differ from line to line and vary as a function of stellar properties, we cannot exclude possible degeneracies between the parameters in some parts of the parameter space. We conclude that the wind properties probed here do not show convincing evidence for a dependence on metallicity.

Fig. 8 Distributions of the wind parameters of LMC and SMC XShootU samples (this work and Brands et al. in prep.). The solid lines indicate the approximate KDE based on the uncertainty on the parameter values. The thin dashed lines indicate the contribution of the individual objects to the total KDE. Objects for which the parameters are unconstrained are excluded.

6 Conclusion

We have analyzed ultraviolet ULLYSES and optical X-Shooting ULLYSES spectra of 13 massive O4–09.5 stars in the SMC using the model atmosphere code FASTWIND and the genetic algorithm-fitting approach KIWI-GA. The KIWI-GA approach is the only way in which as much as 15 stellar and wind properties, of which 8 characterize the outflow, may be constrained simultaneously. The study accounted in detail for clumping, porosity, and velocity-porosity properties of the wind. These processes may impact the mass-loss-sensitive spectral-line diagnostics.

Assuming single-star evolution, we found that the sample stars originate from ~20–60 M_⊙ stars. They are somewhat evolved, with ages of about 2–7 Myr, but all still reside on the main sequence.

Although the ULLYSES sample was selected on the basis of an absence of indications of binarity (in existing spectra), we suspect three targets in our sample (2dFS 163, AV 83, and AV 70) to have composite spectra. Additionally, two stars, AV 372 and AV 469, show enhanced mass-loss rates compared to the rest of the sample. These two stars also show enhanced He abundances. These stars may be physically different from the rest of the sample. They might also be post-interaction binaries.

The main focus of the study lay on the stellar wind properties. Excluding the three possible binaries, we found mass-loss rates ranging from log $\dot{M}$ [M_⊙ yr⁻¹] = −7.65 to −6.02 and terminal velocities of v_∞ = 1650 to 2625 km s⁻¹. We accounted for inhomogeneities in the outflows, which we characterized with a medium that partly consisted of a spectrum of clumps that may be either optically thick or thin and that may have varying degrees of an internal velocity dispersion. Our main findings are listed below.

• The mass-loss rates of our SMC stars are consistent with the predictions of Björklund et al. (2021) to within 0.1–0.2 dex. The theoretical rates of Vink et al. (2001), Krtička & Kubát (2018), and Gormaz-Matamala et al. (2022b) produce a higher mass loss, especially at a relatively low luminosity. Care should be taken with this result. Brands et al. in prep., using the same modeling approach as we did for a set of LMC stars, find mass-loss rates that lie between the predictions of Vink et al. (2001) and Krtička & Kubát (2018). This underlines that the good agreement between empirical and predicted mass loss in one part of the parameter space is not necessarily valid for the full luminosity and metallicity range.

• Incorporating the results from Hawcroft et al. (2021) and Hawcroft et al. (2024a, Brands et al. in prep.) for Galactic and LMC stars, we found that in the luminosity interval log L/L_⊙ = 5.0–6.0, the modified wind momentum D ∝ $\dot{M}$v_∞R^1/2 is not merely a function of metallicity, but also of luminosity. Relatively low-luminosity stars show a stronger dependence on Z than high-luminosity stars. This contradicts the recent finding of Marcolino et al. (2022), who reported the opposite behavior.

• Eliminating the dependence of the terminal velocity on metallicity, we found at a luminosity log L/L_⊙ = 5.75 the same $\dot{M}$ ∝ Z^m as in the pioneering study by Mokiem et al. (2007): m = 1.02 ± 0.30 in our study versus m = 0.83 ± 0.16 in theirs.

• We used a mathematical formalism to describe D(L, Z) for O-stars in both the Galaxy, LMC and SMC. This recipe does not capture the intricate dependences fully and may need revision in the future, but it reveals that in particular, effort is needed to better populate the sample of low-luminosity (log L/L_⊙ ≲ 5.6) Milky Way stars.

• We did not identify trends of wind inhomogeneity properties with metallicity, except for a tentative finding that the clumping factor f_cl and the contrast between the inter-clump density and the mean density f_ic show higher values with higher metallicity. Although this is interesting, firm conclusions require larger samples that are studied in identical ways.

This study is part of the analysis of the ULLYSES and X-Shooting ULLYSES datasets. A full analysis of the ultraviolet and optical spectra of the ~250 targets in this program will greatly further our understanding of the outflowing atmospheres of massive stars. In parallel to these studies, time-dependent hydrodynamical simulations that resolve the clumping properties are needed to assess whether our analytical prescriptions of the wind inhomogeneity properties capture the essential characteristics of the wind structure. Such simulations are now becoming available in 2D (Debnath et al. 2024) and 3D (Moens et al. 2022).

Data availability

Supplementary material is available on https://doi.org/10.5281/zenodo.13929105

Acknowledgements

This publication is part of the project ‘Massive stars in low-metallicity environments: the progenitors of massive black holes’ with project number OND1362707 of the research TOP-programme, which is (partly) financed by the Dutch Research Council (NWO). F.B. and J.S. acknowledge the support of the European Research Council (ERC) Horizon Europe grant under grant agreement number 101044048. Based on observations obtained with the NASA/ESA Hubble Space Telescope, retrieved from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. We thank SURF (www.surf.nl) for the support in using the National Supercomputer Snellius. D.P. acknowledges financial support by the Deutsches Zentrum für Luft und Raumfahrt (DLR) grant FKZ 50 OR 2005. A.A.C.S. acknowledges support 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). A.A.C.S. further acknowledge support 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. OV acknowledges support from the KU Leuven Research Council (grant C16/17/007: MAESTRO). F.N. acknowledges grant PID2022-137779OB-C41 funded by the Spanish MCIN/AEI/ 10.13039/501100011033.

Appendix A Comparison with other works

Here we compare the best fit parameter values of individual stars to those found by Bouret et al. (2021); Martins et al. (2024); Evans et al. (2004a); Heap et al. (2006); Penny & Gies (2009). Fig. A.1 shows the parameter values they found against the ones found here. All studies, except for Penny & Gies (2009), use optical and UV spectroscopy. Martins et al. (2024, XShootU V) study a large sample of ULLYSES targets in order to determine their surface abundances. Their sample and our sample overlap with 8 stars, AV 80, AV 15, AV 95, AV 207, AV 69, AV 469, AV 307 and AV 327. Since the focus lies on stellar abundances their study does not include mass-loss properties, but does include effective temperature, surface gravity, luminosity, rotation, and abundances. Additionally, Bouret et al. (2021) studies the fundamental stellar parameters of a sample of O-type giants and super giants including mass-loss rate. In a similar study Evans et al. (2004a) determine the atmosphere properties OB supergiants in the SMC including winds. Both take clumping into account adopting a volume filling fraction of 0.1. Heap et al. (2006) derive stellar parameters, but exclude mass loss. Penny & Gies (2009) determined rotation velocities of O-type stars using FUSE data. In general Fig. A.1 shows that we find slightly higher temperatures than the other studies. This is consistent with the results of Sander et al. (2024) who compare different analysis approaches to the same data. They find the same method as applied here to result in higher temperatures than other combined optical and UV analyses. These higher temperatures then also gives rise to a higher luminosity, however, a larger scatter is visible. The derived surface gravity values are mostly consistent between the different works. The projected rotation velocities are mostly consistent, however a slight trend with higher values in this work. This could be the result of excluding a separate macro turbulent velocity in our broadening. The mass-loss rates determined here are consistent with those found by Bouret et al. (2021) and Evans et al. (2004a), except for AV 70 and AV 307. The latter has a very low mass-loss rate, possibly hindering accurate determination. AV 70, for which we find the highest mass-loss rate in the sample, is considered a poor fit and likely does not represent the physical properties of the star. We find terminal velocities consistent with or higher than the other studies. The most significant outlier, AV 307, has a low mass-loss rate and therefore more challenging terminal velocity determination.

A.1 SMC mass-loss rate

Here we compare the mass-loss rates of our sample of SMC stars and compare them to other optical and UV mass-loss rate determinations. Fig. A.2 shows the mass-loss rate as function of luminosity from Bouret et al. (2013); Ramachandran et al. (2019); Bouret et al. (2021) and Rickard et al. (2022), as well as this work. Ramachandran et al. (2019) and Rickard et al. (2022) studied the mass-loss rate of stars in the SMC wing and in SMC cluster NGC 346. Their samples are diverse, but mostly dominated by dwarf stars and stars with lower luminosities than in the sample studied here. Bouret et al. (2013) studied a sample of exclusively O-type dwarf stars. In their sample some stars with low luminosity have no significant wind features. For those stars, the values are to be considered upper limits. Similar to this work, Bouret et al. (2021) studied a sample of O-type giants and supergiants. Ramachandran et al. finds a significantly more shallow mass-loss luminosity relation, which is mainly determined by the low-luminosity stars in their sample. Rickard et al. found a very similar slope of the mass-loss rate luminosity relation, but at lower mass-loss rate for at any luminosity. They found two stars with significantly lower mass-loss rates than expected given their luminosity. The bulk of their sample lies just below the mass-loss rates inferred here. Bouret et al. find a similar slopes in both works, which are steeper than the one found here. In the luminosity range studied in this work, the mass-loss rates are comparable to the ones studied in this sample. Overall the mass-loss rates appear consistent, but with a scatter of ~1 dex, at log L/L_⊙ > 5, at lower luminosity mainly upper limits are found, with a larger scatter.

Fig. A.1 Best fit parameter values compared to Bouret et al. (2021); Martins et al. (2024); Evans et al. (2004a); Heap et al. (2006); Penny & Gies (2009). The horizontal axis shows the values found in this work, and the vertical axis the parameter values found in the literature.

Fig. A.2 Mass-loss rate plotted against luminosity for SMC stars from Ramachandran et al. (2019), Rickard et al. (2022), and this work. The dotted line indicates a linear fit to the data points. The upper limits from Rickard et al. (2022) are not considered in the linear fit, neither are stars with a luminosity log L/L⊙ < 5 from Bouret et al. (2013). Stars with poor fits, outlined in red here, are not included in the linear fits.

Appendix B CNO abundances

During the evolution of a star the surface abundances of carbon, nitrogen, and oxygen may change due to the CNO-cycle and internal mixing. However, in this process the total amount of C, N, and O is expected to remain constant at the total initial value. Fig. B.1 shows the total C + N + O surface abundance relative to hydrogen, broken down in their separate components. Additionally, the derived abundances of Martins et al. (2024, XShootU V) are included, as well as the baseline abundances from Dopita et al. (2019). Here, Martins et al. used the same dataset from XshootU. Generally, our total abundance lies above the baseline abundance. This is mostly dominated by the high and uncertain oxygen abundances we find. Additionally, we assume a fixed micro turbulent velocity of v_mic = 15 km s⁻¹, which may affect the abundance determination, including the He abundance. In all stars we find a nitrogen enhancement, except for AV 69 which is consistent with the baseline abundance. Within uncertainties our total abundance agrees with the values found by Martins et al. (2024), however, we do find larger scatter and uncertainties. This is likely due to our analysis not being tailored to abundance determination and the possible correlation between abundances and other stellar parameters.

Fig. B.1 Total carbon, nitrogen, and oxygen abundances for each star, broken down in the separate components. In each bar the uncertainties of the carbon (left), nitrogen (middle), and oxygen abundances are shown. Stars from Martins et al. (2024) that overlap with our sample are shown in the hatched bars. The horizontal dashed line shows the total baseline abundance from Dopita et al. (2019), which is also broken down in the separate components on the right.

Appendix C Spectroscopic and evolutionary mass

The mass of a star can be determined based on their surface gravity and radius. This is referred to as the spectroscopic mass, M_spec. Another way of determining the mass is though their position on the HRD and comparing that to evolutionary models. This gives the evolutionary mass, M_evo. A discrepancy between these two types of masses has been observed (e.g., Herrero et al. 1992). Fig. C.1 shows these masses for our sample and the sample of Brands et al. in prep. A clear discrepancy between the two masses remains visible, with lower spectroscopic masses below ~50 M_⊙, and higher spectroscopic masses above that threshold. This is consistent with the discrepancy found by Markova et al. (2018), who report a similar trend, but a lower “transition” mass of ~40 M_⊙. We do note that our evolutionary masses, determined with BONNSAI, see Section 3.4, are limited to M < 60 M_⊙, and that for stars with higher masses evolutionary tracks of LMC metallicity are used. This could potentially affect the evolutionary masses, increasing the discrepancy.

Fig. C.1 Spectroscopic mass plotted against the evolutionary mass. The solid black line indicates the where the two masses are equal. The LMC stars are from Brands et al. in prep. Stars with unreliable stellar parameters are marked with a red border.

Appendix D Challenges in the spectral fitting

Our spectral fitting relies on 55 spectral lines and line complexes, an overview of which is given in Table 3. To add to our understanding of quantitative spectroscopy of hot massive stars, we discuss some notable complications.

N V 1240. The N V 1240 resonance line, located in the wing of the Lyα absorption, is strong and shows a P-Cygni profile for most stars. However, it requires X-ray emission from wind-embedded shocks to have N V abundant enough to reproduce the observed profiles in our sample. X-ray emission is implemented in FASTWIND (Carneiro et al. 2016), however the exact amount of emergent emission as well as the shock temperature remain uncertain (Crowther et al. 2022). As N V is a trace ion (and the diagnostic a resonance line), the line profile is highly sensitive to the X-ray radiation field and, hence, difficult to predict reliably. Test calculations in which we included the line in the analysis resulted in systematically higher nitrogen abundances (for more details see Brands et al. in prep.). We therefore opted to exclude it.

He I singlets. In our modeling, we consider both singlet and triplet transitions of He I. Najarro et al. (2006) point out the importance of accurately predicting the radiation field near the reasonance transition 1 s^{2 1}S−2p ¹P° (at 584 Å) in order to correctly reproduce He I 4387, 4922 and 6678 Å, lines that have 2p ¹P° as their lower level. In modeling efforts, however, the radiation field near the He I resonance transition may be uncertain due to assumptions on line-blanketing and turbulent velocity. Specifically, the Fe IV ion has two transitions near to 584 Å that may affect the pumping efficiency of the resonance line, hence may impact the diagnostic He I singlet lines. Despite this, we find that in almost all cases the singlet lines are well reproduced. For AV 80 and AV 83 the singlet fits are of relatively poor quality; the triplet lines in these sources fit somewhat better.

The P V problem. Fullerton et al. (2006) signal a discrepancy between mass-loss rates determined from the P V 1118,1128 resonance doublet and those determined from Hα, termed the ‘P V problem’. The P V doublet would systematically result in lower mass-loss rates. The authors suggest the discrepancy can be solved by introducing clumping in the outflow. Such clumps increase the strength of the Hα recombination line while the P V line is not affected by changes in density directly, though it may be through changes in ionization as a result of the higher density in the clumps (see Oskinova et al. 2007). Šurlan et al. (2013) successfully reproduce both the P V lines and Hα line with the same mass-loss rate by including micro clumping (analogous to f_cl) and macro clumping, the latter being a measure for the porosity of the wind.

Despite our relatively sophisticated and flexible prescription of wind inhomogeneities and leaving the phosphorus abundance a free parameter, we do not get satisfactory fits to the P V lines for part of our sample. This may suggest the applied clumping prescription fails to capture intricacies of wind structure key for the formation of P V or some other cause. Though, often we find wind features in the model, while for all stars we observe photospheric line profiles. Issues with photospheric line profiles are unlikely to be resolved with wind clumping properties. Possibly the problem is connected to the poor constraints on the phosphorous abundance.

Appendix E Optical results

Here we briefly discuss the optical-only fits and how they compare to the optical + UV fits. Table E.1 lists the best fit values of the free parameters used in the optical-only fits. Here we fixed β = 1, f_cl = 10, f_ic = 0.1, f_vel = 0.5, v_cl,start = 0.05, v_windturb = 0.1, and v_∞ from Hawcroft et al. (2024b). Generally, the optical-only fits agree well with the fits using the full set of diagnostics. The most significant deviations are found in the lower mass-loss rates, with the optical-only fits resulting in significantly higher values at the low mass loss end. The most extreme example of this is AV 307, where the mass-loss rate is 1.3 dex higher in the optical only fit. The available diagnostics in the optical are not well suited for constraining low rates, as the sensitivity of the line profiles to the mass-loss rate drops significantly. The mass-loss rate needs to increase a significant amount to compensate a small filling in of the Hα line profile. The wind diagnostics in the UV remain sensitive to lower values of the mass-loss rate. This indicates the importance of the UV diagnostic features in constraining low mass-loss rates.

Appendix F Mass loss relations

Here we show the same analysis as Section 5.3, but now fitting directly the mass-loss rates, rather than the modified wind momentum.

Figure F.1 shows the mass-loss rate as function of the luminosity for each of the samples. The dotted lines indicate linear fits to the individual samples. The slope and offsets of these fits are listed in Table F.1. A general trend of increased mass-loss rate for increased metallicity is visible. We fitted the linear fits at log L/L_⊙ = 5.75 to determine the mass loss metallicity relation, and find d log $\dot{M}$ / d log Z = m = 1.01 ± 0.26, consistent with the wind momentum results.

Fig. F.1 Mass-loss rates as function of luminosity for samples of SMC (this work), LMC (Hawcroft et al. 2024a, Brands et al. in prep.), and MW stars (Hawcroft et al. 2021). The dotted lines indicate the linear fits to the data for each of the metallicities, with the shaded region the 1 σ confidence interval on the fit. Stars with poor spectral fits, marked with a red border, are left out the linear fit here. The parameters of the fit are listed in Table F.1.

Appendix G Combined metallicity and luminosity dependence

Here we supply additional information on the combined fit to the modified wind momentum using both a metallicity and luminosity dependence simultaneously. Fig. G.1 shows the fit as displayed in Fig. 5 together with the theoretical predictions of (Krtička & Kubát 2018) and (Björklund et al. 2021). The coefficients of wind momentum relations are listed in Table 8. Most notably the theoretical predictions lie much closer together at low luminosity. The Solar metallicity wind momentum we find at low metalicity is significantly higher than the predictions. This is likely due to a lack of low luminosity stars in our sample. Towards higher luminosities of 10⁶ L_⊙ the predictions and observations converge.

The covariance matrix of the fit of Equation 9 to the empirically determined modified wind momentum rates is $$Cov(a,b,c,d)=\left[\begin{array}{llll}0.18& 0.34& 0.06& 0.03\\ 0.34& 0.90& 0.17& 0.06\\ 0.06& 0.17& 0.06& 0.02\\ 0.03& 0.06& 0.02& 0.01\end{array}\right].$$(G.1)

Fig. G.1 Empirical wind momentum relation with luminosity and metallicity compared to the theoretical predictions. The displayed metallicities are Solar, half Solar, and one fifth Solar, from top to bottom for each relation.

Appendix H Fit overview

Fig. H.1 Overview of the best fits of the UV and wind sensitive optical lines. The normalized observed data is plotted as black vertical bars with the size of the bar indicating the 1σ uncertainty on the data. The green line indicates the best-fit model and the shaded region is the 1σ confidence interval of the models. The thin dashed grey lines are at unity. The major ticks on the x-axis indicate steps of 5 Å, and the minor ticks indicate steps of 1 Å. On the y-axis the major ticks indicate steps of 10% of the continuum and the minor ticks indicate 5%. A more detailed figure of the line profiles can be found online https://doi.org/10.5281/zenodo.13929105

Fig. H.2 As Figure H.1, but for optical lines of mostly photospheric origin.

References

All Tables

All Figures

Fig. 1 Hertzsprung-Russell diagram of our SMC sample of stars with temperatures and luminosities resulting from the fit. Points with a red border indicate stars with unreliable parameters. We overplot the evolutionary tracks (thin solid gray lines) and isochrones (thin dashed gray lines) of Brott et al. (2011). The solid black line indicates the zero-age main sequence. The points are labeled with the name of the object in a small font that is visible when zoomed in.

In the text

Fig. 2 Distribution of the surface mass fractions of He, C, N, O, Si, and P from the FASTWIND/GA fitting of SMC stars. The solid green line indicates an approximate kernel density estimate of the distribution taking the uncertainties on the best-fit values into account. The thin dashed lines show the contribution of the individual sources. The vertical black line shows the baseline mass fraction from Dopita et al. (2019) for all elements except phosphorus, for which we show the scaled (by 0.2) solar abundance from Asplund et al. (2009). The area marked in red shows the surface depletion or enrichment that can be expected due to mixing based on the evolution of a 40 M⊙ star with an initial rotation of 389 km s−1 (Brott et al. 2011). For phosphorus, only stars with a reasonable fit are included. The helium mass fractions are taken from the optical-only fit. All other fractions are from the optical and UV fit.

In the text

Fig. 3 Mass-loss rate as a function of luminosity. The fit to the empirical mass-loss rates and its 1σ uncertainty are indicated with the dotted green line and the shaded region, respectively. The mass-loss rate – luminosity relations of Krtička & Kubát (2018) and Björklund et al. (2021) are shown with the dashed blue and yellow lines, respectively. For the mass-loss predictions of Vink et al. (2001), we used their equation (24) to obtain the predicted rates for each star based on their properties. The magenta V indicate these rates, and the dashed magenta line corresponds to the fit to these points. The light blue stars indicate the mass-loss rates of equation 7 from Gormaz-Matamala et al. (2022b). Points with a red border are unreliable fits that are likely not representative of their physical properties; these are not included in the fit.

In the text

Fig. 4 Modified wind momentum as function of luminosity for samples of SMC (this work), LMC (Hawcroft et al. 2024a, Brands et al., in prep.), and Galactic stars (Hawcroft et al. 2021). The dotted lines indicate the linear fits to each subsample, and the shaded regions show their 1σ uncertainties. The points with a red border represent unreliable fits that are likely not representative of their physical properties; these are not included in the fit.

In the text

Fig. 5 Modified wind momentum as function of luminosity and metallicity, as in Figure 4, but now fit with Equation (9) to the three samples, instead of separate linear fits for each sample. The dashed magenta lines and pink shaded regions indicate the best fit and the 1σ confidence interval, respectively. Equation (9) has been evaluated at Z = 1.0, 0.5, and 0.2 Z⊙. The points with a red border are unreliable fits that are likely not representative of their physical properties; these are not included in the fit.

In the text

	Fig. 6 Fits to the empirically determined modified wind momentum as a function of luminosity. The solid colored lines show the fits to the MW, LMC, and SMC samples presented here. The dashed and dash-dotted black and gray lines indicate the modified wind momentum relations from Mokiem et al. (2007) and Marcolino et al. (2022), respectively. The values representing the fits are listed in Table 7.
In the text

Fig. 7 Slope of the mass-loss rate as a function of metallicity, m, for different values of the luminosity. The top panel shows the empirical determinations of Mokiem et al. (2007), Ramachandran et al. (2019), and Marcolino et al. (2022), along with our determination based on the linear fits as shown in Figs. 4 and 6. The empirical values of m from Rickard et al. (2022) lie above the vertical range of this figure. The bottom panel shows the predicted values of m of Krtička & Kubát (2018) and Björklund et al. (2021), along with our empirical determination. The value of m = 0.83 from Mokiem et al. (2007) is also shown for reference. In all cases, d log v∞/d log Z = n = 0.13 is assumed (Leitherer et al. 1992).

In the text

	Fig. 8 Distributions of the wind parameters of LMC and SMC XShootU samples (this work and Brands et al. in prep.). The solid lines indicate the approximate KDE based on the uncertainty on the parameter values. The thin dashed lines indicate the contribution of the individual objects to the total KDE. Objects for which the parameters are unconstrained are excluded.
In the text

	Fig. A.1 Best fit parameter values compared to Bouret et al. (2021); Martins et al. (2024); Evans et al. (2004a); Heap et al. (2006); Penny & Gies (2009). The horizontal axis shows the values found in this work, and the vertical axis the parameter values found in the literature.
In the text

Fig. A.2 Mass-loss rate plotted against luminosity for SMC stars from Ramachandran et al. (2019), Rickard et al. (2022), and this work. The dotted line indicates a linear fit to the data points. The upper limits from Rickard et al. (2022) are not considered in the linear fit, neither are stars with a luminosity log L/L⊙ < 5 from Bouret et al. (2013). Stars with poor fits, outlined in red here, are not included in the linear fits.

In the text

Fig. B.1 Total carbon, nitrogen, and oxygen abundances for each star, broken down in the separate components. In each bar the uncertainties of the carbon (left), nitrogen (middle), and oxygen abundances are shown. Stars from Martins et al. (2024) that overlap with our sample are shown in the hatched bars. The horizontal dashed line shows the total baseline abundance from Dopita et al. (2019), which is also broken down in the separate components on the right.

In the text

	Fig. C.1 Spectroscopic mass plotted against the evolutionary mass. The solid black line indicates the where the two masses are equal. The LMC stars are from Brands et al. in prep. Stars with unreliable stellar parameters are marked with a red border.
In the text

Fig. F.1 Mass-loss rates as function of luminosity for samples of SMC (this work), LMC (Hawcroft et al. 2024a, Brands et al. in prep.), and MW stars (Hawcroft et al. 2021). The dotted lines indicate the linear fits to the data for each of the metallicities, with the shaded region the 1 σ confidence interval on the fit. Stars with poor spectral fits, marked with a red border, are left out the linear fit here. The parameters of the fit are listed in Table F.1.

In the text

	Fig. G.1 Empirical wind momentum relation with luminosity and metallicity compared to the theoretical predictions. The displayed metallicities are Solar, half Solar, and one fifth Solar, from top to bottom for each relation.
In the text

Fig. H.1 Overview of the best fits of the UV and wind sensitive optical lines. The normalized observed data is plotted as black vertical bars with the size of the bar indicating the 1σ uncertainty on the data. The green line indicates the best-fit model and the shaded region is the 1σ confidence interval of the models. The thin dashed grey lines are at unity. The major ticks on the x-axis indicate steps of 5 Å, and the minor ticks indicate steps of 1 Å. On the y-axis the major ticks indicate steps of 10% of the continuum and the minor ticks indicate 5%. A more detailed figure of the line profiles can be found online https://doi.org/10.5281/zenodo.13929105

In the text

	Fig. H.2 As Figure H.1, but for optical lines of mostly photospheric origin.
In the text